UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Integrated field investigation, numerical analysis and hazard assessment of the Portillo Rock Avalanche.. 2008

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

Item Metadata

Download

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

Full Text

INTEGRATED FIELD INVESTIGATION, NUMERICAL ANALYSIS AND HAZARD ASSESSMENT OF THE PORTILLO ROCK AVALANCHE SITE, CENTRAL ANDES, CHILE   by   DANIELA WELKNER  B.Sc., University of Chile, 1996   A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF   MASTER OF SCIENCE   in   THE FACULTY OF GRADUATE STUDIES  (Geological Engineering)    THE UNIVERSITY OF BRITISH COLUMBIA  (Vancouver)     April 2008   © Daniela Welkner, 2008  ii ABSTRACT This thesis reports a detailed rock slope hazard investigation located in the rugged mountains of the Andean Cordillera of Central Chile, close to the border between Chile and Argentina. A number of large lobate-shaped diamicton deposits were dated by 36Cl and the results showed that they corresponded to two significant rock mass wasting events: the Upper Pleistocene Portillo Rock Avalanche (PRA) and a Holocene rock slump and rockslide. The pre-historic landslide deposits underlie a key transportation route between Chile and Argentina, as well as an important ski resort. The purpose of this research was to investigate the likely failure mechanism and characterise the runout path and volume constraint of the PRA utilizing a “Total Slope Analysis” which relates the rock slope failure analysis with the runout analysis. The insights gained on the back analysis of the slope were used in later stages of the study to assess the hazard potential of a recurring major rockslide which eventually could affect the international corridor.  The distinct element code UDEC was used to evaluate the failure mechanism and identify the geological model that best reproduced the observed failed state. Elasto-plastic modelling results showed that a stress-controlled failure at the toe of the slope together with a structurally-controlled failure in its upper part likely had occurred. Runout analyses were carried out using DAN3D, a dynamic analysis code in which combinations of rheologies were tested and ranked based on their ability to represent the current distribution of the debris. To do so, a pre-failure topography reconstruction and volume estimates of the deposits were performed for the two events. Results showed that the best basal rheological combination, as constrained by the outline of the PRA, was frictional during the rockslide and a Voellmy when entrainment became important. In contrast, the rheology that best represented the Holocene rock slump was a constant frictional rheology.  The performance of the present-day state of the slope was tested under several different triggering scenarios. The first test was carried out using the same rock mass and discontinuity properties as those back calculated. Under this condition the slope proved to be stable; an expected result considering that there have not been major failures in historic times, thus pointing towards a stabilized geometrical profile with time. Next, a coupled hydro-mechanical analysis was performed involving increased pore water pressures at the toe of the slope; nevertheless they were not sufficient to overcome the strength of the rock mass resulting in a stable slope. Finally the modelled slope was subjected to a seismic load to simulate the effect of an M = 7.8 earthquake at the base of the slope. Under this highly unlikely dynamic condition for the site, the crest of the slope failed due to an outward rotation of blocks, probably aided by topographic amplification. The hazard assessment proceeded accounting for the estimated failed volume, for which the runout simulations showed that for both a highly frictional and non-saturated path (dry season) or a snow covered path (winter), the leading edge of the flow would override part of the International Santiago-Mendoza Corridor with the debris coming to rest in a flat-lying area in the upper part of the valley. In both cases, the Portillo Ski Resort shows no direct impact; nevertheless the close proximity with the edge of the deposit could represent potential harm towards the structure involving individual boulders. Overall, though, the results from this integrated study suggest that the hazard level is very low.   iii TABLE OF CONTENTS ABSTRACT .............................................................................................................................. ii TABLE OF CONTENTS ........................................................................................................... iii LIST OF TABLES .................................................................................................................... vii LIST OF FIGURES .................................................................................................................. ix ACKNOWLEDGEMENTS ...................................................................................................... xv DEDICATION ......................................................................................................................... xvi 1. INTRODUCTION .............................................................................................................. 1 1.1. Problem Statement ................................................................................................... 1 1.2. Research Objectives ................................................................................................ 3 1.2.1. Mechanisms of failure ....................................................................................... 3 1.2.2. Runout deposits ................................................................................................ 3 1.2.3. Forward modeling ............................................................................................. 3 1.3. Thesis Structure ....................................................................................................... 4 1.4. Literature Review ..................................................................................................... 5 1.4.1. Rock mass classification systems .................................................................... 5 1.4.1.1. Hoek-Brown failure criteria ..................................................................... 10 1.4.1.2. Mohr-Coulomb failure criteria ................................................................. 12 1.4.2. Slope failure initiation ..................................................................................... 14 1.4.2.1. Mechanisms ............................................................................................ 14 1.4.2.2. Numerical Modeling ................................................................................ 17 1.4.2.3. Triggering factors .................................................................................... 18 1.4.3. Runout prediction of failed slopes .................................................................. 19 1.4.3.1. Empirical methods .................................................................................. 19 1.4.3.2. Numerical methods ................................................................................. 21 2. METHODOLOGY ........................................................................................................... 22 2.1. Total Slope Analysis Procedure ............................................................................. 22 2.2. Field Work and Data Collection .............................................................................. 24 2.2.1. Field mapping ................................................................................................. 24 2.2.2. Rock sampling ................................................................................................ 27 2.2.3. Discontinuity mapping .................................................................................... 27 2.3. Terrestrial Laser Scanning (LiDAR) ....................................................................... 29 2.3.1. Data acquisition .............................................................................................. 29  iv 2.3.2. Data processing and analysis ......................................................................... 31 2.4. Computer-Based Techniques ................................................................................. 33 2.4.1. Data processing .............................................................................................. 33 2.4.1.1. Spheristat 2.2 .......................................................................................... 33 2.4.1.2. RocLab 1.0 ............................................................................................. 33 2.4.1.3. Surfer 8.0 ................................................................................................ 34 2.4.1.4. ArcGIS 9.0 .............................................................................................. 34 2.4.2. Distinct-element analysis (UDEC) .................................................................. 34 2.4.3. Runout analysis (DAN3D) .............................................................................. 35 2.5. Numerical Modeling Procedure .............................................................................. 35 2.5.1. Static (non-inertial) solution ............................................................................ 35 2.5.2. Coupled hydro-mechanical analysis – steady state flow ................................ 37 2.5.3. Dynamic solution ............................................................................................ 39 3. ENGINEERING GEOLOGY OF THE STUDY AREA ..................................................... 42 3.1. Introduction ............................................................................................................. 42 3.2. Tectonic Setting & Seismicity ................................................................................. 43 3.3. Geology & Geomorphology .................................................................................... 47 3.3.1. Bedrock & structural geology .......................................................................... 47 3.3.2. Quaternary geology ........................................................................................ 50 3.4. Age Dating & History of Rock Slope Failure Events ............................................... 58 3.4.1. Theory & background ..................................................................................... 58 3.4.2. Cosmogenic nuclide dating by 36Cl ................................................................. 59 3.4.3. Sampling ......................................................................................................... 60 3.4.4. History of events & failure scenarios .............................................................. 62 3.4.4.1. Upper Pleistocene Portillo Rock Avalanche ........................................... 65 3.4.4.2. Holocene rockslide and rock slump ........................................................ 67 3.5. Geotechnical Characterization of the East Slope ................................................... 68 3.5.1. Description of discontinuity network ............................................................... 68 3.5.1.1. Orientation .............................................................................................. 71 3.5.1.2. Spacing & block volume ......................................................................... 75 3.5.2. Laboratory geomechanical testing .................................................................. 75 3.5.2.1. Intact rock ............................................................................................... 76 3.5.2.2. Discontinuities ......................................................................................... 79 3.5.3. Rock mass characterization & properties ....................................................... 80  v 3.5.3.1. Rock Mass Rating (RMR’) & Geological Strength Index (GSI) ............... 81 3.5.3.2. Adjusted rock mass properties for numerical models ............................. 84 3.5.4. Hydrogeology ................................................................................................. 87 4. FAILURE MECHANISM OF THE PREHISTORIC PORTILLO ROCK AVALANCHE .... 88 4.1. Model Development ............................................................................................... 88 4.1.1. Failure initiation .............................................................................................. 88 4.1.2. Pre-failure topography & failed volume estimate ............................................ 89 4.1.3. Model configuration ........................................................................................ 91 4.2. Base Model Results ............................................................................................... 97 4.3. Parametric Analysis .............................................................................................. 102 4.3.1. Variation in geometry .................................................................................... 103 4.3.1.1. Homocline geometry ............................................................................. 103 4.3.1.2. Toe-buttress geometry .......................................................................... 105 4.3.1.3. Block size reduction .............................................................................. 108 4.3.2. Variations in block constitutive models and material behaviour ................... 110 4.3.2.1. Elasto-plastic base ................................................................................ 110 4.3.2.2. Zone of damaged rock as an equivalent continuum ............................. 112 4.3.3. Variation of mechanical properties ............................................................... 114 4.3.3.1. Rock Mass Strength ............................................................................. 114 4.3.3.2. Joint friction ........................................................................................... 116 4.4. Interpretation of Results ....................................................................................... 118 5. RUNOUT ANALYSIS OF THE PREHISTORIC ROCK SLOPE FAILURE EVENTS .... 121 5.1. Introduction ........................................................................................................... 121 5.2. Model Development ............................................................................................. 122 5.2.1. Volume calculations ...................................................................................... 123 5.2.2. Pre-failure topography reconstruction .......................................................... 125 5.2.3. Rheological relationships .............................................................................. 130 5.3. Model Results ....................................................................................................... 131 5.3.1. The Portillo Rock Avalanche ........................................................................ 131 5.3.2. The Holocene rock slump lobe ..................................................................... 138 5.4. Interpretation of results ......................................................................................... 143 6. FORWARD ANALYSIS OF POTENTIAL LARGE-SCALE ROCKSLIDE EVENTS AT PORTILLO ........................................................................................................................... 145 6.1. Potential Failure Initiation & Volume Assessment ................................................ 146  vi 6.1.1. Model development and results .................................................................... 146 6.1.2. Varying water table ....................................................................................... 148 6.1.3. Seismic triggering (seismic loading of slope) ............................................... 151 6.2. Runout Assessment of Potential Failed Volume .................................................. 156 6.3. Hazard Assessment Summary ............................................................................. 161 7. CONCLUSIONS AND RECOMMENDATIONS ............................................................ 163 7.1. Conclusions .......................................................................................................... 163 7.2. Recommendations for Further Work .................................................................... 165 REFERENCES ..................................................................................................................... 167 APPENDIX A: GEOMECHANICAL LABORATORY TEST RESULTS ................................. 175 APPENDIX B: ROCLAB OUTPUTS ..................................................................................... 187 APPENDIX C: NUMERICAL MODELLING CODES ............................................................ 192     vii LIST OF TABLES Table 2.1. Terrestrial laser scan measurement stations and associated parameters. ....... 29 Table 3.1. Historical and recent seismic events with mb or Ms greater than 5.0. Data taken from Earthquake Hazard Program – USGS (http://neic.usgs.gov/). ................. 46 Table 3.2. Data and ages obtained through 36Cl cosmogenic nuclides. ............................ 62 Table 3.3. Summary of discontinuity parameters obtained from the eastern slope. .......... 70 Table 3.4. Average values of intact rock strength parameters for uniaxial tests. ............... 77 Table 3.5. Triaxial compression test results. ...................................................................... 78 Table 3.6. Direct shear tests results for samples LI-7.1 & LI-7.2 containing mate bedding planes. .............................................................................................................. 79 Table 3.7. Rock mass classification at D1 & D2. ............................................................... 83 Table 3.8. Estimated range of intact rock and rock mass parameters for D1 & D2. .......... 85 Table 3.9. Estimated range of discontinuity strength and stiffness parameters. ................ 86 Table 4.1 Basic assumptions made for the development of the numerical models. ......... 91 Table 4.2. Mechanical properties of the rock mass used for the base model after the initial state was achieved. Source of information same as for Table 3.8. .................. 95 Table 4.3. Mechanical properties and geometric parameters of discontinuities for the base model. ............................................................................................................... 96 Table 4.4. Varying strength properties used in the staging process to develop the base model. ............................................................................................................... 97 Table 4.5. Mechanical parameters of the far-field base modeled as a continuum with a Mohr-Coulomb elasto-plastic constitutive criterion. ........................................ 110 Table 4.6. Range of rock mass strength properties tested for the parametric study. ...... 114 Table 4.7. Joint friction range tested for the parametric study. ........................................ 116 Table 4.8. Qualitative estimation on the influence of the varying parameters in the proposed failure mechanism and volume failed. ............................................. 120 Table 5.1. Geometrical parameters of the prehistoric mass movements. ........................ 125 Table 5.2. Basal rheology and parameter values used in the Portillo Rock Avalanche dynamic analysis. ............................................................................................ 132  viii Table 5.3. Basal rheology and parameters that better conform with the Holocene rock slump deposit. ................................................................................................. 141 Table 6.1 Geometric configurations and mechanical properties of the forward analysis for the initial state. ................................................................................................ 146 Table 6.2. Material properties used in the fluid flow models. ........................................... 149 Table 6.3. Input variables for the UDEC dynamic analysis. ............................................. 153 Table 6.4. Geometrical parameters and rheology of the potential failed volume. ............ 157 Table 6.5. Susceptibility to failure of the eastern slope. ................................................... 161 Table 6.6. Hazard assessment for the eastern slope. ..................................................... 162   ix LIST OF FIGURES Figure 1.1. Photos of the Portillo Rock Avalanche deposits which underlie the Santiago- Mendoza International Corridor and the Portillo Ski Resort. ............................... 2 Figure 1.2. Rock Mass Rating classification system (modified from Bieniawski 1989). ........ 6 Figure 1.3. General chart for GSI estimates based on geological observations (modified from Marinos et al. 2005). ................................................................................... 9 Figure 1.4. The Mohr-Coulomb failure criterion (modified from Hudson & Harrison 2005). 12 Figure 1.5. Approximation of M-C failure envelope from the H-B envelope (modified from Cai et al. 2004). ................................................................................................. 13 Figure 1.6. Illustration of the transition zone utilizing the Prandtl´s prism (modified from Kvapil & Clews 1979 ). ...................................................................................... 15 Figure 1.7. a) Distribution of transversal deformation with respect of the position of the Prandtl transition zone; b) Distribution of the degree of fracturing with respect of the position of the Prandtl transition zone (modified from Kvapil & Clews 1979 ).  .......................................................................................................................... 16 Figure 1.8. Diagram of the fahrböschung (or angle of reach), where H is the vertical drop and L is the horizontal projection of the distal margin of the displaced mass. .. 20 Figure 2.1. Diagram of a total slope failure analysis showing the codes that will be used in this thesis for each specific area (modified from Stead et al., 2006). ............... 22 Figure 2.2. Flow chart of the procedure followed by the Total Slope Back Analysis and Total Slope Forward Analysis (modified from Strouth 2006). ........................... 23 Figure 2.3. Pictures obtained from the eastern wall which shows different views of the scarp at different distances. .............................................................................. 25 Figure 2.4. Aerial photographs at 1:20,000 (top) and 1:80,000 (bottom) scales which were used to map the debris in the runout zone. ....................................................... 26 Figure 2.5. LiDAR (Optech ILRIS-3D) measurement stations. ........................................... 28 Figure 2.6. Approximate location and orientation of the survey stations. ............................ 30 Figure 2.7. Sequence of steps followed during the data processing consisting on scanning of the outcrop, point cloud generation and stereonet analysis. ......................... 32 Figure 2.8. Solution procedure for static analysis (modified from Itasca 2004). ................. 36 Figure 2.9. Diagram showing the relationship between fracture conductivity, mechanical deformation, and water pressures in a coupled hydro-mechanical analysis (modified from Zangerl 2003). ........................................................................... 38  x Figure 2.10. Types of dynamic loading and boundary conditions for a flexible base (left) and a rigid base (right). Modified from Itasca 2004. ................................................. 40 Figure 3.1. Shaded relief map showing the location of the study area ............................... 42 Figure 3.2. Major morphostructural units of the central Andes and schematic block diagrams showing the subducted plate segments. ........................................... 44 Figure 3.3. Earthquakes with mb or Ms greater than 4.0 and their epicentral area. Modified from Earthquake Hazard Program – USGS. ..................................................... 45 Figure 3.4. Google Earth picture showing the structures that are recognized in the study area. Note the relief in the terrain. .................................................................... 48 Figure 3.5. Structural map and cross section of the Aconcagua Fault and Thrust Belt (modified from Ramos et al. 2004). ................................................................... 48 Figure 3.6. Geological map of the study area. Modified from Rivano et al. (1993). ............ 49 Figure 3.7. Photos of the Portillo Rock Avalanche deposit along its path.. ......................... 51 Figure 3.8. Photo of a large epidote-rich nodule. ................................................................ 52 Figure 3.9. Aerial photo 1:20,000 which shows the ridge in the upper part of the valley and a detail of the debris. ......................................................................................... 52 Figure 3.10. Other landslide deposits recognized in the area and assigned to the Pleistocene. ...................................................................................................... 54 Figure 3.11. Pictures of the Holocene deposits showing the eastern (top) and western (bottom) deposits. ............................................................................................. 55 Figure 3.12. Pictures of the alluvial plateaus along the stepped topography of the glacial valley. ................................................................................................................ 56 Figure 3.13. Picture of the talus developed in the valley sides. ............................................ 57 Figure 3.14. Location of the 17 samples obtained for surface exposure dating by Cosmogenic Nuclide using 36Cl.. ...................................................................... 60 Figure 3.15. Pictures showing size, shape, and inclination of sampled boulders. All samples were collected with a hammer and a chisel. ..................................................... 62 Figure 3.16. Aerial photograph (1:80,000) showing the distribution of the slope failure events and their approximate ages.. ............................................................................ 64 Figure 3.17. Cross sections of the potential pre-failure scenarios. Only the predominant joint sets are considered.. ......................................................................................... 65 Figure 3.18. Current topographic profile of the eastern slope showing the different degree of fracturing affecting D1 & D2.. ............................................................................ 69  xi Figure 3.19. Contoured stereographic projection of selected automatically generated patches from scans 1-7, 9, 11-12, together with outcrop measurements, for D1.   .................................................................................................................... 71 Figure 3.20. Detailed picture of the main fracture and stereographic projection showing the kinematics between JS1 & F1. .......................................................................... 73 Figure 3.21. Contoured stereographic projection of hand measurements taken in D2 for JS2 & JS3. ............................................................................................................... 74 Figure 3.22. Picture of samples for a) Uniaxial loading, sample LI-10, H/d ratio of 2:1; b) Triaxial loading, sample LI-8.2, H/d ratio of 2:1. ............................................... 76 Figure 3.23. Shear failure of samples after failure subject to uniaxial and triaxial loading. ... 77 Figure 3.24. Mohr stress circles and shear strength failure envelope obtained from triaxial testing of intact rock. ......................................................................................... 78 Figure 3.25. Failure envelopes for peak and residual shear strength (τ, σn space) from direct shear testing of mated discontinuity samples for JS1. ...................................... 79 Figure 3.26. Detailed pictures obtained from the rough surface of the principal joint set (JS1). ................................................................................................................ 80 Figure 3.27. RMR’89 and GSI values for D1 & D2. ................................................................ 82 Figure 3.28. Diagram showing the relationship between the block size and rock strength properties.. ........................................................................................................ 84 Figure 4.1. a) Contour map of the study area with vector indicating steepness of the slope; b) Tin map showing elevation and location of the cross-section. ...................... 89 Figure 4.2. Area covered by the Portillo Rock Avalanche deposits and its relationship to the volume as proposed by Li (1983).. .................................................................... 90 Figure 4.3. Flow chart of the modelling procedure utilised in this thesis. ............................ 92 Figure 4.4. Concave geometry without daylighting structures at the toe of the slope used for the numerical models in UDEC.. .................................................................. 92 Figure 4.5. Horizontal (σxx) and vertical (σyy) principal stress contours for the initial state of the base model. ................................................................................................ 94 Figure 4.6. Evolution of horizontal displacements for damage states 2 and 4. ................... 97 Figure 4.7. Evolution of vertical displacements and total displacement vectors for damage states 2 and 4.. ................................................................................................. 98 Figure 4.8. Plasticity indicators and shear displacements. Arrows indicate directions of shear displacement. .......................................................................................... 99  xii Figure 4.9. a) Plasticity indicators and shear displacements for modelled state 3; b) Northwest facing view of the slope showing the presence of JS2 acting as release surfaces. ............................................................................................... 99 Figure 4.10. Grid plot magnified 30x showing the bulge evolution from state 2 to state 4. . 100 Figure 4.11. Plot of the principal stress tensors. Blue diamonds indicate principal stresses in tension. ........................................................................................................... 101 Figure 4.12. Flow diagram showing the sequence of trial models. ..................................... 102 Figure 4.13. Homocline geometry. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators; d) Grid plot magnified 30 times. .......................................................................... 104 Figure 4.14. Zoomed in plots of vertical stresses (σyy) in domain 2 for the concave model and the homocline model. ............................................................................... 105 Figure 4.15. E-W cross sections across the Inca Lake and location map showing a protruding ridge along the toe of the Caracoles Range. ................................. 106 Figure 4.16. Toe-buttress geometry. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators. ... 107 Figure 4.17. Reduced block size. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators. ... 109 Figure 4.18. Far-field base modeled with a Mohr-Coulomb constitutive criterion. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators; d) Horizontal principal stress contours (σxx). ................................................................................................ 111 Figure 4.19. Domain 2 as an equivalent continuum. a) Horizontal displacements for state 4; b) Vertical displacements for state 4; c) Plasticity and shear displacement indicators. ........................................................................................................ 113 Figure 4.20. Trial models with varying friction angle. a) Plasticity and shear displacement indicators for upper bound friction angle values; b) History plot of horizontal displacements for upper bound values of friction angle. ................................. 115 Figure 4.21. Trial model with varying tension. a) History plot of horizontal displacements for upper bound values of tension; b) History plot of vertical displacements for upper bound values of tension. ....................................................................... 116 Figure 4.22. Trial models with varying joint friction angle. a) Plasticity and shear indicators for lower, average, and high frictional strength values; b) Shear strain plots for lower, average and high frictional strength values. ......................................... 117 Figure 5.1. Digital terrain model (DTM) constructed through the interpolation of a 1:50,000 digital contour map with 50 m contours (PSAD 56, zone 19S). ...................... 122  xiii Figure 5.2. Aerial photograph (1:80,000) showing the distribution of the mass wasting events and approximate ages. H indicates location of the ski resort. ............. 123 Figure 5.3. Picture comparing the shape of the western lobe with that of a wedge, and indication of the values used for the volume calculation. ................................ 124 Figure 5.4. a) Illustration that shows the Sloping Local Base Level (SLBL) concept which corresponds with the lowest position of a sliding surface (taken from Jaboyedoff, 2005); b) Principles of the SLBL concept applied to the reconstruction of the source area of the PRA. ................................................ 126 Figure 5.5. DTM and cross sections across and along the valley used to reconstruct the pre-failure topography. NE-SW-2 cross section shows the profile of the path of the PRA (see Figure 5.7 for details). ............................................................... 127 Figure 5.6. Pre-failure (top) and post-failure (bottom) topography showing the valley reconstruction.. ............................................................................................... 128 Figure 5.7. NE-SW cross section parallel to the path followed by the Portillo Rock Avalanche showing the chronology of events and its related surfaces (see Figure 5.5 for location of profile). .................................................................... 129 Figure 5.8. Cross section showing the change in basal rheologies according to the change in basal material strength along the path. ....................................................... 132 Figure 5.9. Calibrated simulation of the Portillo Rock Avalanche event at 50 s intervals.. 134 Figure 5.10. Terrain slope map showing major zones of mapped deposition. Contouring of slope angle in the range 0° – 7° has been highlighted. .................................. 135 Figure 5.11. Plot of maximum simulated flow velocities. Maximum velocity contours are at 10 m/s and elevation contours are at 100 m intervals. ........................................ 136 Figure 5.12. Calibrated simulation of the Portillo Rock Avalanche event utilising a constant basal frictional rheology with a basal friction angle of 20°.. ............................ 137 Figure 5.13. Simulations for the Holocene rock slump with varying bulk basal friction angle.. .   .................................................................................................................. 140 Figure 5.14. Calibrated simulation of the Holocene rock slump lobe (S2) utilising a constant basal frictional rheology with a basal friction angle of 33°.. ............................ 142 Figure 5.15. Aerial photograph (1:20,000) showing in detail the Holocene Rock Slump deposit. ........................................................................................................... 143 Figure 6.1. Flow diagram showing the procedure followed for the hazard assessment. .. 145 Figure 6.2. Stable condition for the modelling of the current slope using lower bound strength properties.. ........................................................................................ 147 Figure 6.3. Current profile (topography) of the eastern slope showing the range in which the water table is considered to vary. ............................................................. 149  xiv Figure 6.4. Forward modelling and volume assessment for future events considering a fully saturated slope. .............................................................................................. 151 Figure 6.5.  State of the Caracoles Range after the input of an M=5.5 earthquake.. ........ 154 Figure 6.6. State of the Caracoles Range after the input of an M=7.8 earthquake. .......... 155 Figure 6.7. Views of the zone affected by the outward rotation of blocks.. ....................... 156 Figure 6.8. DAN3D runout assessment of a potential outward rotation of blocks at the crest of the slope. .................................................................................................... 158 Figure 6.9. Maximum velocities of the potential failed volume. a) Accounting for a dry path; b) Accounting for a path covered by snow. ..................................................... 160   xv ACKNOWLEDGEMENTS Foremost I would like to express my sincere gratitude to my supervisor Dr. Erik Eberhardt. His mentoring and guidance, together with the financial support during this process was invaluable. I would also like to thank Dr. Oldrich Hungr for his advice and continuous support during my research. Their valuable suggestions and discussions, as well as the numerous hours of individual and classroom instruction were fundamental in the development of this thesis.  I wish to express my gratitude to Dr. Reginald Hermanns, member of my committee, and to my external examiner Dr. Wayne Savigny, for their time and effort on my behalf.  I also would like to thank the Geological and Mining Survey of Chile (SERNAGEOMIN) for funding the field work and providing the logistics. Particular thanks to Ms. Renate Wall, Dr. Luis Lara, Dr. Estanislao Godoy (Pirzio), Dr. Jorge Clavero, and Dr. Jorge Muñoz for helping me from the very beginning of this project.  Special thanks to the Geological Survey of Canada which through the MAP:GAC Project helped financing part of the field work, and contributed with the Cosmogenic Nuclide dating.  Many thanks to my friend Cristián Cáceres from the Rock Mechanics Laboratory of the Mining Department for his assistance and supervision in the preparation and execution of the UCS and direct shear rock lab tests. Dr. Laxmi Chikatamarla from the Rockphysics Laboratory of the Earth and Ocean Sciences Department also kindly contributed with equipment for the triaxial tests.  Finally, this work would not have been possible without the incredible love, patience and support of my husband, Nicolás, and my children, Diego and Martina, and a heartfelt thank- you to my friends for their unconditional love and support.      xvi DEDICATION               ………..To Peter Welkner, my lovely dad who art in Heaven………..                1 1. INTRODUCTION 1.1. Problem Statement Analysis of massive rock slope failures and subsequent motion of rapid rock avalanches is a challenging task given the difficulty in determining the processes and circumstances under which failures occurred. This task is even more complex when dealing with prehistoric rock avalanches due to the need for paleo-landscape reconstruction, paleo-climate, and knowledge of the pre-failure conditions. In most cases, only approximations of reality can be provided based on judgment and experience due to geologic uncertainty and the inherently variable nature of rock.  In the rugged mountains of the Andean Cordillera of central Chile, close to the Chile- Argentina border, a number of large lobate-shaped diamicton deposits were erroneously identified as moraines owing to their geomorphic similarity with rock avalanche deposits, and the detrital characteristics of both types of material. Detailed mapping, including textural, compositional, and stratigraphic relations to the geomorphic setting, yielded no clear evidence of moraines in the study area located in the valley of Portillo, Chile. This condition was also supported by cosmogenic nuclide dating by 36Cl which revealed that the deposits were post-glacial in age. Instead, the Portillo deposits have been interpreted as belonging to two significant rock slope failure events: the Upper Pleistocene Portillo Rock Avalanche (PRA) initiated from the Caracoles Range in the east side of the valley, and a Holocene rock slump and rockslide initiated from the west and east side of the valley, respectively. This new interpretation changed the assessment of the hazard posed by future rock slope failure events at the site from negligible (glacial deposits assumed) to possible (landslide deposits present).  Research was undertaken through this thesis to investigate the likely failure mechanism, volume and runout path characteristics of the Portillo Rock Avalanche. The insights gained through this back-analysis were then used to assess the hazard potential of a reccurring major rockslide posed to the Portillo site, which includes both a key transportation route between Chile and Argentina (International Santiago – Mendoza Corridor) and a popular ski resort (Figure 1.1).  2  Figure 1.1. Photos of the Portillo Rock Avalanche deposits which underlie the Santiago-Mendoza International Corridor and the Portillo Ski Resort.     3 1.2. Research Objectives The purpose of this research is to investigate the likely failure mechanism of the Portillo Rock Avalanche, to characterize the runout of the debris in terms of path, volume constraints, and rheology, and to use these results to calibrate and constrain forward analyses focused on the likelihood of future catastrophic rockslide events. The specific objectives of the research are described below. 1.2.1. Mechanisms of failure    Investigate the importance of the internal structure of the key discontinuity sets and the mechanical properties of the rock mass with respect to the pre-historic Portillo Rock Avalanche.    Develop numerical models which will shed light on the rockslide kinematics.    Investigate the nature of the toe-release of the slide mass. 1.2.2. Runout deposits    Determine a pre-failure topography reconstruction and volume estimates for the landslide deposits.    Develop three-dimensional runout models (DAN3D) to characterize the travel path, volume constraints and rheological behaviour of the different prehistoric rock slope failure events mapped at the study site. 1.2.3. Forward modeling    Develop numerical models to assess the stability of the present-day Caracoles Range at Portillo.    Determine the effects of the application of a varying water table in the stability of the Caracoles Range.    Determine the effects of the application of a seismic load to the stability of the Caracoles Range (i.e. sensitivity of the slope to an earthquake trigger).  4    Determine extent and reach of rockslide debris in event of worst case scenario of catastrophic rock slope hazard. 1.3. Thesis Structure The problem statement and the research objectives are presented in Chapter 1, together with a critical review of the literature related to the research. General topics such as slope failure initiation and runout prediction of failed slopes are covered.  Chapter 2 introduces the different methodologies used in this thesis. The Total Slope Analysis procedure is presented and defined, and the numerical modeling procedures used for this are discussed in detail. The methodologies used for data collection (e.g. terrestrial laser scanning) are explained. The computer-based techniques used to process and analyze the data are described.  The engineering geology of the study area is presented in Chapter 3. It includes the tectonics, geological setting and seismicity of the site together with the chronology of the prehistoric events. A geotechnical characterization of the Portillo Rock Avalanche site is presented, which includes a description of the discontinuity network, laboratory testing results and scaled rock mass properties and hydrogeology.  Chapter 4 describes in detail the failure mechanism of the Portillo Rock Avalanche, starting with the development of a numerical model, followed by the base model results, parametric analyses, and interpretation of results.  The runout analyses of the prehistoric rock slope failure events mapped in the study area are presented in Chapter 5, where the volume calculations are explained and the pre-failure topography reconstruction is discussed. Both the Portillo Rock Avalanche and the Holocene rock slump model results are presented.  Chapter 6 presents the forward analysis of a potential repeat large-scale rockslide event at Portillo, originating from the same slope as the prehistoric Portillo Rock Avalanche. This commences with a study of the potential failure initiation and volume assessment, followed by a runout analysis. This chapter concludes with a hazard assessment summary.  5 A general discussion of the most important results and conclusions of the research is presented in Chapter 7. Limitations and uncertainties involved in the development of the research are reviewed. Finally, recommendations for further research are presented. 1.4. Literature Review 1.4.1. Rock mass classification systems Rock mass classification systems provide a means to characterize, qualitatively and quantitatively, the rock mass and discontinuities in order to estimate their strength and deformation properties (Cai et al. 2004), so as to predict its behaviour under different loading conditions. The most widely used systems are the Rock Mass Rating system (Bieniawski 1973), the Norwegian Geotechnical Institute’s Q system (Barton et al. 1974), and the Geological Strength Index (Hoek 1994; Hoek et al. 1995; Hoek & Brown 1997; Hoek et al. 1998). The RMR and Q systems are both build on the Rock Quality Designation (RQD; Deere 1964), and give a general idea of the quality of the rock for engineering purposes. Since the inception of the RMR, several updates have been forwarded corresponding to empirical calibrations to observed tunnels or mine behaviour. The most popular are RMR76 (Bieniawski 1976) and RMR89 (Bieniawski 1989) which primarily differ with respect to the joint condition description and ratings. The RMR employs five input parameters, each of which is given a rating of importance for a particular situation. The total rating ranges from less than 25% (worst conditions) to 100% (best rock conditions) and indicates the quality of the rock as shown in Figure 1.2. The parameters are:     The strength of the intact rock    Rock Quality Designation (RQD)    Spacing of fractures    Condition of fractures (influencing their shear strength)    Ground water conditions  The value derived may then be subsequently adjusted to account for adverse discontinuity dip orientations with respect to the engineering objective (Table B in Figure 1.2).   6  Figure 1.2. Rock Mass Rating classification system (modified from Bieniawski 1989).  It is widely recommended that when the RMR rating system is to be used to estimate the rock mass properties, only the first four terms listed above should be used in the calculations (referred to as RMR’76 or RMR’89); the fifth term is set to that for a completely dry rock mass. This allows for the groundwater conditions, as well as corrections for unfavourable structural orientation parameters, to be treated explicitly, for example through an effective stress numerical analysis where their significance can be better dealt with (Marinos et al. 2007).  7 Because the first four parameters refer only to the characteristics of the rock mass itself, they can be used to estimate the Geological Strength Index (GSI) value by using either:  RMR’76 = GSI RMR’89 = GSI -5  [1.1]  For Equation 1.1, the procedure followed must consider setting the Adjustment for Joint Orientation = 0 (very favourable) and Groundwater ratingRMR76 = 10 or Groundwater ratingRMR89 = 15, corresponding to dry conditions (Hoek & Brown 1997; Hoek et al. 1998).  These relationships hold if the conditions correspond to good quality rock (RMR > 25). Where the RMR < 25, it is recommend that the Q-system be used (described in the next paragraph) with the following relationship (Equation 1.2):  RMR = 9*LnQ + 44  [1.2]  The Q or NGI (Norwegian Geotechnical Institute) classification system (Barton et al. 1974) is based on approximately 200 tunnel case histories from Scandinavia (whereas the RMR was based on tunnelling experience in South Africa). It is essentially a weighting process in which positive and negative aspects of the rock mass are assessed. The quality of the rock, Q, is expressed through six independent parameters arranged in three quotients (Equation 1.3), and the rating varies on a logarithmic scale from 0.001 (exceptionally poor) to 1,000 (exceptionally good).  r w n a J JRQDQ x x J J SRF =   [1.3]  where,     RQD: rock quality designation    Jn: number of joint sets    Jr: roughness of the most unfavourable joint set    Ja: alteration of the joint surfaces    Jw: water inflow  8    SRF: stress reduction factor  Like RMR89 the last two of these parameters can be dropped to give Q’ and the effects of groundwater and stress can be treated separately numerically. Tzamos & Sofianos (2007) suggested that RQD/Jn represents the block size, Jr/Ja reflects the inter block shear strength (friction angle), and Jw/SRF represents the effective stress conditions.  The Geological Strength Index (GSI) is a rating system which estimates the reduction of the rock mass strength for different geological conditions, and is based on a geologic description of block structure and joint surface conditions within the rock mass (Hoek et al. 1998). One of its advantages is based on its geological reasoning which allows for the covering of a wide range of rock masses and conditions. The chart in Figure 1.3 illustrates how the geologic descriptions are used in calculating the GSI. The most common critique of the GSI system is that its use requires extensive engineering experience and judgment. For this reason, attempts has been made to better quantify the GSI (Sonmez & Ulusay 1999; Cai et al. 2004), the latter using the block volume and joint condition factor as quantitative characterization factors (Cai et al. 2004). However, Marinos et al. (2007) point out that care must be taken when quantifying the GSI, especially when dealing with rock masses whose structural fabric has been highly disturbed (GSI < 35) because in those cases the quantitative characterization factors are difficult to measure.   9  Figure 1.3. General chart for GSI estimates based on geological observations (modified from Marinos et al. 2005).     10 1.4.1.1. Hoek-Brown failure criteria Hoek & Brown 1980 proposed an empirical failure criterion for estimating the strength of jointed rock masses in some specific cases. Later this criterion was extended for most types of rock (Hoek & Brown 1997; Hoek et al. 1998; Marinos & Hoek 2001), and only the latest version (Hoek et al. 2002) estimates static rock mass strength taking into account the specific stress conditions of slopes and disturbance induced to rock masses due to rock blasting (D factor; see Equations 1.5 & 1.6). The Generalised Hoek-Brown Failure Criterion for jointed rock masses is defined by (Hoek & Brown 1997):  a ci bci sm ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ++= σ σσσσ ''' 331   [1.4]  Where σ1’ and σ3’ are the maximum and minimum effective stress at failure, mb is the reduced value for the Hoek-Brown material constant mi, s and a are empirical constants for the rock mass and σci is the uniaxial compressive strength of the intact rock. This criterion does not consider the intermediate effective stress, σ2 but is still applicable as a first- approximation for two-dimensional cases.  To apply the Hoek-Brown failure criteria and thus obtain the strength and deformability of jointed rock masses, three properties of the rock mass must be estimated (Cai et al. 2004):     The uniaxial compressive strength (σci) of the intact rock: obtained from laboratory triaxial or uniaxial compressive tests.    The Hoek-Brown constant mi for the intact rock: obtained from laboratory triaxial compressive tests.    The GSI for the rock mass: field survey.  When laboratory testing is not performed, mi and σci may be obtained from published tables (Hoek et al. 1995; Hoek & Brown 1997). When the GSI is known, the parameters are given as (Hoek et al. 2002):   11 100exp , 28 14b i GSIm m D −⎛ ⎞= ⎜ ⎟−⎝ ⎠   [1.5]  100exp , 9 3 GSIs D −⎛ ⎞= ⎜ ⎟−⎝ ⎠   [1.6]  20 15 310.5 , 6 GSI a e e − −⎛ ⎞= + −⎜ ⎟⎜ ⎟⎝ ⎠   [1.7]  The parameter m is related to the particle interlocking representing the cohesion of the rock (Hudson & Harrison 2005). Thus, brittle rocks will be represented by a large mi value (greater than 15) and 0.5 ≤ a ≤ 0.6, while more ductile rocks are characterised by lower values. The s parameter relates to the degree of fracturing; it is a representation of the internal friction angle of rocks (Hudson & Harrison 2005). Low values of s (close to 0) correspond with heavily jointed rock masses (Kumar 1998). D is a factor that represents the degree of disturbance of the rock mass when blasting has occurred, and varies from 0 (undisturbed rock masses) to 1 (very disturbed rock masses).  The uniaxial compressive strength for the rock mass (σc) can be obtained by setting the minor principal stress σ3 = 0 in Equation 1.4, resulting in:  a c cisσ σ=   [1.8]  and the tensile strength (σt) can be obtained by setting σ1’=σ3’=σt in Equation 1.4, resulting in:  ci t b s m σσ = −   [1.9]   12 1.4.1.2. Mohr-Coulomb failure criteria In the Mohr-Coulomb failure criterion, the rock mass shear strength is defined by the cohesive strength c and the friction angle φ. Failure will take place on any plane when the shear stress acting across that plane reaches critical shear strength (Figure 1.4). This criterion can be expressed by the following equations:  φστ tannc +=   [1.10]  31 sin1 sin1 sin1 cos2 σφ φ φ φσ − ++−= c   [1.11]  Equation 1.11 corresponds with the representation of the Mohr-Coulomb failure criterion in terms of the minor and major principal stresses, so it is used for two-dimensional stress analysis. Since this criterion was developed for compressive stresses, a tensile cut-off is usually utilized to give a more realistic value for the uniaxial tensile strength (Hudson & Harrison 2005).   Figure 1.4. The Mohr-Coulomb failure criterion (modified from Hudson & Harrison 2005).  13 Equivalent Mohr-Coulomb strength parameters can be derived for a range of minor principal stress values. The procedure consists on simulating a set of full-scale triaxial tests using the Hoek-Brown envelope given by Equation 1.4 and its results are fitted to an average linear relationship given by Equation 1.11 (Hoek & Brown 1997). The adjustment must be done because the non-linearity of the H-B criterion results in different values of M-C parameters as normal stress varies (Figure 1.5). The values of c and φ are then determined from Equations 1.12 and 1.13.  1sin 1 k k φ −= +   [1.12]  2 cmc k σ=   [1.13]  where k is the slope of the line given by Equation 1.11 and σcm is the unconfined compressive strength of the rock mass.   Figure 1.5. Approximation of M-C failure envelope from the H-B envelope. The Mohr-Coulomb strength parameters are obtained from simulated triaxial values derived from the Hoek-Brown failure criterion. Different Mohr-Coulomb envelopes are determined for different levels of σ3 (modified from Cai et al. 2004).  14 Hoek et al. (2002) indicates that in some cases it is important to consider the overall behaviour of the rock mass rather than the detailed fracture propagation process, and introduces the concept of a ‘global rock mass strength’ which could be estimated from the M- C relationship as follows:  φ φσ sin1 cos2 −= c cm   [1.14]  A key issue is to determine the proper σ3max value in order to correctly obtain the M-C strength parameters. Hoek et al. (2002) investigates the σ3 range for tunnels and slopes and their results show that for the case of shallow tunnels the maximum value is given by Equation 1.15, and for the case of slopes the maximum value is given by Equation 1.16:  0.94 3max' '0.47 ' cm cm H σ σ σ γ −⎛ ⎞= ⎜ ⎟⎝ ⎠   [1.15]  0.91 3max' '0.72 ' cm cm H σ σ σ γ −⎛ ⎞= ⎜ ⎟⎝ ⎠   [1.16]  1.4.2. Slope failure initiation A significant amount of work has been published on the failure of large rock slopes. This section presents the up-to-date available work on the mechanisms of failure of large slopes in jointed rock masses, and discusses the most suitable tools to approach the problem. 1.4.2.1. Mechanisms Rock slope failure initiation mechanisms that involve sliding movements along discrete planes with release along persistent joints are traditionally analyzed through kinematic and limit equilibrium techniques (Stead et al. 2004; 2006). These techniques assume translation of a rigid block along one or more surfaces by means of a static analysis (balance between resisting and driving forces), with minimal internal rock mass deformation.  15 In massive natural slopes, where the discontinuity network is not fully persistent, failure may occur through a combination of sliding along discrete planes and brittle fracture propagation through intact rock where extensive internal rock mass deformation occurs (Terzaghi 1962; Einstein et al. 1983; Eberhardt et al. 2004). In this case, instead of dealing with a planar failure, the type of slope failure determined by the shear surface may correspond to a bi- planar or step-path failure mode. As such, the limit equilibrium approach is not suitable for the analysis, owing to its inability to account for slope deformation and intact rock fracturing (Eberhardt et al. 2004; Stead et al. 2004; 2006).  In the case of bi-planar failure, which is the mode of failure applicable to this study, Mencl (1966) and later Kvapil & Clews (1979 ) explained that failure along the bi-planar slide surfaces would require the development of a damage or transition zone between the active upper slope and the resistive lower part of the slope (Figure 1.6). The transition zone is a modified solution based on the Prandtl´s prism theory initially elaborated for problems of foundation construction. According to Kvapil & Clews (1979 ) the transition zone is the area where fracture propagation through the rock continuum occurs (Figure 1.7) and is further characterized by:     Change in direction of the stresses.    Very large transversal deformation (heaving).    Very large degree of fracturing.   Figure 1.6. Illustration of the transition zone utilizing the Prandtl´s prism (modified from Kvapil & Clews 1979 ).  16  Figure 1.7. a) Distribution of transversal deformation with respect of the position of the Prandtl transition zone; b) Distribution of the degree of fracturing with respect of the position of the Prandtl transition zone (modified from Kvapil & Clews 1979 ).  Another mechanism that may affect steep and high dip-slopes in hard rock is buckling. This type of failure mode is closely associated with the internal structure of the rock mass and the geometry of the slope and generally occurs as buckling of the slab near the toe of the slope followed by translational failure above the buckle (Cavers 1981). In the literature available regarding buckling mechanisms (Cavers 1981; Nilsen 1987; Li & Zhang 1990; Froldi & Lunardi 1995; Stead & Eberhardt 1997), it is of general agreement that some geostructural and geomechanical conditions must be achieved, namely:     Presence of rocks with a slab-shaped structure in which discontinuities isolate masses with one dimension clearly smaller than the other two.    Inclination of the main discontinuity is generally higher than their angle of friction.    The surface of key persistent structures (e.g. bedding, foliation, etc.) dips sub- parallel to the surface of the slope (i.e. dip-slope).  In a broad sense, most of the case studies present in the literature (Cavers 1981; Nilsen 1987; Li & Zhang 1990; Eberhardt & Stead 1998) occur in thinly bedded and weak sedimentary lithologies such as interbedded shales and sandstones, or their metamorphic equivalents. Those that have occurred in massive crystalline rocks are related to the development of exfoliation joints (sheeting).   17 Many of the authors cited above have analyzed their cases using force and momentum balance calculations (i.e. Limit Equilibrium approach), basing their solutions on mostly simplified assumptions such as rigid individual blocks with no contact force between the blocks and slope (Cavers 1981). As Stead & Eberhardt (1997) and Eberhardt & Stead (1998) show in their work, numerical modeling through a discontinuum approach provides a better means to assess whether buckling is a possible mechanism or not. 1.4.2.2. Numerical Modeling Several authors (Coggan et al. 1998; Kimber et al. 1998; Costa et al. 1999; Eberhardt et al. 2004; Stead et al. 2004; 2006) conclude that numerical modeling techniques (continuum and discontinuum codes) are a suitable tool to conduct analyses where internal slope deformation occurs because they can best reproduce complex slope failure mechanisms.  The numerical method that is frequently used for modeling failures in natural jointed systems is based on a discontinuum approach where the problem domain corresponds to an array of distinct interactive blocks that can be subjected to external loads, and large strain displacements and block rotations can be accommodated and modelled (Pritchard & Savigny 1990; Eberhardt et al. 2004). The commercial distinct-element code that will be used in this thesis is UDEC (Itasca 2004).  Starfield & Cundall (1988) wisely advise that the ability to develop successful models is in close relation with the capacity of extracting only the geological features that will define the boundary conditions for different types of rock-failure mechanisms. In general, key factors that help bring understanding to the mechanisms that control failures in natural jointed rock slope systems are:     Reasonable representation of the discontinuity geometry, especially persistence and spacing (Terzaghi 1962; Kimber et al. 1998; Eberhardt et al. 2004; Stead et al. 2006).    Reliable estimates of the strength of jointed rock masses based on the internal array of the discontinuities (Singh & Rao 2005) and strength properties of the joints and intact rock (Stead et al. 2006).   18 The most used approach to dynamic analysis in geotechnical design is the pseudo-static analysis where effects of a seismic loading are represented by constant horizontal or vertical accelerations (Wyllie & Mah 2004). The results of this type of analysis are highly dependant on the selected seismic coefficient k, but when the difficulty in assigning the proper k is considered, the reliability on the pseudo-static analysis starts to decrease. Nevertheless, the advance of numerical models have provided an alternative to the dynamic analyses because they allow permanent slope deformations resulting from seismic loading (Wyllie & Mah 2004). 1.4.2.3. Triggering factors The two most common trigger mechanisms for landslides are heavy precipitation, which can induce changes in effective stress that may overcome the material strength due to build up of pore water pressures, and earthquakes which may cause a slope failure due to the extra load induced by ground shaking during a determined period of time.  When a soil or rock slope is affected by dynamic loading, the material strength, slope configuration, pore water pressure, and number of cycles of the ground motion (or duration of the event) are the key aspects which will determine if a slope fails during an earthquake (Keefer 1984; Bommer & Martinez-Pereira 1999). In general it can be said that the number of landslides caused by earthquakes increases with increasing magnitude, however, local geologic conditions and seismic parameters such as depth of the rupture zone and distance of the site from the source may also constrain the number of triggered landslides (Keefer 1984; Rodriguez et al. 1999).  A study from Keefer 1984 reveals that 27 of the rock avalanches reported originated from intensely fractured slopes (centimetric to metric block size), presenting also at least one of the following signs of potential instability:     Distinctive planes of weakness, such as bedding planes, master joints, or fault planes, that dip out of the slope,    Significant weathering of the rock mass, and    Geologic or historical evidence of previous activity.   19 The risk presented by precipitation or earthquake-induced landslides depends on the population distribution and constructed works relative to potential landslide sources. Moreover, a large earthquake in a mountainous region can trigger multiple landslides in a short period of time having an overall impact on human life and economic infrastructure far beyond that of any single landslide (Keefer 1984). 1.4.3.  Runout prediction of failed slopes Efforts to better comprehend and predict catastrophic slope failures have steadily increased since the early works of Terzaghi (1950), mainly because of the increase in frequency of impact of such natural phenomena due to the rapid growth of population and infrastructure and development into marginal areas (Hungr et al. 2005). In step, the runout analysis must fit into the framework of quantitative landslide risk assessment providing the spatial distribution of landslide intensity (Hungr et al. 2005).  Moreover, considering the need to conduct hazard assessments, the use of a common terminology is an essential starting point. As such, the typological classification described by Hungr et al. (2001) and the terminology used by Cruden & Varnes (1996) will be used in this thesis. 1.4.3.1. Empirical methods Empirical methods involve the collection of historical data regarding the characteristic parameters of landslides (volume, fall height, distance traveled) and the control exerted by topography (path). The most common approach is the fahrböschung (Heim 1932) or reach angle (Corominas 1996), which is defined as the inclination of a projected line that connects the crest of the scarp with the distal margin of the displaced mass (Figure 1.8). Heim (1932) observed and described an inverse relation between the fahrböschung and the volume of the deposit; therefore the angle of reach is widely used as a measure of the relative mobility of rock avalanches, and can be interpreted as an expression of the efficiency of energy dissipation (McDougall 2006).   20  Figure 1.8. Diagram of the fahrböschung (or angle of reach), where H is the vertical drop and L is the horizontal projection of the distal margin of the displaced mass.  Several linear regression equations have been proposed which correlate the volume of landslides and the distance traveled by the debris (Scheidegger 1973; Hsu 1975; Davies 1982; Li 1983; Corominas 1996), but the results show a large amount of scatter in the data indicating that the relationship not only depends on the volume of the landslide but on both the type of movement and topography of the runout path (Corominas 1996). Likewise, the mode of failure (e.g. translational vs. rotational) is also an important controlling factor that would contribute to this scatter.  Similar to the volume-fahrböschung relationships, correlation between the volume and the aerial exposure of the deposits, rather than its distal extent, have been proposed by Li (1983) and Iverson et al. (1988).  The difficulty in obtaining a complete dataset of descriptions of the initial volume, conditions of the slope, topography of the path, etc., means that reliable relationships are hard to achieve (Crosta et al. 2003). Consequently, the next step is numerical modeling, which has the ability to better simulate the complex dynamic behaviour of landslide runout.   21 1.4.3.2. Numerical methods In contrast to empirical methods, numerical methods are based on solid and fluid dynamics with an aim to model moving landslides, which often assume complex behaviour and show a continuum passage from sliding to flowing (Hungr et al. 2005). Furthermore, landslide runout have the ability to entrain and deposit material while they are moving resulting in important changes in mass and volume (Crosta et al. 2003). Therefore, the most common solutions are based on continuum methods because of their capacity to accommodate internal deformation.  The dynamic model that will be used in this thesis bases its interpolation on the smooth particle hydrodynamic theory that allows unlimited distortion as well as separation and joining of flowing masses (McDougall & Hungr 2004). The resulting code, DAN3D (McDougall & Hungr 2004), is a three-dimensional dynamic model based on DAN, the two-dimensional “Equivalent Fluid Method” of Hungr (1995).  DAN3D includes methods that account for specific features of landslide such as internal strength, material entrainment and rheology variations. This continuum dynamic model is governed by internal and basal rheological relationships. The internal rheology is assumed to be frictional and is ruled by the internal friction angle (φi). On the other hand, various basal rheologies can be implemented including: laminar, turbulent, plastic, Bingham, frictional, and Voellmy. The selection of a basal rheology together with its associated parameters is done based on an empirical calibration procedure where the occurring landslide is subject to a trial-and-error back-analysis. The calibrated parameters are considered apparent rather than actual material properties (McDougall 2006).  22 2. METHODOLOGY 2.1. Total Slope Analysis Procedure The methodology applied to this study attempts to use field-based tools such as geological and geotechnical mapping and sampling to provide input data for the development of numerical models that can best capture the likely failure mechanisms and stability state of the Caracoles Range. The distinct element code UDEC (Itasca 2004) was used to back analyze the Portillo Rock Avalanche and assess the potential instability of the current slope (Figure 2.1). The insights gained on the behaviour of the slope combined with state-of-the-art dating techniques for the mass wasting events were used to calibrate the rheological behaviour and distribution of the debris. The post-failure motion was analyzed with the dynamic/rheological code DAN3D (McDougall & Hungr 2004). In view of the importance on the assessment for future rock mass failures in the area, a ‘Total Slope Analysis’ approach (Stead et al. 2006) was taken to link the properties, kinematics and mechanisms important to understanding both failure initiation and runout (Figure 2.2).   Figure 2.1. Diagram of a total slope failure analysis showing the codes that will be used in this thesis for each specific area (modified from Stead et al., 2006).  23  Figure 2.2. Flow chart of the procedure followed by the Total Slope Back Analysis and Total Slope Forward Analysis (modified from Strouth 2006).        24 2.2. Field Work and Data Collection The field work carried out included 15 days of mapping of the rockslide source area and runout deposits, and three days of terrestrial laser scanning (LiDAR) to obtain data on the spatial characteristics of the sliding surface and associated structures, key input parameters for the numerical models. The objectives of this work were two fold:     Examine the nature of the runout deposits and possible slide kinematics.    Collect data to constrain the numerical models developed to answer questions regarding rockslide kinematics, triggering mechanism and potential for future failure events. 2.2.1. Field mapping The mapping and characterisation of the eastern wall was done with the aid of several pictures taken of the scarp from diverse angles of sight and distances (Figure 2.3). This was necessary because the steepness of the wall impeded a close survey of the initiation zone. During the visual inspection of the eastern slope, key geometrical characteristics were extracted in order to build the numerical models, and two domains were identified based on distinct orientations of the discontinuities and differential damage to the rock mass generated by a thrust fault that runs along the toe of the Caracoles Range (see section 3.5.1).   25  Figure 2.3. Pictures obtained from the eastern wall which shows different views of the scarp at different distances.  Mapping of the slide deposits was carried out using aerial photographs at 1:80,000 and 1:20,000 (Figure 2.4), taking special care of the stratigraphic relationships between the different lobes, and the spatial relation between distinct lithologies present in the lobes and their potential source area. The mapping was conducted to help answer key issues for the future dynamic analysis, such as whether the deposits corresponded to a single event or were generated over a period of time due to multiple events, and if the deposits originate from a single source or different sources participated in the process.  26  Figure 2.4. Aerial photographs at 1:20,000 (top) and 1:80,000 (bottom) scales which were used to map the debris in the runout zone.  For the bedrock geology and structures a 1:250,000 geological map was used as a base map (Rivano et al. 1993), with some modifications being introduced following updated tectonic interpretations (Cegarra & Ramos 1996; Godoy 2005).  27 2.2.2. Rock sampling Part of the field work involved the collection of rock samples, both intact rock and mated joint surfaces, which were tested in the UBC Mining Department Rock Mechanics Laboratory and the Rockphysics Laboratory of the UBC Earth & Ocean Sciences. A series of uniaxial, triaxial, and direct shear tests were performed and the results were later used to compute elastic modulus, Poisson’s ratio, and the strength of both the intact rock and the discontinuities in order to constrain the numerical models. Tests were carried out in accordance with methods recommended by Brown (1981).  The sampling strategy consisted of collecting bulk samples of ~0.01 m3 from the most representative lithologies present in both, the scarp and the rockslide debris. The chosen size of the blocks in general allowed for the exclusion of major discontinuities that could affect the performance of the intact rock during the triaxial and uniaxial tests. 2.2.3. Discontinuity mapping Key to understanding the kinematics of the failure was the information regarding large-scale structures and persistent discontinuities. However, due to the steepness of the slope walls, as well as rockfall hazard from above, accessibility issues made traditional methods of rock mass characterization and joint mapping restrictive. This problem was identified prior to the field visit and plans were made to carry out a detailed 3-D terrestrial laser scanning (LiDAR) campaign using the UBC/SFU Optech ILRIS-3D laser scanner (Figure 2.5). The joint orientation data was complemented with random measurements around the site (see section 3.5.1).   28  Figure 2.5. LiDAR (Optech ILRIS-3D) measurement stations.   29 2.3. Terrestrial Laser Scanning (LiDAR) 2.3.1. Data acquisition Five measurement stations were established with a number of scans produced ranging from 150 m to 1000 m distance to the slope face (Table 2.1 & Figure 2.6). Once processed and analyzed, these data sets were used to characterize the rock mass, derive rock mass properties, and provide spatial input and constraints for the following numerical models (distinct-element modeling).  Table 2.1. Terrestrial laser scan measurement stations and associated parameters. UTM N UTM E trend plunge2 1 59 21 failure scarp - N 600 fair 3 80 22 failure scarp - center 800 fair 4 111 19 failure scarp - S 1000 fair B 6,366,708 395,518 5 17 23 JS1 150 good 6 43 37 JS1 - JS2 - F1 300 - 400 good 7 43 37 JS1 - JS2 - F1 300 - 400 fair 8 43 37 F1 300 - 400 poor 9 43 37 JS1 & F1 300 - 400 fair 10 81 31 base of slope 400 poor 11 13 33 JS1 - F1 250 good 12 13 33 JS1 - F1 250 fair 13 105 32 JS1 600 poor 14 38 27 JS1 - JS2 - F1 500 fair (1) Estimated distance from the survay station to the center of the target (2) Direction of inclination towards SW (3) JS1: bedding planes; JS2: cross-joints; F1: main fracture E 6,366,608 395,559 6,367,003 395,527 6,366,500 395,300C D A 6,367,204 395,168 Target3 Relative Quality Position (m) Orientation (°)Survey Station Scan # Approx. distance (m)1    30  Figure 2.6. Approximate location and orientation of the survey stations.  Possible source of errors in the final dataset are directly related to the number of laser points that reflect off the discontinuity surface, which depend on the laser resolution, size of the discontinuity, distance from the target, and orientation of the discontinuity relative to the orientation of the laser scanner (Kemeny & Donovan 2005 ).  The better scans were the ones obtained from stations where the line of sight of the scanner was orthogonal to the target and the distance to the target was within the operating range of the instrument (< 800 m for this model, Sturzenegger et al. 2007). The rationale behind the relative quality classification of the scans was based on the combination of both possible sources of error, distance and orientation. The worst combination was that of far distance and oblique line of sight with respect of the target which resulted in few spot reflections thus lowering the resolution of the point cloud (Kemeny & Donovan 2005; Sturzenegger et al. 2007).  31 2.3.2. Data processing and analysis The data was processed using Split-FXTM software (Split_Engineering 2005). A 3-D point cloud of the scanned surface was generated for each station classified between fair and good; poor scans were discarded to avoid the uncertainty associated with them. From the point cloud, discontinuity orientation and spacing (when possible) was extracted (Figure 2.7). The sequence of steps followed during the data processing was the same for each scan and was similar to that suggested by Kemeny & Donovan (2005 ) and Strouth & Eberhardt (2005):     Point cloud orientation: input the trend and plunge of the scanner previously measured in the field so as to make the point cloud match up with the true spatial orientation.     Edit the point cloud: because the laser scanner is a line of sight device, the first object in the path of the light is recorded and thus the point cloud generally must be ‘cleaned’ before it is used (e.g. get rid of anomalous data points associated with talus, vegetation, etc.).     Mesh generation: the software creates a triangulated surface mesh which incorporates data smoothing by varying the mesh density (points per mesh grid cell). The mesh density recommended by the Split-FXTM software manual was in the range of 30 points per mesh grid cell, which proved to work fine for this study.     Patch generation: patches correspond to flat surfaces in the mesh that are automatically generated and interpreted as discontinuities found in the point cloud. The flat surfaces are found by calculating the normal to each triangle and then searching for more triangles that fulfill the flatness criteria corresponding to the maximum degrees that the neighbouring triangles may deviate in the patch. Another criterion corresponds with the size of the patch which refers to the minimum number of triangles that form the patch. For this study, the neighbour angle and the patch size ranges were set to 4º-8º and 10-40 grid cells per patch, respectively, as suggested by Split Engineering (2005).   32    Patch editing: depending on the density of the mesh, sometimes the neighbour angle or the patch size were too small generating few patches or noise in the stereonet. For that reason, a manual inspection and editing after the automatic patch generation was always performed.     Stereonet plotting: after the editing, the average orientations were plotted on a stereonet and the principal joint sets were extracted. Each patch was plotted as one point in the stereonet, and its size was adjusted relative to the patch area giving an idea of the relative importance of the joint set.   Figure 2.7. Sequence of steps followed during the data processing consisting on scanning of the outcrop, point cloud generation and stereonet analysis.  For some of the scans (i.e. 5, 6, 7 & 14) it was also possible to estimate the spacing of the discontinuities, following the methodology presented by Strouth & Eberhardt (2005) and Strouth (2006). This proved valuable for building the numerical models used to investigate the likely trigger and failure mechanisms for the main rockslide and the potential for future failures.         33 2.4. Computer-Based Techniques Several computer techniques were utilised during this thesis which can be divided into two main groups.     Data processing tools involving software used to process, validate, and interpret the field and laboratory data collected.    Numerical modelling software tools, which were used for the stability and initiation analysis (UDEC) and for the runout analysis (DAN3D).  Both the data processing and numerical modeling tools are described in more detail in the following sections. 2.4.1. Data processing 2.4.1.1. Spheristat 2.2 After identifying the principal joint sets by means of field mapping and terrestrial laser scanning, the sets were exported from Split-FXTM and analyzed using Spheristat 2.2 created by Pangea Scientific (Stesky & Pearce 2007). A single stereonet was produced for each domain and its point density distribution, which included all the joint sets of the domain, was calculated using the Schmidt (1%) counting method. The average orientation obtained from each joint set was then considered for subsequent UDEC numerical models. 2.4.1.2. RocLab 1.0 The data was processed in RocLab 1.0 (Rocscience 2002), a programme used to determine rock mass strength parameters based on the Hoek-Brown constitutive criterion (mb, s and a) and equivalent Mohr-Coulomb strength parameters (c, φ). For each domain (see section 3.5.1 for explanation) a lower and upper bound value of σci and GSI was entered and the corresponding Hoek-Brown and equivalent Mohr-Coulomb strength parameters were obtained. The final range of scaled rock mass strength parameters was used as input for the numerical models.  34 2.4.1.3. Surfer 8.0 Surfer 8.0 (Golden_Software 2002) is a grid-based software which interpolates irregularly spaced XYZ data into a regularly spaced grid. The grids are then used to produce contour, vector and surface maps among other applications. For this thesis, Surfer 8.0 was used to create the pre-failure topography and to estimate volumes of the detached masses. The pre- failure topography was created using the current topography and modifying its regularly spaced grid. The resultant grid file corresponded to a surface that simulated a deeper valley (no rock avalanche debris), with no lobes and lake. Whenever a pre-failure topography was available (path or initiation zone), its grid was subtracted from the current one and the failed (or exceeding) volumes were calculated. The values obtained were always compared to the volumes calculated by other methods (see section 5.2) and if more than 20% error existed between the methods, the grid was adjusted until an error less than 20% was obtained. 2.4.1.4. ArcGIS 9.0 The spatial distribution of the data was analysed in a GIS platform ArcGIS 9.0 (ESRI 2004). The input data included the geology and topography which was used to produce the digital geological map (see Figure 3.6), locate the cosmogenic nuclide and lab test samples (see Figure 3.14) and to generate a DTM which gave information of the form of the slopes such as elevation, shape, dimensions, angle and aspect (see Figure 5.1). The EZ profile extension was used to create the cross-sections later used to help with the pre-failure topography reconstruction and the 2-D numerical models. 2.4.2. Distinct-element analysis (UDEC) In the distinct element code UDEC, the rock mass is represented as an assemblage of discrete interacting blocks which are subdivided into a deformable finite-difference mesh that follows linear or non-linear stress-strain laws (Itasca 2004). These interacting blocks are subjected to external loads, which may then undergo significant motion with time (Itasca 2004).   35 The discontinuous nature of the rock mass in the eastern slope at Portillo together with the most likely mode of failure hypothesized for the prehistoric rockslide made UDEC the best technique for developing the numerical models relative to other numerical tools (e.g. finite element using an equivalent continuum approach). 2.4.3. Runout analysis (DAN3D) DAN3D (McDougall & Hungr 2004) provides the user with a number of different basal rheologies that can be separately or combined. For this particular study two rheologies were selected and their combinations tested and ranked based on their ability to reproduce the mapped distribution of the debris observed during the field investigation. The results were further tested against temporal constraints corresponding with the sequence of events established by cosmogenic nuclide dating.  The runout analysis performed required two key inputs in addition to the material properties, a volume estimate and a pre-failure topography. The main constraints on model parameters were the actual distribution of the debris in the valley and the results provided by UDEC in terms of likely volume and position of failed mass. 2.5. Numerical Modeling Procedure 2.5.1. Static (non-inertial) solution The procedure followed for the development of the numerical models was obtained from the UDEC user’s manual (Itasca 2004). The main components that need to be specified in a problem so as to set-up and run a numerical model:     Problem geometry;    Constitutive behaviour and material properties; and    Boundary and initial conditions  The constitutive behaviour and material properties define the type of response the model have once disturbed, while the boundary and initial conditions simulate the in-situ state  36 before disturbance is applied. Once the components are defined, an alteration is made and the resulting response is calculated. An example of this procedure is illustrated in Figure 2.8.   Figure 2.8. Solution procedure for static analysis (modified from Itasca 2004).  37 A static solution is reached in UDEC when the rate of change of kinetic energy in a model approaches a negligible value by means of damping the equations of motion. When a model is under a static solution stage it means that it is either at a state of force equilibrium or at a state of steady flow of material (Itasca 2004). In the procedure of a static (non-inertial) analysis, UDEC uses a mechanical damping algorithm known as local damping to reach a force equilibrium state under the initial and boundary conditions applied. With local damping, the amount of damping on a node is proportional to the magnitude of the unbalanced force. Another form of damping used in the static solutions is the so called adaptive global damping, in which the damping constant is adjusted automatically in direct proportion to the rate of change of kinetic energy. Both forms give similar results when applied to static analyses (Itasca 2004). 2.5.2. Coupled hydro-mechanical analysis – steady state flow To assess the potential destabilizing effect of increased pore pressure developing within the discontinuity network in the slope, a coupled hydro-mechanical analysis can be performed using UDEC’s steady-state flow algorithm developed by Itasca (2004). For this type of analysis the basics are:     Fluid only flows through the fractures as the blocks are assumed to be impermeable. This results in total stresses being calculated inside the impermeable blocks, and effective normal stresses along the block contacts;    The coupling includes mechanical deformation occurring in the form of normal joint displacements (i.e. joint closure or opening) affecting the joint aperture and joint hydraulic conductivity, and conversely, joint water pressures affect the mechanical computations (Olsson & Barton 2001) (Figure 2.9).   38  Figure 2.9. Diagram showing the relationship between fracture conductivity, mechanical deformation, and water pressures in a coupled hydro-mechanical analysis (modified from Zangerl 2003).  Fluid flow through fractures in UDEC is assumed (idealized) as being a laminar flow between two smooth parallel plates with a mean velocity given by Equation 2.1:  fk iυ = ⋅   [2.1]  where i is the hydraulic gradient, and kf the fracture hydraulic conductivity which can be calculated using Equation 2.2 as follows:  2 12f a gk ν=   [2.2]  where,  a  is the fracture width; ν is the kinematic viscosity of the fluid; and g is the acceleration of gravity.    39 The flow rate per unit width can be represented by Equation 2.3, which is also known as the “cubic flow law”:  3 12 a gQ iν= ⋅   [2.3]  In reality, joints are rough and their aperture varies, limiting in part the validity in the application of Equation 2.3 as representative of flow behaviour in natural jointed rock systems. Mechanical joint aperture is defined as the average perpendicular distance between two rock joint surfaces, while hydraulic aperture is a theoretical smooth wall aperture measured by analysis of the fluid flow (Olsson & Barton 2001).  In the analysis carried out in this study (see forward analysis in section 6.1.2), a steady-state condition for the flow will be applied thus several simplifications are introduced in the algorithm. One of the most important is the elimination of the influence of the fluid stiffness in the mechanical timestep resulting in a much more efficient algorithm, and with the approaching of the steady-state condition, the pressure variation in each fluid step becomes smaller which permits the implementation of several fluid steps for each mechanical step without loss of accuracy (Itasca 2004). 2.5.3. Dynamic solution Dynamic solutions are required for problems involving high frequency and short duration loads such as seismic loadings. The calculations used for a dynamic analysis in UDEC solve the complete equations of motion, including inertial terms, using the real mass of the rigid blocks rather than scaled inertial masses as in the static solution. The generation and dissipation of kinetic energy is directly affected by the solution (Itasca 2004). A static equilibrium calculation always precedes a dynamic analysis. Similar to the case of the static solution, there are three main aspects to be considered when performing a dynamic analysis:     Dynamic loading and boundary conditions    Mechanical damping (Raleigh)    Wave transmission through the model   40 In UDEC, the external dynamic loading is applied as a stress (force) or as a velocity either at the model boundary or to internal blocks. Wave propagation at model boundaries are reduced by setting either quiet or free-field boundary conditions (Figure 2.10; Itasca 2004). Mechanical damping is applied to account for energy losses due to internal friction in the intact material and slippage along the discontinuities (Itasca 2004; Wyllie & Mah 2004).   Figure 2.10. Types of dynamic loading and boundary conditions for a flexible base (left) and a rigid base (right). Modified from Itasca 2004.  To assess the potential for failures in the eastern slope under an earthquake trigger (see section 6.1.3), the dynamic input was applied as a velocity history to the boundaries of the deformable blocks. The velocity history, which is treated as a multiplier by UDEC, was represented by a harmonic sine function. Because in UDEC the velocity input cannot be directly applied to quiet boundaries (see Figure 2.10), the velocity record was transformed into a stress record and then applied to a quiet boundary following the recommended procedure specified in the UDEC users manual (Itasca 2004). To do so, the following equations were used which assume plane-wave condition:  npn C νρσ )(2 ⋅=   [2.4]  2( )s s sCτ ρ υ= ⋅   [2.5]  ρ 3 4GK Cp +=   [2.6]   41 ρ GCs =   [2.7]  Where σn = applied normal stress;   τs = applied shear stress;   ρ = mass density;   Cp = speed of p-wave propagation through medium;   Cs = speed of s-wave propagation through medium;   νn = input normal particle velocity;   νs = input shear particle velocity.  The damping in the dynamic simulation should reproduce the energy losses occurring either in the joints or rock blocks of the natural system. In the forward analysis, Raleigh damping was used because it provides damping that is approximately frequency independent, similar to the hysteretic behaviour of natural damping (Itasca 2004).    42 3. ENGINEERING GEOLOGY OF THE STUDY AREA 3.1. Introduction The study area is located in Portillo, Chile, on the western flank of the central Andean Cordillera between 32°S and 33°S and approximately 70°05’W (Figure 3.1). In this region the average altitude is around 4,000 m.a.s.l, with summit altitudes greater than 6,000 m.a.s.l. The highest peak, Mount Aconcagua, has an elevation of 6,960 m and is the summit of the Americas (Figure 3.4).   Figure 3.1. Shaded relief map showing the location of the study area (red box). Source: http://www.lib.utexas.edu/maps/americas.html  43 3.2. Tectonic Setting & Seismicity The tectonic and seismic setting of central Chile is characterised by the subduction of the oceanic Nazca Plate beneath the continental South American lithosphere forming an east- dipping seismic zone. The collision of the plates takes place in the N78°E direction at a convergent rate of approximately 7.89 cm/yr (DeMets et al. 1994). The subduction zone along the Chilean trench is seismically active, where interplate underthrusting earthquakes (M>8) have occurred along the entire coast between latitudes 18°S and 46°S (Comte et al. 1986; Beck et al. 1998) with return periods in the order of 80 to 130 years for central Chile (Barrientos et al. 2004).  Several studies point to the segmentation of the subduction zone based on the maximum observed rupture length of major earthquakes together with along-strike variations in the dip angle of the Nazca Plate (Lavenu & Cembrano 1999; Pardo et al. 2002; Yañez et al. 2002; Barrientos et al. 2004) (Figure 3.2). Beneath central Chile and western Argentina (28°S – 33°S), where the study area is located, the subducted plate is almost horizontal, with a dip angle of less than 10° that extends east for hundreds of kilometres at a depth of approximately 100 km before resuming its downwards trend (Cahill & Isacks 1992). It is believed that in this segment, volcanism terminated approximately 9 – 10 Ma ago because of the geometry of the subducted oceanic slab (Kay et al. 1988), while there is active volcanism south of the 33°S latitude.   44  Figure 3.2. Major morphostructural units of the central Andes and schematic block diagrams showing the subducted plate segments (modified from Giambiagi & Ramos 2002; Fock et al. 2005). Note that the study area lies on top of the flat subduction zone (red box).  In general the seismic energy release in Chile originates from two sources: the subduction zone where large (M > 8), shallow (0-50 km) thrust earthquakes occur; and from compressional (and tensional) events of moderate to large magnitude (4 < M ≤ 7) and intermediate depth focus (70–100 km) that take place in the subducting Nazca Plate (Pardo et al. 2002; Barrientos et al. 2004). In the Andean region of central Chile, along the flat slab segment between 28°S–33°S, most of the earthquakes are compressional with moderate magnitude and depth focus of ~80 km (Pardo et al. 2002; Barrientos et al. 2004). Nevertheless, very shallow (0-20 km) damaging earthquakes (5.9 ≤ M ≤ 6.9) have also been recorded in the area (Lomnitz 1961; Barrientos & Eisenberg 1988; Barrientos et al. 2004). It is believed that the shallow events are produced by the deformation of the overriding continental plate in response to the differential coupling that exists in the region due to along- strike variations in the dip angle of the Nazca Plate. Figure 3.3 shows the distribution of intermediate and shallow earthquakes in central Chile (32°S – 34°S) and Table 3.1 lists the earthquakes that are contained within the map area. For brevity, the list only includes M > 5 intraplate earthquakes.   45  Figure 3.3. Earthquakes with mb or Ms greater than 4.0 and their epicentral area. Modified from Earthquake Hazard Program – USGS (http://neic.usgs.gov/). Rectangle indicates the study area.  Based on frequency-magnitude relationships in central Chile for the period between 1986 and 2001, the return period for events with magnitudes 3.5 ≤ M ≤ 5.5 is represented by Equation 3.1:  ( )log 7.8 0.1 1.4( 0.02)N M= ± − ±   [3.1]  In Equation 3.1, N is the number of events with magnitude equal or greater than M, therefore a magnitude 5 event should occur every 3 years and a magnitude 5.5 every 16 years (Barrientos et al. 2004). Considering that for higher magnitudes this equation is no longer valid due to a breakdown of the f-M relationship (Okuda et al. 1992), estimates on the return period for large thrust earthquakes (M>8) along the coast are mainly based on historical records since 1575 (Comte et al. 1986).  46 Table 3.1. Historical and recent seismic events with mb or Ms greater than 5.0. Data taken and modified from Earthquake Hazard Program – USGS (http://neic.usgs.gov/). YEAR MO DAY LAT (°S) LONG (°W) DEPTH (km) mb & Ms YEAR MO DAY LAT (°S) LONG (°W) DEPTH (km) mb & Ms 1575 3 17 33.4 70.6 7.3 1973 4 23 34.0 70.6 85 5.2 1687 7 12 32.8 70.7 7.3 1974 1 23 32.2 69.8 115 5.2 1850 12 6 33.8 70.2 7.3 1974 3 24 33.0 70.3 104 5.3 1927 4 14 32.0 69.5 110 7.1 1974 11 12 33.2 70.6 90 5.5 1928 1 22 33.5 69.5 5.5 1974 12 29 33.0 70.0 99 5.5 1928 10 17 34.0 69.5 5 1975 1 2 33.1 70.0 108 5.2 1928 12 6 33.5 69.5 5.4 1975 6 14 32.5 70.7 94 5.7 1929 7 27 34.0 69.5 5.7 1975 9 14 33.7 70.5 37 5.2 1930 1 9 33.0 69.5 5.8 1979 12 30 32.6 70.5 78 5.2 1930 3 27 33.5 70.0 5.1 1980 7 13 33.5 70.2 103 5.6 1931 8 17 32.5 69.5 120 5.8 1981 7 11 32.2 71.4 56 6.1 1931 11 24 33.0 69.5 5.3 1983 12 15 33.1 70.1 100 6.0 1932 5 8 32.5 69.5 110 5.5 1984 10 30 33.6 70.5 92 5.4 1932 5 10 32.0 70.0 110 5.5 1985 2 24 33.1 69.7 33 5.0 1932 7 9 32.5 69.5 5.1 1985 11 14 32.4 69.6 120 5.2 1932 11 29 32.0 71.0 110 6.8 1985 11 14 33.2 70.1 23 5.0 1933 6 5 33.0 70.0 5.2 1986 1 5 32.2 70.8 97 5.0 1933 10 19 32.0 70.0 5.7 1986 11 23 32.0 70.3 105 5.1 1933 11 14 32.0 69.5 110 6.5 1989 4 1 32.8 69.9 109 5.5 1934 6 9 32.7 69.7 5.8 1989 11 9 33.9 70.6 96 5.0 1942 6 29 32.0 71.0 100 6.9 1990 7 16 32.5 70.0 102 5.7 1945 9 13 33.3 70.5 100 7.1 1990 12 21 32.7 69.6 116 5.3 1955 11 4 33.5 69.5 100 6.8 1991 12 1 32.1 69.5 115 5.2 1958 9 4 33.5 69.5 6.7 1993 8 6 33.2 70.0 109 5.0 1958 9 4 33.8 70.2 10 6.7 1996 10 18 32.4 70.0 121 5.2 1963 10 6 33.9 70.0 101 5.1 1997 3 25 33.5 70.5 84 5.5 1964 9 10 33.0 69.8 95 5.2 1997 6 19 33.2 70.1 106 5.4 1965 5 3 32.4 70.3 81 5.6 1997 6 19 32.8 69.9 117 5.1 1966 1 15 33.5 70.0 17 5.2 1999 8 2 33.0 70.2 95 5.2 1966 1 15 33.5 69.8 50 5.5 2000 6 16 33.9 70.1 120 6.2 1967 9 26 33.5 70.5 78 5.7 2001 5 6 32.5 70.9 52 5.1 1969 12 13 32.8 70.1 103 5.4 2003 1 7 33.8 70.1 110 6.0 1970 4 9 33.9 70.1 120        5.2 2004 4 30 33.5 70.5 93 5.2 1972 10 2 33.9 70.9 80          5.3 2007 7 11 32.7 70.3 93 5.0    47 3.3. Geology & Geomorphology 3.3.1. Bedrock & structural geology The study area is located on the eastern border of an extensional volcano-tectonic basin that developed between the middle to late Eocene and Oligocene (Godoy et al. 1999; Charrier et al. 2002) on the western (Chilean) side of the Andean Cordillera. This basin was developed under extensional conditions and underwent subsequent tectonic inversion (Charrier et al. 2002) in response to the extreme compression that initiated around 26±1 Ma with the beginning of the orthogonal convergence of the Nazca and South American plates (Tebbens & Candie 1997).  At the toe of the Caracoles Range (Figure 3.4), the primary source area of the Portillo Rock Avalanche, a back-thrust was recognized which structurally separates W-dipping rocks of the hanging-wall from E-dipping rocks of the foot-wall. It is proposed  that the thrust fault that outcrops in the study area originated during a compressional episode ca. 10-8 Ma which is believed responsible for the uplift of the Aconcagua and its associated large scale thrusting (Aconcagua Fault and Thrust Belt – AFTB; Ramos et al. 2004) observable to the east of the study area (Godoy 2002; Yañez et al. 2002) (Figure 3.5).  Figure 3.6 shows the geological map of the study area. A description of the main rock units is given in the following paragraphs and an explanation of the unconsolidated units is given in section 3.3.2.  In the study area the bedrock corresponds to late Eocene - early Miocene andesitic lavas, volcaniclastic and sedimentary bedded rocks that dip between 50°W and 60ºW and belong to the Abanico Formation (Charrier et al. 2002) (EOla in Figure 3.6). Towards the base of the Caracoles Range, and separated by the thrust fault, the bedding planes dip approximately 30ºE. On the opposite side of the valley, the volcaniclastic rocks dip roughly 35ºW, gradually reaching a subhorizontal position towards the summit of the mountain (Figure 3.6).   48  Figure 3.4. Google Earth picture showing the structures that are recognized in the study area. Note the relief in the terrain.    Figure 3.5. Structural map and cross section of the Aconcagua Fault and Thrust Belt (modified from Ramos et al. 2004).    49  Figure 3.6. Geological map of the study area. Modified from Rivano et al. (1993).     50 The bedrock is occasionally intruded by minor stocks and sills (Mqd and Mdp in Figure 3.6, respectively) which commonly leave an imprint of altered and fractured bedrock around the intrusive contact. Localized relationships between these weakness zones and initiation zones for minor rockfalls have been recognized in the study area and its surroundings. 3.3.2. Quaternary geology The Last Glacial Maxima (LGM) in the Andes of northern-central Chile covered the area ~14,500 14C yr BP (Heusser 2003) and left the valleys with abundant evidence of recent glacial shaping (chiselled mountain peaks, cirques, and U-shaped valleys). Other than the morphology, indications of these glacial events are the moraines left in the high Andes, which in some occasions have been mapped as extending far beyond their likely true limits given the misidentification of landslide deposits as being glacial in origin (Heusser 2003). This is the case for most of the landslide deposits that crop out in the study area, which were initially interpreted as moraine deposits (Caviedes & Paskoff 1975; Rivano et al. 1993; Godoy 1994) due to its texture and distinctive lobate morphology.  The valley is marked by barren, boulder-covered lobes or ridges that belong to the Portillo Rock Avalanche debris (Plpra in Figure 3.6). The deposit can be divided into different granulometric domains (proximal, medial and distal) according to their relative position to the scarp (Figure 3.7).    51  Figure 3.7. Photos of the Portillo Rock Avalanche deposit along its path. a) Coarse proximal phase with scarce matrix; b) Switch-back deposit with predominance of metric cobbles immersed in a silty sand matrix; c) Hummocky terrain made up of fine sand, silt, and some clay containing some metric cobbles.  The proximal deposit has negligible to no matrix and is mainly composed of boulders and blocks as large as ~750 m3 (Figure 3.7-a). They correspond to angular self-supported blocks and boulders of similar lithology to the Abanico Formation, characterised by distinctive reddish volcanic breccias with large (greater than 5 cm) chlorite and epidote rich nodules (Figure 3.8). The fragments often show fault striations on their free surfaces. A ridge is recognized in the upper part of the valley which presents a finer granulometry characterised by cobbles of red volcanic breccias with green-epidote-rich nodules immersed in a fine gravel and sand matrix, the latter made up of pale biotite flakes, pyroxene and altered olivine crystals (Figure 3.9). Due to the granulometric differences and lithologic similarities between the deposits that out crop in the upper part of the valley, it is likely that the ridge deposit corresponds to a mix of stripped alluvial material and debris from the Portillo Rock Avalanche event.  52  Figure 3.8. Photo of a large epidote-rich nodule.   Figure 3.9. Aerial photo 1:20,000 which shows the ridge in the upper part of the valley and a detail of the debris.  Down the valley, the deposit appears to follow the paleo-topography which could be compared to that of a glacially carved steep undulating valley floor. Abrupt faces of coarse debris where the valley steepens and narrows, or plateaus more recently covered by alluvial sediments (Ha in Figure 3.6) where the valley floor flattens are common features. In the  53 switch-back zone (~2500 m.a.s.l) the medial (or transitional) deposits are finer and better graded, with the coarse portion made up of cobbles and coarse gravel. The lithology of the fragments is similar to the one described for the upper part of the valley, while the matrix, made up of sand and silt, is rich in feldspar, quartz, pyroxene, and epidote (Figure 3.7-b).  When the bottom of the valley is reached (~2200 m.a.s.l), some cobbles appear with coarse sand and gravel forming the coarse portion of the deposit, with a high content of fine sands, silt, and some clay as part of the matrix. At the distal edge of the rockslide deposit, the material takes the form of hummocky terrain (Figure 3.7-c). The base of the deposit is unexposed.  It is likely that the major impact of the Portillo Rock Avalanche on the prehistoric landscape was the obstruction of the NS valley drainage and the formation of Inca Lake (see Figure 3.7) located at the toe of the Caracoles Range (Godoy 1994). Field observations of the deposit in the upper segment of the valley strongly support this hypothesis.  In the study area, other deposits have been recognized as originating from landslide events (Pll in Figure 3.6). Of these, only the northern deposit has a recognized source scarp, while for the southern deposit the scarp is unclear (Figure 3.10). The northernmost debris is made up of blocks and cobbles of reddish volcanic breccias and light gray tuffs immersed in a sandy gravel matrix. The southern deposit is made up of blocks and cobbles with negligible matrix. Based on their stratigraphic position and relative relation with the other landslide deposits, they were assigned to the Pleistocene, being most likely older than the Portillo Rock Avalanche. Because of the uncertainty of their origin, age and volume, they were not included in the back analysis.   54  Figure 3.10. Other landslide deposits recognized in the area and assigned to the Pleistocene.  A set of Holocene events was also identified, which are characterised by two distinct deposits found along the east and west side of the valley and correspond to a rockslide and rock slump respectively (Hsl in Figure 3.6). The eastern deposits are made up of cobbles and boulders of andesitic lavas and breccias with some gravel and sand-rich matrix (Figure 3.11). Even though its scarp has not been clearly identified, it is evident that it came from the eastern wall aided by the steeply westward dip of the beds. These deposits overlie the PRA deposits and are overlaid by recent ocoitic-rich talus material (Godoy 1994).  The west lobe deposit is made up of self-supported angular boulders and blocks of porphyritic andesitic lavas, lapilli tuffs, fine grained andesites and characteristic boulders of white dacitic porphyry (Figure 3.11). The dacitic blocks come from the sill that out crops only  55 in the western wall and intrudes the volcaniclastic rocks of the Abanico Formation. Even though at the surface of the lobe there is a predominance of blocks, boulders and cobbles, deeper into the deposit finer material (gravel and sand) starts to appear. The true granulometric distribution of the deposit is indeterminable, but considering its short runout and the degree of jointing in the wall, it is most likely that the coarse fraction predominates.   Figure 3.11. Pictures of the Holocene deposits showing the eastern (top) and western (bottom) deposits.  Fine to very fine sediments (sand, silt and some montmorillonitic clay) are commonly observed in flat areas or plateaus that highlight the stepped topography of the valley at ~2930, ~2800, ~2525, and ~2150 m.a.s.l (Ha in Figure 3.6). Most of them are believed to have accumulated and filled open spaces left by the previously deposited PRA debris once drainage was resumed after the catastrophic event (Figure 3.12).  56  Figure 3.12. Pictures of the alluvial plateaus along the stepped topography of the glacial valley.    57 The colluvium (Hc in Figure 3.6) mainly corresponds to thick talus developed in the steep slopes of the valley sides (Figure 3.13). It is interpreted as forming due to the intense frost shattering of the rock exposed to the daily alterations of freezing and thawing (Caviedes & Paskoff 1975).   Figure 3.13. Picture of the talus developed in the valley sides.      58 3.4. Age Dating & History of Rock Slope Failure Events In the study area, a number of different deposits were identified as belonging to multiple pre- historic landslide events based mainly on the identification of different initiation zones, local relative relationships and, when possible, distinct lithology.  The surface exposure dating by Terrestrial Cosmogenic Nuclide (TCN) was used to help constrain the chronology of events identified from the geological/geomorphological investigations. This novel method helped place the deposits in an absolute timescale by taking samples from the break away wall and debris and dating them for 36Cl. 3.4.1. Theory & background This section provides a brief overview of the TCN theory and applications. A more comprehensive review is provided by Gosse & Phillips (2001).  The chlorine-36 nuclide is created primarily through two production avenues, spallation of 40Ca and 39K under cosmic ray bombardment, and neutron activation of 35Cl. The rates of accumulation are proportional to the cosmic ray flux and concentration of target nuclides in the surface material. Evidence of cosmic ray spallation is evidence that the material in question has been exposed as a surface of the body of which it is part, and gives a means of measuring the length of time of exposure (Gosse & Phillips 2001; Dunne & Elmore 2003).  Both spallation of Ca and K and neutron activation of Cl are affected by site-specific conditions. Elevation plays an important role in the production because the cosmic ray flux is attenuated by passage through the atmosphere, so the higher the altitude the higher the production rate. Also, snow and vegetation cover can significantly modify production rates, as can geometric shielding from nearby obstructions (Dunne et al. 1999; Dunne & Elmore 2003). Among the systematic errors associated with this analytic method, the most common are:    59    Poorly known exposure history,    Significant erosion of the exposed surface, and    Invalid assumptions of isotope production rates.  Depending on the surface preservation and exposure history, the dating technique has an effective range from the Pliocene (2.65 Ma) to the late Holocene (<10 ka). The advantage of this analytical method is that only the landslide material is required and that it is possible to obtain an exact date due to the fact that catastrophic events are instantaneous (Ballantyne et al. 1998 ). 3.4.2. Cosmogenic nuclide dating by 36Cl Lobes with overlying rims and levees formed by rock avalanche deposits have been mapped in piedmont areas and broad valleys of the NW Argentinean Andes as belonging to distinct events separated in time (Hermanns & Strecker 1999). Surface exposures dating of deposits and/or break away walls have confirmed this interpretation (Hermanns et al. 2001; Hermanns et al. 2004 ).  The geomorphology at Portillo, however, is more complicated due to its narrow valley setting, which was glaciated during the last glacial maximum (LGM), and thus, deposits were interpreted as being partially deposited by glaciers and partially by landslides (Caviedes & Paskoff 1975; Rivano et al. 1993; Godoy 1994).  Different lobes were mapped and sampled systematically for cosmogenic nuclide dating by 36Cl. A total of 17 samples were taken from the various lobes and the break away wall. As a preliminary step, 6 samples were dated to answer the following:     Whether the lobes are glacially derived or correspond to post-glacial slope failure events;    Whether the deposits in the valley belong to the same event or represent various landslide events separated in time.    60 The samples were prepared and concentrated at the PRIME Lab at Purdue University. The ages derived from CN for 36Cl were calculated using the program CHLOE developed by Phillips & Plummer (1996). 3.4.3. Sampling The location of the 17 samples is presented in Figure 3.14. For each location the guidelines presented by Gosse & Phillips (2001) were used and the following parameters were systematically recorded:   Figure 3.14. Location of the 17 samples obtained for surface exposure dating by Cosmogenic Nuclide using 36Cl. Stars indicate the samples that were dated.  61    Rationale for sample selection. In the runout zone, the most distinct lobes and those deposits clearly outlined in the aerial photointerpretation (API) and geological mapping were targeted. For the break away wall, the sampling was done where it was clearly determined that the area corresponded to the slip surface.     Description of object sampled. For every sampled lobe, large boulders with flat surfaces which showed no evidence of movement after the failure were selected (Figure 3.15). Faces corresponding to previously unexposed fractures sampled from the center of the boulder rather than close to corners or edges were preferred.     Description of sample material. For every sample, a lithological characterisation of the sampled material was carried out, including a hand-sample and thin-section description.     Location, orientation, and sample thickness. Spatial parameters such as location and orientation of the samples were recorded using a hand-held GPS compass and an inclinometer. The Dip-Dip direction (DDR) of the sample was measured only when the inclination was greater than 5°. The thickness of the sample was variable between 2 and 4 cm. In some cases, due to the characteristics of the rock it was not possible to obtain such a thickness, so small chips were sampled.     Shielding geometry. Sample shielding above the horizon was done by drawing a circle with azimuthal intervals (e.g. 0°-30°, 30°-60°, etc.). Major features were localized and the corresponding average angle between the horizon and the skyline was recorded for every interval. For the case of shielding due to snow cover, reasonable estimates were made related to average snow accumulation in Portillo. The samples were then corrected accordingly.      62  Figure 3.15. Pictures showing size, shape, and inclination of sampled boulders. All samples were collected with a hammer and a chisel.  3.4.4. History of events & failure scenarios The summary data displayed in Table 3.2 presents two different ages obtained by the superficial exposure dating which consider an erosion rate of 0.56 mm/ka and an erosion rate of 2.22 mm/ka. The choice of erosion rates was made following the work done by Costa & Gonzalez-Diaz (2007 ) in the Argentinian Patagonia, which assumed that erosion rates in hard volcanic rocks are quite low (typically 0.5–2 mm/ka).  Table 3.2. Data and ages obtained through 36Cl cosmogenic nuclides. P0402 P0404 P0405 P0410 P0507 P0512 Rock type Andesite Andesitic  Tuff Andesitic  Tuff Andesitic  Porphyry Andesitic  Tuff Andesite Latitude (S) 1 32.8 32.8 32.8 32.9 32.9 32.9 Longitude (W) 1 70.1 70.1 70.1 70.1 70.1 70.2 Elevation (m) 2990 3120 2952 2953 2574 2215 Orientation sample surface (º) 0 40 36 20 0 18 Block geometry (m) 2 2.5x1.0x1.2 wall 1.0x0.5x0.5 4.0x3.0x2.5 5.0x4.0x2.0 2.5x2.0x1.5 Maximum angle of shielding geometry (°) 33 40 34 30 40 44 Calculated age (ka) for erosion rate: 0.56 mm/1000  yr 4.1 ± 0.3 12.7 ± 1.0 4.2 ± 0.5 12.1 ± 0.5 14.2 ± 1.0 13.6 ± 0.6 Calculated age (ka) for erosion rate: 2.22 mm/1000 yr 4.1 ± 0.3 12.3 ± 1.0 4.1 ± 0.5 11.5 ± 0.5 13.2 ± 0.9 12.9 ± 0.5 (1) Datum WGS 84 (2) Geometry corresponds with LxWxH (3) Ages are reported with an error at the level of 1σ  63 All data point towards a post glacial origin for the deposits and the occurrence of at least two separate prehistoric events in the Portillo valley. Due to the overlapping uncertainties within the Upper Pleistocene and Holocene values, precise intervals of landslide events may be determined with ages clustering in the range of 15.6 to 11.6 ka with a mean of ca. 13.2 ka, and a range of 4.7 to 3.7 ka with a mean of 4.2 ka, respectively.  The older deposit corresponds to the Portillo Rock Avalanche event and is dated at about 13.2 ka. The younger event(s) is dated to 4.2 ka and is associated with two different lobate deposits located at the bottom of both the east and west slopes in the upper part of the glacial valley (S1 and S2 in Figure 3.16, respectively). Considering the significant reworking that has occurred as part of the operations related to the ski resort located in the upper parts of the glacial valley, it is difficult to establish whether the 4.2 ka event derives entirely from the Caracoles Range along the east side of the valley and/or possibly in part from the west side of the valley as well.  Even though the largest boulders possible were always sampled to minimize the effects of erosion which may lead to underestimated ages (Gosse & Phillips 2001; Hermanns et al. 2004 ), some variations due to erosion processes cannot be excluded, and thus all ages represent lower bound limits of the absolute age.   64  Figure 3.16. Aerial photograph (1:80,000) showing the distribution of the slope failure events, their respective scarps and approximate ages. The star indicates the approximate location of the Portillo Ski Resort.    65 3.4.4.1. Upper Pleistocene Portillo Rock Avalanche Comparing the two events and related deposits, the Portillo Rock Avalanche debris is the only one that can be identified as having a common and characteristic lithology and well- defined source area. Its sliding surface is largely planar, rising up to a peak of 4050 m.a.s.l with a vertical relief of ca. 1000 m.  Considering the uncertainty inherent in reconstructing a prehistoric slope and its associated pre-failure geometry, two different scenarios were accounted for: a failed slab with slope parallel daylighting structures and a failed slab without slope parallel daylighting structures (Figure 3.17-a,b). The key constraint for the pre-failure topography was the estimated thickness of the slab (see Chapter 4, point 4.1.2 for more details) and the fact that the sliding surface corresponded with the current dip-slope. Based on the likelihood of occurrence and its potential contribution to the understanding of the failure mechanism, the scenarios were ranked and only one was considered for further analysis.   Figure 3.17. Cross sections of the potential pre-failure scenarios. Only the predominant joint sets are considered. a) Failed slab with daylighting structures; b) failed slab without daylighting structures. (For location of cross-sections see Figure 4.1-b).  The control for the presence or absence of slope parallel structures daylighting at the toe of the failed slab is given by the position of the thrust fault. Even though the trace of the fault is completely covered by talus, its existence is evident because in the upper part of the eastern slope the volcaniclastic rocks dip to the west, and at the toe of the slope the same rocks dip towards the east. Two separate outcrops of east-dipping rocks help constrain the lowermost position of the thrust fault at the toe of the Caracoles Range, while the exposure of west  66 dipping rocks above the talus constrain the uppermost position (see geological map of Figure 3.6).  In order to have daylighting structures, the fault must appear at the lowest elevation permitted by the outcrop of volcaniclastic rocks that dip to the East. This scenario was discarded for two reasons:     Its appearance would coincide with the lowermost position possible, which does not fit well with field observations,    The daylighting scenario would have had a clear behaviour in terms of its mechanism of failure, which would have been primarily sliding along the bedding planes that would likely have been unstable long before glaciation.  Therefore, based on the reasons abovementioned, it was considered unnecessary to replicate the geometry with daylighting structures at the toe of the slope. It also must be mentioned that the position of the fault increased or decreased the amount of failed volume by approximately 10%, depending on whether it was placed in the lowermost or uppermost position, respectively.  Given the steepness of the bedding/slope, one that likely exceeds the frictional strength along the bedding planes, and in absence of any more recent rock-slide events, kinematic release may have been triggered by a large earthquake, perhaps aided in part by oversteepening of the valley walls during glacial advance and debuttressing during retreat (i.e. after the LGM, approximately 14.5 ka). A progressive strength degradation mechanism may also have contributed to failure, both in the form of asperity breakdown along the bedding planes and loss of coherence in the stronger units at the toe of the slope through brittle fracture processes.     67 3.4.4.2. Holocene rockslide and rock slump The NNE-SSW trending lobe at the toe of the W slope (see Figure 3.11, bottom) was initially interpreted and mapped as a lateral moraine due to its distinctive morphology (Caviedes & Paskoff 1975; Rivano et al. 1993; Godoy 1994). This interpretation was discounted once the cosmogenic nuclide dating by 36Cl was available and indicated that all the sampled lobes were post-glacial in age. These results were also supported by the field observations and detailed inspection of the deposits.  In the upper parts of the valley, significant reworking of the material has occurred as part of the operations related to a ski resort (see Figure 3.16 for location). This alteration of the terrain makes it difficult to establish whether the two younger deposits are derived entirely from the eastern Caracoles Range or correspond to two simultaneous events derived from the east (S1 deposit) and west sides of the valley (S2 deposit), most likely being triggered by an earthquake (given their temporal and spatial proximity).  In the case of the first hypothesis, the leading edge of the rockslide must have had to travel approximately 2 km westward across the valley overriding the previously deposited Portillo rockslide debris, run up the base of the west slope and land at its toe as it fell back forming a fall-back ridge.  For the second hypothesis, the source of the lobate deposit found at the foot of the western slope could be explained as originating from a rock slump/collapse whose short runout could be attributed to unfavourable geometrical conditions for motion created by the dip of the beds on that side of the valley (i.e. into the slope). In this case, the lack of continuity between the eastern and western deposits would then need to be explained by the reworking of the terrain for the development of the ski resort.  Comparing the arguments supporting each hypothesis, it was decided that the best alternative was that of multiple initiation sources and simultaneous events. This scenario was later supported by the dynamic runout analyses performed, which are explained in detail in Chapter 5.  68 3.5. Geotechnical Characterization of the East Slope 3.5.1. Description of discontinuity network To assess the geotechnical characteristics of the rock slopes in the study area, a geological engineering survey was performed in the field which included geological description, structural measurements, rock mass characterisation, and general discontinuity surveys. Due to the steepness of the slope, only the base of the wall was surveyed in detail. The characterization was complemented by the use of a terrestrial-based laser scanner (LiDAR).  The slope was divided into two distinct domains to account for the different degree of rock mass damage observed in the slope due to the presence of the thrust fault: domain 1 (D1) which corresponds with the hanging wall and includes the upper slope, and domain 2 (D2) which corresponds with the footwall and includes the zone of damage rock (Figure 3.18). The rationale behind this was to better assess the rock mass strength and have a closer representation of reality in the numerical models.   69  Figure 3.18. Current topographic profile of the eastern slope showing the different degree of fracturing affecting D1 & D2. In the topographic profile, solid line indicates rock exposure and broken line indicates unconsolidated deposits (talus and rockslide debris).   70 Several terrestrial laser scans were performed in D1, enabling the extraction of the discontinuity patterns, including orientation and spacing. For D2, only random mapping was done due to the limited access to the outcrop.  Two dominant discontinuity sets were recognized in D1: a persistent bedding plane (JS1) and a cross-cutting joint set (JS2). It is believed that the generation of JS2 is syntectonic to the generation of the thrust fault recognized at the toe of the eastern slope, which explains its continuity into D2. In the highly fractured D2, the layering of the beds was considered as another joint set (JS3) in addition to JS2.  The bedding planes and cross-joints were assumed to be fully persistent based on their geological origin, namely stratification in the volcaniclastic rocks together with their tectonic fabric respectively. The uncertainty associated with this assumption is higher for the case of the cross-joints. The results of the discontinuity assessment for the eastern slope are shown in Table 3.3. These values were used in the numerical models presented in Chapters 4 and 6.  Table 3.3. Summary of discontinuity parameters obtained from the eastern slope. min (-σ) max (+σ) avg min (-σ) max (+σ) avg min (-σ) max (+σ) avg outcrop LiDAR Bedding (JS1) N=131 Cross-joints (JS2) N=38 Bedding (JS3) N=6 Cross-joints (JS2) N=4 D om ai n 1 D om ai n 2 Dip-direction (°) 48/225 62/255 55/240 22/19 40/45 31/032 30/117 40/129 35/123 33/008 37/014 35/011 Outcrop Spacing (m) LiDAR Spacing (m) 3 7 5 20.2 14.45 4 8.68 12.165 9 7 2 4 3 Both joint sets also develop a centimetric spacing generating a densely fractured rock mass      JS3/JS2 = 1/1 N/A avg spacing ratio (JS i /JS j ) 5/7 = 0.71 3/4 = 0.75 29.08 20.42 3 5 ⇒     71 3.5.1.1. Orientation The joint orientation data obtained from the analysis of the 3-D point cloud was exported and plotted in lower hemisphere stereographic projections and analysed using the software SpheriStat 2.2 (Stesky & Pearce 2007). The plotted data was grouped by density of discontinuities using contour plots, and the principal joint sets were defined (Figure 3.19).   Figure 3.19. Contoured stereographic projection of selected automatically generated patches from scans 1-7, 9, 11-12, together with outcrop measurements, for D1. Joint set orientations are measured as dip/dip direction.  In domain 1, the most important joint set in terms of its persistence and spacing is the one generated by the bedding planes (JS1), whose principal orientation is 55/240. Most of the scans were able to pick the bedding planes and generate patches large enough to give  72 confidence to the remote data. The coherence between the hand measurements and the remote data also add to the reliability of the orientation determined for JS1. The large surface extent of the patches was also used as an indicator of the smoothness and persistence of JS1, while the scatter in its azimuth was interpreted as a large-amplitude undulation.  A slope-scale feature which corresponds to a persistent fracture (F1) with a mean orientation of 62/177 (± 9/008) and a very low fracture frequency was remotely recorded by scans 6, 7, & 9. Its orientation with respect to JS1 suggests that the failed slab could have been divided into large tabular blocks due to the cross-cutting of this major fracture, which would have helped in the release of the slide mass with sliding along the bedding planes (Figure 3.19). This fracture was not considered in the numerical models because its participation in the slope failure was mainly related to the lateral release of the slab in the out of plane direction to the 2-D cross-section used in the UDEC modelling (see Figure 4.1).  The orientation of the cross joints (JS2) was mainly estimated from the automatically generated patches of scans 6 & 7 with a mean orientation of 31/032. Even though this joint set can be very clearly seen in the main wall (Figure 3.20), they are hard to detect in the processed LiDAR data, only showing up in small patches with a high uncertainty. The low reflectivity of the surfaces due to the over-hanging condition of these joints was one factor that contributed to making them more difficult to recognize in the point clouds. Field measurements were therefore used to help constrain the orientation of these joint set and give information on the joint profile, which corresponds to undulating rough surfaces following Barton (1973). The variation of the strike of JS2 between the two domains could be explained either by an error associated with the automatic patch generation, or by a large scale folding of the structures very common in this structural setting.  73  Figure 3.20. Detailed picture of the main fracture and stereographic projection showing the kinematics between JS1 & F1.  In domain 2, the orientation of JS2 & JS3 was obtained solely through hand measurements due to the densely fractured characteristic of the outcrop which compromises the effectiveness of the LiDAR when used in this type of rock mass. The high persistence and close spacing (< 1 m) of the joint sets was used as a valid argument to consider the few localized measurements collected as representative of the whole rock mass.   74 The mean orientation of JS2 and JS3 in D2 is 35/011 and 35/123, respectively (Figure 3.21). Both joint sets present a characteristic dense fracturing, which depending on the scale, ranges from centimetres to meters (Figure 3.18 bottom right and left, respectively). In this domain the fracturing is a clear representation of the rock damage produced by the thrust fault.   Figure 3.21. Contoured stereographic projection of hand measurements taken in D2 for JS2 & JS3. Joint set orientations are measured as dip/dip direction.   75 3.5.1.2. Spacing & block volume The average in situ block size for D1 obtained from direct observation and characterisation of the base of the slope ranges between 225 to 945 m3. These calculated block volumes were made by using an estimated spacing of 15 m for F1, 5 times the spacing of the main joint set as suggested by Palmstrom (2005) for large joint spacing. When the influence of the thrust fault at the base of D1 is considered, a closer spacing ranging between 0.5 and 1.5 m is observed for JS1 and JS2.  Accounting for the size of the slope to be modelled (i.e. approx. 1000 m high), scaling of the in situ block size was deemed necessary for computational efficiency. For this, the scaled block size used for the 2-D models was constrained by the average ratio JS1/JS2 (Table 3.3). In D1, the spacing ratio JS1/JS2 varies between 3/5 and 7/9, with a mean of JS1/JS2 = 5/7. When the LiDAR data is used to estimate the joint spacing the results show larger spacing for the discontinuity sets, with an average spacing ratio JS1/JS2 = 3/4.  The lack of LiDAR data in D2 required that digital photos be used to estimate the block size scaling ratio. A similar analysis was followed for D2 and showed that at the meter scale the JS3/JS2 ratio varied between 2/3 and 4/5, while at the centimetre scale the spacing ratio became JS3/JS2 = 1/1. In the numerical models the latter was the ratio used because it better represent the relation between discontinuities in the fractured rock mass. 3.5.2. Laboratory geomechanical testing Representative rock samples from the eastern slope and the Portillo Rock Avalanche deposits were cored and tested using the facilities of the Rock Mechanics Laboratory of the Mining Department and the Rockphysics Laboratory of the Earth and Ocean Sciences Department, both from the University of British Columbia. Triaxial and uniaxial compression tests and direct shear tests on rock discontinuities were performed. Basic index properties like density were also measured. The results were used to compute elastic modulus, Poisson´s ratio, and intact rock strength. The rock mass and discontinuity properties were then calculated using the empirical scaling relationships described in section 3.5.3.  76 3.5.2.1. Intact rock All tested samples were taken from volcanic and volcaniclastic rocks of the Abanico Formation collected in the field. Detailed test results are included in Appendix A. Rock samples were cored following the standard D 4543–01 (ASTM 2001) and their specific gravity was measured following the standard D 5779–95a (ASTM 1996a). Rock strength was determined using the uniaxial compressive test method described in the standard D 2938–95 (ASTM 1995) and the triaxial compression tests followed the standard D 2664–95a (ASTM 1996b).  Uniaxial compression tests were carried out using a MTS servo-controlled stiff testing machine. Axial and circumferential strain measurements were obtained using high precision extensometers attached to the centre of the specimen (Figure 3.22-a). Axial stress and strain measurements were automatically logged at evenly spaced time intervals. The triaxial compressive strength tests were conducted in a standard triaxial pressure vessel within a servo-controlled GCTS 1000 kN (225 kips) loading frame (Figure 3.22-b). All tests were performed in dry conditions at room temperature.   Figure 3.22. Picture of samples for a) Uniaxial loading, sample LI-10, H/d ratio of 2:1; b) Triaxial loading, sample LI-8.2, H/d ratio of 2:1.    77 Two andesite and two andesitic tuff samples representative of the lithologies at the site were tested under uniaxial compression for the unconfined intact strength (UCS), Young’s modulus and Poisson’s ratio (Table 3.4). The Young’s modulus was obtained from the slope of the axial stress strain curve at 50% UCS. The Poisson’s ratio was determined from the ratio of the circumferential and axial strain. All samples failed through a combination of axial fractures and shear (Figure 3.23).  Table 3.4. Average values of intact rock strength parameters for uniaxial tests. Sample Rx Type Diameter: d (mm) Height: H (mm) Ratio: H/d Loading Rate (kN/min) UCS (MPa) Young Modulus: E ( GPa ) Poisson Ratio: ν LI-10.1 Andesite 46.42 92.64 2.00 18.00 97.2  42.31 0.17 LI-10.2 Andesite 46.43 82.71 1.78 16.00 97.4  42.29 0.18 LI-7 Andesitic Tuff 46.47 88.53 1.91 18.00 241.7  42.08 0.30 LI-7.1 Andesitic Tuff 46.42 91.39 1.97 16.00 258.2  55.61 0.34     Figure 3.23. Shear failure of samples after failure subject to uniaxial and triaxial loading.  Three volcaniclastic samples were investigated for their static triaxial compressive strength properties. The tests were performed using confining pressures of 5, 15 and 25 MPa. The results, shown in Table 3.5, were used to define the Mohr-Coulomb failure envelope and obtain the intact rock cohesion and friction strength properties (Figure 3.24).  78 Table 3.5. Triaxial compression test results. Sample Lithology Confining     Stress (MPa) Axial Stress at Failure (MPa) E (GPa) ν Cohesion (MPa) Friction (°) LI-8.3 5.4 333.7 40.8 0.18 LI-8.2 15.2 425.8 41.39 0.22 LI-8.1 24.4 436.2 45.7 0.2 Andesitic lapilli tuff 56 49     Figure 3.24. Mohr stress circles and shear strength failure envelope obtained from triaxial testing of intact rock.  79 3.5.2.2. Discontinuities Two direct shear strength tests were performed (standard D 5607–02; ASTM 2002) on samples obtained from the eastern slope corresponding to a mated bedding plane (JS1 from D1; see section 3.5.1). Three normal stresses were applied to each sample with approximate values of 270, 515 & 750 kPa, starting from the lower load and ending with the higher one. The results of the tests are presented in Table 3.6 and plotted in Figure 3.25.  Table 3.6. Direct shear tests results for samples LI-7.1 & LI-7.2 containing mate bedding planes. σ n  levels (kPa) 269.19 507.05 741.26 275.06 518.12 757.44 size (cm 2 ) φ  peak (º) φ  residual (º) LI-7.1 10.3 x 6.5 LI-7.2 10.4 x 6.3 51 47 46 45     Figure 3.25. Failure envelopes for peak and residual shear strength (τ, σn space) from direct shear testing of mated discontinuity samples for JS1.  Peak & residual joint friction angles for sample LI-7-1 τ = 1.22σn R2 = 0.15 τ = 1.04σn R2 = 0.94 0 200 400 600 800 1000 1200 0 100 200 300 400 500 600 700 800 900 1000 normal stress (kPa) sh ea r s tr es s  (k Pa )      (k Pa ) Φresidual = 46º Φpeak = 51º Peak & residual joint friction angles for sample LI-7-2 τ = 1.06σn R2 = 0.93 τ = 0.99σn R2 = 0.97 0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200 normal stress (kPa)     s he ar  s tr es s   (k Pa )   (k Pa ) Φresidual = 45º Φpeak = 47º  80 Even though at the scale of the slope the sliding surface, which corresponds with the principal joint set (JS1), looks smooth and slightly stepped, at a more detailed scale (i.e. centimetres) the high peak friction angles obtained, between 47° and 51°, may be explained by the degree of roughness and asperities encountered along the JS1 surfaces (Figure 3.26).   Figure 3.26. Detailed pictures obtained from the rough surface of the principal joint set (JS1).  3.5.3. Rock mass characterization & properties Rock mass characteristics were assessed in the field using the Geological Strength Index (GSI, Hoek & Brown 1997), and the Rock Mass Rating (RMR, Bieniawski 1973). Considering the rationale behind the GSI, which provides a means to estimate the rock mass strength for different geological conditions, it was evaluated as a range of values rather than as a unique value for each domain. Similarly, as discussed in the previous chapter, a modified RMR was evaluated for each domain for which dry conditions and favourable joint orientation with respect to the slope were applied (i.e. RMR’).  81 3.5.3.1. Rock Mass Rating (RMR’) & Geological Strength Index (GSI) To obtain the RMR of the slope, the RQD was estimated utilising the volumetric joint count Jv (Palmstrom 1985) defined as the number of joints intersecting a volume of 1 m3 as follows:  n v SSSS J 1.......111 321 ++++=   [3.2]  where Si (with i = 1, 2……n) is the average spacing for the joint sets. Once Jv was calculated, it was used as input in Equation 3.3 which relates RQD with Jv as follows (Palmstrom 2005):  vJRQD 5.2110 −=   [3.3]  Figure 3.27 shows the RMR’89 and GSI tables with the respective values chosen for domains 1 and 2. The common parameters in both the GSI and RMR systems are those related to the rock structure and joint surface conditions where the rock structure may be quantified by the block size (Tzamos & Sofianos 2007).   82  Figure 3.27. RMR’89 and GSI values for D1 & D2.  83 The rock mass classification values are summarised in Table 3.7. The RMR’ and equivalent GSI values indicate a good quality rock mass for D1 and fair for D2. Following the RMR classification, the rock mass in D1 is in general GOOD to VERY GOOD, with rough to smooth discontinuity conditions, centimetre to metre spacing and moderately weathered surfaces. In D2, the rock mass is clearly affected by the presence of the thrust fault leaving a FAIR and locally disturbed rock mass. The discontinuity conditions are smooth to slickensided, with centimetre spacing and moderately to highly weathered surfaces. In both domains the discontinuities are considered highly persistent (> 20 m), with 1 – 5 mm of separation and no infilling material.  Table 3.7. Rock mass classification at D1 & D2. GSI RMR' 89 GSI RMR89 (1) GSI classification RMR classification  (1) Equivalent GSI obtained from: GSI = RMR89 - 5 (see section 1.4.1) DOMAIN 1 DOMAIN 2 70 - 50 50 - 35 GOOD to VERY GOOD FAIR 56 - 4482 - 67 77 - 62 51 - 39 BLOCKY to VERY BLOCKY VERY BLOCKY to DISTURBED   Table 3.7 also shows that the GSI values calculated using RMR’ are comparable to the GSI values assessed directly in the field, although slightly lower. These varied from 70 to 50 for D1, corresponding with a BLOCKY to VERY BLOCKY rock mass with discontinuities having good to fair surface conditions, to a GSI of 50 to 35 for D2, corresponding to a VERY BLOCKY to DISTURBED rock mass with surface conditions varying from poor to fair.   84 3.5.3.2. Adjusted rock mass properties for numerical models As previously discussed, in order to make the numerical models efficient and run in a reasonable period of time, it is not possible to use the exact in situ block size and its corresponding intact rock properties. Instead, the intact rock properties obtained from laboratory testing can be scaled to account for the discontinuities incorporated into the equivalent continuum when adopting a larger block size using the GSI (Figure 3.28).   Figure 3.28. Diagram showing the relationship between the block size and rock strength properties.  The scaling of the intact rock properties was performed using RocLab 1.0 (Rocscience 2005), which is based on the Hoek-Brown strength criterion. A table of rock mass strength properties (Table 3.8) was constructed and guided by consulting specific literature. Details of the RocLab 1.0 outputs are presented in Appendix B. The Mohr-Coulomb strength parameters obtained with RocLab 1.0 were calculated considering a σ3max for slopes given by Equation 1.16 in section 1.4.1.2.     85 Table 3.8. Estimated range of intact rock and rock mass parameters for D1 & D2. DOMAIN 1 DOMAIN 2 SOURCE OF INFORMATION UCS (MPa) 90 - 250 30 - 60 Lab testing; Portillo Report1 GSI 50 - 70 35 - 50 Field Estimate Intact rock parameter: mi Estimated with RocLab 1.0 using the triaxial test results Intact Young's Modulus: Ei (GPa) 40 - 55   9 - 30 Lab testing; RocLab 1.0 using the modulus ratio MR (Hoek & Diedrichs, 2006) Poisson's Ratio: ν Lab testing Density (Kg/m 3 ) 2700 2600 Lab testing; Portillo Report1 Hoek-Brown parameter: m b  2 - 4  0.9 - 1.8 Estimated with RocLab 1.0 Hoek-Brown parameter: s 0.004 - 0.04 0.0007 - 0.004 Estimated with RocLab 1.0 Hoek-Brown parameter: a Estimated with RocLab 1.0 Cohesion (MPa)  4 - 10  2 - 4 Estimated with RocLab 1.0 Friction angle (°) 33 - 46 20 - 30 Estimated with RocLab 1.0 Tensile Strength:     T 0 (MPa) 0.2 - 2  0.02 - 0.1 Estimated with RocLab 1.0 Rock Mass Deformation Modulus: Em (GPa)  12 - 40  1 - 9 Estimated with RocLab 1.0 Bulk Modulus:          K (GPa)  10 - 22 0.9 - 5 Shear Modulus:       G (GPa)  5 - 16  0.4 - 4 (1) Unpublished report done by ASC Ingenieria Ltda. for the Portillo Ski Resort 11.3 0.5 0.3 - 0.2In pu t: in ta ct  ro ck  p ro pe rti es O ut pu t: ro ck  m as s pr op er tie s )1(2 υ+⋅= EmG )21(3 υ−⋅= Em K   Table 3.9 describes the properties for numerical modeling of rock discontinuities. The range of joint friction angles was obtained both from the lab tests and published data by Kulhawy (1975) for selected andesites and sandstones (equivalent texture of volcaniclastic rocks). The latter represent the lower bound of the adopted joint friction angles.  The UDEC models also require a joint normal and shear stiffness (kn & ks respectively) which relates the normal and shear stresses on the joint with the normal and shear displacement  86 (Huang et al. 1995). Joint stiffness is not an easily measured or well known parameter; therefore some methods of estimating joint stiffness have been derived such as the one developed by Barton (1972). This method is based on the deformation properties of the rock mass and the intact rock, and assumes a single joint set with an average spacing L oriented perpendicular to the direction of loading. For this case, the average spacing between JS1- JS2 in D1 and JS3-JS2 in D2 was considered. The relation between the joint stiffness and intact/rock mass modulus can be explained by the following equations:  ( ) i m n i m E Ek L E E ⋅= ⋅ −   [3.4]  ( ) i m s i m G Gk L G G ⋅= ⋅ −   [3.5]  Table 3.9. Estimated range of discontinuity strength and stiffness parameters. Bedding JS1 Cross-joints JS2 Bedding JS3 Cross-joints JS2 Joint friction angle (º) Lab testing;                  Kulhawy, 1975 Joint cohesion (kPa) Kulhawy, 1975 Joint normal stiffness:  Kn (GPa/m) Joint shear stiffness:  Ks (GPa/m) L: joint spacing 30 - 50 20 - 30 0 - 1000 0 - 50   3 - 25  1 - 12  1 - 10 Discontinuity property DOMAIN 1 SOURCE OF INFORMATION DOMAIN 2  0.5 - 6      87 3.5.4. Hydrogeology Little information could be derived for the hydrogeological conditions at the site. No seeping water was observed and the slope contains no vegetation that would suggest the presence of water. Based on this information it is inferred that the slope is mostly dry with the water level near the toe of the slope coinciding with Inca Lake at 2800 m.a.s.l.  Wyllie & Mah (2004) point out that sometimes the seepage rate is lower than the evaporation rate resulting in a slope that is dry in appearance but with significant pore water pressure within the rock mass. Assuming that groundwater flow for the eastern slope would be controlled by fracture permeability, the discontinuity network derived for this study, although simplified, suggests that it would flow parallel to bedding in the upper parts of the slope until it reaches the flow barrier created by the fault at the toe of the slope. Here pore water pressures and flow may dissipate through the more heavily fractured rock mass. Again, it should be noted that no indication of seeping was observed here.  Because it is not possible to completely discard the presence of groundwater, even though the slope appears to be dry, due to the possibility of seasonal variations a varying water table was considered in the forward modelling to quantify the effects of elevated water pressures on the stability of the present-day slope.    88 4. FAILURE MECHANISM OF THE PREHISTORIC PORTILLO ROCK AVALANCHE 4.1. Model Development Numerical modeling techniques based on discontinuum mechanics present a valuable means to investigate the role of, and interaction between, geological structures and internal rock mass deformation of a jointed rock slope (Eberhardt et al. 2004; Stead et al. 2006). Based on the site investigation and assumed characteristics for the prehistoric Portillo Rock Avalanche, the distinct element code UDEC (Itasca 2004) was used to determine the most likely failure mechanism with varying geometries and strength parameters conforming to the estimated failed volume of the rock mass. 4.1.1. Failure initiation Initial assessments based on visual observations suggest that failure occurred primarily along slope parallel bedding planes in a translational manner. However, analysis of the post- failure structural and kinematic conditions of the eastern slope wall reveals no daylighting of persistent structures. Given that the steepness of the bedding/slope (65º-50º) exceeds the frictional strength along the bedding planes (the upper limit of which was determined to be 50º), and in absence of any recent rock-slide events, kinematic release is assumed to have occurred through yielding of the weaker fault material at the toe of the slope due to the loading imposed by the upper slope. With failure of the toe material, which would have acted as a buttress, kinematic released of the upper slope by sliding along the highly persistent bedding planes would be enabled.     89 4.1.2. Pre-failure topography & failed volume estimate The cross-section to be modelled was selected based on the assumption that the mass slid in the direction of the steepest dip of the slope face (Figure 4.1-a,b). The pre-failure topography of the Caracoles slope was built accounting for the current concave shape in cross section and the absence of slope parallel daylighting structures at the toe. To do so, the dip of the bedding planes varied from 65°W in the upper part of the slope to 50°W in its lower section (Figure 3.17-b). The fault was placed on an average location according to its uppermost and lowermost position as previously discussed in Chapter 3, section 3.4.4.1.   Figure 4.1. a) Contour map of the study area with vector indicating steepness of the slope; b) Tin map showing elevation and location of the cross-section.  90 Volume estimates for the Portillo Rock Avalanche were conducted using the empirical correlation for major rockslides proposed by Li (1983), that relates the logarithm of the aerial exposure (A) in m2 with the logarithm of the volume (V) in the form (Equation 4.1):  Log(A) = 1.9 + 0.57*Log(V)  [4.1]  With the aid of ArcGIS 9.0 (ESRI 2004) a planar area for the Portillo Rock Avalanche deposits was approximated as being 2.1x106 m2, the calculated volume of the debris was 68x106 m3. Figure 4.2 shows the location of the rock slide in the mean empirical relationship plot obtained by Li (1983), with upper and lower bounds calculated based on the standard deviation of the data. A fragmentation volume increase of 25% was considered due to break- up and transport of the detached mass, following the recommendations suggested by Hungr & Evans (2004). Neglecting any volume contributed through entrainment, this would give a pre-failure (in place) volume of the rockslide of approximately 51x106 m3 and an estimated thickness of the tabular slide mass of 65 m ± 10 m. The thickness was calculated based on the surface estimate of the detached zone of approximately 88x104 m2 ± 10x104 m2, which was obtained from a DTM using the software Surfer 8.0 (Golden_Software 2002).   Figure 4.2. Area covered by the Portillo Rock Avalanche deposits and its relationship to the volume as proposed by Li (1983). Broken lines represent the standard deviation of the data from the mean.  91 4.1.3. Model configuration Based on the field observations, a number of assumptions were made and used to build the geometry of the prehistoric slope and constrain the numerical models (Table 4.1). A modelling procedure was created and followed for all the numerical analyses performed and is shown as a flow diagram in Figure 4.3.  Table 4.1 Basic assumptions made for the development of the numerical models. ASSUMPTIONS MODEL CONSTRAINS Bedding dip is parallel to topography No daylighting structures Movement is in the direction of steepest dip 2-D analysis (plane strain) Discontinuities are persistent Slope and failed mass divided into tabular blocks The slope is dry No pore pressure   The model development started with the generation of the geometry. The 2-D cross-section was divided into three regions according to the slope characterisation:     Domain 1 (D1): made up of steeply dipping volcaniclastic beds,    Domain 2 (D2): corresponds with the toe of the slope, and    Base: represents the far-field portion of the slope.  The modelling procedure continued with the addition of the key controlling discontinuity sets using orientations and spacings obtained directly by outcrop mapping and remotely aided by a terrestrial-based LiDAR scanner (Chapter 3, section 3.5.1). A simplified model of the slope was constructed using a bedding plane spacing of 30 m (JS1) and a spacing of 40 m for the east-dipping cross-joints (JS2) in domain 1. For domain 2, a 20 m spacing was adopted for both the east-dipping bedding planes and cross-joints (JS3 & JS2, respectively) to more accurately represent the fractured nature of the zone of damaged rock (Figure 4.4). This concave surface profile was similar to that used by Benko & Stead (1998) for their analysis of the Frank slide.   92  Figure 4.3. Flow chart of the modelling procedure utilised in this thesis.   Figure 4.4. Concave geometry without daylighting structures at the toe of the slope used for the numerical models in UDEC. a) Location of the history points; b) Automatic mesh generation made up of quadratic and triangular-shaped finite difference zones. X-axis is meters from the origin; y-axis is meters above sea level.  93 Based on the table of reasonable values established during the data collection (see Chapter 3, Table 3.8 & Table 3.9), mechanical properties were selected as input parameters for constructing the model that best represented the pre-historic state in terms of proposed failure mechanism and estimated failed volume. A first trial was done using the average values of the mechanical properties but the proposed failure mechanism was not clearly represented (self-stabilized; see section 4.3.3). The best representation was found with the use of close-to-lower bound strength and stiffness properties, therefore that trial was selected as the base model (BM). A parametric analysis was carried out afterwards to test the influence of different geometries and rock mass strength properties in the failure mechanism and failed volume.  The source area of the Portillo Rock Avalanche was treated as a discontinuous rock mass and a Mohr-Coulomb elasto-plastic constitutive criterion was adopted to represent block deformation. Movement along discontinuities was modeled using an elasto-plastic joint area contact with Coulomb slip failure, based upon elastic stiffness, frictional, cohesive and tensile strength properties. The base below the slope was treated as an elastic, continuous and isotropic material to represent a strong, relatively rigid buttress against and underlying the zone of damaged rock. The base was subdivided to enable a better element size gradation between the elasto-plastic zones used to model the fractured nature of the zone of damaged rock and the linear-elastic base (Figure 4.4).  The boundary conditions were set for which the velocities on the left and right boundaries were fixed so that only movements in the y-direction were possible (i.e. rollers). Along the bottom of the model, the y-velocities were fixed (i.e. pins). Then the loading conditions were entered which included the in situ stress state. Stresses were initialized assuming a compressive regime (i.e. horizontal to vertical stress ratio of 2). This was based on the regional stress tensor obtained by the analysis of the focal mechanisms of a number of earthquakes recorded in the central Andes (Pardo et al. 2002; Barrientos et al. 2004).  To solve for the initial equilibrium state, the rock mass and discontinuity strength parameters were set to artificially high values to prevent yielding during the initial loading (Figure 4.5). History points were set near the crest of the slope, close to the free-face and at the toe to track the behaviour of the slope with time-stepping (see Figure 4.4 for location of points).   94  Figure 4.5. Horizontal (σxx) and vertical (σyy) principal stress contours for the initial state of the base model.  Once the initial equilibrium state was achieved, the rock mass strength properties were set to lower, more realistic values (Table 4.2 and Table 4.3). A staged analysis was then run until the ratio (R) of the maximum unbalanced force to the representative internal force reached 0.001%. The idea behind this was to simulate the progressive strength degradation mechanism that must have operated in the slope. Rock mass cohesion was lowered to simulate the loss of coherence in the stronger units at the toe of the slope through brittle fracture processes, while joint cohesion was reduced to represent the destruction of asperities and intact rock bridges between non-persistent joints.                95 Table 4.2. Mechanical properties of the rock mass used for the base model after the initial state was achieved. Source of information same as for Table 3.8. Rock Mass Property Domain 1 Domain 2 Base UCS (MPa) 100 60 250 GSI 55 40 70 Intact rock parameter: mi Density (Kg/m 3 ) 2700 2600 2700 Cohesion (MPa)  10 - 8 - 6 - 4  8* - 6 - 4 - 2 10 Friction angle (º) 30* 25 46 Tensile Strength:     T 0 (MPa) 0.20 0.05 2.00 Intact Young's Modulus: Ei (GPa) 40 20 50 Poisson's Ratio: ν 0.25 0.35 0.20 Rock Mass Deformation Modulus: Em (GPa) 12 3 36 Bulk Modulus:          K (GPa) 8 3 20 Shear Modulus:       G (GPa) 5 1 15 Constitutive Model Mohr - Coulomb      elasto-plastic Mohr - Coulomb elasto-plastic Elastic isotropic (*) Strength properties obtained from Kulhawy (1975). 11.3      96 Table 4.3. Mechanical properties and geometric parameters of discontinuities for the base model. Bedding JS1 Cross-joints JS2 Bedding JS3 Cross-joints JS2 Geometry Spacing: 30 m Dip: 65º - 50º Spacing: 40 m Dip: 35º Spacing: 20 m Dip: 25º* Spacing: 20 m Dip: 35º Joint friction angle (º) 30 30 25 25 Joint cohesion (kPa) 1000 - 10 - 0.1 - 0 50.0 0.0 0.0 Joint normal stiffness:       Kn (GPa/m) 10.0 10.0 5.0 5.0 Joint shear stiffness:       Ks (GPa/m) 1.0 1.0 1.0 1.0 Joint constitutive model (*) Difference in dip with respect to the value obtained from the stereonet analysis corresponds with the correction from true-dip to apparent-dip. Domain 2 Joint area contact with Coulomb slip failure Discontinuity Property Domain 1   The analysis was performed without considering pore-water pressures given the uncertainty on the prehistoric conditions and the general absence of springs or seeps at the toe of the slope. However, to partially account for the uncertainty in the amount and distribution of pore pressures in the prehistoric slope, a sensitivity analysis with varying joint shear strength (i.e. joint friction) was performed (see section 4.3.3.2). The rationale behind this was that by finding the limit equilibrium joint friction angle, the influence of water in reducing the effective strength is implicitly accounted for.    97 4.2. Base Model Results When comparing the results of state 2 with state 4 (Table 4.4), the models indicate that failure initiates when the rock mass cohesion of the slope and zone of damaged rock is set to 8 and 6 MPa respectively and joint cohesion approaches 10 kPa. When the rock mass and joint cohesion is reduced to the values corresponding to state 4, the amount of displacement in the x-direction increases and migrates upwards from the bottom of the zone of damaged rock to the upper slope leading to failure (Figure 4.6).  Table 4.4. Varying strength properties used in the staging process to develop the base model. High Strength STATE 1 Moderately High Strength STATE 2 Moderately Low Strength STATE 3 Low Strength STATE 4 Rock mass cohesion in D1* (MPa) 10 8 6 4 Rock mass cohesion in D2* (MPa) 8 6 4 2 Joint cohesion of JS1* in D1 (kPa) 1000 10 0.1 0 (*) D1: domain 1, D2: domain 2, JS1: bedding planes; sensu section 3.5.1    Figure 4.6. Evolution of horizontal displacements for damage states 2 and 4.   98 With yielding of the passive toe support, y-displacements for stage 4 show an increase downwards (parallel to bedding) in the upper slope (domain 1), with more bedding slabs failing and therefore an increased slide volume compared to that for stage 2 (Figure 4.7). At the base of the zone of damaged rock (domain 2) an outward movement of the toe of the slope is observed. This movement is interpreted as the response of the weaker material being squeezed between more competent rock units within the zone of compression induced at the toe of the slope.   Figure 4.7. Evolution of vertical displacements and total displacement vectors for damage states 2 and 4. Negative values indicate downwards displacements.  The bi-linear geometry of the slip surface that develops for both states is indicated by the displacement vectors which go parallel to the bedding planes in the upper slope and rotate out through continuum failure at the toe of the slope. For the two states, the greater displacements are concentrated mainly around the base of the zone of damaged rock and the hinge (see Figure 4.4-a for location) of the slope.  Plasticity indicators together with shear displacements are shown in Figure 4.8. The overloading of the weakened damaged zone leads to yielding mostly in shear of the rock mass at the toe of the slope indicating the initiation of the rupture surface. Tensile failure above domain 2 could be indicative of a staged failure, controlled in part by a “Prandtl”-like wedge between active and passive zones (e.g. as shown by Kvapil & Clews 1979 ).   99  Figure 4.8. Plasticity indicators and shear displacements. Arrows indicate directions of shear displacement.  Major shearing is concentrated at the base of domain 2 (i.e. passive zone), probably augmented by the active driving compressive forces exerted by the upper slope. Planar sliding begins at state 3, where the stepped pattern shown by the shear indicators along the persistent bedding planes in domain 1 may explain the blocky aspect of the rock mass observed in the upper part of the mapped sliding surface (Figure 4.9-a,b), which also coincide with the change in angle of the bedding planes. Here the cross joints (JS2) act as rear kinematic release surfaces which permit the slab motion.   Figure 4.9. a) Plasticity indicators and shear displacements for modelled state 3; b) Northwest facing view of the slope showing the presence of JS2 acting as release surfaces.  100 Based on the plasticity, shear, displacement, and velocity indicators, the predicted rupture surface extends deeper into the slope compared to the mapped one, therefore adding extra volume. The failed volume obtained from the base model exceeds the average in place field- calculated volume (see section 4.1.2) by approximately 20% (or ca. 12 million m3). The biggest difference between the predicted failure surface and the observed one is in domain 2, where the greatest model deformations occur.  The unstable part of the slope is captured in the magnified (30x) grid plot of Figure 4.10 which shows the en echelon array of the failed bedding planes which increase with the yielding of intact rock at the toe of the slope. The evolution of the bulge at the toe is constrained by the interaction between the normal stresses applied to the hangingwall and the outward movement of domain 2. The bulging zone is characterised by the clustering of tensile failure indicators, probably aided by the negligible to low confinement that exists close to the free-slope face (Figure 4.11). A similar situation, but to a lesser degree, is a cluster of tensile damage observed around the hinge of the slope which can be interpreted as a subsidiary Prandtl wedge. Comparable to the propagation of the shear failure, tensile damage initiates at the toe of the slope and propagates upwards to the hinge. With continuing time-stepping, a partial buckling mechanism around the hinge may be observed as there is some opening and separation (outward displacement) of the bedding planes.   Figure 4.10. Grid plot magnified 30 times that shows the bulge evolution from state 2 to state 4.  101  Figure 4.11. Plot of the principal stress tensors. Blue diamonds indicate principal stresses in tension.  Thus, based on these modelling results, development of the Portillo Rock Avalanche likely initiated through yielding and shear failure of the rock mass at the toe of the slope. Upward propagation of tensile damage and brittle fracturing, aided by the cross cutting joints, would then enable kinematic release at the toe, which in the absence of a cross-cutting daylighting discontinuity dipping out of the slope face would not be otherwise kinematically possible. The failure of the lower portion of the slope, which provides support for the upper slope, subsequently allows the rest of the slope to fail through sliding and shearing along the bedding planes, aided in part by the existence of the cross-joints (Welkner et al. 2007a).   102 4.3. Parametric Analysis Given the uncertainty in several of the input parameters, a parametric study was carried out to gain confidence in the assumptions made. The range of the varied model conditions (including both geometrical factors and material properties) was taken from the table of reasonable values (upper and lower bound) obtained from the scaled lab test data collection and from consulting the literature (see Chapter 3, Table 3.8 and Table 3.9), and constrained by the geological observations made in the field. Every model was run until it reached R=0.001% and then compared to the base model. The influence of the varying parameters in the failure mechanism and failed volume was investigated following the flow diagram shown in Figure 4.12. Only significant changes in the model results are discussed below.   Figure 4.12. Flow diagram showing the sequence of trial models.  The sequence followed for the development of the trial models was taken from the modelling methodology in which the geometry was set first, then the constitutive models used for the deformable blocks, and then the input of the mechanical properties. For every trial model, changes were made one at a time while the rest of the parameters were kept constant (equal to those of the base model – State 4). For example, when the geometry was changed to that of a homocline, only the topography was modified and the remaining parameters were kept the same as the ones presented in Table 4.4 for damage state 4.  103 4.3.1. Variation in geometry Different geometries were tested based on several variants judged possible due to the complex geological history of the Caracoles Range (see Chapter 3 for more details). The first two trial models (sections 4.3.1.1 & 4.3.1.2) gave insights into the influence of the topography on the stress distribution and therefore its influence on the failure mechanism, while the other studied the influence of changing block size (section 4.3.1.3). 4.3.1.1. Homocline geometry The Caracoles Range is a concave slope in cross section which was simplified to that of a homocline following Starfield & Cundall (1988) recommendations to start with a simple conceptual model in which the important mechanisms are represented and then add complexity as required. For this case, the absence of daylighting structures was kept as the primarily constraint in the geometry.  Contour plots of horizontal and vertical displacements, displacement vectors and joint shear indicators were used to examine the behaviour of the failed mass for the homocline geometry (Figure 4.13-a). The distribution of the horizontal displacements is concentrated more along the zone of damage rock compared to the one observed in the base model where the driving forces are partly distributed within the upper slope. The resultant effect was that the lower boundary of the proposed slide surface did not match the mapped one which corresponds with the current profile of the slope. On the contrary, the distribution of the vertical displacements was similar in both the base and homocline models (Figure 4.13-b).  Plasticity indicators (Figure 4.13-c) show the development of elements that yield in shear in domain 2, crowned by elements that yield in tension; few tensile indicators are observed close to the crest of the slope in domain 1. This distribution together with the presence of joint shear displacement indicators along the bedding planes in domains 1 and 2 are evidence of yielding of intact rock at the toe of the slope and sliding through the bedding planes. When analyzing the grid plot magnified 30 times, it is clear that the bulging at the toe is developing in the same way as in the base model. Differences are encountered in the slip of the layers along the bedding planes which is less developed, as evidenced by Figure 4.13- c, and the absence in the upper slope of the additional yielding at the hinge point.  104  Figure 4.13. Homocline geometry. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators; d) Grid plot magnified 30 times.  In conclusion, the homocline geometry was able to reproduce the proposed failure mechanism, but there was a significant difference in the distribution of the horizontal displacements at the toe of the slope. A possible explanation for the differences between the geometries could be the development of higher normal stresses for the homocline model in domain 2 compared to those in the concave model. These stresses limit the movements and reduce the displacements in the upper slope (in the normal direction, parallel to bedding). This can be seen indirectly in Figure 4.14 which compares the vertical stresses for both models.   105  Figure 4.14. Zoomed in plots of vertical stresses (σyy) in domain 2 for the concave model and the homocline model.  4.3.1.2. Toe-buttress geometry A cross-section representative of a toe-buttress geometry was tested based on a series of E- W profiles across Inca Lake which showed the existence of a protruding ridge along the toe of the Caracoles Range (Figure 4.15). These cross-sections revealed that for the current topography the protrusion was missing at the location of the rock avalanche, indicating that it may have existed prior to failure and was part of the slide mass.   106  Figure 4.15. E-W cross sections across the Inca Lake and location map showing a protruding ridge along the toe of the Caracoles Range.  The trial model is presented in Figure 4.16 and the results from it show that the addition of the buttress at the toe of the slope provides a more stable configuration. The larger horizontal displacements are observed close to the bottom of domain 2 resulting in an undefined base for the proposed slip surface, while the vertical displacements are notable only below the hinge.  When the plasticity and shear indicators are analyzed, the elements that fail in shear and tension in domain 2 together with the joint shear displacements along the bedding planes of domain 1 are not completely developed so as to create a continuous shear surface comparable to the actual sliding surface. Moreover, in domain 2, close to the surface, there are no elements that yield in shear as in the other two models (base and homocline). With more time-stepping, the shear surface that develops penetrates deeper into the slope.   107 Displacements at the base of domain 2 (Figure 4.16-c) are the result of more joints that fail in shear, which is also a by-product of the extra load exerted in the passive zone. Considering that joint shear and normal stiffness values (Ks and Kn respectively) are not constant but raise with increasing normal stress, enabling the joint to carry excess load before shearing (Bhasin & Kaynia 2004), these displacements may also been artificially high due to the combination of the greater normal stress (addition of the ridge) and the simplifying assumption of a constant Ks and Kn, using the same values as those for the base model - State 4.  The presence of a ridge at the toe of the slope helps limit yielding of elements in domain 2 by means of the increase in load exerted by the added material placed close to the passive zone. This geometrical configuration is similar to what is commonly used as a stabilization measure in most slope stability earthworks which consists of placing fill at the toe of the slope (Hutchinson 1984).   Figure 4.16. Toe-buttress geometry. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators.  108 4.3.1.3. Block size reduction The following trial model tested the influence of the block size in the failure mechanism and failed volume. In domain 1 the JS1 spacing (bedding planes) was reduced to 15 m and the JS2 spacing was reduced to 20 m, following the results presented in Table 3.3. In domain 2 the spacing of all discontinuities was reduced to 10 m.  On the other hand, when analyzing the plasticity indicators, yielding in shear of intact rock in domain 2 bounded on top by tensile failure, compares well with the damage zone observed in the base model. Nevertheless, a difference is encountered in the joint shear displacements which seem to be inhibited due to the smaller block size in the upper slope, considering that they are well developed below the hinge (Figure 4.17-c).  The results presented in Figure 4.17 show that the distribution of the horizontal and vertical displacements is comparable to the ones obtained for the base model however, the slide surface tends to be shallower resulting in a smaller failed volume. The situation is clearly shown by Figure 4.17-b, where contouring every 75 cm needs to be used to observe displacements.   109  Figure 4.17. Reduced block size. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators.        110 4.3.2. Variations in block constitutive models and material behaviour 4.3.2.1. Elasto-plastic base In the base model the far-field base below the slope was modeled as an elastic, continuous and isotropic material to represent a strong relatively rigid buttress underlying domain 2. A trial model was developed with the base modeled using a Mohr-Coulomb elasto-plastic constitutive criterion. The strength parameters used for the base of the slope are presented in Table 4.5, and correspond with the higher strength values to coincide with the assumption of stronger rock below the fault. The rest of the strength parameters correspond to those representing State 4 (Table 4.4).  Table 4.5. Mechanical parameters of the far-field base modeled as a continuum with a Mohr-Coulomb elasto-plastic constitutive criterion. Parameter Value UCS (MPa) 250 Young's Modulus (GPa) 70 Poisson's ratio 0.2 Density (kg/m 3 ) 2700 Rock mass cohesion (MPa) 10 Rock mass friction (º) 46 Rock mass tensile strength (MPa) 2   Analysis of the plasticity indicators show shearing along the right boundary of the model (Figure 4.18-c). However, these and the high horizontal stresses (σxx) that occur here (Figure 4.18-d) are only boundary effects that have no impact on the distribution of stresses  111 or yield in the region of interest. This situation only arises when the constitutive model used for the base is changed to a non-linear Mohr-Coulomb material.  The trial model results are basically equivalent to the base model results, with the horizontal and vertical displacements showing a fairly good agreement with the actual slide surface (Figure 4.18-a,b).   Figure 4.18. Far-field base modeled with a Mohr-Coulomb constitutive criterion. a) Horizontal displacements for state 4; b)Vertical displacements for state 4; c) Plasticity and shear displacement indicators; d) Horizontal principal stress contours (σxx).  112 4.3.2.2. Zone of damaged rock as an equivalent continuum For this series, the densely fractured and weak domain 2 was modeled as an equivalent continuum for which the mechanical properties were invariant with respect to the base model – State 4 (Table 4.4). This effectively resulted in a higher total rock mass strength for the zone of damaged rock, relative to the other models due to the explicit absence of the discontinuities in the equivalent continuum.  The results are shown in Figure 4.19, and are basically equivalent to the base model results. The horizontal and vertical displacements show a good agreement with the proposed slide surface. On the other hand, the plasticity indicators show a smaller amount of elements yielding in shear close to the free-face at domain 2 resulting in much better defined spoon- shape geometry for the toe of the sliding surface. This situation probably is in response to the absence of small-sized blocks which are prone to yield due to the loading exerted from the upper slope. The elements that yield in tension are located in the same place compared to those of the base model. Only additional yield is observed at the crest of the slope probably representing the tension cracks that commonly form at the crest of failing slopes. Importantly, representing the zone of damaged rock as an equivalent continuum significantly reduced the model runtime.  Considering that this model showed to be capable of effectively reproduce the failure mechanism and estimated failed volume, and accounting on the performance of the slope with time-stepping, the best representation for domain 2 is that of an equivalent continuum material.   113  Figure 4.19. Zone of damaged rock as an equivalent continuum. a) Horizontal displacements for state 4; b) Vertical displacements for state 4; c) Plasticity and shear displacement indicators.        114 4.3.3. Variation of mechanical properties 4.3.3.1. Rock Mass Strength A set of trial models was performed in which the rock mass strength properties were varied from those used in the base model – State 4, following the table of reasonable values presented in Chapter 3 (Table 3.8). While for most of the parameters the base model included the lowermost strength values, this set of trial models used the rest of the range in order to constrain the values for which the proposed failure mechanism and estimated volume were possible. Given that a staged procedure with decreasing cohesion (rock mass and discontinuity) was used for the base model, only the block friction and tension were varied in this analysis.  The procedure involved varying one property at a time rather than simultaneously to test the influence of the changing parameter instead of determining the effects of a weaker or stronger rock mass. When one parameter was changed, the rest were kept constant and equal to the values corresponding with State 4 of the base model (Table 4.6).  Table 4.6. Range of rock mass strength properties tested for the parametric study. D1 D2 Friction (º) 35 - 46 25 - 30 30 25 Tension (MPa)  1 - 2 0.1 0.2 0.05 (*) Friction & tension used in Base Model - State 4 BM-State 4* Strength Properties Domain 1 Domain 2   The first set of trial models show the performance of the slope under a varying friction angle (Figure 4.20). The plot of plasticity indicators together with the joint shear displacements show that even with the use of the highest friction angles the toe fails in shear. However, the history plot of the horizontal displacements shows they eventually reach a steady state. When these two elements are considered together, it is likely that the model self-stabilizes rather than catastrophically fails.  115 The other set of models test the influence of the tensile strength in the stability and failed volume of the slope. As was expected, even with the use of the highest values of tensile strength, failure of the slope still occurs (Figure 4.21). The reason was because the most important mode of failure of the slope was shear of the toe rather than tensile failure. Thus, the variation of tensile strength was not going to produce major changes in the modelled failure mechanism or failed volume. In contrast to those models where higher values resulted in stable slope conditions, as indicated by their history plots, in this model an accelerating condition was observed.   Figure 4.20. Trial models with varying friction angle. a) Plasticity and shear displacement indicators for upper bound friction angle values; b) History plot of horizontal displacements for upper bound values of friction angle.    116  Figure 4.21. Trial model with varying tension. a) History plot of horizontal displacements for upper bound values of tension; b) History plot of vertical displacements for upper bound values of tension.  In summary, the lack of sensitivity of the proposed failure mechanism to the varying block properties was taken as an element of confidence in the mechanical properties used as it helped reduce some of the uncertainty inherent in back-analyzing a prehistoric rockslide. 4.3.3.2. Joint friction The effective joint friction angle incorporates the influence of pore water pressure and surface roughness. The net effect of pore pressure is to allow the joint to shear at a lower angle of applied total stress. For the back-analysis these values were completely unknown therefore a series of trial models were developed to test the influence of the joint friction in the failure mechanism and failure mode of the slope (Table 4.7).  Table 4.7. Joint friction range tested for the parametric study. D1 D2 Joint Friction (º) 30 - 50 20 - 30 30 25 Discontinuity Property  Domain 1          (upper slope) Domain 2 (zone of damaged rock) BM-State 4    117 Figure 4.22 shows the results of the modeling with varying joint friction angle in both domains. Plots of plasticity and joint shear indicators illustrate that with high joint friction values no important shear displacement along bedding planes was observed compared to lower values, but more elements yield in shear in domain 2 (particularly in its lower portion). This behaviour could be explained by the response of the zone of damaged rock to deformation which is triggered by the shearing along discontinuities. The higher the joint friction, the lower the shearing along the discontinuities so more elements representing the intact rock must yield in shear to cope with the strain imposed.   Figure 4.22. Trial models with varying joint friction angle. a) Plasticity and shear indicators for lower, average, and high frictional strength values; b) Shear strain plots for lower, average and high frictional strength values.  When the shear strain plots are analyzed, deformable blocks show that with progressive increase of frictional strength the zone with higher shear strain migrates from the top of domain 2 downwards to the bottom of it. At low frictional strength, the blocks at the top of domain 2 are deformed triggered by the high rate of shearing along the bedding planes in domain 1 which are the ones that transfer the active driving forces from the upper slope to  118 the passive blocks in the toe. With increasing frictional strength, shearing along the bedding planes is no longer important so deformation was transferred to the base of the zone of damaged rock, along the cross-joints, where shear is most effective due to the alignment of these structures with the local induced stresses. This situation is also in direct relation with the localization of the elements yielding through the continuum, as explained in the paragraph above.  The range of values used did not change the failure mechanism. However, when jfricFZ was set to 20º (jfricFZ = 25º in the BM) shear in the cross-joints was highly exaggerated and when jfricSLOPE = 45º was used (jfricSLOPE = 30º in the BM) shear along the bedding planes was negligible. The latter situation is not realistic given the geometry of the slope and the most acceptable kinematic solution for motion.  Thus, because the joint friction angles considered here are effective values, it is most likely that if an analysis was performed that explicitly included water pressures, then the shear strength properties would be in the order of the upper bound values, closer to the values obtained in the lab testing. However, the failure mechanism would not change and therefore the insights and understanding gained can be transferred to the forward analysis of the present day stability state. 4.4. Interpretation of Results In interpreting the results, consideration must be given to the explicit time-stepping solution procedure employed by UDEC and the diagnostic indicators available (Itasca 2004). For the modeling of the Portillo Rock Avalanche the primary indicators used included the evolution of the unbalanced forces, gridpoint velocities, plasticity indicators and the displacement histories at the crest and toe of the slope.  The use of discontinuum numerical modeling techniques (UDEC, Itasca 2004) applied in the analysis of the Portillo Rock Avalanche helped support the field observations with respect to the most probable failure mechanism. Shear displacements and plasticity indicators show a mechanism that involves toe failure through brittle fracture processes and rock mass yielding. Loading and shearing of the bedded rock at the toe of the slope, enabled translational kinematic release of the upper beds along bedding up to the crest of the slope  119 (Welkner et al. 2007a). This failure mode was accommodated and highly controlled by the adverse geometry of the slope and bedding.  This failure mechanism was common to all models, independent of the strength parameters, geometry and constitutive models used. This provided a certain degree of confidence with respect to model uncertainty and focus was shifted to parameter uncertainty and the model that best reproduced the proposed failure surface along the Caracoles Range.  Table 4.8 contains the results of all the trial models and summarizes the effect of the varying parameters in the proposed failure mechanism and failed volume. It was concluded that the geometrical variations did not dramatically influence the proposed manner in which the slope failed, however, in most of the cases the failed volume was not well constrained. By means of the parametric analysis, the results showed that the best model in terms of estimated volume and mapped failure surface was the one having domain 2 represented as a continuum, which also had the added benefit of being the most computationally efficient. It should be noted though, that if domain 2 is represented as a continuum then problems arise if its required to include a water table because of the need of a network of discontinuities to let the water flow through the model (given the coupled hydro-mechanical formulation used in UDEC).                 120 Table 4.8. Qualitative estimation on the influence of the varying parameters in the proposed failure mechanism and volume failed. FAILURE MECHANISM FAILED VOLUME OBSERVATIONS Homocline Same as proposed Not clearly reproduced Distribution of horizontal displacements did not clearly outlined the base of the slip surface. Toe-buttress Not fully developed Not reproduced Plasticity indicators were not able to produce a linear trend coinciding with the proposed slip surface. Reduced block size Same as proposed Less than proposed Shallower sliding surface Base as Mohr- Coulomb Same as proposed Fairly reproduced Some boundary effects. Domain 2 as equivalent continuum Same as proposed Clearly reproduced Best way to represent the zone of damaged rock. Tension Same as proposed Same as proposed No effect in the model because tensile failure was not a predominant failure mode. Failure mechanism reproduced within range: D1: 30° – 40° D2: 25° – 30° Failure mechanism reproduced within range: D1: 30° – 37° D2: 25° – 30° Fairly reproduced within new range Fairly reproduced within new range Discontinuity strength Joint friction angle Same as proposed within new range Friction Same as proposed within new range Geometry Constitutive model Rock mass strength   121 5. RUNOUT ANALYSIS OF THE PREHISTORIC ROCK SLOPE FAILURE EVENTS 5.1. Introduction This chapter presents the pre-failure terrain reconstruction and 3-D dynamic flow modelling used to characterize the runout path and volume constraints of the prehistoric rock slope failure events recognized in the study area. The main objective of this chapter was to understand the flow behaviour and constrain the volumes involved in order to unravel the complex history and sequence of rock avalanche events in the valley. Some of the results of this chapter were previously presented by Welkner et al. (2007b).  The importance of this back-analysis is that it provides insights to the mobility (runout) of the flow, which will be used in later stages of the study to predict the potential impact of any future major rockslides on the transportation route between Chile and Argentina (Santiago – Mendoza international corridor) and on an important ski resort.  Continuum dynamic analyses were carried out using the numerical code DAN3D (McDougall & Hungr 2004), where combinations of rheologies were tested and ranked based on their ability to reproduce the mapped distribution of the debris, and further tested against temporal constraints corresponding with the sequence of events established by cosmogenic nuclide dating. Specifically, the back-analysis of the post-failure motion of the pre-historic slope failure events tried to answer the following questions:     Is the origin of the Inca Lake caused by the obstruction of the valley during the PRA event (i.e. rockslide dam)?    Do the simulations of the pre-historic events agree with the chronology of events established by the cosmogenic nuclide dating?    Does the Holocene event have a multiple (east and west slope) or single (east slope) initiation zone?  122 5.2. Model Development The key inputs required for the runout modeling were the pre-failure (in place) volume estimate and the pre-failure topography. The model results were highly sensitive to the location of the pre-failure volume and the amount of detail considered in the estimated pre- failure surface. For this reason, the identification of the detached areas was carefully surveyed utilizing air photographs at different scales and ground-based photos taken of the observed detachment zones. For the reconstruction of the pre-failure topography, a digital terrain model (DTM) of the current topography was used as a base map. This map was constructed by interpolating the only digital map contours available for the zone at 1:50,000 scale, with 50 m contours. The detail of the terrain model was limited to the amount of detail existing in the original data (Figure 5.1).   Figure 5.1. Digital terrain model (DTM) constructed through the interpolation of a 1:50,000 digital contour map with 50 m contours (PSAD 56, zone 19S).   123 5.2.1. Volume calculations The volume of debris left behind by the Holocene rock slump (S2 in Figure 5.2) was estimated by comparing the geometry of the lobe that crops out at the foot of the western slope to that of a truncated pyramid (i.e. frustum) (Figure 5.3). As before for the Portillo Rock Avalanche, a bulking factor of 25% was adopted to account for the break-up and transport of the detached mass (Hungr & Evans 2004). On account of the short runout (Fahrböschung = 35º) the volume contribution due to entrainment was assumed to be negligible. The resulting pre-failure (in place) volume of the rock slump was calculated to be approximately 16x106 m3 with an estimated average thickness of 30 m ± 5 m (Table 5.1). The latter was derived by means of the surface estimate of the detached zone of approximately 50x104 m2 calculated from the DTM using Surfer 8.0. Details on the estimated failed volume of the Portillo Rock Avalanche were already presented in section 4.1.2.   Figure 5.2. Aerial photograph (1:80,000) showing the distribution of the mass wasting events and approximate ages.  124  Figure 5.3. Picture comparing the shape of the western lobe with that of a wedge, and indication of the values used for the volume calculation.  A summary of the geometrical parameters of the Upper Pleistocene Portillo Rock Avalanche (PRA) and the Holocene rock slump (S2) is given in Table 5.1. All values are approximated with an error estimate no greater than 30% based on a consideration of the possible limits of variations in the calculation of the exposed areas. Only the upper bound value (see Figure 4.2) was considered in the error estimates for the volume calculations of the PRA due to its better agreement with field observations; the lower bound value was ruled out due to inconsistencies with the field observations (e.g. a predicted thickness of the deposit of only 15 m).        125 Table 5.1. Geometrical parameters of the prehistoric mass movements. PRA Rock Slump Detachment zone Scarp surface (m2) 88x104 50x104 Thickness (m) 60-65 30* Detached volume (m3) 51x106 14x106 Vertical drop (m) 900 1200 Slope angle (º) 50-65 40-45 Debris Area (m2) 2.1x106 4x105 Volume (m3) 68x106 18x106 Fahrböschung (º) 14 35 * Average thickness  5.2.2. Pre-failure topography reconstruction The detachment surfaces and debris were identified during field mapping and delineated on aerial photographs. After constraining the volumes, the procedure used to best estimate the pre-failure topography of both prehistoric mass movements was the same, and was conducted sequentially following the chronology of events.  The source of the rock mass events was reconstructed using a DTM constructed from the current 1:50,000 scale topography. The previously calculated volume of the detached mass was added to the source area following the same principles in which the ‘Sloping Local Base Level (SLBL)’ concept is based (Jaboyedoff et al. 2004; Jaboyedoff et al. 2005). This method consists of defining the maximum volume that can be affected by gravitational movements, which is found by searching for a surface that is the locus of all the points that correspond with the lowest level of erosion also known as base level. The Sloping Local Base Level is  126 then defined as the lowest position of a sliding surface and corresponds with the surface joining all local minima of altitude of a valley, such as rivers, which are considered fixed during the processing of the SLBL (Jaboyedoff et al. 2004) (Figure 5.4-a).  It is clear that this GIS methodology was developed to identify potential, large, unstable rock/soil slopes, i.e. slopes that have not failed yet and show no or minor evidence of instability. For the case of the prehistoric mass movements, the SLBL criterion was used but in a reverse mode, i.e. the debris of failed mass in the valley was compacted based on a bulking factor of 25% and then added to the area outlined by the scarp. The mapped sliding surface (ie. the current slope profile) was approximated as corresponding with the SLBL, and the addition of the in-place volume was derived considering the intersection between the outline of the scarp and the map contours as the fixed points (Figure 5.4-b; see Figure 5.6 for location in the study area).   Figure 5.4. a) Illustration that shows the Sloping Local Base Level (SLBL) concept which corresponds with the lowest position of a sliding surface (modified from Jaboyedoff, 2005); b) Principles of the SLBL concept applied to the reconstruction of the source area of the PRA.  To reconstruct the path of the rock mass events, the DTM plus a series of E-W and N-S cross sections, across and along the valley respectively, was used to better constrain the thickness of the deposits (Figure 5.5). The reconstruction of the topography before the Portillo Rock Avalanche event followed the same methodology and principles as the ones used in the source areas, but this time instead of adding volume, the debris was subtracted  127 from the valley floor. In order to follow the chronology of events, the volume of the rock slump was subtracted from the depositional zone and Inca Lake was removed by redrawing the N-S valley based on a series of E-W cross sections (see Figure 5.5). The methodology followed involved extending the valley walls where the rock outcrops to the lake and the current valley bottom. Also, the thalweg of the probable pre-existing N-S river was projected from N to S beneath the lake (Figure 5.6).   Figure 5.5. DTM and cross sections across and along the valley used to reconstruct the pre-failure topography. NE-SW-2 cross section shows the profile of the path of the PRA (see Figure 5.7 for details).  128  Figure 5.6. Pre-failure (top) and post-failure (bottom) topography showing the valley reconstruction. Note around the edges of the PRA’s scarp the addition of material (V+), and in the location where the rock slump debris should crop out, the zone where material was subtracted (V-). In the pre-failure topography, drainage has been kept for guidance purposes.     129 A simple back-calculation of the depth of the deposit based on its volume and exposed area yielded a constant thickness of approximately 40 m ± 7 m. On the other hand, the profile of the path of the Portillo Rock Avalanche estimated from the 1:50,000 topographic map of the area through the projection of the valley side beneath the debris yields a thickness of 60 m ± 10 m for the deposit. The difference in thickness obtained from the two methods was attributed to pre-existing glacio-fluvial sediments in the valley floor (Figure 5.7).   Figure 5.7. NE-SW cross section parallel to the path followed by the Portillo Rock Avalanche showing the chronology of events and its related surfaces (see Figure 5.5 for location of profile).      130 5.2.3. Rheological relationships The main constraints for the dynamic analysis were the volume, position of the failed mass, and the actual distribution of the debris in the valley as mapped during the field investigation. Based on the type of landslide and geological materials involved, preliminary analyses were tested varying the basal rheological relationships that can be implemented in DAN3D (McDougall & Hungr 2004). The selected rheologies with their respective equations were:     Frictional:  (1 )tan ;  where  tan (1 )tanu b ur rτ σ φ φ φ= − − = −  [5.1]     Voellmy:  2 ;   where   tanx b gf fρ ντ σ φξ −⎛ ⎞= − + =⎜ ⎟⎝ ⎠  [5.2]  The frictional Equation 5.1 shows that the shear strength component represented by a frictional basal resistance is proportional to the effective bed-normal stress at the base. Equation 5.1 can be simplified to a one-independent variable by replacing the dynamic basal friction angle (φ) by a bulk basal friction angle (φb) with constant pore pressure ratio (ru). In this case, the loading response of the overridden material lies between purely drained, where little entrainment is expected, and undrained, where failure and mobilization of the bed material is expected due to excess pore water pressures (McDougall 2006). Because ru is kept constant, Equation 5.1 still has a frictional character even though it represents material friction and pore-fluid pressures (Hungr et al. 2005).  The Voellmy rheology is a two-parameter frictional (f) - turbulent (ξ) resistance relationship. In Equation 5.2, the first component is equivalent to Equation 5.1, where f = tanφb, and the second component accounts implicitly for all possible sources of velocity-dependent resistance (McDougall 2006). As a result of the inclusion of the second term, the use of a Voellmy basal rheology in landslide simulations predicts lower maximum velocities for a  131 given overall displacement compared to a constant frictional rheology and more uniform distribution of debris in the deposition area (Hungr et al. 2005). 5.3. Model Results 5.3.1. The Portillo Rock Avalanche In Chapter 4, it was demonstrated that the most likely mechanism of failure of the Portillo Rock Avalanche was shearing of intact rock at the toe of the slope followed by sliding along the bedding planes. Once the failed mass reached the valley floor, fragmentation and entrainment of saturated glacio-fluvial material corresponding to the valley fill (see Figure 5.7) would have become important. The mass then continued flowing down the valley until it stopped 7 km away from the source.  It is likely that a major impact of the Portillo Rock Avalanche on the prehistoric landscape was the obstruction of the drainage of the N-S valley and the formation of Inca Lake located at the toe of the Caracoles Range (Figure 5.2). Field observations of the deposit in the upper segment of the valley strongly support this hypothesis.  To account for the mobility of the rock avalanche due to the character of the superficial material encountered, a combined basal rheology along the path was tested and proved to be the one that best reproduced the physical characteristics of the PRA event. The two phases of motion simulated involved using a frictional basal rheology in the proximal path and a Voellmy basal rheology in the distal path; these are in accord with the probable change in behaviour of the moving mass in terms of response to the shear resistance encountered at the basal interface (Figure 5.8). The parameters listed in Table 5.2 were systematically adjusted until they matched with the field constraints aforementioned.   132  Figure 5.8. Cross section showing the change in basal rheologies according to the change in basal material strength along the path.  Table 5.2. Basal rheology and parameter values used in the Portillo Rock Avalanche dynamic analysis. Frictional Voellmy Rheology elevation (m) 4050 - 2800 2800-2050 Bulk basal friction angle (º): φ b 30 N/A Pore-fluid pressure coefficient: r u 0.18 N/A Internal friction angle (º): φ i 35 N/A Basal friction coefficient: f N/A 0.1 Turbulence parameter (m/s 2 ): ξ N/A 500   133 In the proximal path, under a basal frictional rheology, the selected basal bulk friction angle was set to φb=30º, which is in the range of angles proposed by Hsu (1975) for rock slides. In the distal path, and under a basal Voellmy rheology, the basal friction coefficient and the turbulence parameter were set to f=0.1 (equivalent φb = 5.7º) and ξ=500 m/s2 respectively. The Voellmy parameters were adjusted so as to obtain the best fit in terms of the total displacement, and correspond with the values reported by Hungr & Evans (1996) after back- analysis of 23 rock avalanches.  It is assumed that some entrainment must have occurred given the runout of the event over saturated post glacial materials (e.g. Abele 1997) and the characteristics of the deposit in its distal portion (i.e. hummocky terrain; see Figure 5.2). Nevertheless, the entrained volume was likely small in comparison to the enormous volume of rock produced by the rockslide. So, no change of volume was considered due to the high uncertainty involved in the estimation of the amount and location of the entrained material.  The results of the PRA dynamic analysis are presented in Figure 5.9, and were seen to provide good agreement with the prehistoric event in a number of ways. First, the model provides a close match with the general extent and distribution of the deposit, especially with respect to the reach of the leading edge of the debris. It is interesting to note that there are three distinct zones where deposition of the mapped and modelled debris occurs that fall between the upper valley and the area where the flow mass came to rest: Z1, Z2 & Z3 in Figure 5.10. This is not an arbitrary situation, and its explanation is based on the pre-failure topography reconstruction which was done using an average valley slope of 5.7º, thus “conditioning” these zones to make them more favourable for deposition. The average simulated deposit depth is approximately 45 m, which falls within the range of the thickness obtained from back-calculating the areal exposure and its associated volume, as well as the range of the thickness estimated in the mapped cross-sections.    134  Figure 5.9. Calibrated simulation of the Portillo Rock Avalanche event at 50 s intervals. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. Simulated and estimated trimlines are shown in dashed and dotted lines, respectively.   135  Figure 5.10. Terrain slope map showing major zones of mapped deposition. Contouring of slope angle in the range 0° – 7° has been highlighted.  In addition, the runout analysis supports the hypothesis of Inca Lake having been formed by a rockslide dam, as a significant amount of the modelled runout material is deposited at the south end of where the lake currently sits (Figure 5.9). The interpretations based on the cosmogenic nuclide dating are also indirectly supported as there is no indication of material being deposited in the area of the S2 lobe, confirming a different source for the Holocene age-dated deposit.  Given the age of the PRA event, no solid evidence of a trimline could be discerned, limiting modelled constraints of the travel path to that of the outline of the deposit. As such, when comparing the simulated trimline with the outline of the deposits (i.e. estimated trimline) small but noticeable differences are observed around the margins. These differences may also be related to redistribution of the original deposits, specifically at the margins, due to reworking of the landscape mainly attributed to erosion once drainage was resumed. Other minor  136 differences may be attributed to the omission of small-scale pre-slide terrain features in the reconstruction of the prehistoric path.  A plot of maximum simulated flow velocities is presented in Figure 5.11. Even though no field evidence exists to constrain these values (e.g. superelevation-based estimates), according to the literature available the velocities are in the typical range for catastrophic landslides (Hungr et al. 2001). The highest velocities, up to 85 m/s, are located in the initiation zone, where the mass started its descent down the slope. A second peak of velocities is shown, between 60 m/s and 70 m/s, around the area where it is assumed that entrainment started thus coinciding with the change from a basal frictional to a Voellmy rheology for the simulation.   Figure 5.11. Plot of maximum simulated flow velocities. Maximum velocity contours are at 10 m/s and elevation contours are at 100 m intervals.   137 Similar results as the ones presented in Figure 5.9 were obtained by modeling the Portillo Rock Avalanche with a constant Voellmy model, with a tanφb = 0.1 and ξ = 500 m/s2, showing slight differences in the distribution of the deposits along the path around the Z2 deposition zone. This result is supported by Hungr & Evans (1996) which reported that for most (70%) of the 23 cases analyzed; a constant basal Voellmy rheology was needed to fully represent the dynamics of the rock avalanches in which entrainment was believed to have occurred.  The use of a constant frictional rheology was also attempted, and after several calibrations a basal friction angle of 20º with ru coefficient of 0.48 showed to be the parameters which more closely conformed to the current deposit. Nevertheless, these results did not represent properly the longitudinal debris distribution and morphology of the lobes (Figure 5.12). The maximum simulated flow velocities were in the range of 110 – 120 m/s, also recorded in the initiation zone, with velocities at the toe of the slope, ranging between 90 - 100 m/s. As previously noted, Hungr et al. (2005) state that high velocities are a well-known characteristic of the frictional model.   Figure 5.12. Calibrated simulation of the Portillo Rock Avalanche event utilising a constant basal frictional rheology with a basal friction angle of 20°. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. Estimated trimline is shown in dotted line.  138 5.3.2. The Holocene rock slump lobe In the upper parts of the valley significant reworking of the material has occurred as part of the operations related to a ski resort (see Figure 5.2 for location). This alteration of the terrain makes it difficult to establish which of the following two scenarios better represents the initiation of the Holocene event (see Figure 5.2):     The two deposits (S1 & S2) are derived entirely from the Caracoles Range or,    S1 & S2 correspond to two simultaneous events derived from the east (Caracoles Range; S1 deposit) and west sides (S2 deposit) of the valley.  For the first case, the leading edge of the rockslide would have had to travel approximately 2 km west across the valley overriding the previously deposited Portillo rockslide debris, run up the base of the west slope, and land at its toe as it fell back forming a fall-back ridge. In this case, the lack in continuity between the eastern and western deposits would then need to be explained by the reworking of the terrain.  For the second scenario, the source of the lobate deposit found at the foot of the western slope could be explained as a dry rock slump/collapse whose short runout could be attributed to unfavourable geometrical conditions for motion created by the dip of the beds on that side of the valley (i.e. into the slope). For this situation, given the temporal and spatial proximity of both S1 and S2 deposits it is likely that their failure could have been triggered by an earthquake.  A series of dynamic runout analyses were performed using the code DAN3D (McDougall & Hungr 2004), which helped constrain the rheology of the flow and evaluate the likelihood of the two scenarios. There was no need to conduct a special analysis for the first scenario as it was indirectly ruled out by the results of the modelling of the PRA and taking advantage of the starting zone on the eastern slope. As explained in section 5.3.1 (Figure 5.12), the modelling of the PRA showed that not enough material was left by the flow at the toe of the western slope comparable to the lobe-shaped deposit. Moreover, if the greater volume of the PRA compared to that of the Holocene rockslide is accounted for, the smaller momentum of the latter event together with the highly frictional characteristics that should have had its path  139 (overriding the PRA debris) makes the common origin scenario even more difficult to be accepted.  The methodology used for the dynamic analysis of the Holocene rock slump was the same as that presented for the Portillo Rock Avalanche. To reconstruct the path, a volume equivalent to 14 million m3 was removed from the base of the west slope (see Table 5.1), and the volume of debris corresponding to the earlier Portillo Rock Avalanche was left in place, following the chronology of events (see Figure 5.6).  Considering the failure kinematics, the distribution/morphology of the deposits and the short runout, a constant basal frictional rheology was used along the path to develop the dynamic analysis, utilising as input parameters an internal friction angle and a bulk basal friction angle with its associated pore pressure ratio (ru) to account for any pore-fluid pressure along the path. Based on the same argument of a short runout, no entrainment was considered. A series of model runs were performed varying the bulk basal friction angle (φb) and its related pore-pressure ratio (ru) using Equation 5.1 (Figure 5.13).   140  Figure 5.13. Trial runs for the West lobe rock slump with varying bulk basal friction angle. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. Estimated trimline is shown in dotted red line and simulated trimline in dashed-blue line.  The overall motion of the mass was best represented when the bulk basal friction angle was close to the internal friction angle, ranging between 33° and 35°, resulting in an almost dry- frictional rheology (Table 5.3). Such characteristics are typical of small rockslides, or those that fail in a sequence of partial detachments (e.g. Randa; Eberhardt et al. 2004).         141 Table 5.3. Basal rheology and parameters that better conform to the Holocene rock slump deposit. Bulk basal friction angle (º): φ b Pore-fluid pressure coefficient: r u Internal friction angle (º): φ i Frictional 33 0.07 35   The results of the runout analysis are shown in Figure 5.14, which corresponds with φb = 33° because it better reproduces the features of the current deposit. The bulk of the simulated flow is deposited proximally close to the base of the west slope. The rest of the flow continues moving several meters towards the SE and finally comes to rest generating a small ridge adjacent to the NNE margin of the deposit. This feature is in close agreement with a similar observable ridge in the actual deposit (Figure 5.15). However, some terrain features were not reproduced by the simulation such as a small depression that exists between the talus of the west slope and the west margin of the lobe/debris (Figure 5.15). Again, this situation is most likely attributed to the loss of topographic detail in the reconstruction of the pre-failure topography, especially in the geometry of the source volume.  Good agreement was observed between the modelled and estimated trimlines, losing conformity towards the margins. Also the model predicted deposition of debris in the upper portion of the runout path probably related to a change in slope, but no evidence has been found that may support or negate this result, with terrain artefacts always being a possible explanation.   142  Figure 5.14. Calibrated simulation of the Holocene rock slump lobe (S2) utilising a constant basal frictional rheology with a basal friction angle of 33°. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. Simulated and estimated trimlines are shown in dashed blue and dotted red lines, respectively.   143  Figure 5.15. Aerial photograph (1:20,000) showing in detail the Holocene Rock Slump deposit.  5.4. Interpretation of results The use of DAN3D in the dynamic analysis proved to be a powerful tool, with insights being provided that agreed with field observations in a number of ways including the following:     The suggested formation of Inca Lake being caused by the blocking of the valley and drainage due to the debris from the Portillo Rock Avalanche.    Model results appear to discount the possibility of the two Holocene deposits (S1 & S2) as originating from a single source (eastern slope). Instead, it is more likely that they were formed through two rock slope failure events originating from both the  144 east and west sides of the valley, probably temporally linked to the same triggering event (i.e. a large earthquake).    The runout models supported the interpretations of the chronological distribution of the deposits based on the cosmogenic nuclide dating by 36Cl.  Even though there are significant sources of uncertainty in the reconstruction of the pre- failure topography, a critical issue for 3D dynamic analyses, DAN3D was able to simulate with good accuracy the general extent and distribution of the prehistoric deposits.  Data shows that the best basal rheological combination for the Portillo Rock Avalanche is that of a frictional rockslide in the proximal path and a Voellmy rheology when entrainment becomes important. The rheology that provided the best fit to field observations for the Holocene rock slump was a constant basal frictional rheology. These results showed that, in general, the frictional basal resistance model is best suited for dry granular behaviour while the Voellmy basal model is better used when entrainment of saturated deposits becomes an important factor.  Because the Voellmy rheology is largely empirical, justification for the use of the selected parameters is mainly based on previous modelling experience following an iterative trial-and- error method (McDougall 2006). Moreover, the use of two independent parameters in the Voellmy rheology indicates that there is no unique solution for the simulated displacement of a landslide. For this reason, a unique pair of φb and ξ can be selected only when estimated flow velocities are available and comparisons with the simulated velocities can be done (Hungr et al. 2005). Keeping a complete calibrated database supported by back-analysis of real case histories is crucial for proper forward analysis used to assess landslide hazards.  A conclusive statement regarding the Upper Pleistocene Portillo Rock Avalanche and the Holocene rock slump from the opposite valley wall is unlikely due to the number of uncertainties inherent with back analyzing prehistoric slope failure events. However, the use of numerical simulation tools such as the dynamic code DAN3D combined with geological mapping and engineering judgment, enables significant insights to be gained to unravel the complex geological history in rockslide-prone valleys and to later extend these results to hazard assessments of future rockslide events.  145 6. FORWARD ANALYSIS OF POTENTIAL LARGE-SCALE ROCKSLIDE EVENTS AT PORTILLO Chapters 4 and 5 presented the back analysis of the Portillo Rock Avalanche with a detailed look at its failure initiation and runout path characteristics. Following the ‘Total Slope Analysis’ procedure (section 2.1), the insights gained in relation to the failure mechanism, geometrical factors, and material properties and rheology were used to develop and constrain a forward analysis in order to assess the hazard potential of a reoccurring major rockslide (Figure 6.1).   Figure 6.1. Flow diagram showing the procedure followed for the hazard assessment.  UDEC modelling was used to identify the potential mode of failure and estimate a potential failed volume. The major constraint for the initial stage of the forward model was that the slope has been stable for ca. 4.2 ka, which means that any attempt in modeling the current state must consider an equilibrium state. Once this state was achieved in the numerical model, an increased water table was introduced and the stability of the slope was tested. Also the stability of the slope was tested under different dynamic conditions for an earthquake trigger. The worst case scenario was defined as the one that produced the largest failed volume, quantity that was later used to develop the runout analysis using DAN3D. The simulations helped quantify, in terms of deposit depth and maximum velocity, the impact area of a modern rock slope failure at the Portillo site.  146 6.1. Potential Failure Initiation & Volume Assessment 6.1.1. Model development and results The procedure followed for the forward modelling of the initial state was the same as for the back analysis (see Figure 4.3). The development of the numerical model considered each of the insights gained during the back analysis with the major model characteristics presented in Table 6.1.  Table 6.1 Geometric configurations and mechanical properties of the forward analysis for the initial state. FORWARD MODEL CHARACTERISTICS PARAMETERS USED RATIONALE Geometry Current slope topography Forward analysis. Controlling Structures Slope: 40 x 30 m block size Zone of damaged rock: 20 x 20 m block size More computationally efficient compared to a smaller block size. Zone of Damaged Rock Densely fractured - discontinuum analysis Despite the fact that a continuum representation of domain 2 proved to be marginally better, the application of a varying water table in the forward analysis requires a discontinuity network for groundwater flow. Base Elastic and jointed 100 x 300 m block size An elastic material is more computationally efficient and has no influence on the results. A discontinuity network is required within the base to enable groundwater flow. Strength Properties Low strength, equivalent to BM - State 4 Simulate the worst case scenario by means of using lower bound strength properties (average properties proved to be stable).   The geometry was built using the current slope topography and the average material properties provided in Table 3.8 and Table 3.9. After time-stepping (R = 0.001%, see section 4.1.3 for explanation), the slope proved to be stable so the modelling conditions were shifted to a worst scenario corresponding to the lower bound material properties used in the base model - State 4. The numerical model was cycled (until R = 0.001%) and the results showed that even with the use of low strength material properties the current slope was still stable (Figure 6.2), meeting the key constraint for the initial stage (slope has not failed in the last  147 13.2 ka). The rest of the models (varying water table & dynamic loading) were built on this case.  Figure 6.2-a shows that even though some elements yield in tension and shear they are not enough to align and form a continuous yielding surface that could be interpreted as a slide surface. This stability condition is reinforced when the history plots are analyzed and a steady state condition is reached by the horizontal and vertical displacements (Figure 6.2-b). The y-velocity plot also shows equilibrium (v = 0 m/s) after a series of erratic movements. When both, plasticity indicators and history plots are considered, the constant evolution of the displacements is clearly representing an equilibrium state.   Figure 6.2. Stable condition for the modelling of the current slope using lower bound strength properties. a) Plasticity indicators (up); b) History plots, down-left: x-displacements at the toe, down- middle: y-displacements at the crest, down-right: y-velocity at the crest.    148 6.1.2. Varying water table Even though it was stated in section 3.5.4 that there was no field evidence for the presence of significant groundwater in the slope, in order to assess the potential of future major rock slope failures triggered by high pore pressures due to a series of heavy prolonged precipitation events, numerical models with varying water tables were conducted using upper and lower bounds corresponding to:      Fully saturated slope (water table at 4,000 m.a.s.l)    Partly saturated slope (water table at 3,500 m.a.s.l)  The modelling procedure consisted of using the stable model presented in section 6.1.1 to perform a coupled hydro-mechanical analysis with a steady-state fluid flow. A minor geometrical modification was introduced at the crest of the slope (top-left) just to simplify the hydromechanical analysis. The models were run independently, using the properties presented in Table 6.2 and a range of water tables as showed in Figure 6.3.  Fluid flow in UDEC is calculated based on the assumption of fracture permeability where the intact blocks are impermeable. For this a cubic law is used to correlate fracture aperture to a hydraulic conductivity. For this analysis a joint aperture of 5 mm (under zero normal stress) was assumed, decreasing to a minimum of 2 mm. These values were based on field observations. The cross-joints of the upper slope (JS2) were the only ones with different (more impermeable) hydraulic properties to better represent the field conditions in which no seepage is observed.            149 Table 6.2. Material properties used in the fluid flow models. MATERIAL PROPERTIES PARAMETERS USED REFERENCES permeability factor (k j )*: 0.83x10 2 Pa-1 sec-1 residual hydraulic aperture (a r ): 2x10 -3 m aperture at zero normal stress (a 0 ): 5x10 -3 m permeability factor (k j )*: 0.83x10 2 Pa-1 sec-1 residual hydraulic aperture (a r ): 2x10 -4 m aperture at zero normal stress (a 0 ): 5x10 -4 m (*) k j  = 1/12μ; μ = 10-3 Pa*sec is the dynamic viscosity of water Fluid properties UDEC Users Manualdensity (ρw): 1000 kg/m3 Same as BM-Stage 4 Joint hydraulic properties (bedding planes) Initial state - Tables 4.3 & 4.4; sections 4.1 & 4.2 Field observations & UDEC Users Manual Joint hydraulic properties (cross joints) Field observations & UDEC Users Manual Rock Properties & Joint Mechanical properties     Figure 6.3. Current profile (topography) of the eastern slope showing the range in which the water table is considered to vary.  150 Results for the partially saturated slope condition proved to be completely stable, therefore only the results of the fully saturated model with the water table at 4000 m.a.s.l will be presented here as a worst case scenario.  Figure 6.4 shows the results of the numerical model considering a fully saturated slope. When analyzing the plasticity indicators it is clear that there are more tensile damage indicators compared to the dry slope case presented in section 6.1.1 (Figure 6.2). These elements are directly related to the opening of the joints due to the effect of the increased pore water pressures (Figure 6.4-a). However, no significant increase above that for the dry case was seen in terms of elements at the toe of the slope failing in shear. When the history plots for points at the toe of the slope are studied, they show an almost constant horizontal and vertical displacement rate (Figure 6.4-b top), likely pointing towards a steady state where self-stabilization apparently develops. The reduced permeability of the cross joints, and therefore the reduced drainage at the face of the slope face results in that the pore pressures in the bedding joints corresponds to an hydraulic head of approximately 100 m. History plots containing the evolution of the pore water pressure is showed in Figure 6.4-b (bottom). Oscillating pressure waves are observed as the solution converges followed by stable or steady state conditions towards the end.      151  Figure 6.4. Forward modelling and volume assessment for future events considering a fully saturated slope; a) Plasticity and shear indicators, b) History plots. Displacements are measured in meters.  In summary it can be stated that even when modelling a fully saturated slope, pore water pressures are not enough to overcome the shear strength of the slope provided by the stresses that act parallel to the bedding planes and normal to the cross-joints. Nevertheless, it must be considered that this coupled hydro-mechanical analysis was done under partially drained conditions, which means that pore water pressures are not able to build-up at the toe of the slope.  6.1.3. Seismic triggering (seismic loading of slope) To assess the stability of the slope during a seismic event, a dynamic load was applied to the base of the slope to simulate an earthquake trigger. The sequential approach followed involved starting with the numerical model presented in section 6.1.1 (dry slope, rock mass and joint strength properties from BM-State 4), followed by the application of the seismic loading.   152 The dynamic modelling procedure consisted of changing the boundary conditions to those of a “quiet boundary”, which will simulate the free-field earthquake motion. This means that plane waves propagating upwards are properly absorbed in the boundaries and suffer no distortion (Itasca 2004).  A critical component of a seismic analysis is the earthquake ground motion which is well defined by an acceleration time-history. Parameters such as peak ground acceleration (PGA), mean period (Tm), and effective duration (DE; Bommer & Martinez-Pereira 1999) may be used to characterize the intensity, dominant frequency, and duration of ground motion, respectively.  The PGA is the maximum acceleration experienced by a particle in the ground during the course of the earthquake motion, and is measured in % of g (g = 9.8 m/s2). The dominant frequency depends upon the source characteristics and propagation medium (Bhasin & Kaynia 2004), and describes the distribution of the amplitude of a ground motion among different frequencies containing most of the vibration energy. The effective duration is defined as the time (in seconds) in which the strong motion phase of an earthquake occurs, and represents a continuous time window which contains the motion of engineering importance. At a rock site in the near-field of an earthquake (close to the rupture zone), DE may be related to the moment magnitude (Mw) of the earthquake using Equation 6.1 (Bommer & Martinez-Pereira 1999):  log( ) 0.69 3.70E WD M= −   [6.1]  The input variables for the dynamic analysis are presented in Table 6.3 and its calculations are explained below.         153 Table 6.3. Input variables for the UDEC dynamic analysis. PARAMETERS USED RANGE OF VALUES REFERENCES Peak Ground Acceleration: PGA %g (m/s 2 ) 0.8 - 1.2g Tanner & Shedlock (2004) Frequency: f (Hz)  2 - 5 Bhasin & Kaynia (2004) Ground Particle Velocity: v s  (m/s)  0.23 - 0.94 Bhasin & Kaynia (2004) Shear Wave Velocity: Cs (m/s) 620 - 1360 Equation 2.7 Duration: T (s)  2 - 40 Bommer & Martinez-Pereira (1999) Shear Stress: τ  (MPa)  0.8 - 7 Equation 2.5   The earthquake event was simulated using a simple harmonic sinusoidal wave applied to the base of the slope. The frequency of the wave was obtained following Bhasin & Kaynia (2004), who noted that measurements of seismic motions on rock sites are normally in the range of 2-5 Hz. The effective duration (period) was calculated using Equation 6.1, accounting for an earthquake ranging between 5 ≤ Mw ≤ 8.0. Even though the upper limit is higher than the stronger damaging earthquakes recorded close to the site (see section 3.2 for details), once again the worst case scenario was simulated for this forward analysis.  The maximum ground particle velocity was calculated from Equation 6.2 (Bhasin & Kaynia 2004). The PGA range was obtained from Tanner & Shedlock (2004), and accounts for a 2% chance of exceedance in 50 years (2%/50) for rock sites, which corresponds to a 2475-year return period. The authors consider this time frame period conservative for hazard assessments because it includes very rare large (Mw ≥ 8) subduction earthquakes.  2s PGAV f = ⋅ Π ⋅   [6.2]  Finally the excitation was applied in the form of a shear stress using equations 2.5 & 2.7 from section 2.5.3.   154 The first model series started considering the lower bound characteristics of the input seismic wave, which are a good representation of the magnitude and duration of a number of events that have occurred in the region, close to the Portillo site. Based on the absence of important failures during historic times, as expected the dynamic load applied did not result in the triggering of a catastrophic slope failure in the model. The most pronounced effect was an increased number of blocks in the upper slope yielding in tension (Figure 6.5), which may correspond with small volume of rock fall episodes during and after the earthquake.   Figure 6.5.  State of the Caracoles Range after the input of an M=5.5 earthquake. a) Plasticity indicators showing some elements yielding in tension; b) Velocity histories showing damping of the shear wave (top: x-velocity; bottom: y-velocity).  Figure 6.6 shows the performance of the Caracoles Range slope after the input of a seismic wave equivalent to an M=7.8 earthquake. Analysis of the plasticity indicators show that unlike the back analysis case, failure in tension predominates at the slope face. Moreover, when velocity vectors are overlaid, both indicators reveal the development of an outward rotation of blocks resulting in a rock slump failure for the crest of the slope. The elements that yield in tension below the hinge probably represent unravelling of the lower portion of slope associated with some degree of buckling. Close to the base of the slope, tensile failure and  155 large displacements (> 9 m) are associated to heaving due to the dilation suffered by the weaker zone of damaged rock.   Figure 6.6. State of the Caracoles Range after the input of an M=7.8 EQ. a) Plasticity indicators showing elements yielding in tension; b) Zoom to the crest of the slope showing the outward rotation of blocks (top-right) and shear strain contours (bottom); positive values indicate extension and negative indicate compression; c) Vertical displacements (bottom-left), horizontal displacements (bottom-right).     156 Even though no observable shear develops at the base of the toppled wedge to clearly indicate the possible detachment of the mass, the assumption was that important microfracturing was occurring within the zones where widespread tensile failure was observed in the model and that eventually shear was going to take place.  The rock slump yields an estimated in-place (source) failed volume ranging between 0.5 and 1 million m3, considering a failed area of ~1x104 m2 and assuming a maximum length for the scarp between 50 and 100 m based on topographic constrains (Figure 6.7).   Figure 6.7. Views of the zone affected by the rock slump. a) Google Earth picture showing the modelled cross-section and the failed crest of the slope; b) Picture of the eastern slope indicating the proposed failed volume.  6.2. Runout Assessment of Potential Failed Volume The failed volume obtained for the earthquake triggered reverse topple modelled below the crest of the slope was increased by 25% to account for bulking of the failed mass due to fragmentation. This volume was then used as an input condition for the hazard assessment. For the forward modelling there were two main scenarios:     Failure occurs in a dry season, so dry conditions for the path are expected; i.e. failed material must override previously deposited and non saturated highly frictional debris belonging to the PRA, Holocene rockslide, and recent talus.  157    Failure occurs in a wet season, where the highly frictional material in the path is overlaid by approximately 6 m of snow (IGM 2002).  Using the values obtained from the PRA back analysis, and accounting for the material path, a constant basal frictional rheology was used for the case of failure occurring under dry conditions, while a basal Voellmy rheology was used for the case where the flow would be expected to override and entrain snow. Table 6.4 shows the geometrical parameters of the proposed failure and the values used for each rheology.  Table 6.4. Geometrical parameters and rheology of the potential failed volume. Scarp surface (m2) 1x104 Out-of-plane scarp length (m) 50 - 100 Detached volume (106 m3) 0.5 - 1 Thickness* (m) 45 Vertical drop (m) 900 Slope angle (º) 60 * Average thickness REVERSE TOPPLING INPUT PARAMETERS Voellmy parameters 0.1 Turbulence parameter (m/s2): ξ 500 Detachment zone geometry Frictional rheology parameters Bulk basal friction angle (º): φb Pore-fluid pressure coefficient: ru 30 0.18 Internal friction angle (º): φi 35 Basal friction coefficient: f   For the frictional case occurring in dry conditions, the leading edge of the flow runs over part of the International Santiago-Mendoza Corridor and stops in a flat area in the upper part of the valley where between 1 and 12 m of debris will be deposited (Figure 6.8-a). For this simulation there is no direct indication of material covering the facilities of the Portillo Ski Resort, as most of the debris accumulates in the southern part of the deposit, close to the stream channel and local drainage basin. Nevertheless, given the close proximity of the leading edge of the runout debris to the resort hotel, a minor amount of debris may impact in the form of outlying boulders and rock debris.  158  Figure 6.8. DAN3D runout assessment of a potential rock slump failure at the crest of the slope. a) Top: Maximum depth of the deposit considering a highly frictional path (no snow); b) Bottom: Maximum depth of the deposit considering the path covered with snow.  159 For the second case, as previously noted, a constant Voellmy rheology was used to simulate the influence of snow along the runout path (Figure 6.8-b). The results show a slightly different runout pattern from the ones obtained for the frictional rheology, specifically in the distribution of debris. For this case the flow tends to separate in two streams: a southern stream which carries most of the volume and is spatially related to the drainage network (i.e. topographic lows), and a northern stream which appears to follow another channel like feature in the topography. The leading edge in this case reaches approximately 90 m ahead of the constant frictional case, an expected behaviour which is explained by its more mobilised behaviour. In general terms, this simulation clearly highlights the sensitivity of the Voellmy rheology to the path topography.  The velocity profiles of the frictional (Figure 6.9-a) and Voellmy (Figure 6.9-b) rheologies also show several differences, the most noticeable one being the lower maximum particle velocity reached by the Voellmy rheology. Once again, this result is expected due to the well known behaviour of frictional flows that tend to overestimate maximum velocities (Hungr et al. 2005). Nevertheless, both rheologies show that the maximum velocities are reached by the southern limit of the flow, which coincides with the local stream channel. This result is important because if a more mobilised flow or more material is considered, the potential reach and impact of the runout increases down the valley.   160  Figure 6.9. Maximum velocities of the potential failed volume. a) Accounting for a dry path; b) Accounting for a path covered by snow.  161 6.3. Hazard Assessment Summary The performance of the present-day slope state was tested under different scenarios. If unstable, potential failed volumes were estimated using the interpreted slip surface by means of UDEC and extrapolating in the out-of-plane direction (constrained by topographic and geomorphic features). The volume was then bulked, to account for fragmentation and transport of the slide mass, and then used as input for a runout analysis using DAN3D. Two different rheologies were utilised and a range of parameters were tested. The results were then presented in the form of deposit depth, impacted area, and maximum velocity contour maps. The outline of the velocity maps, which coincides with the velocity = 0 m/s contour, was cut-off to exactly match the corresponding depth = 1 m contour in each depth map.  Table 6.5 briefly summarizes the susceptibility to failure of the eastern slope after different conditions were applied, and Table 6.6 shows the hazard assessment analysis presented in this chapter for the present-day state of the slope.  Table 6.5. Susceptibility to failure of the eastern slope. SUSCEPTIBILITY TO FAILURE OF THE PRESENT DAY SLOPE Static Stable N/A State Condition Estimated in-place failed volume (m 3 ) N/A Seismic loading (M=7.8) Failed 1x106 Fully saturated Stable Seismic loading (M=5.5) Stable N/A          162 Table 6.6. Hazard assessment for the eastern slope. Seismic loading (M=7.8) PRESENT DAY SLOPE HAZARD ASSESSMENT State Runout behaviour Approximate maximum runout distance (m) Frictional (dry season) 2180 Voellmy (wet season) 2270  The first trial, using the same rock mass and discontinuity properties as those back calculated in Chapter 4, proved to be stable, a result when combined with the fact that no further major failures have occurred in recorded times, points towards a stabilized geometrical profile with time. Results from the second trial that build on the first by simulating a fully saturated slope condition in which pore water pressures preferentially build up at the toe of the slope along bedding planes (due to the nature of the tighter cross joints) like wise were not able to sufficiently overcome the strength of the discontinuities and indirectly the strength of the rock mass.  The slope failed only in the third trial, when the slope was subjected to a seismic load equivalent to a M=7.8 earthquake, a magnitude significantly greater than the typical for the vicinity of the site (see section 3.2). In this case, a different failure mode was observed than that for the back analysis. The outward rotation of blocks at the crest of the slope resulted in a potential rock slump failure, probably aided by topographic amplification which in general is produced at convex topographies due to seismic waves focusing effects (Assimaki & Gazetas 2004). Below the hinge of the slope, buckling is likely to occur probably responding as unravelling producing rockfall material.  Based on this worst case scenario and calculated failed volumes, the hazard assessment proceeded with runout simulations for both a highly frictional path and a path covered by deep snow. The results suggest that the leading edge of the flow would override part of the International Santiago-Mendoza Corridor and the debris would come to rest in a flat-lying area in the upper part of the valley. In both cases the Portillo Ski Resort does not appear to be directly impacted, nevertheless the proximity of the distal edge of the runout indicates probable damage of the resort’s structures due to impact by rockslide boulders at the margins of the flow.  163 7. CONCLUSIONS AND RECOMMENDATIONS 7.1. Conclusions This study investigates the prehistoric rock mass wasting events occurred in the Portillo site, and the potential for another catastrophic rockslide that may put in risk the structures in that area. For this, an integrated approach containing detailed field mapping and investigations combined with numerical modelling were employed.  The use of cosmogenic nuclide dating by 36Cl in conjunction with detailed geologic and geomorphologic mapping helped constrain the origin of the deposits at the toe of the Caracoles Range. The chronology of events set them in a post glacial period. Late Pleistocene and Holocene ages clustered around 4.2 ka and 13.2 ka, thus ruling out a glacial origin for any of the lobes.  The use of discontinuum numerical modeling techniques applied in the analysis of the Portillo Rock Avalanche supported the field observations with respect to the most probable failure mechanism. A stress-controlled failure at the toe of the slope represented by shearing of the rock mass continuum enabled a structurally controlled failure in the upper part of the slope by means of sliding and shearing along bedding planes. The resulting model also provided new insights into the way the slide mass failed, suggesting the possibility of a staged failure through an active-passive Prandtl wedge. The lack of daylighting discontinuities dipping shallowly out of the slope and the unchanging failure mechanism to variations in rock mass and discontinuity strength parameters, strongly suggests a failure mode controlled by the adverse geometry, with high and steep dipping bedding planes combined with back thrust structures at the toe of the slope. Several possibilities could be mentioned that might have helped or triggered the slope failure, the most important being:     Glacial rock mass damage which led to progressive failure mechanisms, both in the form of asperity breakdown along the bedding planes and loss of coherence in the stronger units at the toe of the slope through brittle fracture processes, attributed to shearing of the valley sides by the ice-sheet that advanced along the valley floor;  164    Glacial oversteepening of the valley walls during glacial advance and debuttressing of the paleoslope during glacial retreat (i.e. after the LGM, approximately 15 ka);    An exceptionally large magnitude earthquake.  The use of DAN3D in the dynamic analysis proved to be a powerful tool for modelling rockslide runouts, with insights being provided that agreed with field observations in a number of ways:     The formation of Inca Lake arising through the blocking of the natural valley drainage system due to the debris of the Portillo Rock Avalanche is almost certain;    Model results appear to discount the possibility of the two Holocene deposits, found at the bottom of the east and west facing slopes, as originating from a single source. Instead, it is more likely that they were formed through two rock slope failure events, one from the eastern slope and one from the western slope, perhaps temporally linked to the same triggering event (i.e. a large earthquake);    The runout models support the interpretation of the chronological distribution of the deposits, based on the cosmogenic nuclide dating by 36Cl.  The forward modelling helped characterise the potential hazard posed by the eastern slope to the Portillo site. The static and dynamic conditions under which the slope was tested revealed that:     The present-day stability of the eastern slope under static conditions and low strength mechanical properties is probably the product of a more stable slope profile following the prehistoric Portillo Rock Avalanche event;    No failure of the current slope occurred when using a coupled hydro-mechanical analysis representing a fully saturated slope. The nature of the bedding and fracture permeability network for which the water pressures would be concentrated at the toe of the slope, still did not produce an unstable condition exceeding the rock mass effective strength;    A potential instability may develop for the present-day slope in the event of a high magnitude (M=7.8) earthquake generated near or below the slope. The modelled effects of the seismic loading resulted in the instability of the crest of the slope, probably augmented by topographic amplification;  165    The proposed mode of failure for the crest of the slope in the event of an earthquake trigger scenario is a rock slump, where the failed wedge yielded an estimated volume of 1 million m3;    The results of the hazard assessment suggested that the leading edge of the flow could override part of the International Santiago-Mendoza Corridor considering both, a path covered by snow and a dry path, and that for the first case the maximum runout distance would be approximately 90 m ahead of the dry case;    The Portillo Ski Resort does not appear to be directly impacted by the flow debris; however its close proximity with the edge of the deposit may result in some damage to the structure.  It must be stressed though that an earthquake of this magnitude has not been recorded in the Portillo region, thus the nature of the hazard is very low.  The ‘Total Slope Analysis’ proved to be an efficient methodology towards the overall understanding of rock slope behaviour because of its advantageous coupling of the back analysis and the forward modelling. A conclusive statement regarding the Upper Pleistocene Portillo Rock Avalanche and the Holocene rock slump is limited due to the number of uncertainties inherent in back analyzing any prehistoric slope failure event. However, the use of numerical simulation tools such as the distinct element code UDEC and the rheological runout code DAN3D combined with geological mapping, field observations and engineering judgment, enables significant insights to be gained to unravel the complex geological history in rockslide prone valleys and to later extend these results to hazard assessments of future rockslide events.  7.2. Recommendations for Further Work Although the back analysis proved to be consistent with field observations and dating tools, some limitations were identified for which recommendations for further studies can be made. These include:  1. The use of a strain softening constitutive model in the UDEC analysis that will better capture the degradation of rock strength with increasing shear strain.  166 2. A staged analysis in the UDEC modelling that includes the glacial loading, unloading and oversteepening of the valley walls.  The results obtained in the forward analysis are by no means considered as conclusive, and their limitations must be accounted for in future works. The most important aspects that were left out of this research either because they did not directly correspond to the research objectives are:  1. Assess the likelihood of potential runout flows channelized by the valley drainage network by means of changing Voellmy parameters that make the flow more mobile when snow/water is entrained. 2. Assess the rockfall potential of the western slope, source of the Holocene rock slump, which currently is a very active talus slope with permanent minor rock detachments from the crest of the slope.      167 REFERENCES Abele, G. 1997. Rockslide movement supported by the mobilization of groundwater- saturated valley floor sediments. Zeitschrift fur Geomorphologie, 41: 1-20. Assimaki, D., and Gazetas, G. 2004. Soil and topographic amplification on canyon banks and the 1999 Athens earthquake. Journal of Earthquake Engineering, 8(1): 1-43. ASTM 1995. Standard Test Method for Unconfined Compressive Strength of Intact Rock Core Specimens. American Society for Testing and Materials (ASTM), Designation: D 2938 – 95 (Reapproved 2002). ASTM 1996a. Standard Test Method for Field Determination of Apparent Specific Gravity of Rock and Manmade Materials for Erosion Control. American Society for Testing and Materials (ASTM), Designation: D 5779 – 95a (Reapproved 2001). ASTM 1996b. Standard Test Method for Triaxial Compressive Strength of Undrained Rock Core Specimens Without Pore Pressure Measurements. American Society for Testing and Materials (ASTM), Designation: D 2664 – 95a. ASTM 2001. Standard Practices for Preparing Rock Core Specimens and Determining Dimensional and Shape Tolerances. American Society for Testing and Materials (ASTM), Designation: D 4543 – 01. ASTM 2002. Standard Test Method for Performing Laboratory Direct Shear Strength Tests of Rock Specimens Under Constant Normal Force. American Society for Testing and Materials (ASTM), Designation: D 5607 – 02. Ballantyne, C.K., Stone, J.O., and Fifield, L.K. 1998. Cosmogenic Cl-36 dating of postglacial landsliding at The Storr, Isle of Skye, Scotland. The Holocene, 8(3): 347–351. Barrientos, S., and Eisenberg, A. 1988. Secuencia sismica en la zona cordillerana al interior de Rancagua. In 5th Chilean Geological Conference. Edited by. Santiago, Chile, pp. F121-F132. Barrientos, S., Vera, E., Alvarado, P., and Monfret, T. 2004. Crustal seismicity in central Chile. Journal of South American Earth Sciences, 16(8): 759 – 768. Barton, N. 1972. A model study of rock-joint deformation. International Journal of Rock Mechanics and Mining Sciences, 9: 579-602. Barton, N. 1973. Review of a new shear strength criterion for rock joints. Engineering Geology, 7: 287-322. Barton, N.R., Lien, R., and Lunde, J. 1974. Engineering classification of rock masses for the design of tunnel support. Rock Mechanics, 6(4): 189 – 236. Beck, S., Barrientos, S., Kausel, E., and Reyes, M. 1998. Source characteristics of historical earthquakes along the central Chile subduction zone. Journal of South American Earth Sciences, 11(2): 115-129. Benko, B., and Stead, D. 1998. The Frank slide: a reexamination of the failure mechanism. Canadian Geotechnical Journal, 35: 299-311. Bhasin, R., and Kaynia, A.M. 2004. Static and dynamic simulation of a 700-m high rock slope in western Norway. Engineering  Geology, 71: 213-226.  168 Bieniawski, Z.T. 1973. Engineering classifications of jointed rock masses. Transactions of South African Institute of Civil Engineering, 15(12): 335 – 344. Bieniawski, Z.T. 1976. Rock mass classification in rock engineering. In Exploration for rock engineering. Edited by. Cape Town. Vol. 1, p. 97–106. Rotterdam: Balkema, Vol.1, pp. 97-106. Bieniawski, Z.T. 1989. Engineering Rock Mass Classifications. John Wiley, New York. Bommer, J., and Martinez-Pereira, A. 1999. The effective duration of earthquake strong motion. Journal of Earthquake Engineering, 3(2): 127-172. Brown, E.T. 1981. Rock Characterization Testing and Monitoring – ISRM Suggested Methods. Oxford: Published for the Commission on Testing Methods, International Society for Rock Mechanics. Pergamon. Cahill, T., and Isacks, B. 1992. Seismicity and shape of the subducted Nazca Plate. Journal of Geophysical Research, 97(B12): 17503 - 17529. Cai, M., Kaiser, P.K., Uno, H., Tasaka, Y., and Minami, M. 2004. Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the GSI System. International Journal of Rock Mechanics & Mining Sciences, 41(1): 3 – 19. Cavers, D.S. 1981. Simple methods to analyze buckling of rock slopes. Rock Mechanics, 14(2): 87-104. Caviedes, C., and Paskoff, R. 1975. Quaternary Glaciations in the Andes of North-Central Chile. Journal of Glaciology, 14(70): 155-170. Cegarra, M., and Ramos, V. 1996. The Aconcagua fold and thrust belt. In Geology of the Aconcagua region; San Juan and Mendoza provinces, Argentina. Dirección Nacional del Servicio Geológico. pp. 387 - 422. Charrier, R., Baeza, O., Elgueta, S., Flynn, J.J., Gana, P., Kay, S.M., Muñoz, N., Wyss, A.R., and Zurita, E. 2002. Evidence for Cenozoic extensional basin development and tectonic inversion south of the flat-slab segment, southern Central Andes, Chile (33º- 36º S.L.). Journal of South American Earth Sciences, 15: 117-139. Coggan, J.S., Stead, D., and Eyre, J. 1998. Evaluation of techniques for quarry slope stability assessment. Transactions of the Institute of Mining and Metallurgy, Sect. B: Applied Earth Sciences, 107: B139 - B147. Comte, D., Eisenberg, A., Lorca, E., Pardo, M., Ponce, L., Saragoni, R., Singh, S.K., and Suarez, G. 1986. The 1985 central Chile earthquake: a repeat of previous great earthquakes in the region? Science, 233: 449-453. Corominas, J. 1996. The angle of reach as a mobility index for small and large landslides. Canadian Geotechnical Journal, 33: 260-271. Costa, C.H., and Gonzalez-Diaz, E.F. 2007. Age constraints and paleoseismic implication of rock avalanches in the northern Patagonian Andes, Argentina. Journal of South American Earth Sciences, 24: 48–57. Costa, M., Coggan, J.S., and Eyre, J.M. 1999. Numerical modelling of slope behaviour at Delabole slate quarry. International Journal of Surface Mining, Reclamation & Environment, 13: 11 – 18. Crosta, G.B., Imposimato, S., and Roddeman, D.G. 2003. Numerical Modeling of Large Landslides Stability and Runout. Natural Hazards and Earth System Sciences, 3: 523 - 538.  169 Cruden, D.M., and Varnes, D.J. 1996. Landslide types and processes, National Academy Press, Washington D.C. Davies, T.R. 1982. Spreading Of Rock Avalanche Debris by Mechanical Fluidization. Rock Mechanics, 15: 9 - 24. Deere, D.U. 1964. Technical description of rock cores for engineering purposes. Rock Mechanics & Engineering Geology, 1(1): 17 – 22. DeMets, C., Gordon, R.G., Argus, D.F., and Stein, S. 1994. Effect of recent revisions to the geomagnetic reversal time scale on estimate of current plate motions. Geophysical Research Letters, 21(20): 2191-2194. Dunne, J., Elmore, D., and Muzikar, P. 1999. Scaling factors for the rates of production of cosmogenic nuclides for geometric shielding and attenuation at depth on sloped surfaces. Geomorphology, 27: 3-11. Dunne, J.A., and Elmore, D. 2003. Monte Carlo simulations of low-energy cosmogenic neutron fluxes near the bottom of cliff faces. Earth and Planetary Science Letters, 206: 43-49. Eberhardt, E., and Stead, D. 1998. Mechanisms of slope instability in thinly bedded surface mine slopes. In 8th Congress International Association of  Engineering Geology and the Environment. Edited by. Vancouver. Balkema, p. 3011 – 3018. Eberhardt, E., Stead, D., and Coggan, J.S. 2004. Numerical analysis of initiation and progressive failure in natural rock slopes - the 1991 Randa rockslide. International Journal of Rock Mechanics & Mining Sciences, 41(1): 69-87. Einstein, H., Veneziano, D., Baecher, G.B., and O'Reilly, K.J. 1983. The effect of discontinuity persistence on rock slope stability. International Journal of Rock Mechanics & Mining Sciences, 20(5): 227 – 236. ESRI 2004. www.esri.com/. Fock, A., Charrier, R., Farias, M., Maksaev, V., Fanning, M., and Alvarez, P. 2005. Exhumation and uplift of the western Main Cordillera between 33º and 34ºS. In 6th International Symposium on Andean Geodynamics (ISAG 2005), Extended Abstracts. Edited by. Barcelona, Spain, pp. 273-276. Froldi, P., and Lunardi, P. 1995. Buckling Failure Phenomena and Their Analysis. In Mechanics of Jointed and Faulted Rocks. Balkema, Rotterdam. pp. 595-604. Giambiagi, L.B., and Ramos, V.A. 2002. Structural evolution of the Andes in a transitional zone between flat and normal subduction (33°30' - 33°45'), Argentina and Chile. Journal of South American Earth Sciences, 15: 101-116. Godoy, E. 1994. Derrumbes de Cerros Holocenos en los Andes Centrales de Chile. In Proceedings of the 7th Chilean Geological Conference. Edited by. Concepcion, Chile, Vol.1, pp. 310-314. Godoy, E. 2002. Remociones en masa de aspecto Cuaternario en el frente tectonico cordillerano de Chile central, los depositos de Colo-Coya y Gualtatas, Chile. In International Symposium of Environmental Geology for Land Use Planning. Edited by. Puerto Varas, Chile, Vol.1, pp. 58-61. Godoy, E. 2005. Notes for a Viña del Mar-Portillo geologic transect (33ºS). Chilean Geological and Mining Survey - Dalhousie University, Halifax, p. 8.  170 Godoy, E., Yañez, G., and Vera, E. 1999. Inversion of an Oligocene volcano-tectonic basin and uplifting of its superimposed Miocene magmatic arc in the Chilean Central Andes: first seismic and gravity evidences. Tectonophysics, 306: 217–236. Golden_Software 2002. Surfer 8.0. Golden Software, Inc., Golden, Colorado. Gosse, J.C., and Phillips, F.M. 2001. Terrestrial in situ cosmogenic nuclides: theory and application. Quaternary Science Reviews, 20: 1475-1560. Heim, A. 1932. Landslides and Human Lives. BiTech Publishers Ltd. Hermanns, R.L., and Strecker, M.R. 1999. Structural and lithological controls on large Quaternary rock avalanches (sturzstroms) in arid northwestern Argentina. Geological Society of America Bulletin, 111(6): 934-948. Hermanns, R.L., Niedermann, S., Ivy-Ochs, S., and Kubik, P.W. 2004. Rock avalanching into a landslide-dammed lake caus-ing multiple dam failure in Las Conchas valley (NW Argen-tina) - evidence from surface exposure dating and strati-graphic analyses. Landslides, 1: 113-122. Hermanns, R.L., Niedermann, S., Villanueva-Garcia, A., Sosa-Gomez, J., and Strecker, M.R. 2001. Neotectonics and catstro-phic failure of mountain fronts in the southern intra- Andean Puna Plateau, Argentina. Geology, 29(7): 619-623. Heusser, C.J. 2003. Ice age southern andes: a chronicle of paleoecological events. Elsevier, Tuxedo, USA. Hoek, E. 1994. Strength of rock and rock masses. News J ISRM, 2(2): 4 – 16. Hoek, E., and Brown, E.H. 1980. Underground excavations in rock. London: Institution of Mining and Metallurgy. Hoek, E., and Brown, E.T. 1997. Practical estimates of rock mass strength. International Journal of Rock Mechanics, Mining Sciences & Geomechanics Abstracts, 34(8): 1165 – 1186. Hoek, E., Kaiser, P.K., and Bawden, W.F. 1995. Support of Underground Excavations in Hard Rock. Rotterdam: Balkema. Hoek, E., Marinos, P., and Benissi, M. 1998. Applicability of the geological strength index (GSI) classification for weak and sheared rock masses: the case of the Athens schist formation. Bulletin of Engineering Geology and the Environment, 57(2): 151-160. Hoek, E., Carranza-Torres, C.T., and Corkum, B. 2002. Hoek-Brown failure criterion - 2002 edition. In Proceedings of the North American Rock Mechanics Society (NARMS- TAC), Mining Innovation and Technology. Edited by. Toronto, Ontario, Canada. July, 2002, pp. 267-273. Hsu, K.J. 1975. Catastrophic debris streams (Sturzstroms) generated by rockfalls. Geological Society of America Bulletin, 86: 129-140. Huang, T.H., Chang, C.S., and Yang, Z.Y. 1995. Elastic Moduli for Fractured Rock Mass. Rock Mechanics and Rock Engineering, 28(3): 135-144. Hudson, J.A., and Harrison, J.P. 2005. Engineering Rock Mechanics. An Introduction to Principles. Pergamon. Hungr, O. 1995. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Canadian Geotechnical Journal, 32: 610-623.  171 Hungr, O., and Evans, S.G. 1996. Rock avalanche runout prediction using a dynamic model. In Proceedings of the 7th International Symposium on Landslides. Edited by. Trondheim, Norway. A.A. Balkema, pp. 233-238. Hungr, O., and Evans, S.G. 2004. Entrainment of debris in rock avalanches; an analysis of a long run-out mechanism. Geological Society of America Bulletin, 116(9/10): 1240- 1252. Hungr, O., Corominas, J., and Eberhardt, E. 2005. Estimating landslide motion mechanism, travel distance, and velocity. In International Conference on Landslide Risk Management. Edited by. Vancouver, Canada. A.A. Balkema, p. 764. Hungr, O., Evans, S.G., Bovis, M.J., and Hutchinson, J.N. 2001. A Review of the Classification of Landslides of the Flow Type. Environmental & Engineering Geoscience, 7(3): 221 - 238. Hutchinson, J.N. 1984. An Influence line approach to the stabilization of slopes by cuts and fills. Canadian Geotechnical Journal, 21(2): 363-370. Itasca 2004. Itasca Software Products—FLAC, FLAC3D, UDEC, 3DEC, PFC, PFC3D. Iverson, R.M., Schilling, S.P., and Vallance, J.W. 1988. Objective delineation of lahar inundation zones. Geological Society of America Bulletin, 110: 972-984. Jaboyedoff, M., Baillifard, F., Couture, R., Locat, J., and Locat, P. 2004. Toward preliminary hazard assessment using DEM topographic analysis and simple mechanical modeling by means of sloping local base level. In Proceedings of the Ninth International Symposium on Landslides. Edited by. Rio de Janeiro, Brazil. Balkema, Vol.1, pp. 199-205. Jaboyedoff, M., Baillifard, F., Couture, R., Derron, M.H., Locat, J., and Locat, P. 2005. Coupling kinematic analysis and sloping local base level criterion for large slope instabilities hazard assessment - a GIS approach. In Proceedings of the International Conference on Landslide Risk Managment. Edited by. Vancouver, Canada. Balkema, Vol.1, pp. 615-622. Kay, S., Maksaev, V., Moscoso, R., Mpodozis, C., Nasi, C., and Gordillo, C. 1988. Tertiary magmatism in Chile and Argentina between 28°S and 33°S: correlation of magmatic chemistry with changing Benioff zone. Journal of South American Earth Sciences, 1(1): 21-38. Keefer, D.K. 1984. Landslides caused by earthquakes. Geological Society of America Bulletin, 95: 406-421. Kemeny, J., and Donovan, J. 2005. Rock mass characterization using LiDAR and automated point cloud processing. Ground Engineering, 38(11): 26-29. Kimber, O.G., Allison, R.J., and Nicholas, J.C. 1998. Mechanisms of Failure and Slope Development in Rock Masses. Transactions of the Institute of British Geographers, New Series, 23(3): 353 - 370. Kulhawy, F.H. 1975. Stress deformation properties of rock and rock discontinuities. Engineering Geology, 9: 327-350. Kumar, P. 1998. Shear failure envelope of Hoek-Brown criterion for rock mass. Tunnelling and Underground Space Technology, 13(4): 453-458.  172 Kvapil, R., and Clews, K.M. 1979. An examination of the Prandtl mechanism in large- dimension slope failures. Transactions of the Institute of Mining and Metallurgy, Sect. A: Mining Industry, 88: A1 – A5. Lavenu, A., and Cembrano, J. 1999. Compressional- and transpressional-stress pattern for Pliocene and Quaternary brittle deformation in fore arc and intra-arc zones (Andes of Central and Southern Chile). Journal of Structural Geology, 21: 1669-1691. Li, Q., and Zhang, Z. 1990. Mechanism of buckling and creepbuckling failure of the bedded rock masses on the consequent slopes. Journal of Chengdu College of Geology, 17(4): 97-103. Li, T. 1983. A mathematical model for predicting the extent of a major rockfall. Zeitschrift fur Geomorphologie, 27(4): 473-482. Lomnitz, C. 1961. Los terremotos del 4 de Septiembre en el Cajon del Maipo. Anales de la Facultad de Ciencias Fisicas y Matematicas, 3(18): 279-306. Marinos, P., and Hoek, E. 2001. Estimating the geotechnical properties of heterogeneous rock masses such as flysch. Bulletin of Engineering Geology and the Environment(60): 85-92. Marinos, P., Marinos, V., and Hoek, E. 2007. The Geological Strength Index (GSI): A characterization tool for assessing engineering properties for rock masses. In Proceedings of the International Workshop on Rock Mass Classification in Underground Mining. Edited by. Vancouver, Canada. Department of Health and Human Services. Marinos, V., Marinos, P., and Hoek, E. 2005. The geological strength index: applications and limitations. Bulletin of Engineering Geology and the Environment, 64(1): 55-65. McDougall, S. 2006. A new continuum dynamic model for the analysis of extremely rapid landslide motion across complex 3D terrain, University of British Columbia, Vancouver. McDougall, S., and Hungr, O. 2004. A model for the analysis of rapid landslide motion across three-dimensional terrain. Canadian Geotechnical Journal, 41: 1084-1097. Mencl, V. 1966. Mechanics of landslides with non-circular slip surfaces with special reference to the Vaiont slide. Geotechnique, 16(4): 329 – 337. Nilsen, B. 1987. Flexural buckling of hard rock; a potential failure mode in high rock slopes? In Sixth International Congress on Rock Mechanics. Edited by. International Society for Rock Mechanics, Vol.1, pp. 457-461. Okuda, S., Ouichi, T., and Terashima, T. 1992. Deviation of magnitude frequency distribution of earthquakes from the Gutenberg-Richter law: detection of precursory anomalies prior to large earthquakes. Physics of The Earth and Planetary Interiors, 73(3-4): 229- 238. Olsson, R., and Barton, N. 2001. An improved model for hydromechanical coupling during shearing of rock joints. International Journal of Rock Mechanics and Mining Sciences, 38(3): 317-329. Palmstrom, A. 1985. The volumetric joint count - A useful and simple measure of the degree of rock mass jointing. In International Association of Engineering Geologists Congress. Edited by. New Delhi, pp. V221-V228.  173 Palmstrom, A. 2005. Measurements of and correlations between block size and rock quality designation (RQD). Tunnelling and Underground Space Technology, 20: 362–377. Pardo, M., Comte, D., and Monfret, T. 2002. Seismotectonic and stress distribution in the central Chile subduction zone. Journal of South American Earth Sciences, 15(1): 11 - 22. Phillips, F.M., and Plummer, M.A. 1996. CHLOE: A program for interpreting in situ cosmogenic nuclide data for surface exposure dating and erosion studies. Radiocarbon, 38: 98-99. Pritchard, M.A., and Savigny, K.W. 1990. Numerical Modeling of Toppling. Canadian Geotechnical Journal, 27(6): 823 - 834. Ramos, V.A., Zapata, T., Cristallini, E., and Introcaso, A. 2004. The Andean Thrust System - Latitudinal variations in structural styles and orogenic shortening. In Thrust tectonics and hydrocarbon systems: AAPG Memoir 82. pp. 30-50. Rivano, G.S., Sepúlveda, H.P., Boric, P.R., and Espiñeira, T.D. 1993. Hojas Quillota y Portillo. In Carta Geológica de Chile, No. 73. Rocscience 2002. RocLab User's Guide. http://www.rocscience.com, Toronto. Rocscience 2005. http://www.rocscience.com/. Rodriguez, C.E., Bommer, J.J., and Chandler, R.J. 1999. Earthquake-induced landslides: 1980-1997. Soil Dynamics and Earthquake Engineering, 18: 325-346. Scheidegger, A. 1973. On the prediction of the reach and velocity of catastrophic landslides. Rock Mechanics, 5: 231-236. Singh, M., and Rao, K.S. 2005. Empirical methods to estimate the strength of jointed rock masses. Engineering Geology, 77: 127-137. Sonmez, H., and Ulusay, R. 1999. Modifications to the geological strength index (GSI) and their applicability to stability of slopes. International Journal of Rock Mechanics & Mining Sciences, 36(6): 743 – 760. Split_Engineering 2005. Split-FX User's Manual; Beta version 1.0. Split Engineering LLC, p. www.spliteng.com. Starfield, A.M., and Cundall, P.A. 1988. Towards a methodology for rock mechanics modelling. International Journal of Rock Mechanics, Mining Sciences & Geomechanics Abstracts, 25(3): 99-106. Stead, D., and Eberhardt, E. 1997. Developments in the analysis of footwall slopes in surface coal mining. Engineering Geology, 46: 41-61. Stead, D., Coggan, J.S., and Eberhardt, E. 2004. Realistic simulation of rock slope failure mechanisms; the need to incorporate principles of fracture mechanics. International Journal of Rock Mechanics & Mining Sciences, 41(1): 563-568. Stead, D., Eberhardt, E., and Coggan, J.S. 2006. Developments in the characterization of complex rock slope deformation and failure using numerical modelling techniques. Engineering Geology, 83(1-3): 217-235. Stesky, R., and Pearce, W. 2007. SpheriStat 2.2. http://www.pangaeasci.com/.  174 Strouth, A. 2006. Integrated Use of Terrestrial Laser Scanning and Advanced Numerical Methods for a Total Slope Analysis of Afternoon Creek, Washington, UBC, Vancouver. Strouth, A., and Eberhardt, E. 2005. The use of LiDAR to overcome rock slope hazard data collection challenges at Afternoon Creek, Washington. In 41st US Symposium on Rock Mechanics, Golden. Edited by. American Rock Mechanics Association, Vol.CD:06-993. Sturzenegger, M., Yan, M., Stead, D., and Elmo, D. 2007. Applications and limitations of ground-based laser scanning in rock slope characterization. In Proceedings of the 1st Canada-US Rock Mechanics Symposium,Vancouver, Canada, 27-31 May. Edited by. Eberhardt, Stead & Morrison. Taylor & Francis Group, Vol.1, pp. 29-36. Tanner, J.G., and Shedlock, K.M. 2004. Seismic hazard maps of Mexico, the Caribbean, and Central and South America. Tectonophysics, 390(1-4): 159-175. Tebbens, S.F., and Candie, S.C. 1997. Southeast Pacific tectonic evolution from early Oligocene to present. Journal of Geophysical Research, 102(6): 12061-12084. Terzaghi, K. 1950. Mechanism of Landslides. Engineering Geology: 83 - 123. Terzaghi, K. 1962. Stability of steep slopes on hard unweathered rock. Geotechnique, 12: 251-270. Tzamos, S., and Sofianos, A.I. 2007. A correlation of four rock mass classification systems through their fabric indices. International Journal of Rock Mechanics & Mining Sciences, 44(4): 477 – 495. Welkner, D., Eberhardt, E., and Hermanns, R.L. 2007a. Investigation of possible failure mechanisms of the Portillo Rock Avalanche, Central Andes, Chile. In Proceedings of the 1st Canada-US Rock Mechanics Symposium. Edited by. Vancouver, Canada. Balkema, Vol.2, pp. 909-916. Welkner, D., Eberhardt, E., and Hungr, O. 2007b. A Dynamic Analysis of Prehistoric Rock Slope Failure Events in the Central Andes of Chile. In OttawaGeo2007 - The Diamond Jubilee. Edited by. CGS, pp. 748-755. Wyllie, D.C., and Mah, C.W. 2004. Rock Slope Engineering: Civil and Mining. Spon Press, New York. Yañez, G., Cembrano, J., Pardo, M., Ranero, C., and Selles, D. 2002. The Challenger-Juan Fernandez-Maipo major tectonic transition of the Nazca-Andean subduction system at 33-34°S: geodynamic evidence and implications. Journal of South American Earth Sciences, 15: 23-38. Zangerl, C.J. 2003. Analysis of surface subsidence in crystalline rocks above the Gotthard Highway Tunnel, Switzerland, Swiss Federal Institute of Technology (ETH), Zurich.      175               APPENDIX A: GEOMECHANICAL LABORATORY TEST RESULTS               176 A.1 Description and Data Analysis for Triaxial Compressive Strength Test  Triaxial compressive strength tests were conducted in a standard triaxial pressure vessel within a servo-controlled 1000 kN (225 kips) loading frame (see Figure 3.22-b). The axial stress and confining pressure were independently controlled using electro-hydraulic servo- controlled feedback systems. Two linear variable differential transducers (LVDT), attached to the end caps recorded axial displacement. Circumferential displacement measurements were recorded using diametric deformation LVDT mounted on chain length, which in turn wrapped around circumference of the sample. Axial load was measured with an internal load cell and confining pressure was measured throughout the test with a hydraulic pressure transducer. The general procedure for the triaxial compressive test is summarized below:     A cylindrical core sample is prepared for testing and their ends ground parallel according to International Society for Rock Mechanics (ISRM) and American Society for Testing and Materials (ASTM) standards. A length to diameter ratio of 2:1 is recommended by ASTM and ISRM to obtain representative mechanical properties of the sample. However, it is difficult to get such ideal samples from the field and physical dimensions of the specimen are recorded after it is prepared.    The specimen is then placed between two end caps and a heat-shrink Teflon jacket is placed over the specimen to secure to the end caps for triaxial test.    Two LVDTs for axial strain and a radial LVDT are mounted in the end caps and on the lateral surface of the specimen, respectively.    The pressure vessel is lowered on to the specimen assembly and the loading frame is brought in contact with a loading piston that allows for the application of the axial load.    Confining pressure is increased to the desired hydrostatic testing pressure for triaxial test.    The axial load is programmed to increase with a controlled axial strain rate of 0.05% to 0.1% per minute until the specimen fails or axial strain reaches a 5% of strain while the confining pressure is held constant.    The axial stress is reduced to the initial hydrostatic condition after the sample fails or reaches a desired axial strain.  177    The confining pressure is reduced to zero and the sample stack is disassembled and photographed.  The data analysis commences with the determination of the axial stress by dividing the measured load by the initial cross-sectional area of the specimen. Deviator or differential axial stresses are plotted against both axial strain εL = ΔL/Lo (where Lo is the initial length and ΔL is the length change) and radial strain εR = ΔD/Do (where Do is the initial diameter and ΔD is the diameter change) (see Figure A.1). Differential/Deviator stress (σd) is defined as the difference between the total axial stress (σ1) and the confining pressure (Pc). For the sign conventions, compressive stress and contraction (shortening) are considered positive. Therefore, a positive axial strain indicates a shortening of the specimen length and a negative radial strain indicates an increase of the specimen diameter during the test.  The ultimate compressive strength of the specimen is defined as the maximum total stress at failure (= σd + Pc). Static Young’s modulus (E) was determined by the initial linear-least- square slope of the differential stress versus the axial strain curve at the interval of 0 to 50% of the ultimate strength value of a given sample, i.e. secant modulus as per ASTM D5407. Static Poisson’s ratio (ν) was determined by the linear-least-square slope of the radial strain versus the axial strain curve at the same interval where the Young’s modulus was determined.                178  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa -0.2 0 0.2 0.4 0.6 0.8 1 0 100 200 300 Axial Radial  Axial stress vs Axial and Radial strains (a)  Axial strain (%) R ad ia l s tra in  (% ) 0 0.2 0.4 0.6 0.8 0.9 -0.2 -0.15 -0.1 -0.05 0 y=0.01335-0.1775x  Axial strain vs Radial strain (c)  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa 0 0.2 0.4 0.6 0.8 0.9 0 100 200 300 400 Peak Sd=333.7 MPa y=0.1751+408.3x  Axial stress vs Axial strain (b)  Time in seconds St ra in  (% ) 0 200 400 600 800 -0.4 0.0 0.4 0.8 1.2 Axial Radial  Axial and Radial strain vs Time (d) Figure A. 1. Triaxial test results for sample LI-8.3.          179  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa -4 -3 -2 -1 0 1 2 0 100 200 300 400 Axial Radial  Axial stress vs Axial and Radial strains (a)  Axial strain (%) R ad ia l s tra in  (% ) 0 0.2 0.4 0.6 0.8 1 1.2 -0.5 -0.3 -0.1 0.1 y=0.053-0.1895x  Axial strain vs Radial strain (c)  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa 0 0.4 0.8 1.2 1.6 2 0 100 200 300 400 500 600 Peak Sd=425.8 MPa y=-8.112+413.9x  Axial stress vs Axial strain (b)  Time in seconds St ra in  (% ) 0 400 800 1200 1600 -1.2 -0.4 0.4 1.2 AxialRadial  Axial and Radial strain vs Time (d) Figure A. 2. Triaxial test results for sample LI-8.2.          180  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 100 200 300 400 Axial Radial  Axial stress vs Axial and Radial strains (a)  Axial strain (%) R ad ia l s tra in  (% ) 0 0.2 0.4 0.6 0.8 1 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 y=0.01652-0.1783x  Axial strain vs Radial strain (c)  Strain (%) A xi al  D ev ia to r s tre ss , S d, M Pa 0 0.2 0.4 0.6 0.8 1 0 100 200 300 400 500 Peak Sd=436.2 MPa y=14.66+457.3x  Axial stress vs Axial strain (b)  Time in seconds St ra in  (% ) 0 400 800 1200 -0.5 -0.1 0.3 0.7 1.1 Axial Radial  Axial and Radial strain vs Time (d) Figure A. 3. Triaxial test results for sample LI-8.1.        181 A.2 Uniaxial Compression Strength (UCS) Test Results  Sample ID: LI-10.1 Test Date: May 17th, 2007 Tested by: D. Welkner Failure Mode: Shear Young Modulus, E (GPa)  42.31 Poisson Ratio, ν 0.17 Diameter, (φ) Area, (A) Height, (h) Ratio Peak Load (mm) (mm2) (mm) h/φ (kN) (MPa) (psi) 46.42 1692.4 92.64 2.0 164.5 97.2 14,096.2 σUCS Stress-Strain Curves 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 -0.0002 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 0.0016 0.0018 0.0020 Strain ( ) St re ss  (M Pa ) Axial Circumferential Unconfined Compressive Strength Test 0.0 20.0 40.0 60.0 80.0 100.0 120.0 0 100 200 300 400 500 600 time (s) A xi al  S tr es s (M Pa ) Max Axial Stress = 97.2MPa (14,096.2psi)  Figure A. 4. Uniaxial Compressive Test Results for Sample LI-10.1.     182 Sample ID: LI-10.2 Test Date: May 17th, 2007 Tested by: D. Welkner Failure Mode: Shear Young Modulus, E (GPa)  42.29 Poisson Ratio, ν 0.18 Diameter, (φ) Area, (A) Height, (h) Ratio Peak Load (mm) (mm2) (mm) h/φ (kN) (MPa) (psi) 46.43 1693.1 82.71 1.8 165.0 97.4 14,127.8 σUCS Stress-Strain Curves 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 0.0016 0.0018 0.0020 Strain ( ) St re ss  (M Pa ) Axial Circumferential Unconfined Compressive Strength Test 0.0 20.0 40.0 60.0 80.0 100.0 120.0 0 100 200 300 400 500 600 time (s) A xi al  S tr es s (M Pa ) Max Axial Stress = 97.4MPa (14,127.8psi)  Figure A. 5. Uniaxial Compressive Test Results for Sample LI-10.2.         183 Sample ID: LI-7 Test Date: May 17th, 2007 Tested by: D. Welkner Failure Mode: Shear Young Modulus, E (GPa)  42.08 Poisson Ratio, ν 0.30 Diameter, (φ) Area, (A) Height, (h) Ratio Peak Load (mm) (mm2) (mm) h/φ (kN) (MPa) (psi) 46.47 1696.0 88.53 1.9 409.9 241.7 35,043.5 σUCS Stress-Strain Curves 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 Strain ( ) St re ss  (M Pa ) Axial Circumferential Unconfined Compressive Strength Test 0.0 50.0 100.0 150.0 200.0 250.0 300.0 0 200 400 600 800 1000 1200 1400 1600 time (s) A xi al  S tr es s (M Pa ) Max Axial Stress = 241.7MPa (35,043.5psi)  Figure A. 6. Uniaxial Compressive Test Results for Sample LI-7.         184 Sample ID: LI-7.1 Test Date: May 17th, 2007 Tested by: D. Welkner Failure Mode: Shear Young Modulus, E (GPa)  55.61 Poisson Ratio, ν 0.34 Diameter, (φ) Area, (A) Height, (h) Ratio Peak Load (mm) (mm2) (mm) h/φ (kN) (MPa) (psi) 46.42 1692.4 91.39 2.0 436.9 258.2 37,435.7 σUCS Stress-Strain Curves 0.0 20.0 40.0 60.0 80.0 100.0 120.0 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 0.0016 0.0018 0.0020 Strain ( ) St re ss  (M Pa ) Axial Circumferential Unconfined Compressive Strength Test 0.0 50.0 100.0 150.0 200.0 250.0 300.0 0 200 400 600 800 1000 1200 1400 1600 time (s) A xi al  S tr es s (M Pa ) Max Axial Stress = 258.2MPa (37,435.7psi)  Figure A. 7. Uniaxial Compressive Test Results for Sample LI-7.1.         185 A.3 Rock Direct Shear Test Results       Figure A. 8. Direct Shear Test results for sample LI-7.1 at 16.25 Kg (top), 32.5 Kg (middle), and 48.5 Kg (bottom) seating normal load.  186       Figure A. 9. Direct Shear Test results for sample LI-7.2 at 16.25 Kg (top), 32.5 Kg (middle), and 48.5 Kg (bottom) seating normal load.  187                APPENDIX B: ROCLAB OUTPUTS                   188 B.1 Domain 1 Upper and Lower Bound RocLab Outputs  Hoek Brown Classification Hoek Brown Classification sigci 250 MPa sigci 90 MPa GSI 70 GSI 50 mi 11.312 mi 11.312 D 0 D 0 Ei 55000 Ei 40000 Hoek Brown Criterion Hoek Brown Criterion mb 3.87457 mb 1.89677 s 0.035674 s 0.003866 a 0.501355 a 0.505734 Failure Envelope Range Failure Envelope Range Application Slopes Application Slopes sig3max 21.2777 MPa sig3max 18.6111 MPa Unit Weight 0.027 MN/m3 Unit Weight 0.027 MN/m3 Slope Height 1000 m Slope Height 1000 m Mohr-Coulomb Fit Mohr-Coulomb Fit c 10.4798 MPa c 4.11619 MPa phi 46.0029 degrees phi 33.1332 degrees Rock Mass Parameters Rock Mass Parameters sigt -2.3018 MPa sigt -0.183435 MPa sigc 47.0061 MPa sigc 5.42045 MPa sigcm 73.662 MPa sigcm 16.638 MPa Erm 40304.8 MPa Erm 12287.4 MPa K 22391.56 MPa K 10239.5 MPa G 16793.67 MPa G 4725.923 MPa STRUCTURAL DOMAIN 1 Lower BoundUpper Bound             189   Figure B. 1. Principal stress and Normal vs. Shear Stress plots for D1 upper bound properties (top) and D1 lower bound properties (bottom).  190 Hoek Brown Classification Hoek Brown Classification sigci 60 MPa sigci 30 MPa GSI 50 GSI 35 mi 11 mi 11.3 D 0 D 0 Ei 30000 Ei 9000 Hoek Brown Criterion Hoek Brown Criterion mb 1.84445 mb 1.10891 s 0.003866 s 0.00073 a 0.505734 a 0.51595 Failure Envelope Range Failure Envelope Range Application Slopes Application Slopes sig3max 17.3176 MPa sig3max 15.8057 MPa Unit Weight 0.026 MN/m3 Unit Weight 0.026 MN/m3 Slope Height 1000 m Slope Height 1000 m Mohr-Coulomb Fit Mohr-Coulomb Fit c 3.37411 MPa c 1.96924 MPa phi 30.1525 degrees phi 21.6196 degrees Rock Mass Parameters Rock Mass Parameters sigt -0.125758 MPa sigt -0.019754 MPa sigc 3.61363 MPa sigc 0.72245 MPa sigcm 10.9449 MPa sigcm 3.9665 MPa Erm 9215.58 MPa Erm 1020.66 MPa K 5119.767 MPa K 850.55 MPa G 3839.825 MPa G 392.5615 MPa Lower BoundUpper Bound STRUCTURAL DOMAIN 2                191    Figure B. 2. Principal stress and Normal vs. Shear Stress plots for D2 upper bound properties (top) and D2 lower bound properties (bottom).  192               APPENDIX C: NUMERICAL MODELLING CODES               193 C.1 UDEC code for Base Model – Back analysis  ;;Input File for UDEC Base Model ;;Author: Daniela Welkner  ***************************************************** ****************problem geometry**************** *****************************************************  new ro 0.3 block 0,2092 0,3539 799,4048 1055,3500 1519,2942 1875,2740 2900,2500 2900,2092 crack 0,2259 1365,3128          ;fault crack 483.26,2092 1676.725,2852.487        ;closure zone 3 crack 850,2092 1875,2740          ;grad mesh reduction  ***************************************************** ****************region definitions****************** *****************************************************  jregion id=1 0,2828 0,3539 799,4048 1055,3500 delete     ;upper slope 65 jregion id=2 0,2259 0,2828 1055,3500 1365,3128 delete     ;upper slope 50 jregion id=3 0,1782.543 0,2259 1365,3128 1676.725,2852.487 delete  ;fault zone jregion id=4 483.26,2092 1676.725,2852.487 2900,2500 2900,2092 delete ;base jregion id=5 0,2259 0,3539 799,4048 1365,3128 delete     ;zone 1&2  ***************************************************** ****************joint sets*************************** *****************************************************  jset -65,0 650,0 0,0 30,0 1055,3500 range jreg=1 jset -50,0 550,0 0,0 30,0 1055,3500 range jreg=2 jset 32.5,0 1750,0 0,0 20,0 1365,3128 range jreg=3 jset 32.5,0 1650,0 0,0 40,0 1365,3128 range jreg=5     ;fault jset 25,0 1750,0 0,0 20,0 1365,3128 range jreg=3 jdelete delete range area 10         ;delete small blocks  ***************************************************** ****************zoning****************************** *****************************************************  gen quad 50 gen edge 100    194 ***************************************************** ****************assigning materials*************** *****************************************************  change mat=2 cons=3 range jreg=1       ;upper slope 65 (M-C) change mat=2 cons=3 range jreg=2       ;upper slope 50 (M-C) change mat=3 cons=3 range jreg=3       ;fault zone (M-C)  change jmat=3 jcons=2 range jreg=1       ;cross joints change jmat=3 jcons=2 range jreg=2       ;cross joints change jmat=2 jcons=2 range angle -68 -48      ;bedding 65-50 change jmat=4 jcons=2 range jreg=3             ;fault 2 change jmat=1 jcons=2 range jreg=4       ;artificial cracks  ***************************************************** ****************material properties**************** *****************************************************  ; upper slope (mohr-coulomb) prop mat=2 dens=2700 k=20e9 g=12e9 c=10e6 f=40 ten=10e6 prop jmat=2 jkn=1e10 jks=1e9 jfric=40.0 jcoh=5e6 prop jmat=3 jkn=1e10 jks=1e9 jfric=40.0 jcoh=5e6  ; fault zone (mohr-coulomb) prop mat=3 dens=2600 k=20e9 g=12e9 c=10e6 f=40 ten=10e6 prop jmat=4 jkn=1e10 jks=1e9 jfric=40.0 jcoh=5e6  ; base (elastic) prop mat=1 dens=2700 k=20e9 g=12e9 prop jmat=1 jkn=1e10 jks=1e9 jfric=40.0 jcoh=5e6  ***************************************************** **************** boundary conditions************* *****************************************************  bound xvel=0 range -1,1 2091,3540              ;left border bound xvel=0 range 2899,2901 2091,2501         ;right border bound yvel=0 range -1,2901 2091,2093           ;bottom insitu stress -214e6 0 -107e6 ygrad 5.28e4 0 2.64e4 & szz -107e6 zgrad 0 2.64e4 grav 0 -9.8  ***************************************************** **************** Initial equilibrium***************** ***************************************************** his unbal            ;his 1 solve save banalysis.sav     195 **************************************************************************************** ****************Lower strength prop for mat and jsets (Stage 1)****************** ****************************************************************************************  new restore banalysis.sav reset disp jdisp  ;slope prop mat=2 dens=2700 k=8e9 g=5e9 c=10e6 f=30 ten=2e5 prop jmat=2 jkn=1e10 jks=1e9 jfric=30 jcoh=1e6 prop jmat=3 jkn=1e10 jks=1e9 jfric=30 jcoh=5e4  ;fault zone prop mat=3 dens=2600 k=3e9 g=1e9 c=8e6 f=25 ten=5e4 prop jmat=4 jkn=5e9 jks=1e9 jfric=25 jcoh=0  ;base prop mat=1 dens=2700 k=2e10 g=1.5e10  solve save BMicn.sav  **************************************************************************************** ****************Lower strength prop for mat and jsets (Stage 2)****************** ****************************************************************************************  new restore BMicn.sav  ;lower 2 orders of magnitude jcoh of bedding  ;slope prop mat=2 dens=2700 k=8e9 g=5e9 c=8e6 f=30 ten=1e5 prop jmat=2 jkn=1e10 jks=1e9 jfric=30 jcoh=1e4  ;fault zone prop mat=3 dens=2600 k=3e9 g=1e9 c=6e6 f=25 ten=5e4  solve save BM1cn.sav  **************************************************************************************** ****************Lower strength prop for mat and jsets (Stage 3)****************** ****************************************************************************************  new restore BM1cn.sav  ;lower 2 orders of magnitude jcoh of bedding   196 ;slope prop mat=2 dens=2700 k=8e9 g=5e9 c=6e6 f=30 ten=1e5 prop jmat=2 jkn=1e10 jks=1e9 jfric=30 jcoh=1e2  ;fault zone prop mat=3 dens=2600 k=3e9 g=1e9 c=4e6 f=25 ten=5e4  solve save BM2cn.sav  **************************************************************************************** ****************Lower strength prop for mat and jsets (Stage 4)****************** ****************************************************************************************  new restore BM2cn.sav ;no jcoh of bedding  ;slope prop mat=2 dens=2700 k=8e9 g=5e9 c=4e6 f=30 ten=1e5 prop jmat=2 jkn=1e10 jks=1e9 jfric=30 jcoh=0  ;fault zone prop mat=3 dens=2600 k=3e9 g=1e9 c=2e6 f=25 ten=5e4  ; histories his unbal              ;his 2 his xdisp 1500,2950   ;his 3 (toe) his ydisp 860,3870    ;his 4 (crest) step 200000 save BM3cn.sav                197 C.2 UDEC code for coupled hydro-mechanical analysis – Forward Modelling  ;;Input File for UDEC Coupled hydro-mechanical analysis ;;Author: Daniela Welkner  ********************************************************************************************* ****************restore stable present-day slope and input water table**************** **********************************************************************************************  new restore fw_bmprop_new.sav reset disp jdisp hist set delc off fluid dens=1000 set flow steady  ***************************************************** ****************joint flow properties ************** *****************************************************  prop jmat=1 jperm=1e2 azero=5e-3 ares=2e-3 prop jmat=2 jperm=1e2 azero=5e-3 ares=2e-3 prop jmat=3 jperm=1e2 azero=5e-4 ares=2e-4 prop jmat=4 jperm=1e2 azero=5e-3 ares=2e-3  ***************************************************** **************** fluid flow boundary conditions ** *****************************************************  bound imperm range -1,3101 2199,2201     ;impermeable base bound pp=40e6 pygrad=-1e4 range -1,1 2199,4000 bound pp=28e6 pygrad=-1e4 range 3099,3101 2199,2800  *********************************************************************************** **************** groundwater table - fully saturated (100%)****************** ************************************************************************************  table 1 0,4000 250,3980 500,3940 750,3900 974.747,3758.647 1040,3545 & 1453.523,3085.555 1767.178,2993.203 1850,2950 2160,2870 2557.213,2800 3100,2800 *********************************************************************************** **************** initialization of fluid pressure ****************** ************************************************************************************  insitu stress -214e6 0 -107e6 ygrad 5.28e4 0 2.64e4 szz -107e6 zgrad 0 2.64e4 & ywtable table 1  ;histories  198 hist unbal             ;his 1 hist xdisp 1540,3050 ;his 2 (toe mid) hist pp  1540,3050 ;his 3 (toe mid) hist ydisp 1540,3050 ;his 4 (toe mid) hist xdisp 900,3840   ;his 5 (crest top) hist pp 900,3840    ;his 6 (crest top) hist ydisp 900,3840   ;his 7 (crest top) step 200000 save fw_100water200k_impxjts.sav  C.3 UDEC code for dynamic loading analysis – Forward Modelling  ********************************************************************************************************* ;;Input File for UDEC Base Model – present day geometry ;;Author: Daniela Welkner *********************************************************************************************************  new ro 0.5  ; problem geometry block 0,2200 0,4000 830,4000 974.747,3758.647 1040,3545 1453.523,3085.555 1767.178,2993.203 1850,2950 3100,2950 3100,2200  crack 0,2816.784 1040,3545        ;hinge crack 0,2243.944 1355.985,3193.440      ;zone 2&3 crack 634.552,2200 1767.178,2993.203      ;zone 3&4  jregion id=1 0,2816.784 0,4000 906.054,4000.1616 1040,3545 delete   ;US 73 jregion id=2 0,2243.944 0,2816.784 1040,3545 1355.985,3193.440 delete  ;US 48 jregion id=3 634.552,2200 -62.759,2200 1355.985,3193.440 1767.178,2993.203 delete ;FZ jregion id=4 634.552,2200 1767.178,2993.203 3100.0,2950.0 3100.0,2200.0  delete;base jregion id=5 0,2243.944 0,4000 830,4000 1355.985,3193.440 delete  ;zone 1&2  jset -73,0 1000,0 0,0 30,0 1040,3545 range jreg=1 jset -48,0 1000,0 0,0 30,0 1040,3545 range jreg=2 jset 35,0 2100,0 0,0 20,0 1355.985,3193.440 range jreg=3 jset 35,0 2000,0 0,0 40,0 1040,3545 range jreg=5 jset 20,0 2100,0 0,0 20,0 1355.985,3193.440 range jreg=3 jset 35,0 1000,0 0,0 100,0 1767.178,2993.203 range jreg=4 jset -24,0 1200,0 0,0 300,0 1767.178,2993.203 range jreg=4 jdelete delete range area 15     ;delete small blocks  ; meshing gen quad 50  199 gen edge 100         ;triang gen  ; assigning material change mat=2 cons=3 range jreg=1       ;US 73 (M-C) change mat=2 cons=3 range jreg=2       ;US 48 (M-C) change mat=3 cons=3 range jreg=3       ;FZ (M-C)  ; assigning joints change jmat=3 jcons=2 range jreg=1      ;US xjts change jmat=3 jcons=2 range jreg=2      ;US xjts change jmat=2 jcons=2 range angle -75 -45     ;US bedding change jmat=4 jcons=2 range jreg=3             ;FZ change jmat=1 jcons=2 range jreg=4             ;base  ; material properties  ; upper slope (mohr-coulomb) prop mat=2 dens=2700 k=15e9 g=9e9 c=50e6 f=40 ten=10e6 prop jmat=2 jkn=1e10 jks=1e9 jfric=50.0 jcoh=5e6    ;US bedding prop jmat=3 jkn=1e10 jks=1e9 jfric=50.0 jcoh=5e6    ;US xjts  ; fault zone (mohr-coulomb) prop mat=3 dens=2700 k=1.6e9 g=0.7e9 c=50e6 f=40 ten=10e6 prop jmat=4 jkn=1e10 jks=1e9 jfric=50.0 jcoh=5e6    ;FZ xjts & bedding  ; base (elastic) prop mat=1 dens=2700 k=15e9 g=9e9 prop jmat=1 jkn=2.5e10 jks=1e10 jfric=50.0 jcoh=1e6   ;base jts  ; boundary conditions (rollers) insitu stress -214e6 0 -107e6 ygrad 5.28e4 0 2.64e4 & szz -107e6 zgrad 0 2.64e4 grav 0 -9.8 bound xvel=0 range -1,1 2199,4001             ;left border bound xvel=0 range 3099,3101 2199,2951        ;right border bound yvel=0 range -1,3101 2199,2201          ;bottom  ; histories his unbal                     ;his 1  solve elastic save fw_ini_EQ.sav           200 ********************************************************************************************************* ;;Change to Base Model properties – present day geometry ;;Author: Daniela Welkner *********************************************************************************************************  new restore fw_ini_EQ.sav reset disp jdisp his  ;base model properties ;slope prop mat=2 dens=2700 k=8e9 g=5e9 c=4e6 f=30 ten=2e5 prop jmat=2 jkn=1e10 jks=1e9 jfric=30 jcoh=0 prop jmat=3 jkn=1e10 jks=1e9 jfric=30 jcoh=5e4  ;fault zone prop mat=3 dens=2600 k=3e9 g=1e9 c=2e6 f=25 ten=5e4 prop jmat=4 jkn=5e9 jks=1e9 jfric=25 jcoh=0  ; histories his unbal              ;his 2 his xdisp 1380,3150  ;his 3 (toe up) his xdisp 1540,3050  ;his 4 (toe mid) his xdisp 1710,3000   ;his 5 (toe down) his ydisp 900,3840    ;his 6 (crest top) his ydisp 1170,3370    ;his 7 (below hinge) his yvel 900,3840    ;his 8 (crest top) his xvel 1380,3150  ;his 9 (toe up) step 300000 save fw_base_EQ.sav  ********************************************************************************************************* ;;Input free-field boundary conditions ;;Author: Daniela Welkner *********************************************************************************************************  new restore fw_base_EQ.sav  ;generate free-field both lateral boundaries & fixed bottom ffield gen left y 2200,4000 np=1800 ffield gen right y 2200,2950 np=750 ffield change mat=2 cons=3  ;initialize FF stresses ffield ini sxx -214e6 5.28e4 ffield ini syy -107e6 2.64e4 ffield ini szz -107e6 2.64e4  ;cycle with FF not attached to model for equilibrium ffield base xvel=0  201 ffield base yvel=0 reset time hist hist ffyd 4000 1 hist ffsxx 4000 1  step 100000 save ffield_ini.sav  ********************************************************************************************************* ;;Input a seismic wave of 0.8 MPa, 2 Hz, T=2 s & Raleigh damping ;;Author: Daniela Welkner *********************************************************************************************************  new restore ffield_ini.sav  ;apply dynamic boundary condition bound mat=2 bound ff range -1,1 2199,4001 bound ff range 3099,3101 2199,2951 bound xvisc range -1,3101 2199,2201  ;amplitude of shear wave: 0.8 MPa, freq:2 Hz, T=2 sec bound stress 0 0.8e6 0 hist sine(2,2) range -1,3101 2199,2201  ;fix yvel at bottom bound yvel=0 range -1,3101 2199,2201  ;apply free-field boundary conditions ffield base sxy=0.8e6 hist sine (2,2) ffield base yvel=0 ffield base xvisc  reset hist time disp hist xvel 950,3760  yvel 900,3840    ;concave edge of US hist xdis 950,3760  ydis 900,3840    ;concave edge of US hist xvel 1022,3556 yvel 1022,3556   ;hinge hist xdis 950,3760  ydis 950,3760 ;hinge hist ffxd 4000,1 ffxd 2950,2 hist ffxvel 4000,1 ffyvel 4000,1 ffxvel 2200,2 hist xdisp 1380,3151 ydisp 1380,3151  ;damping conditions damp 0.02 10 cycle time 10 save fw_EQ_ff8.sav       202 ********************************************************************************************************* ;;Input a seismic wave of 7 MPa, 5 Hz, T=40 s & Raleigh damping ;;Author: Daniela Welkner *********************************************************************************************************  new restore ffield_ini.sav  ;apply dynamic boundary condition bound mat=2 bound ff range -1,1 2199,4001 bound ff range 3099,3101 2199,2951 bound xvisc range -1,3101 2199,2201  ;amplitude of shear wave: 7 MPa, freq:5 Hz, T=40 sec bound stress 0 7e6 0 hist sine(5,40) range -1,3101 2199,2201  ;fix yvel at bottom bound yvel=0 range -1,3101 2199,2201  ;apply free-field boundary conditions ffield base sxy=7e6 hist sine (5,40) ffield base yvel=0 ffield base xvisc  reset hist time disp hist xvel 950,3760  yvel 900,3840 ;concave edge of US hist xdis 950,3760  ydis 900,3840 ;concave edge of US hist xvel 1022,3556 yvel 1022,3556 ;hinge hist xdis 950,3760  ydis 950,3760 ;hinge hist ffxd 4000,1 ffxd 2950,2 hist ffxvel 4000,1 ffyvel 4000,1 ffxvel 2200,2 hist xdisp 1380,3151 ydisp 1380,3151  ;damping conditions damp 0.02 15 cycle time 45 save fw_EQ_ff15.sav

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 8 0
United States 3 0
Japan 2 0
United Arab Emirates 2 0
Australia 1 0
Slovenia 1 0
Serbia 1 0
New Zealand 1 0
City Views Downloads
Beijing 6 0
Unknown 4 2
Mountain View 2 0
Tokyo 2 0
Shenzhen 1 0
Maribor 1 0
Wantirna South 1 0
Changsha 1 0
Ashburn 1 0

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

Share

Share to:

Comment

Related Items