"Science, Faculty of"@en . "Earth, Ocean and Atmospheric Sciences, Department of"@en . "DSpace"@en . "UBCV"@en . "McDougall, Scott"@en . "2010-01-16T21:22:35Z"@en . "2006"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "Runout analysis, the prediction of landslide motion and its effects, is an essential component\r\nof landslide risk assessment. A new continuum dynamic model has been developed for the\r\nrunout analysis of extremely rapid, flow-like landslides, including rock avalanches, debris\r\navalanches, debris flows and flow slides. The new model, DAN3D, is a 3D extension of the\r\nexisting 2D model DAN. It uses a meshless, Lagrangian numerical method adapted from\r\nSmoothed Particle Hydrodynamics to discretize and solve the depth-averaged equations of\r\nmotion for an \"equivalent fluid\", a hypothetical material governed by simple rheological\r\nrelationships. The required rheological parameters, rather than measured, are calibrated\r\nthrough back-analysis of real landslide case studies. The key capabilities of the model\r\ninclude: 1) the ability to simulate motion across complex 3D terrain without the need to input\r\na pre-defined path direction or width and without introducing problems due to mesh\r\ndistortion; 2) the ability to simulate strain-dependent, non-hydrostatic, anisotropic internal\r\nstresses due to 3D deformation of material with internal shear strength; 3) the ability to\r\nsimulate mass and momentum transfer due to entrainment of path material; 4) the ability to\r\nsimulate variations in rheology along the path and within the landslide; and 5) efficient and\r\nsimple operation. The model outputs the simulated spatial distribution of hazard intensity\r\nparameters, including flow velocity and depth, which are required for delineating the\r\npotential impact area, estimating the vulnerability of elements within this area and designing\r\nprotective measures. It has been tested using both analytical and experimental methods, and\r\nits general behaviour has been demonstrated using a series of simple parametric analyses.\r\nThe model has also been applied at full-scale to the simulation of a wide variety of landslide\r\ntypes. These back-analyses form the basis for a more thorough calibration, but some useful\r\npatterns have already emerged. This experience has been applied in practice to landslide\r\nrunout prediction, although so far only in a parametric way. With continued back-analysis of\r\nreal cases and the development of a probabilistic approach to forward-analysis, true landslide\r\nrunout prediction, including quantification of uncertainty, should eventually be possible\r\nusing the new model."@en . "https://circle.library.ubc.ca/rest/handle/2429/18500?expand=metadata"@en . "A N E W CONTINUUM D Y N A M I C M O D E L FOR THE A N A L Y S I S OF E X T R E M E L Y RAPID LANDSLIDE MOTION ACROSS C O M P L E X 3D TERRAIN by SCOTT M C D O U G A L L B A . S c , University of Toronto, 1998 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE F A C U L T Y OF G R A D U A T E STUDIES (Geological Engineering) THE UNIVERSITY OF BRITISH C O L U M B I A August 2006 \u00C2\u00A9 Scott McDougall, 2006 ABSTRACT Runout analysis, the prediction of landslide motion and its effects, is an essential component of landslide risk assessment. A new continuum dynamic model has been developed for the runout analysis of extremely rapid, flow-like landslides, including rock avalanches, debris avalanches, debris flows and flow slides. The new model, DAN3D, is a 3D extension of the existing 2D model D A N . It uses a meshless, Lagrangian numerical method adapted from Smoothed Particle Hydrodynamics to discretize and solve the depth-averaged equations of motion for an \"equivalent fluid\", a hypothetical material governed by simple rheological relationships. The required rheological parameters, rather than measured, are calibrated through back-analysis of real landslide case studies. The key capabilities of the model include: 1) the ability to simulate motion across complex 3D terrain without the need to input a pre-defined path direction or width and without introducing problems due to mesh distortion; 2) the ability to simulate strain-dependent, non-hydrostatic, anisotropic internal stresses due to 3D deformation of material with internal shear strength; 3) the ability to simulate mass and momentum transfer due to entrainment of path material; 4) the ability to simulate variations in rheology along the path and within the landslide; and 5) efficient and simple operation. The model outputs the simulated spatial distribution of hazard intensity parameters, including flow velocity and depth, which are required for delineating the potential impact area, estimating the vulnerability of elements within this area and designing protective measures. It has been tested using both analytical and experimental methods, and its general behaviour has been demonstrated using a series of simple parametric analyses. The model has also been applied at full-scale to the simulation of a wide variety of landslide types. These back-analyses form the basis for a more thorough calibration, but some useful patterns have already emerged. This experience has been applied in practice to landslide runout prediction, although so far only in a parametric way. With continued back-analysis of real cases and the development of a probabilistic approach to forward-analysis, true landslide runout prediction, including quantification of uncertainty, should eventually be possible using the new model. i i TABLE OF CONTENTS A B S T R A C T ii T A B L E OF CONTENTS i i i LIST OF T A B L E S viii LIST OF FIGURES ix LIST OF S Y M B O L S xii A C K N O W L E D G E M E N T S xiv DEDICATION xv 1 INTRODUCTION 1 1.1 Introduction 1 1.2 Purpose, scope and structure of this thesis 4 1.3 Runout analysis in the framework of landslide risk assessment 6 1.3.1 Estimating the spatial probability of impact 7 1.3.2 Estimating vulnerability 8 1.3.3 Designing protective measures 9 1.4 Continuum dynamic modelling as a runout analysis method 10 1.4.1 Empirical methods 11 1.4.2 Analytical methods 13 1.5 Discussion 16 1.6 Conclusion 17 2 OBJECTIVES FOR M O D E L D E V E L O P M E N T 18 2.1 Introduction 18 2.2 Objective #1: the ability to simulate motion across complex 3D terrain 19 2.3 Objective #2: the ability to simulate the influence of internal strength 21 2.3.1 Rankine earth pressure theory 21 2.3.2 Strain-dependent internal stresses in landslides 25 2.4 Objective #3: the ability to simulate entrainment 28 2.4.1 Mechanisms of path material mobilization 29 2.4.2 Volume change 33 2.4.3 Momentum transfer 35 2.5 Objective #4: the ability to simulate variations in rheology 38 2.5.1 Variations along the path 38 2.5.2 Variations within the landslide 41 2.6 Objective #5: efficiency 42 2.7 Discussion 43 2.8 Conclusion 44 ii i 3 PREVIOUS W O R K 45 3.1 Introduction 45 3.2 Model classification 45 3.2.1 2D vs. 3D 45 3.2.2 Eulerian vs. Lagrangian 46 3.3.3 Measurement vs. calibration.. 47 3.3 Existing continuum dynamic models 49 3.3.1 Extension of hydrodynamic methods 49 3.3.2 Incorporation of path-dependent rheology 50 3.3.3 Incorporation of earth pressure theory 51 3.3.4 Incorporation of entrainment capabilities 53 3.3.5 Towards a synthesis 54 3.4 Discussion 57 3.5 Conclusion 58 4 GOVERNING EQUATIONS 59 4.1 Introduction 59 4.2 Equivalent fluid concept 59 4.3 General conservation laws... 61 4.4 Boundary conditions.. 64 4.4.1 Stress conditions 65 4.4.2 Kinematic conditions 65 4.5 Depth-averaging 66 4.5.1 Depth-averaging the mass balance equation 67 4.5.2 Depth-averaging the momentum balance equations 67 4.6 Further simplifications 71 4.6.1 Lagrangian reference frame 71 4.6.2 Shallow flow assumption 73 4.6.3 Stress state normalization 74 4.7 Basal shear resistance 76 4.7.1 Laminar 77 4.7.2 Turbulent 77 4.7.3 Plastic 78 4.7.4 Bingham 78 4.7.5 Frictional 78 4.7.6 Voellmy 80 4.8 Discussion 81 4.9 Conclusion 82 5 N U M E R I C A L SOLUTION METHOD 84 5.1 Introduction 84 5.2 Smoothed Particle Hydrodynamics 85 5.2.1 Background 85 5.2.2 Theory 86 5.3 SPH-based mass balance 88 5.3.1 Depth and depth gradient interpolation 88 iv 5.3.2 Interpolating kernel 91 5.4 Strain interpolation 95 5.4.1 Incremental tangential strain state 95 5.4.2 Strain rosette method 96 5.4.3 Simplified strain rosette method 99 5.5 Stress update 100 5.5.1 Assumed orientation of principle axes 100 5.5.2 Internal yield criterion 101 5.5.3 Stress incrementation 104 5.6 Entrainment 106 5.6.1 Momentum transfer based on erosion velocity and growth rate 106 5.6.2 Mass transfer 109 5.7 Intensity parameter interpolation 111 5.8 Implementation \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 112 5.8.1 Initial conditions 114 5.8.2 Time stepping 115 5.9 Discussion 118 5.10 Conclusion 119 6 TESTS.-. .' 120 6.1 Introduction 120 6.2 ID dam-break 120 6.2.1 Description 120 6.2.2 Methodology 121 6.2.3 Results and discussion 123 6.3 Slump test \u00E2\u0080\u00A2 125 6.3.1 Description 125 6.3.2 Methodology 125 6.3.3 Results and discussion 127 6.4 Parametric analyses 129 6.4.1 Description 129 6.4.2 Methodology 129 6.4.3 Results and discussion: # particles and particle smoothing 131 6.4.4 Results and discussion: internal friction 133 6.4.5 Results and discussion: internal stiffness 134 6.4.6 Results and discussion: velocity smoothing 136 6.4.7 Results and discussion: basal rheology 137 6.4.8 Results and discussion: entrainment 138 6.5 Channelization 155 6.5.1 Description :......: 155 6.5.2 Methodology 155 6.5.3 Results and discussion 157 6.6 Chute experiments 158 6.6.1 Description 158 6.6.2 Methodology 158 6.6.3 Results and discussion 159 v 6.7 Deflection experiments 161 6.7.1 Description 161 6.7.2 Methodology 163 6.7.3 Results and discussion 163 6.8 Discussion 166 6.9 Conclusion 167 7 C A S E STUDIES 169 7.1 Introduction 169 7.2 Frank 169 7.2.1 Description 169 7.2.2 Methodology 172 7.2.3 Results and discussion 173 7.3 Val Pola 175 7.3.1 Description 175 7.3.2 Methodology 176 7.3.3 Results and discussion 178 7.4 Cervinara 181 7.4.1 Description 181 7.4.2 Methodology 181 7.4.3 Results and discussion 183 7.5 Quindici 185 7.5.1 Description 185 7.5.2 Methodology 185 7.5.3 Results and discussion 186 7.6 Nomash River 188 7.6.1 Description \u00E2\u0080\u00A2 188 7.6.2 Methodology 189 7.6.3 Results and discussion 191 7.7 Zymoetz River 195 7.7.1 Description 195 7.7.2 Methodology 205 7.7.3 Results and discussion 207 7.8 McAuley Creek 212 7.8.1 Description 212 7.8.2 Methodology 214 7.8.3 Results and discussion 215 7.9Kolka 217 7.9.1 Description 217 7.9.2 Methodology 219 7.9.3 Results and discussion 220 7.10 Other cases 225 7.11 Discussion 227 7.12 Conclusion 229 vi 8 C O N C L U S I O N S 230 8.1 Introduction 230 8.2 Summary of work completed 230 8.3 Recommendations for future work 234 8.4 Discussion 237 8.5 Conclusion 237 R E F E R E N C E S 238 vii LIST OF TABLES Table 1.1 Correlation coefficients for proposed volume-fahrboschung relationships 13 Table 6.1 Parameter values used in each parametric analysis 130 Table 6.2 Initial and final smoothing lengths corresponding to analyses (a) to (h) 133 Table 6.3 Limiting values of the stress coefficients in run (c) 134 Table 6.4 Configurations used in the deflected flow experiments 162 Table 6.5 Observed vs. predicted runup distances for Experiments 3, 5 and 7 166 Table 7.1 Parameter values used in each case study 228 viii LIST OF FIGURES Figure 1.1 Screen capture from a runout analysis using D A N - W 3 Figure 1.2 Flooding of Lake Saniba, Russia caused by a landslide dam 8 Figure 1.3 Hazard intensity map for Cheekeye Fan, British Columbia 9 Figure 1.4 Hazard intensity modified by deflection berms, Cheekye Fan 10 Figure 1.5 Schematic definition of the fahrboschung or angle of reach 12 Figure 1.6 The classical sliding block model 14 Figure 2.1 Branching and merging of debris avalanches/flows 20 Figure 2.2 Idealized loading conditions on soil adjacent to an embedded retaining wall.... 22 Figure 2.3 Schematic illustrations of idealized (a) active and (b) passive local failure 22 Figure 2.4 Mohr circles showing the two limiting stress conditions for the idealized loading shown in Figure 2.4 23 Figure 2.5 Theoretical range of horizontal earth pressure coefficients for a range of internal friction angles 25 Figure 2.6 Total stress state on an element of material within a landslide 27 Figure 2.7 Hypothetical Mohr circle representation of the 3D stress state shown in Figure 2.6 27 Figure 2.8 Schematic illustration showing entrainment of path material by both plowing at the margins and erosion/scour at the base 29 Figure 2.9 Aerial view of the 1965 Hope rock avalanche 30 Figure 2.10 Depth of bed instability predicted using Equations [2-4] and [2-5] 32 Figure 2.11 The Tsing Shan debris flow 34 Figure 2.12 Energy grade line analysis of a single sliding block 35 Figure 2.13 Energy grade line analysis of two sliding blocks 37 Figure 2.14 Sliding block model demonstrating the influence of a change in basal material strength along the path 39 Figure 2.15 Aerial view of McAuley Creek rock avalanche 40 Figure 2.16 The distal end of a debris flow deposit near Prince Rupert, British Columbia... 41 Figure 3.1 Eulerian vs. Lagrangian reference frames 47 Figure 4.1 Schematic illustration of the \"equivalent fluid\" approach 61 Figure 4.2 Generic orientation of the reference coordinate system 65 Figure 4.3 Orientation of the Lagrangian reference coordinate system used in DAN3D 72 Figure 4.4 Normalized total stress state on an element of material within a landslide 75 Figure 4.5 Relationships between frictional parameters 79 Figure 5.1 A physical interpretation of SPH 90 Figure 5.2 Comparison of different interpolating kernels 93 Figure 5.3 Plane strain measurement using strain gauge rosettes 96 Figure 5.4 Example of a DAN3D incremental tangential strain interpolation 98 Figure 5.5 Example of a DAN3D incremental tangential strain interpolation using the simplified method 100 Figure 5.6 Hypothetical Mohr circle representation of the normalized 3D stress state shown in Figure 4.4 102 Figure 5.7 Flowchart showing the main functions in DAN3D 113 Figure 5.8 Example of a DAN3D initial particle distribution 115 ix Figure 6.1 Analyses of the dam-break problem 124 Figure 6.2 Isometric view of initial geometry used in the slump test 126 Figure 6.3 Initial particle locations used in the slump test 126 Figure 6.4 Isometric view of interpolated initial geometry 127 Figure 6.5 Results of the slump test 128 Figure 6.6 Section through the final deposit 128 Figure 6.7 Initial configuration of the parametric analyses 131 Figure 6.8 Parametric analyses varying the number of particles and the particle smoothing coefficient 141 Figure 6.9 Parametric analyses varying the internal friction angle 145 Figure 6.10 Prevailing longitudinal and lateral stress states recorded at each time interval during run (c) 147 Figure 6.11 Parametric analyses varying the dimensionless stiffness coefficient 148 Figure 6.12 Parametric analyses varying the velocity smoothing coefficient 150 Figure 6.13 Parametric analyses varying the basal rheology and its parameters 152 Figure 6.14 Parametric analyses neglecting centripetal acceleration and varying the entrainment growth rate 154 Figure 6.15 Erosion map corresponding to run (b) shown in Figure 6.14b 155 Figure 6.16 Initial configuration of the hypothetical particle channelization tests 156 Figure 6.17 Particle positions in channel at time / = 30 s 157 Figure 6.18 Initial configuration of the tests based on experiments presented by Gray et al. (1999) and Wieland et al. (1999) 159 Figure 6.19 Simulation of the straight runout chute experiment presented by Gray et al. (1999) 160 Figure 6.20 Photograph of the experimental apparatus 162 Figure 6.21 Calibrated simulation of Experiment 1 164 Figure 6.22 Comparisons of observed and predicted maximum runup distances 165 Figure 7.1 Oblique aerial view of the Frank Slide 170 Figure 7.2 Calibrated simulation of the Frank Slide 174 Figure 7.3 Oblique aerial view of the Val Pola landslide deposit 176 Figure 7.4 Calibrated simulation of the Val Pola landslide 179 Figure 7.5 Plot of the maximum simulated flow velocities recorded along the path 180 Figure 7.6 View of the source slope of the Cervinara landslide 182 Figure 7.7 Simulation of the Cervinara landslide 184 Figure 7.8 Two simulations of the Quindici landslide 187 Figure 7.9 View of the fan apex of the Quindici landslide 188 Figure 7.10 Oblique aerial views of the Nomash River rock slide -debris avalanche 189 Figure 7.11 Simulation of the Nomash River landslide 192 Figure 7.12 Comparison of simulated and measured trimlines 193 Figure 7.13 Erosion map corresponding to the simulation shown in Figure 7.11 193 Figure 7.14 View of the source slope of the Nomash River landslide 194 Figure 7.15 Simulation of the Nomash River landslide neglecting entrainment 195 Figure 7.16 Location of the June 8, 2002 Zymoetz River rock slide - debris flow 196 Figure 7.17 Aerial view of the Zymoetz River rock slide - debris flow 196 Figure 7.18 Plot of landslide volume vs. the tangent of the fahrboschung 197 Figure 7.19 Enlarged aerial view of the proximal part of the path 198 x Figure 7.20 View from toe of source slope toward source area 200 Figure 7.21 Rock avalanche deposits on top of snow in the upper part of the cirque basin. 200 Figure 7.22 Oblique aerial view of the cirque basin 202 Figure 7.23 Enlarged aerial view of the distal part of the path 202 Figure 7.24 Oblique aerial view up the channel towards the source area 204 Figure 7.25 Oblique aerial view of the large fan created by the landslide 205 Figure 7.26 Calibrated simulation of the Zymoetz River landslide 208 Figure 7.27 Plot of the maximum simulated flow velocities recorded along the path 209 Figure 7.28 Simulated trimline and distribution of material after 600 s, superimposed on the aerial photograph 210 Figure 7.29 Aerial view of the McAuley Creek landslide 213 Figure 7.30 Simulation of the McAuley Creek landslide 216 Figure 7.31 View of the main deposit of the Kolka Glacier event 218 Figure 7.32 Calibrated simulation of the Kolka event 221 Figure 7.33 Matrix of maximum velocity plots 223 Figure 7.34 Comparison of arrival times along the path 224 xi LIST OF SYMBOLS A area b bed depth; maximum erosion depth; depth of instability B particle smoothing coefficient c constant shear strength C consequence; probability of loss given an occurrence C Chezy coefficient C velocity smoothing coefficient D stiffness coefficient E Young's elastic modulus Es entrainment growth rate; displacement-dependent rate of material entrainment Et erosion velocity; time-dependent rate of material entrainment at the base / Voellmy friction coefficient / value of a generic function g gravitational acceleration g gravitational acceleration vector h flow depth H difference in elevation between crest of source and toe of deposit / central column/particle identification number j neighbouring column/particle identification number k stress coefficient (kx, ky , k , kx_, k , kv_, k.x and k.v) k node identification number ka active earth pressure coefficient kp passive earth pressure coefficient I smoothing length L length of the horizontal projection of the streamline of motion m mass n Manning roughness coefficient n number of neighbouring particles within a local radius of influence N total number of particles P(L) annual frequency of landslide occurrence P(L0L) annual probability that a specific person will lose their life P(S:L) spatial probability of impact given a landslide occurrence P(T.S) temporal probability of impact given spatial impact R(prop) risk to property; annual loss of property value R bed-normal radius of curvature of the path in the direction of motion r position vector r2 correlation coefficient ru pore pressure ratio; ratio of pore fluid pressure to total normal stress xii s distance t time T stress tensor u pore fluid pressure V velocity V depth-averaged velocity V velocity vector V volume V (DT) vulnerability; the probability of death given impact V v(prop:T) vulnerability; the degree of loss given impact, on an increasing scale of 0 to 1 w interpolating kernel a inclination/angle of the bed from horizontal ad deflection plane true dip angle J3 angle between original and corrected particle velocity vectors 8 Dirac delta function 8 deflection angle A increment operator e normal strain

b bulk basal friction angle nternal friction angle r unit weight r engineering shear strain M dynamic viscosity /^Bingham Bingham viscosity e angle in local x - y plane counter-clockwise from positive x -axis P material density a total normal stress effective normal stress effective horizontal normal stress effective vertical normal stress X shear stress Tf shear strength tyield Bingham yield stress V Poisson's ratio Voellmy turbulence parameter momentum correction coefficient V gradient operator \u00E2\u0080\u00A2 dot product \u00C2\u00AE tensor product xiii ACKNOWLEDGEMENTS First, I have to thank my lovely wife and best friend, Samantha Ward, for her amazing support and patience. Thanks also to our families, who have been incredibly supportive of both of us and helped make this little adventure possible. I am very grateful to have had the opportunity to work closely with Oldrich Hungr, who has been a fantastic advisor. His technical yet practical approach to complex subjects has been particularly influential, and I think that influence is quite obvious in this thesis. He also provided generous funding for a number of conference trips, including to exotic places like Brazil, Nepal and Saskatchewan. Thanks also to my thesis and examining committee members, Dave McClung, Roger Beckie, Garry Clarke, Roland Stull and Scott Dunbar for their encouragement and for all of the tough questions and stimulating discussions that led to many improvements. Special thanks to my external examiner, Barry Voight, for a thorough and truly constructive review. I would also like to thank another great mentor, Erik Eberhardt, and my officemates, Jordan Severin, Lara Fletcher, Mike Scholz, Jonathan Mackin, Andrea Cade and Alex Strouth, for their input, friendship and general tolerance. Several other friends, colleagues and collaborators contributed to this work in a number of ways, including Steve Evans, Rejean Couture, Bob Gerath, Marina Pirulli, Claudio Scavia, Paola Revellino, Francesco Guadagno, Giovanni Crosta, Heike Willenburg, Suzanne Chalindar, Nikolai Hungr, Nichole Boultbee, Marc-Andre Brideau, Doug Stead, Jim Schwab, Marten Geertsema, Kevin Turner, Rick Guthrie, Dana Ayotte, Olga Tutubalina, Sergey Chernomorets, Dmitry Petrakov and Valeriy Drobyshev. Their contributions are gratefully acknowledged. Finally, this work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC). The Department of Earth and Ocean Sciences was instrumental in securing this and other funding, and also provided a terrific work environment. I would especially like to thank Alex Allen for all of her help over the past few years. xiv To Samantha xv 1 INTRODUCTION 1.1 Introduction Flow-like landslides such as rock avalanches, debris avalanches, debris flows and flow slides travel at extremely rapid velocities and can impact large areas, often far from their source (cf, Hungr et al. 2001 for landslide classification). When the source of a potentially mobile landslide is identified, stabilization is not always a practical option and the consequences of failure must be considered. The extent of the potential impact area, the distribution of intensity (i.e., destructiveness) within it and the possibility of secondary effects, such as landslide-generated displacement waves, must all be estimated. The prediction of landslide motion and its effects is called runout analysis. As an essential component of landslide risk assessment, it can provide valuable guidance for land use planning decisions and the development of mitigation strategies, including the design of protective structures. There is an increasing need for effective methods of performing landslide runout analysis, and this is the general motivation for this thesis. Worldwide, the annual cost of landslides, in terms of both life and economic losses, is rising, at least in part due to increasing exposure of people and property as development continues to push into marginally-stable land (Petley et al. 2005). These methods must be quantitative and as accurate, objective and accessible as possible. Quantification facilitates communication and rational decision making within a risk-based framework. Accuracy is required to minimize underestimation but also overestimation of risk, as decisions that are too conservative, particularly involving risk avoidance through land sterilization, can be extremely costly (e.g., Hungr 2004). Objectivity helps produce results that are repeatable and defensible, factors that promote widespread use and development of the methods in practice. In terms of accessibility, the methods must be made available to practitioners and should be as simple and as easy to understand as possible. Simplicity facilitates application of the methods as well as interpretation of the results, a step in the runout analysis process when subjective judgment plays an important role. 1 Continuum dynamic modelling of landslide motion has emerged as a method that may be able to satisfy all of these criteria. The spatial distribution of quantitative landslide hazard intensity parameters, such as flow velocity and depth (Hungr 1997), which are required to perform the calculations at each time step, can be easily output. Dynamic models can be developed into user-friendly software and experiences can be shared among users to aid in further development. Models with a minimal number of adjustable input parameters are simple and can be reasonably objective. Finally, dynamic modelling is potentially more accurate than traditional empirical methods because the unique geometry and mechanics of each case can be accounted for explicitly. As a result, a large number of dynamic models have been proposed or are currently in development. Most are based on established hydrodynamic theory, with some landslide-specific modifications. The model proposed by Hungr (1995), D A N , simulates motion along a user-prescribed 2D path and includes several unique features to account for the most important characteristics of extremely rapid landslides, including the effects of internal strength, material entrainment and variations in rheology both within the slide mass and along the path (Figure 1.1). It is based on a highly efficient Lagrangian solution of the depth-averaged equations of motion for an \"equivalent fluid\" (Hungr 1995), a hypothetical material governed by simple internal and basal rheologies. The equivalent fluid represents an idealization of the actual landslide material, which in reality may be very complex and difficult to model. D A N is a calibration-based model, which means that the appropriate rheological parameters must be constrained by trial-and-error back-analysis of previous real landslides. The trials are judged in terms of their ability to simulate the bulk characteristics of the prototype event, including the total travel distance, the distribution of deposits and the velocities estimated along the path. A large number of case.studies have been analyzed and a valuable database of calibrated parameters has been created (e.g., Hungr and Evans 1996; Ayotte and Hungr 2000; Hungr et al. 2002; Revellino et al. 2004; Pirulli et al. 2004). For landslides of similar type and scale, the calibrated parameters are generally well-constrained, suggesting that accurate, first-order runout prediction is possible. 2 g\u00C2\u00BB. DAH-W , Ocledie 4 P:\iimcdougan\Mp Doi:unrenU\ScnmPA\u00C2\u00BB3D\C 0.95 generally indicates a strong correlation). The sample stratification used by Corominas (1996), based on landslide type and path morphology, improved the strength of the correlations, but a significant amount of scatter remained. Similar sample stratification was proposed by Nicoletti and Sorriso-Valvo (1991). One of the problems with sample stratification is that it divides an already small database into even smaller parts. Empirical/statistical models for snow avalanches are generally much better constrained because a much larger database of case histories is typically available. With improved record keeping and an expansion of the landslide case history database, the reliability of empirical landslide runout analysis methods may improve. Other empirical methods already developed for snow avalanches (e.g., Lied and Bakkebeii 1980; McClung and Mears 1991; McClung 2000) could also be investigated. 12 Reference Proposed Relationship r2 Scheidegger (1973) log 1 0 (H/L) = -0.15666log 1 0 V + 0.62419 0.82 L i (1983) log l 0 (H/L) = -0.1529log l 0 V + 0.6640 0.7784 Corominas (1996) all landslides log 1 0 (H/L) = -0.085log l 0 V - 0.047 0.625 Corominas (1996) rock falls log 1 0 (H/L) = -0.109 log 1 0 V + 0.210 0.759 Corominas (1996) obstructed rock falls log 1 0 (H/L) = -0.091 log 1 0 V + 0.231 0.834 Corominas (1996) unobstructed rock falls log, 0 (H/L) = -0.119log 1 0 V + 0.167 0.924 Table 1.1. Correlation coefficients for some proposed volume-fahrboschung relationships (an example of weak correlation in empirical methods). H/L is the tangent of the fahrboschung (as shown in Figure 1.5) and V is the landslide volume. The proposed relationships and correlation coefficients are shown as originally published in the cited references, with no modifications to the number of significant figures. 1.4.2 Analytical methods In contrast to empirical methods, analytical methods are based on mechanics and involve the solution of a governing system of equations of motion, either closed-form or numerically. Models that perform a time-wise numerical solution of the equations of motion and advance the location of the simulated landslide incrementally are known as dynamic models. The use of analytical methods is somewhat motivated by the limitations of purely empirical methods, as the unique geometry and materials involved in each case can be accounted for explicitly and a statistically-significant database of previous events is not necessarily required. The simplest analytical model is the classical sliding block model (Figure 1.6), which is based on work-energy theory (e.g., Muller-Bernet in Heim 1932; Sassa 1988). Internal deformation and its associated energy dissipation are neglected and the landslide is treated as a lumped mass. At any point along the path, the sum of the potential energy, kinetic energy and net energy loss equals the initial potential energy. This energy balance can be visualized using the concept of energy grade lines, as shown in Figure 1.6. 13 Figure 1.6. The classical sliding block model, based on work-energy theory. The concept of energy grade lines is useful for visualizing the energy balance, v is the velocity ofthe block, g is the vertical acceleration due to gravity and v2/lg is known as the velocity head, which is the kinetic energy of the block normalized by the product of its mass and g . The same normalization of net energy loss is known as head loss. Note that the positions of the energy lines are referenced to the centre of mass of the block and that the true energy line and mean energy line do not necessarily coincide. Given the initial position of the centre of mass and a suitable relationship to approximate the energy losses, the position and velocity of the block can be determined at any given time. Depending on the complexity of the energy loss relationship, which may include basal resistance as well as additional influences, including momentum transfer to entrained material and momentum corrections due to abrupt slope changes (e.g., Perla et al. 1980), the model can be solved in closed-form or using a numerical method. The sliding block model has been used with various energy loss relationships by a number of workers to model both landslides and snow avalanches (e.g., Voellmy 1955; Salm 1966; Korner 1975; Pariseau 1980; Perla et al. 1980; McClellan and Kaiser 1984; McEwen and Malin 1989; Salm et al. 1990). Strictly speaking, the sliding block model should only apply to the motion of the centre of mass of a rigid body, but in practice it is often applied to the motion of the leading edge of a deforming body, for example, to estimate runup height against an adverse slope (e.g., Mears 1981), or to back-calculate velocities based on observed runup heights (e.g., Costa 1991). More complex multi-block dynamic models have been developed to account explicitly for deformation during motion, including across 3D terrain (e.g., Kobayashi and Kagawa 1987). 14 Distinct element methods (Cundall and Strack 1979) have also been used to model mass movements in both 2D and 3D. Similar to multiple sliding block models, the distinct element method is a discontinuum method in which the main mass is modelled as an assemblage of discrete masses (distinct elements) that interact to simulate large-scale deformations. The distinct elements can be a variety of shapes and sizes and a variety of inter-particle and particle-bed contact relationships can be implemented, although this can require the input of certain parameters with limited physical basis. Cleary and Prakash (2004) demonstrated the potential of the method in landslide modelling, and some work has been done using the commercially-available Particle Flow Code (PFC) developed by HCItasca (e.g., Calvetti et al. 2000; Gonzalez et al. 2003; Pirulli et al. 2003). The method has been limited to the simulation of dry flows, but a recently-implemented fluid analysis feature in PFC may expand its applicability (cf, www.itascagc.com). Nevertheless, relatively long computing times remain a weakness (Gonzalez et al. 2003; Pirulli et al. 2003). Much more attention has been given to the development of continuum methods, which forego the simulation of computationally-expensive small-scale interactions by concentrating on large-scale field effects. In their common depth-averaged form (cf, Chapter 4), continuum methods combine the simplicity and familiarity of sliding block models with a method to account for internal deformation that is based on established hydrodynamic theory. The same basal resistance and momentum transfer relationships that have been developed for sliding block models can be implemented in the depth-averaged continuum framework and, as in shallow water flow models, spreading is driven simply by internal stress gradients. At the same time, the relationships governing the internal stresses can be modified to account for the strength and saturation of the material using established earth pressure theory. But the method is not too complex; in its simplest form, a landslide continuum dynamic model can be governed by just two physically-based parameters: the basal and internal friction angles. A further advantage is that established continuum solution methods can be employed. Closed-form and simple numerical continuum models have been developed for estimating debris flow runout onto fans (e.g., Takahashi and Yoshida 1979) and debris flow and snow avalanche superelevation around bends and runup against adverse slopes (e.g., Hungr et al. 15 1984; Hungr and McClung 1987; McClung and Mears 1995; McClung 2001). Advances in both theory and computational capacity have led to the development of a large number of continuum dynamic models, some of which, including D A N , are commercially-available. But despite this effort, the most important characteristics of extremely rapid landslide motion have not yet been captured by a single model (cf, Chapters 2 and 3). The model proposed in this thesis is based on D A N because, despite its inherent limitations as a 2D model, it includes all of the basic features required to simulate these important characteristics, including methods to account for internal strength, material entrainment and rheology variations. It is also possible that the experience already gained using D A N , including the existing database of calibrated cases, can be transferred to the new model. 1.5 Discussion Continuum dynamic modelling appears to have potential as a quantitative, accurate, objective and accessible method for analyzing extremely rapid landslides. The relatively large amount of attention focused on it suggests that it is widely viewed as the ultimate method, but this view ignores the relative strengths of the other approaches. Empirical/statistical methods may currently lack accuracy but are relatively easy to apply and, most importantly, are shaped by actual events. Simple analytical methods like the sliding block model may not be able to account for deformation but are also easy to apply and can provide spatial information, as well as a simple basis for a better understanding of the complex dynamics that are involved. Discontinuum methods may be relatively new and untested, but are already capable of fully 3D analysis, unlike depth-averaged continuum models, and are potentially better suited to the simulation of coarse-grained rock slides and rock avalanches, which may not strictly satisfy the continuum assumption. As a result, more than one method is likely to be useful in any given case and a multi-tool approach is recommended, which should always be supplemented by expert judgment. This thesis is not a quest for the ultimate method, just an effort to improve the state-of-the-art. 16 To do so, a semi-empirical approach has been adopted, in order to exploit the advantages of both empirical and analytical methods. This approach accounts for the unique geometrical and material characteristics of each case, but is shaped by the behaviour of previous events through model calibration. The micro-mechanics, which are considered intractable, are not directly considered and attention is instead focused on the large-scale characteristics of an event, which can be modelled using relatively simple rheological relationships. Still, the debate continues over whether to use the historical record to build dynamic models or simply to test them. From a purely scientific perspective, it may be a matter of what is physically correct (e.g., Iverson 2003), but from an engineering perspective, it is arguably more a matter of what works, and the utility of the semi-empirical approach has already been demonstrated using D A N and other similar models. The possibility that this approach could contribute to more informed decisions that can potentially reduce future life and property losses warrants further investigation. Besides, the two perspectives are not necessarily mutually exclusive, and work on both fronts could turn out to be mutually beneficial. 1.6 Conclusion With continuing development of marginal areas and consequent increasing exposure of people and property to landslides, there is an increasing need for better methods of performing landslide runout analysis within the landslide risk assessment framework. Continuum dynamic modelling is one alternative that appears promising but still requires improvement. To address these needs, this thesis is focused on the development of a new continuum dynamic model that includes unique features to account for the most important characteristics of extremely rapid landslides. The objectives guiding these features are the subject of the next chapter. 17 2 OBJECTIVES FOR MODEL DEVELOPMENT 2.1 Introduction Continuum dynamic models are based largely on the principles of hydrodynamics. However, while rock avalanches, debris flows, debris avalanches and flow slides all exhibit flow-like movement similar to that of fluids, landslide motion is fundamentally different from open channel water flow. Landslides travel across steep and irregular terrain, they contain earth materials that can resist internal strain, they can entrain substantial volumes of path material, their rheology can vary internally and change in the course of movement and their motion is highly unsteady and nonuniform. Modifications to the basic shallow water equations and powerful computational methods are required in order to capture these unique characteristics. Specifically, a comprehensive model should have the following capabilities: 1) It should be able to simulate motion across complex 3D terrain without the need to input a pre-defined path direction or width. At the same time, it should be able to handle large displacements and deformations, including bifurcation. 2) It should be able to simulate strain-dependent, non-hydrostatic, anisotropic internal stresses due to 3D deformation of material with internal shear strength. 3) It should be able to simulate mass and momentum transfer due to entrainment of path material. 4) It should be able to simulate variations in rheology along the path and within the landslide. 5) It should be efficient and easy to use. 18 The first four criteria are strictly technical, as they are essential for the accurate simulation of real landslides. The last criterion is mainly practical. While it is possible to incorporate the four technical criteria into a complex model that runs on a supercomputer or for several days on a microcomputer, such a model would not be very useful in practice. The last criterion necessitates innovative approaches to satisfying the first four. DAN3D has been designed specifically to satisfy all five criteria, which is unique among existing continuum dynamic models. Other models are reviewed in Chapter 3 in the context of these objectives. The purpose of this chapter is to describe these objectives in detail. 2.2 Objective #1: the ability to simulate motion across complex 3D terrain Topography has a significant influence on landslide dynamics. It controls the orientation of the gravity driving stress, which dominates over the internal stresses on steep slopes. Sudden redirections due to abrupt changes in slope angle or orientation, even by minor surface irregularities, result in momentum and energy losses. Gradual slope changes also influence the dynamics, as centripetal acceleration due to local curvature of the path alters the total bed-normal stress, which in turn can affect both the internal stresses and the basal shear resistance. Topography may or may not provide confinement. On unconfined slopes the path direction and width can be difficult to predict. A significant amount of spreading is possible and, in response to local topography, bifurcation (branching of the flow) can occur (Figure 2.1). Mergers are also possible. Lateral spreading also dissipates energy, which can reduce the travel distance. In contrast, confinement typically decreases the slope angles at which landslides begin to deposit and increases their travel distance (Hungr et al. 1984). However, sudden constrictions along the path can also cause momentum and energy losses. In confined reaches, lateral and transverse internal stresses can influence the flow depth, width and discharge. The shape of the channel bed cross-section plays an important role in the development of these stresses and, as in 19 fluvial hydraulics (e.g., Chow 1959), also controls the size of the basal surface area and therefore the magnitude of the basal resistance. Figure 2.1. Branching and merging of debris avalanches/flows triggered by heavy rainfall in May 1998, Campania Region, Italy. (Photograph taken one day after the events, courtesy of Prof. Francesco M . Guadagno, University of Sannio, Benevento, Italy) Momentum and energy losses can also occur in channel bends. At the same time, lateral and transverse internal stress changes caused by deformation can be significant and can contribute to superelevation or runup (Hungr 1995; McClung and Mears 1995; McClung 2001). Avulsion can occur when runup exceeds the height of the channel on the outside of a bend and bifurcation results when some of the flow remains in the channel. Although the path direction of a confined landslide is generally well-constrained, the locations of potential channel avulsions can be difficult to predict. 20 Whether confined or unconfined, landslides are sensitive to topography. In either case, they can travel extremely long distances and undergo large deformations. A model should be able to handle this behaviour without introducing numerical problems, for example, due to mesh distortion. In addition, it should be able to predict not only the distal limit of motion, but also the path direction and width and the spatial distribution of intensity parameters (cf, Chapter 1) within this area. Finally, simulation of an event should be based on a reasonably accurate representation of the actual topography. 2.3 Objective #2: the ability to simulate the influence of internal strength Flow-like landslides (cf, Hungr et al. 2001) exhibit motion similar to fluids. In contrast to fluids, however, landslides are composed predominantly of earth materials that have internal shear strength. As a result, they can sustain internal shear stresses and resist internal strain. Complex 3D deformations may give rise to non-hydrostatic, anisotropic internal stresses that can have a significant influence on landslide dynamics. Sassa (1988) and Hutter and Savage (1988) were the first to incorporate this idea into their models, a major departure from the basic shallow water assumptions of hydrostatic, isotropic internal stresses. The concept is based on established Rankine earth pressure theory (e.g., Terzaghi and Peck 1967). 2.3.1 Rankine earth pressure theory Rankine (1857) studied the limiting stress conditions for a vertical column of dry, cohesionless soil undergoing 1D horizontal strain (plane strain), an idealization of the loading conditions on soil adjacent to a retaining wall (Figures 2.2 and 2.3). Note that the distinction between dry and wet soil is important, as pore fluid pressures can significantly influence the stress conditions. The effective normal stress, cr' (the stress that is actually transmitted to the soil fabric), is the difference between the total normal stress, o~, and the pore fluid pressure, u . For cohesionless material, the shear strength, Tf, is a function of cr' and the internal friction angle of the material, ty.: 21 [2-1] r y = ( < 7 - \u00C2\u00AB ) t a n $ =cr ' tan$ The two idealized limiting stress conditions are illustrated by the Mohr circles in Figure 2.4. Since it is assumed that no shearing occurs in the horizontal and vertical planes, the effective horizontal and vertical normal stresses, . . A 4 5 \u00C2\u00B0 H / 2 4 ^ _ Figure 2.2. Idealized loading conditions on soil adjacent to an embedded retaining wall (wall friction is neglected and plane horizontal strain is assumed). The zones of horizontal expansion and contraction on either side of the wall are known as the active and passive zones, respectively. In each case, the dashed lines indicate the inclination of the local failure planes from horizontal. (a) \"active\" failure (b) \"passive\" failure Figure 2.3. Schematic illustrations of idealized (a) active and (b) passive local failure. The arrows indicate the directions of the prevailing horizontal normal strain (not to be confused with the horizontal normal stress, which is always compressional). 22 normal stress o \"passive\" circle (horizontal contraction) Figure 2.4. Mohr circles showing the two limiting stress conditions for the idealized loading shown in Figure 2.3. The limiting value of ah' is lower than o~v' during horizontal expansion and higher than av' during horizontal contraction. The dashed lines indicate the inclination of the local failure planes from horizontal. The prevailing stress conditions are strain-dependent. During horizontal expansion, failure occurs on a set of planes inclined at 45\u00C2\u00B0 + $/2 from the horizontal. This is known as the \"active\" Rankine state and is illustrated by the small Mohr circle in Figure 2.4. In this state, the shear stresses that are mobilized on the failure planes act against gravity to resist expansion. The horizontal stress is lower than the vertical stress, and so the ratio o~h'/cr,', which in this limiting condition is known as the active earth pressure coefficient, ka, is less than one. From the geometry of the small Mohr circle in Figure 2.4: [2-2] \u00C2\u00A3 f l = t a n 2 ( 4 5 \u00C2\u00B0 - # / 2 ) 23 The response is different during horizontal contraction. Failure occurs on a set of planes inclined at a lower angle of 45\u00C2\u00B0 - ) r 1- tancf y tand). y 31 The only difference between the two expressions is the first term in the numerator on the right side, but the results are significantly different (Figure 2.10). Assuming drained response, Equation [2-4] suggests that a certain amount of entrainment is possible for any value of y less than y / ( l - tan al tan > az a r . . a r da. \u00E2\u0080\u0094\u00E2\u0080\u0094+- -dx dy dz + PS; Equations [4-6] to [4-9] are the most general form of the incompressible continuum mass and momentum balance equations. 4.4 Boundary conditions Extremely rapid landslides are generally well-bounded by the landslide-atmosphere interface at the top and the landslide-bed interface at the base. The stress state and kinematic conditions at these two interfaces comprise the boundary conditions. It is useful to orient the reference coordinate system so that the depth of the flow, h, and the depth of the bed, b, are measured in the z direction, with z = 0 corresponding to the bottom of the bed (Figure 4.2). For the moment, this direction can still be somewhat arbitrary and does not necessarily have 64 to correspond with either the vertical or bed-normal directions (the directions in which subsequent depth-averaging is typically performed). Figure 4.2. Generic orientation of the reference coordinate system. A more specific orientation is employed in DAN3D (cf, Figure 4.3). 4.4.1 Stress conditions The top surface is assumed to be stress free (relative to the atmospheric pressure). The bed-normal stress at the base balances the bed-normal component of the total weight of material above, including the influence of centripetal acceleration due to bed curvature. The basal shear stress can depend on many factors, including the composition of the material near the base and its strength and drainage characteristics. A description of alternative constitutive relationships for the basal shear strength, consistent with the equivalent fluid approach, is deferred to Section 4.7 in order to keep this derivation as general as possible. 4.4.2 Kinematic conditions It is assumed that material does not enter or leave the landslide at the free surface. This approach neglects entrainment of material from bank failures that may occur as the landslide passes and also neglects possible ejection of material from the free surface. The kinematic boundary condition at the free surface is then: 65 d(b + h) d(b + h) d(b + h) [4-10] - i '- + v. - i >- + v - i J\u00E2\u0080\u0094v = 0 3r A,--='+\"\u00C2\u00BB dx y{z-M,) By It is assumed that material may only enter the landslide at the bed, due to scour of material in the path. Although plowing is technically a free surface phenomenon, which occurs at the flow front, it is treated here mathematically as a component of basal erosion. Consistent with the previous assumption of constant density, it is assumed that the bulk density of the entrained material is the same as that of the landslide. This assumption is valid, for example, for debris avalanches that initiate in the same surficial deposits that are encountered along the path. It is also a reasonable assumption for rock avalanches that entrain saturated material. It must be acknowledged, however, that some common path materials, including ice and snow, can have bulk densities significantly different from those of an overriding landslide, which limits the validity of the constant density assumption in certain cases. The time-dependent rate at which bed material enters the landslide, E[, is also known as the \"erosion velocity\" (Takahashi 1991). The erosion velocity is considered to be positive during erosion, consistent with Takahashi's (1991) definition. The kinematic boundary condition at the bed is: VA , n d b db db r [4-11] \u00E2\u0080\u0094 + v r \u00E2\u0080\u0094 + v\u00E2\u0080\u009E v =-E, dt \u00E2\u0080\u00A2x'\"\"> dx *<--\" dy ''--='> 4.5 Depth-averaging It would be computationally-expensive to perform a numerical solution of the full system of governing equations presented in Section 4.3. Integrating these equations between the bed and the free surface, imposing the boundary conditions, substituting appropriate depth-averaged values and making appropriate assumptions based on scaling arguments produces a much more useful system of equations. The classical St. Venant shallow water equations are based on the same depth-averaging procedure. 66 Two important manipulations are used frequently during the depth-averaging process. The first is Leibniz's rule for interchanging the order of differentiation and integration (e.g., z=b+h dx z=b+h dz \u00E2\u0080\u00A2\u00E2\u0080\u00A2 dx d(b + h) + \u00E2\u0080\u0094), The second defines a depth-averaged value of a function as the depth-wise integral of the function divided by the depth (e.g \u00E2\u0080\u0094 I f vv = \u00E2\u0080\u0094 vxdz , where the overbar denotes a depth-averaged value). z=b+h z=b 4.5.1 Depth-averaging the mass balance equation Integrating Equation [4-6] in the z direction from z = b to z = b + h: =b+h [4-12] dvx dvv dv. \u00E2\u0080\u0094 ^ + \u00E2\u0080\u0094\u00E2\u0080\u00A2- + \u00E2\u0080\u0094 ^ dx dy dz dz = 0 d(v~h) d{Vyh) - + -dx dy d(b + h) d(b + h) v \u00E2\u0080\u0094i '- + V \u00E2\u0080\u0094 - - V db db V,. h V\u00E2\u0080\u009E V. Substituting the kinematic boundary conditions from Equations [4-10] and [4-11] into Equation [4-12], the depth-averaged mass balance equation is: dh d(v 0 (a negative value of signifies deposition). This constraint must be imposed explicitly on the momentum balance equations because the process of deposition is external to the chosen reference frame. Material that leaves simply removes its own share of the momentum, without influencing the momentum of the material that remains within the reference frame (Hungr 1995). A negative erosion velocity in the momentum balance equations would imply, erroneously, that the process is internal to the flow and would produce a thrust in the direction of motion similar to that of a rocket (cf, Cannon and Savage 1988; Hungr 1990; Cannon and Savage 1990; Erlichson 1991). 68 Because the depth-wise velocity profile may be nonuniform, resulting in differential advection of momentum, momentum corrections must be applied to relate v v 2 to vv and yvv to v t v v , which allows the derivation to proceed with a further separation of terms (e.g., vx2 = gvx , where g is a momentum correction coefficient). When the velocity in the y direction is zero, for example, a value of g = 1.2 applies i f the depth-wise velocity profile is parabolic (i.e., no basal sliding) and a value of g = 1 applies i f the depth-wise velocity profile is uniform (i.e., full basal sliding) (Savage and Hutter 1989). Since rapid landslides typically exhibit some degree of basal sliding, the required momentum correction is typically small and it is therefore commonly assumed that v t = vv and vvv =y vv . Substituting these relationships into Equation [4-15], expanding the partial derivatives and collecting terms: z=b+h [4-16] 3vx | 8(v,v 2) | d(vxvy) | 8 ( V J V I ) dt dx dy dz V J dz = P v V V dh d(v,h) a ( v / z ) \u00E2\u0080\u0094 + \u00E2\u0080\u0094 - \u00E2\u0080\u0094 - + \u00E2\u0080\u0094 - \u00E2\u0080\u0094 -dt dx dy + h rdvx -dvx. -dv ^ \u00E2\u0080\u0094 - + v \u00E2\u0080\u0094 dt V dx y dy Substituting the depth-averaged mass balance Equation [4-13] into Equation [4-16]: [4-17] z=b+h z=b P V v dvx | a ( y , 2 ) | a ( v , y , ) | d{vxvz) \u00E2\u0080\u00A2+ dt dx dy dz dz = P ( fdvx -dvx -dvx^ \u00E2\u0080\u0094 - + v v \u00E2\u0080\u0094 - + v v \u00E2\u0080\u0094 -dt dx ' dy Now, integrating the right side of Equation [4-7] in the z direction from z = b to z = b + h: 69 z=b+h z=b d(Tx dtyx 3 T A- + \u00E2\u0080\u0094 \u00E2\u0080\u0094 + \u00E2\u0080\u00A2 - v dx dy dz + Pgx dz [4-18] d((Jxh) | d(Tvxh) dx dy d(b + h) d(b + h) 36 ^ 36 CT hT T x<^ > 3x ,VV(^> 3y Z\",I=4* + Phgx Substituting the stress free boundary condition at the top into Equation [4-18]: =b+li z=b da' drvx 3 r , dx dy dz + Pgx dz [4-19] d{oxh) d(r\u00E2\u0080\u009Eh) r x '\u00E2\u0080\u00A2+\u00E2\u0080\u0094\u00E2\u0080\u0094\u00E2\u0080\u0094- + dx dy db db CT + T T '<--*' dx dy + P\"gx As mentioned, the same series of manipulations can be performed in the y and z directions. The resulting depth-averaged momentum balance equations in the x, y and z directions are, respectively: f [4-20] h V V 8vv \u00E2\u0080\u0094 dvY \u00E2\u0080\u0094 dvr \u00E2\u0080\u0094 - - l - v , . \u00E2\u0080\u0094 - + v , \u00E2\u0080\u0094 dt v dx y dy rd[Vxh) d^h) J dx - + -db \ dy + dx + T db r'\"\" dy -T_. V + phgx dv.. \u00E2\u0080\u00A2 +v.. [4-21] dt A dx f .dv \u00E2\u0080\u0094dv \u00E2\u0080\u0094 ^^y \u00E2\u0080\u0094 dy d[Txyh) | d{Vyh) dx dy db r + db T \u00E2\u0080\u0094 + a T dx dy '\u00E2\u0080\u00A2 \u00E2\u0080\u00A2y (.-=*) + phgv 70 dv. \u00E2\u0080\u0094 dv. \u00E2\u0080\u0094 dv \u00E2\u0080\u0094- + vY\u00E2\u0080\u0094^ + v \u00E2\u0080\u0094-dt x dx y dy [4-22] (d(z~zh) dfch) ( d b dx dy db T,_ + Tv_ CJ v dx dy *<--\" Equations [4-13], [4-20], [4-21] and [4-22] are the most general form of the Eulerian, depth-averaged governing equations. 4.6 Further simplifications 4.6.1 Lagrangian reference frame A number of additional manipulations can be performed to simplify this system. First, the equations can be recast in their Lagrangian form by substituting Lagrangian derivatives for , \u00E2\u0080\u00A2 , \u00E2\u0080\u00A2 / Dh their Eulerian counterparts (e.g., = dh \u00E2\u0080\u0094 dh \u00E2\u0080\u0094 dh \u00E2\u0080\u0094 + v t \u00E2\u0080\u0094 + vv \u00E2\u0080\u0094 dt dx ' dy where the left side is the Lagrangian form and the right side is the Eulerian form). Lagrangian derivatives are taken with respect to a reference frame that moves with the material and are therefore also known as material, substantive, advective or convective derivatives. Expanding the partial derivatives in the mass balance Equation [4-13], collecting terms and substituting the Lagrangian derivative of h for its corresponding Eulerian derivative gives the following Lagrangian form of the depth-averaged mass balance equation: r . \u00E2\u0080\u009E\u00E2\u0080\u009E-, Dh , 4-23 \u00E2\u0080\u0094 + h Dt dx dy V = E. Second, several terms cancel out of the momentum balance equations with an appropriate reorientation of the reference coordinate system. In DAN3D, the z direction is aligned with the local bed-normal direction (Figure 4.3) and the x direction is aligned with the local 71 direction of motion, which eliminates spatial derivatives of b as well as y and z direction components of velocity. Figure 4.3. Orientation of the Lagrangian reference coordinate system used in DAN3D. The z direction is aligned with the local bed-normal direction (the x direction, not shown, is aligned with the local direction of motion). Substituting the Lagrangian derivatives of v t , vv and v. for their corresponding Eulerian counterparts and reorienting the local reference coordinate system, the Lagrangian forms of the depth-averaged momentum balance equations in the x , y and z directions are, respectively: [4-24] d(axh) + phgx dx dy 72 A Lagrangian approach has been adopted for the following three reasons: 1) as shown above, the convective acceleration terms are eliminated from the momentum balance equations, which facilitates subsequent numerical integration and improves efficiency; 2) higher resolution is possible because available computational resources can be concentrated within the simulated slide mass, where they are needed (cf, Chapter 3); and 3) DAN3D becomes a logical extension of D A N , which was also formulated in the Lagrangian framework. 4.6.2 Shallow flow assumption Further simplifications can be made using the following scaling argument, which is well-established (cf, Savage and Hutter 1989; Gray et al. 1999). Assuming that the depth varies gradually and is small relative to the length and width of the landslide, the terms containing shear stress derivatives of rv_ and tv. in Equation [4-26] are correspondingly small relative to the total bed-normal stress at the base, cr , and can therefore be neglected. This is the classical shallow flow assumption (e.g., Chow 1959). When the flow is in contact with the bed, the bed-normal curvature of the flow is equal to the bed-normal curvature of the path. It follows that, unless the material is airborne, the Lagrangian derivative of v. is equal to the centripetal acceleration due to bed-normal curvature of the path in the direction of motion: [ 4-27] ^ = ^ Dt R where R is the bed-normal radius of curvature of the path in the direction of motion, which is considered to be positive for concave curvature. The z direction component of gravity is: [4-28] g. =-gcosa 73 where g is the acceleration due to gravity and a is the inclination of the bed from horizontal. Substituting Equations [4-27] and [4-28] into Equation [4-26], neglecting the relatively small shear stress derivatives and rearranging terms produces the following expression for the total bed-normal stress at the base: \u00E2\u0080\u00942 \ g cos a + R K J \u00E2\u0080\u00942 V Note that the value of a. becomes negative when - i - > gcosa, which implies that the material becomes airborne. Although ballistic/freefall conditions commonly occur in extremely rapid landslides, they are not accounted for explicitly in DAN3D (or in any other existing continuum dynamic model). As a result, net energy losses associated with impact after freefall must be implicitly accounted for in the basal shear strength term. A method to account for such effects explicitly could possibly be implemented in a future version of the model. 4.6.3 Stress state normalization Further simplifications are possible with the use of classical soil mechanics techniques. Assuming that all stresses increase linearly with depth below the free surface, and consistent with the Rankine earth pressure theory described in Chapter 2, it is useful to normalize the stress state with respect to the total bed-normal stress using stress coefficients, denoted by k (e.g., \u00E2\u0080\u0094^r^-+phgy Dt dx dy z X Figure 4.4. Normalized total stress state on an element of material within a landslide. The normalization is with respect to the total bed-normal stress (e.g., kx = crx/cr_). The stress coefficients are considered positive as shown. Compare with Figure 2.6 in Chapter 2. dk Finally, it is assumed that spatial variations in the normalized stress state (e.g., \u00E2\u0080\u0094-) are dx relatively small and can therefore be neglected (Gray et al. 1999). Expanding the partial derivatives in Equations [4-30] and [4-31], neglecting terms containing derivatives of stress coefficients and rearranging terms, the Lagrangian forms of the depth-averaged momentum balance equations in the x and y directions are, respectively: [ 4 .32] ph^L = phgx-kxcr. ^--ka ^- + T.x -p(Vx-vx )Et 75 [4-33] ph = phg\u00E2\u0080\u009E-ka_ ka. dh dx The terms on the left side of Equations [4-32] and [4-33] represent local depth-averaged accelerations of a moving reference column of material (multiplied by the mass of the column per unit basal area). The first three terms on the right side represent depth-averaged gravity, normal and transverse shear stresses, respectively. Additional terms appear in the x direction momentum balance equation due to the assumed orientation of the reference coordinate system. The fourth term on the right side of Equation [4-32] represents the basal shear stress, which is described in more detail in the next section. The fifth term represents momentum flux due to entrainment of path material. Note again that Et > 0 . Equations [4-23], [4-29], [4-32] and [4-33] are the fundamental system of depth-averaged, Lagrangian mass and momentum balance equations used in DAN3D. Additional simplifications due to assumptions specific to the numerical solution method are described in Chapter 5. 4.7 Basal shear resistance The basal shear stress, r . / ( , opposes motion and, due to the chosen reference coordinate system orientation, is always negative. Consistent with the concept of equivalent fluid, r is governed by a basal rheology that may be different from the internal rheology. To allow the simulation of different types of rapid landslides involving different geological materials, a variety of basal rheological relationships can be implemented in DAN3D, including laminar, turbulent, plastic, Bingham, frictional and Voellmy rheologies. The user can change the basal rheology along the path or within the slide mass, in accordance with Objective #4 in Chapter 2. 76 4.7.1 Laminar Laminar flow is characteristically smooth and streamlined and is exhibited by Newtonian fluids (in which shear stresses are linearly proportional to shear rate) at relatively low Reynolds numbers (the ratio of inertial to viscous stresses). Laminar basal shear resistance is proportional to the depth-averaged flow velocity and inversely proportional to the flow depth. In the DAN3D coordinate system: [4-34] r = \u00E2\u0080\u0094 where ju is the dynamic viscosity. Note that, in the equivalent fluid framework, the landslide material does not have to be strictly Newtonian for the laminar resistance v relationship to be applied. 4.7.2 Turbulent At relatively higher ratios of inertial to viscous stresses, flow may transition to a turbulent regime, which is characterized by intense mixing. Turbulent basal shear resistance is proportional to the square of the depth-averaged flow velocity and can be calculated using the Manning equation: __pgn v v [4-35] : where n is the Manning roughness coefficient. A common alternative to Equation [4-35] is the Chezy equation: [4-36] _ pgvx C2 77 where C is the Chezy coefficient, which is related to Manning's n by C = /z 6 /\u00C2\u00AB . Again, turbulent resistance does not necessarily apply only to Newtonian fluids within the equivalent fluid framework. 4.7.3 Plastic Plastic behaviour occurs when the basal shear strength is constant: [4-37] T = -c where c represents a constant shear strength, such as the undrained shear strength of clay. 4.7.4 Bingham The Bingham resistance model combines plastic and viscous behaviour. A so-called Bingham fluid behaves as a rigid material below a threshold yield strength, but as a viscous material above. The following cubic equation must be solved to determine the basal shear resistance: [4-38] T.v 3 + 3 ^yield _^ MBingham ^x ^ T 3 where Tyield is the Bingham yield stress and p B i n g h a m is the Bingham viscosity. 4.7.5 Frictional Frictional behaviour was introduced in Chapter 2 and is often exhibited by granular materials. Frictional basal resistance is proportional to the effective bed-normal stress at the 78 base, cr.' , which is the difference between the total bed-normal stress at the base, c * (--=\u00C2\u00BB) '(.-\u00C2\u00BB/.) and the pore fluid pressure at the base, u : [4-39] r =-(cr_-w) tan0 = -cr. ' tan0 where 0 is the dynamic basal friction angle. Pore fluid pressure can be related to the total bed-normal stress using the pore pressure ratio, ru = u/a., in which case: [4-40] T =-cr_ ( 1 - r ) tan0 Equation [4-40] can be simplified to include only one independent variable, a bulk basal friction angle, , dW ds yu 5.3.2 Interpolating kernel A number of interpolating kernels that satisfy the criteria in Equations [5-2] and [5-3] can be used. However, as recommended by Monaghan (1992), a Gaussian kernel has been used almost exclusively in DAN3D. In 2D form (appropriate for a depth-averaged model for 3D terrain): [5-15] W\u00E2\u0080\u009E = - ^ e u J V (Oaussian) JI ^ [5-16] dW ds 2SiJ I f = \u00E2\u0080\u0094J-e v where the smoothing length at any given time is assumed to be the same everywhere. The Gaussian kernel and kernel gradient are asymptotic to zero at an infinite distance, but in 91 practice they are negligible at distances greater than about It, so only the n particles within this radius need to be included in the summations. Because the kernel never actually reaches zero, a cutoff depth representing the margin of the flow must be specified by the user. Other kernels have explicit cutoffs. Kernels based on linear, quadratic and cubic functions, in S-which W and VW are zero i f -j- > 2, have also been implemented as options in DAN3D. In 2D form they are, respectively: [ 5- 1 7l W ^ , r t K t \\u00C2\u00B0 A i f - ^ < 2 I [5-18] dW ds 'J 'linear) v 8 v i f 2 [5-19] W\u00E2\u0080\u009E y (qttadnn nt \2) f t 4 v V J 2 \ i f - ^< 2 [5-20] dw ds y (ifiimlrtilic) K2J 1 \u00E2\u0080\u0094 v JJ i f - ^ < 2 [5-21] W\u00E2\u0080\u009E =< K t I14 A 2 - ^ J nt v ' A 3 + \u00E2\u0080\u0094 4 i f l< -2-<2 I if ^ < 1 t [5-22] 9 ^ ds K t \"J (cubic) 30 15 V v 1 4 A 2-1 K t <' S ^ JL t 2 \ i f 1 < X L < 2 e i f ^ < i \u00C2\u00A3 92 The four alternative interpolating kernels are compared in Figure 5.2. The linear and quadratic kernels are not strictly smooth, as the resulting interpolated free surface contains peaks at each particle location. With the cone-shaped linear kernel, neighbouring particles \"push\" each other apart with the same force at any distance within two smoothing lengths of one another, which causes instabilities. The Gaussian and cubic kernels are very similar. In fact, with appropriate different smoothing lengths (i.e., when ^ s -y/lO/7 I ^ ) they are virtually identical. As a result, like other SPH-based models (e.g., Monaghan 1992), DAN3D is not very sensitive to the choice between them. Marginal efficiency can be gained with the cubic kernel, but because it is a well-known function with a simple physical interpretation (Monaghan 1992), the Gaussian kernel has been used for all of the analyses presented in this thesis. 0 0.5 1 1.5 2 2.5 3 Normalized Distance s/C Figure 5.2. Comparison of different interpolating kernels that can be implemented in DAN3D. Physically, an interpolating kernel represents a cross-section of a single particle with unit volume. The Gaussian and cubic kernels form similar bell-shaped particles that produce smooth free surfaces. 93 To maintain resolution during large deformations (i.e., to have a relatively constant number of neighbouring particles within the radius of influence), DAN3D uses an adaptive smoothing length, analogous to an adaptive finite element mesh: where B is a dimensionless particle smoothing coefficient, which influences the smoothness which means that the number of particles within a radius of one smoothing length of each particle is, on average, equal to KB1 . For example, a value of B = 4 gives each particle, on average, about 50 neighbours within a radius of one smoothing length. The sensitivity of the model to the value of B is investigated in Chapter 6. This method of adaptive smoothing was recommended by Monaghan (1985) and has been used for all of the analyses presented in this thesis. But because the smoothing length is the same for all of the particles (i.e., spatially-constant), some resolution is sacrificed in regions with local particle number densities that deviate strongly from the overall average. This can occur near the flow margins, within channelized reaches of the path or within secondary flows that become separated from the main mass. A spatially-variable smoothing length is an alternative and can be implemented by limiting the summation in Equation [5-23] to a suitable radius of influence. Preliminary testing of this alternative has suggested that, while resolution can be improved, the ability of the model to handle highly channelized cases may be simultaneously sacrificed. This is because relatively higher local smoothing lengths result in channelized sections of the path, where each particle has a reduced number of neighbours in the lateral direction. The increased smoothing results in weaker inter-particle forces, which are then dominated by gravitational forces that eventually cause the particles to line up in the downstream direction. A directional spatially-variable smoothing method is a possible B of the interpolated flow depth. The term is the average areal particle number density, 94 solution to this problem. However, sensitivity analyses have suggested that both of these issues (resolution and particle alignment during channelization) can be resolved simply with the use of more particles (cf, Chapter 6). Because the capabilities of personal computers continue to improve, this solution may be the easiest to implement in future versions of the model. 5.4 Strain interpolation 5.4.1 Incremental tangential strain state As described in Chapter 2, stresses in earth materials are strain-dependent, and strain within a 3D deforming landslide mass may be extremely complex. Similar to the stress tensor in Equation [4-5], the 3D strain tensor is: [5-24] e = 7^ 2 2 r, xy_ 2 2 2 LL 2 where e and y denote normal and engineering shear strains, respectively. In the framework of a depth-averaged model, estimates of strain tangential to the bed are required in order to determine the tangential stresses using the Rankine theory. Plane strain in the tangential x-y plane is assumed. As described in Section 5.5, the updated tangential stresses are coupled with the updated bed-normal and basal shear stresses to obtain the updated 3D stress state. The incremental tangential strain tensor is: [5-25] Ae = Ae. 95 5.4.2 Strain rosette method The particle-tracking nature of the SPH method provides a framework for approximating the incremental tangential strain state at each reference column location. Tangential flow deformation can be approximated from relative particle motion; converging particles indicate flow compression while diverging particles indicate flow extension. D A N 3 D uses an approach analogous to strain measurement in flat plates using strain gauge rosettes (Figure Figure 5.3. Plane strain measurement using strain gauge rosettes, a common technique in structural and mechanical engineering. The complete plane strain state (e.g., in a flat metal plate) can be determined using only three directional strain measurements. The directions can be arbitrary, as long as they are different, but typical rosette configurations are shown: (a) rectangular and (b) delta. Neglecting the influence of reorientation of the local reference frame during the small time interval, the incremental tangential strain between two particles can be approximated by: where the old separation distance is primed. Convergence of particles produces a positive strain. The relative direction to the neighbouring particle can be approximated by: 5.3). (a) (b) [5-26] A\u00C2\u00A3jj = 96 [5-27] ft, = tan\" zJL where ft., is an angle in the local x-y plane, measured counter-clockwise from the positive x -axis. Using conventional strain compatibility theory, the particle-centred incremental tangential strain state in Equation [5-25] can be represented by the function: [5-28] Ae(ft),. f Asr + As,, } f Aer- As,, } _. (Ay > cos 2ft + V sin 2ft where As (ft) is the incremental tangential strain in the direction ft , ft is a parametric angle in the local x-y plane, measured counter-clockwise from the positive x-axis, Aex is the incremental strain in the x -direction, Ae is the incremental strain in the y -direction and Ayxv is the incremental transverse engineering shear strain. Equation [5-28] has the form of an equation of a plane, with cos 2ft and sin 2ft being the independent variables. The procedure for interpolating the tangential strain state is as follows. At each time step and at each particle, values of As^ and 0.. corresponding to each neighbouring particle within a distance of one smoothing length are calculated. Substituting these values into Equation [5-28] produces a system of equations with three unknowns (Aex, Asy and Ayxv). Analogous to the three gauge rosette method shown in Figure 5.3, when exactly three neighbouring particles are present there are three equations and the system is determined. With less than three neighbours the system is indeterminate and the influence of strain must be neglected; when two or three particles are isolated together the deformations are extremely limited, making this assumption reasonable, and when a single particle is isolated the flow depth gradient is zero and the local strain and stress states become irrelevant (cf, governing Equations [4-32] and [4-33]). With more than three neighbours, which is much more typical 97 when a smoothing coefficient close to 5 = 4 is used, the system is redundant and a least squares approximation is used to fit a plane to the data points, providing estimates of A\u00C2\u00A3 t , Aev and A y v v . A sample interpolation from a DAN3D simulation is shown in Figure 5.4. Figure 5.4. Example of a DAN3D incremental tangential strain interpolation, (a) The relative positions of neighbouring particles are updated at each time step, (b) Ae^ (Equation [5-26]) and 0.} (Equation [5-27]) are calculated for each neighbouring particle, and a least squares fit of Equation [5-28] provides an estimate of the column/particle-centred strain state. This strain interpolation method requires no assumptions regarding the orientation of the principal incremental tangential strains. Denlinger and Iverson's (2004) model is the only other existing continuum dynamic model with this capability. The interpolated strain state can then be used to increment the internal stress state accordingly, but this step is not straightforward, as the process of stress redistribution following yielding within a rapidly deforming landslide is extremely complex and poorly understood. Denlinger and Iverson (2004) proposed an iterative method for redistributing the stresses at each time step to satisfy a plastic yield criterion. A similar iterative method was originally implemented in DAN3D, but such an approach uses ad hoc assumptions regarding stress redistribution while sacrificing considerable computational efficiency. 98 5.4.3 Simplified strain rosette method The transparency and efficiency of the strain interpolation and resulting stress redistribution can be improved significantly by explicitly assuming a suitable orientation of the principal incremental tangential strains. To help determine a reasonable assumption, the interpolated strains associated with randomly selected particles were tracked during a series of DAN3D runs in which the stress state was assumed to be hydrostatic. These tests showed that the direction of maximum straining typically corresponds very closely with the direction of motion, implying predominantly streamlined flow. Exceptions occurred, but did not persist for long, when there were very abrupt changes in direction that caused strong oblique shocks. The strain interpolation shown in Figure 5.4, taken from one of these tests, is a good example of this behaviour. The maximum and minimum points on the curve, which represent the interpolated principal incremental strains, correspond with directions at approximately 90\u00C2\u00B0 increments from the direction of motion. As a result, a simplified approach has been adopted. Assuming that Aex (in the direction of motion) and Ae approximate the principal incremental tangential strains (i.e., assuming Ayxv = 0), Equation [5-28] becomes: f Aer + Ae ^ [5-29] Ae{d\ = y_ Ji f Ae-Ae ^ + cos 26 This approach is analogous to strain measurement using a two gauge rosette, in which the gauges are aligned with the assumed directions of the principal strains. Equation [5-29] has the form of an equation of a straight line. Similar to the two gauge rosette method, when exactly two neighbouring particles are present there are two equations and the system is determined. With less than two neighbours the system is indeterminate and the influence of strain must be neglected. With more than two neighbours, the system is redundant and a much more efficient two-parameter least squares approximation is used to fit a line to the data points, providing estimates of Aex and Aey . An example is shown in Figure 5.5. 99 -1 -0.5 0 0.5 1 cos29 Figure 5.5. Example of a DAN3D incremental tangential strain interpolation using the simplified method. Similar to the rigorous method shown in Figure 5.4 (and using the same data), the relative positions of neighbouring particles are updated at each time step and Aejj and Q.. are calculated for each neighbouring particle. However, a least squares fit of linear Equation [5-29], rather than planar Equation [5-28], provides an estimate of the column/particle-centred strain state. Trial comparative runs using typical slide geometries have indicated that this simplified method produces results comparable to the rigorous method. Similar assumptions regarding the orientation of the principal strains have been made by other workers, as described in the next section. Note that the use of this simplified method does not preclude the implementation of the more rigorous method in the future, when justifiable assumptions about stress redistribution can be made and computational efficiency is not such a significant constraint. 5.5 Stress update 5.5.1 Assumed orientation of principal axes Referring to Figure 2.6 in Chapter 2, with the assumption of Ayxy = 0 and again neglecting the influence of reorientation of the local reference frame during the small time interval, 100 shear stresses r and Tyx equal zero and ay (the cross-stream normal stress) is therefore one of the 3D principal stresses. It is useful to compare this assumption with those used by other workers. Gray et al. (1999) used a similar approach, but instead assumed that the principal axes correspond with the mean downslope direction, as measured in a reference system. When the mean downslope direction equals the actual downslope direction and the direction of motion is exactly downslope, these two assumptions are identical. In contrast, Chen and Lee (2000) assumed that the principal axes are aligned with a global coordinate system, which makes the results dependent on the choice of global coordinate orientation. Similarly, Iverson and Denlinger (2001) assumed that the principal axes are oriented at a fixed angle from arbitrarily-oriented local coordinate axes. Again, in DAN3D the principal axes are aligned with the direction of motion. This is by no means an exact assumption, especially in the case of strongly curving flow or flow around obstacles, but one that appears the most intuitive of those above. As mentioned in the previous section, the only other known alternative is to make assumptions about complex 3D stress redistribution, similar to Denlinger and Iverson (2004), at the expense of model transparency and efficiency. With the assumption that the transverse shear stresses equal zero, the governing x and y direction momentum balance Equations [4-32] and [4-33] become, respectively: [5-30] ph^^ = phg-kxa. ^- + T.r -p(v -vr \E, [5-31] ph\u00E2\u0080\u0094y- = phg-kv(j_ \u00E2\u0080\u0094 Dt y & y y dy 5.5.2 Internal yield criterion In accordance with the equivalent fluid concept (cf, Chapter 4) and the assumptions of other workers, a frictional model, distinct from the basal rheology, has been adopted to simulate 101 the internal rheology of the slide mass. The yield criterion is governed by a single parameter, the internal friction angle, . Similar to the bulk basal friction angle, b (Equation [4-41]), the influence of pore fluid pressure can be accounted for implicitly with (j>.. A hypothetical Mohr circle representation of the normalized 3D stress state is shown in Figure 5.6. Figure 5.6. Hypothetical Mohr circle representation of the normalized 3D stress state shown in Figure 4.4 in Chapter 4. Since the transverse shear stress coefficients, kxv and k , are assumed to be zero, kv is a principal stress coefficient. As in the Rankine theory (cf, Chapter 2), the stress coefficients are limited by the yield criterion. The minimum and maximum limiting values correspond with the active and passive stress states, respectively. Similar to the approach used by Gray et al. (1999), and consistent with the assumption that downstream deformations dominate, the x -direction is given priority. This assumption is also made implicitly in D A N (Hungr 1995) and other 2D IT models that neglect the influence of lateral strain. As long as after Savage and Hutter (1989): < tan$ (Iverson 1997), 102 [5-32] k . . . =2 L J x (min/ max) 1\u00C2\u00B1 l - c o s 2 # 1+ ^ cos $ The limiting \u00C2\u00A3 y values are then a function of the prevailing value of kx [5-33] ky[mi min) fk+l) 2 2 l - s in$ . ^ \u00E2\u0080\u0094: + J + I 2 J V I 2 J ^ ' 7 1 + sin d>. [5-34] k y(inax) v \" y v \" y 2 f r v ^ 1 + sin $ N v l - s i n # j On the other hand, when > tan : l + s in 2 $ [5-35] \u00C2\u00A3 V ( m i l l / m a x ) 2 cos $ [5-36] k ^(min) 1 + sin d>. [5-37] k 1 y(inax) 1 - sin Note that a value of $ = 0\u00C2\u00B0 produces the liquefaction condition of kx -ky=\ and can therefore be used to simulate fluid flow. 103 5.5.3 Stress incrementation Assuming elastic-plastic behaviour with a constant value of the stress response to deformation is implemented as follows. Initially, stress coefficients kx and ky are set to 1 (the tangential normal stresses are assumed to be hydrostatic and isotropic). Use of the classical \"at rest\" earth pressure coefficient (Terzaghi 1920; Jaky 1944), which is less than 1, may seem more appropriate from a theoretical standpoint, but does not have a significant influence on the model results, as the stress conditions change rapidly following initial deformation. At each time step, the change in the local stress state is calculated based on: where old values are primed, E is the Young's elastic modulus and v is the Poisson's ratio. The first term on the right side of Equation [5-38] represents the incremental change in tangential stress due to a change in the total bed-normal stress. The second term is an expression of Hooke's Law and represents the incremental change in tangential stress due to tangential strain, as estimated using the method described in Section 5.4. The matrix in this term is the stiffness matrix for an elastic, isotropic material in plane strain. It must be acknowledged that the assumption of linear elastic response is a major simplification. In reality, the stress-strain relationship within a rapidly deforming landslide may be extremely complex (e.g., non-linear; highly dependent on pore pressure response; spatially, temporally and cyclically variable). Furthermore, assuming for simplicity that the x and y direction responses are de-coupled (i.e., that the Poisson's ratio is zero), Equation [5-38] simplifies to: [5-38] 104 Technically, such de-coupling is inconsistent with the assumption of incompressibility, but it does not have a significant influence on the model results. Expanding Equation [5-39] and substituting appropriate stress coefficients: [5-40] krcr_ -k'cr. A A X i kyaz-ky'az k'CT.-k'CT. A A A i. k'CT.-k'CT. y * y ' Dividing Equation [5-40] by the updated total bed-normal stress, a., and rearranging terms produces the following equation for the updated tangential stress coefficients: In order to satisfy the assumption that all of the internal stresses increase linearly with depth from zero at the free surface (cf, Chapter 4), the normalized elastic modulus, E/ar, is assumed to be constant: where D is a dimensionless stiffness coefficient, which is assumed to be spatially and temporally constant. This assumption is consistent with the observation that the stiffness of cohesionless material increases with increasing confining pressure. Stiffness is also assumed to be independent of strain. The current default value in DAN3D is D = 200 for both compression and extension. This is within the range of dimensionless stiffness coefficients estimated by Hungr (1995), which were based on measured values of stiffness of granular soil behind retaining walls during active and passive deformation (Terzaghi and Peck 1967). The sensitivity of the model to the value of D is investigated in Chapter 6. The results show that DAN3D is not very sensitive to the specified value of D, as Hungr (1995) also found with D A N , but it is possible that the incremental approach can aid numerical stability. Stress 105 response occurs over several time steps, smoothing anomalous strain interpolations and preventing numerical oscillation between active and passive stress states (Tai and Gray 1998). As described in Chapter 3, with the exception of Denlinger and Iverson (2004), other existing models assume instantaneous stress response, equivalent to using an infinite D value. The apparent insensitivity of models to the assumed stress-strain relationship suggests that deformation is so large and occurs so fast in extremely rapid landslides that the stress response assumptions may be virtually immaterial. 5.6 Entrainment 5.6.1 Momentum transfer based on erosion velocity and growth rate As described in Chapter 4, a momentum flux term due to entrainment of path material appears in the x direction momentum balance Equation [5-30]. Assuming, as in D A N (Hungr (1995), that the entrained material is initially stationary (i.e., v =0), Equation [5-30] becomes: [5-43] p h ^ L = phg-k.a. ^ T + T.x -pvE, Dt ''--^ dx u'h) The momentum flux term (the last term on the right side of Equation [5-43]) arises because momentum must be transferred, by solid collisions, fluid thrust and friction, from the landslide to the path material in order to accelerate it (in this case from rest) to the reference frame velocity, vv . As mentioned in Chapter 2, this process results in a velocity-dependent inertial resistance, which is additional to the basal shear resistance (Perla et al. 1980). A n alternative to this explicit formulation is to account for the momentum transfer effect implicitly during calibration of a bulk basal rheology. However, separating the terms is useful because it facilitates independent analysis of shear and entrainment related resistances, permits simulation of cases involving isolated entrainment zones that can be independently quantified and ensures compatibility between momentum and mass balances. 106 Entrained material is assumed to have the same constant bulk density as the overriding landslide, which is approximately valid in most practical cases (cf., Chapter 4). Spatial and temporal variations in density can be handled by the SPH-based discretization method, which was originally developed for the simulation of compressible flows, but the necessary modifications have not been incorporated into DAN3D. Some of the most challenging assumptions concern the availability of erodible material and the rate at which it is entrained by the landslide. DAN3D requires that maximum erosion depths be estimated independently, for example, using surficial geological maps, aerial photographs or field surveys, and input at nodal locations on a fixed reference grid. Finite or infinite values can be used to simulate supply-limited or supply-unlimited conditions, respectively (cf, Chapter 2). Entrainment occurs when a particle is centred within a fixed grid cell containing erodible material. As the flow velocity increases, the bed erosion velocity, Et, can have an important influence on the momentum balance Equation [5-43]. As described in Chapters 2 and 3, limited effort has been directed at understanding actual entrainment mechanisms and constraining erosion rates, and so it appears that an empirical approach must be adopted at this time. Takahashi (1991) and Hungr (1995) used different approaches to ensure that the predicted maximum erosion depth at a point is reached only after the entire landslide passes. In D A N (Hungr 1995), the erosion rate increases in proportion to the flow depth, resulting in a depth-proportional distribution of entrained material and natural exponential growth of the landslide with displacement. Although largely empirical, the method has a physical basis, as the changes in stress conditions leading to failure within the path material can be related to changes in the total bed-normal stress and therefore the flow depth. DAN3D uses a similar empirical approach based on a user-prescribed parameter, Es, which represents the bed-normal depth eroded per unit flow depth and unit displacement or, equivalently, the displacement-dependent natural exponential growth rate (with dimensions L~[). Note the difference between the erosion velocity, E,, a time-dependent erosion rate 107 (with subscript t) and the growth rate, Es, a displacement-dependent erosion rate (with subscript s). A constant growth rate is independent of flow velocity. For example, a value of Es=0.0\ m\"1 produces a 1% increase in local flow volume per metre travelled, irrespective of flow velocity. The assumption of natural exponential growth with displacement is used for its simplicity and may eventually serve as a baseline for more complex entrainment modelling, including the development of constitutive relationships for growth rates that are based on actual entrainment mechanisms. Similar to the open basal rheological kernel used in DAN3D, the user could potentially select an appropriate entrainment relationship from a menu of proposed models. Such models could incorporate dependence on other factors, including flow velocity, slope angle, path curvature, surface roughness or the strength and drainage characteristics of local path material. In the meantime, different values of Es can be assigned to distinct entrainment zones to empirically account for such factors. The erosion velocity and the growth rate are related by: [5-44] E, = EshVx The x direction momentum balance Equation [5-43] then becomes: [5-45] p h ^ L = phg-kxo. ^r + t.x -phEv* Dt x U:'h) dx The momentum flux term in Equation [5-45], and therefore the inertial resistance produced by entrainment, is then a function of the flow depth and the square of the flow velocity. This velocity-squared term has the potential to dominate the resistance at high flow velocities. However, the momentum transfer effect may be countered by a coinciding reduction in the effective basal shear strength caused by undrained loading or the incorporation of additional water and/or weak path material (cf, Chapter 2). DAN3D allows the user to modify the basal rheology at the onset of entrainment to simulate this effect. Note that Es is 108 automatically set to zero if the grid cell in which the particle is located does not contain any entrainable material (including the condition when available material has been depleted by preceding flow). 5.6.2 Mass transfer SPH provides a simple framework for handling mass transfer between the landslide and an erodible bed. The basal area independently covered by particle / , Aj, can be approximated by (Wang and Shen 1999): [5-46] A,=^-The volume of each particle, and therefore the total volume of the landslide, may increase in a time step due to entrainment. Following from Equation [5-44], the volume of material entrained by particle / , AVn which travels a distance, As,, in a time step is: [5-47] AV, = EshiAjAsj Substituting Equation [5-46] into Equation [5-47] produces the following incremental form of the natural exponential growth equation: [5-48] AV,=EsV;As; As path material is added to each particle it is simultaneously removed from the local fixed reference grid using the following relationship, based on Equation [5-10]: [5-49] A 6 , = - \u00C2\u00A3 A ^ . 109 where Abk is the change in bed depth at node k and j = 1 to n are nearby particles. Equation [5-49] ensures mass balance between the landslide and the bed. As an option, DAN3D can account for the resulting changes in bed elevation, which may cause self-channelization that can influence the lateral spreading of open slope landslides such as debris avalanches. Entrainment proceeds within a fixed grid cell until the available material is exhausted or the entire landslide passes. One of the limitations of this simple approach is the use of an additional user-specified parameter, Es, which must be adjusted by trial-and-error in order to obtain a reasonable distribution of entrained material. For example, using an average growth rate, Es, for a specific entrainment zone produces results similar to Takahashi et al. (1992) and Hungr (1995), in which the entrainable material is distributed throughout the entire length of a passing landslide. A useful preliminary estimate of the average growth rate for a specific entrainment zone can be obtained from the following natural exponential growth equation: [5-50] Vf=Voexp{E~s) where Vf is the estimated total volume of the landslide exiting the zone, VQ is the estimated total volume of the landslide entering the zone and S is the approximate average path length of the zone. Rearranging Equation [5-50] gives: [5-51] Es= X 'J 0 / S Under supply-limited conditions, Vf- can be approximated by summing V0 and the total volume of entrainable material in the entrainment zone, assuming an appropriate path width. In such a case, using a value of Es higher than the estimated average growth rate would distribute more of the available material toward the flow front, which could be used to simulate plowing. Alternatively, a relatively higher rate could be assigned to the margins. 110 5.7 Intensity parameter interpolation SPH provides a framework for estimating the value of intensity parameters anywhere within the domain of interest, which is required for landslide hazard mapping (cf, Chapter 1). Each particle carries information about the local flow properties, including depth and velocity, as it is advected with the flow. This information can be output directly, for use in scatter-type plots that show the flow properties at each particle location. At the same time, variations of Equation [5-8] can be used to interpolate any of these properties from the moving particles back to the fixed reference grid for output in the form of intensity contour maps or 3D surfaces. This eliminates the need to post-process the particle-centred properties using other software. For example, similar to Equation [5-10], the flow/deposit depth at each grid node can be interpolated using: [5-52] hk=JV^Vkj where j = 1 to n are nearby particles. The following variation of Equation [5-8], which eliminates edge effects near the flow margins, is useful for estimating the value of other flow properties at each grid node: For example, the flow velocity at each node can be interpolated using: [5-54] vk=^-\u00C2\u00B1VjVjWkJ K j=i where vk and v. are velocity vectors in the global coordinate system. Equation [5-53] represents a simple weighted average, in which the influence of a single particle is weighted according to its relative contribution to the local flow depth. I l l The landslide depths and velocities calculated using Equations [5-52] and [5-54] can be output at user-specified intervals. Over the course of a simulation, the program also stores maximum recorded values of these parameters at each grid node. The maximum recorded flow depth data is useful for generating a map of the simulated landslide trimline, which delineates the direct impact area. As described in Chapter 1, additional intensity parameters, such as discharge, impact pressure and kinetic energy, can be derived using the depth and velocity data. DAN3D can be configured to perform these calculations automatically and output the results. 5.8 Implementation The model has been coded in C++ and runs on an IBM compatible personal computer. The program reads spatial data in the form of user-created grid files, which contain the following data at nodal locations on a global fixed grid: 1) the bed elevation; 2) the landslide source depth (the distance between the original ground surface and the rupture surface, measured in the vertical direction and increased as appropriate to account for dilation/bulking of the failed mass); and 3) the depth of the erodible bed (measured in the bed-normal direction). The user inputs the required rheological and control parameters, such as the time step interval. A variety of data files can be output at user-prescribed intervals and processed by a 3D surface-modelling program, such as Surfer (Golden Software Inc.) or a suitable GIS package. The final output takes the form of hazard intensity contour maps (cf, Chapter 1) and static or animated isometric views. The duration of each simulation depends on the number of particles used, the size of the global reference grid, the length of the time step interval, the frequency of output required and the nature of the landslide in question. For a typical simulation using 4000 particles, a 2 GHz computer processes approximately one time step per second. A flowchart summarizing the main program functions is shown in Figure 5.7 (for proprietary reasons, the actual program code cannot be published in this thesis). 112 I Start j Input grid file names^ Read input files Input rheological and control parameters^ Distribute particles Output initial diagnostic information Proceed? ~NO^ f~End YES t + At Jf t< maximum t_ ~YES~\" NO End Correct particle positions and calculate bed inclination angles using bicubic interpolation T. Correct particle velocity vectors and calculate radii of curvature of path Calculate smoothing length Transform neighbouring particle positions into local reference frames Interpolate flow depths Interpolate flow depth gradients Calculate total bed-normal stresses Calculate basal shear stresses Interpolate incremental tangential strains Increment internal stresses Calculate local accelerations Output diagnostic information and grid files if / is an output time^ Update particle velocities Smooth particle velocities 31 Update particle positions Increment particle volumes Update bed depths Figure 5.7. Flowchart showing the main functions in DAN3D. The column on the right represents the loop that is performed during each time step. 113 5.8.1 Initial conditions From the source depth data, the program calculates the area covered by the source landslide and the initial volume of slide material contained in each cell of the global reference grid (using the trapezoidal rule). The total initial landslide volume is divided evenly among N particles. The initial average particle density is estimated and Equation [5-23] is used to calculate the initial particle smoothing length. Two different particle distribution methods can then be used, depending on whether the source landslide depths are irregular or uniform. When the initial depths are irregular, the particles are systematically placed in the source area, one at a time, according to the following procedure: 1) the program searches the reference grid for the cell requiring the largest volume of material (to match the input data); 2) a particle is placed at a random location within this cell; 3) Equation [5-52] is used to update the required nodal depths within the radius of influence of the particle, and 4) the required cell volumes are updated. The procedure is repeated until all of the particles have been placed. This process is analogous to filling a number of bins (the reference grid cells) with equal buckets of material (the particles) until each bin contains the correct amount of material. The particles are smooth and therefore typically influence more than one grid cell in their neighbourhood, or analogously, the buckets of material can spill into neighbouring bins. This method is capable of reproducing complex landslide configurations and does not require a matrix inversion. The random placement of the particles within each cell prevents particle stacking and alignment. The input depths can typically be reproduced with an error of less than a few percent, which tends to be slightly larger near the landslide margins, where the smooth particles have a tapering effect. An example is shown in Figure 5.8. When the initial depth is uniform, the following alternative particle distribution method can be used. Rather than being placed randomly within each cell, the particles are placed evenly on an orthogonal grid. The spacing of the grid is a multiple of the fixed global grid spacing and is adjusted by the program to maximize the number of particles used in the simulation (up to the limiting value). This method is capable of reproducing a specified uniform initial depth with negligible uniform errors, except near the margins. 114 (a) given source depth (b) interpolated source depth Figure 5.8. Example of a DAN3D initial particle distribution for a very irregular source comparing (a) the given source depth with (b) the interpolated source depth following distribution of the particles (the images are isometric views of 3D surface plots). In this case, the source area is approximately 1 km x 1 km and the maximum depth is about 185 m. The maximum discrepancy between (a) and (b) is 10 m. 4000 particles were used with a particle smoothing coefficient of 5 = 5. 5.8.2 Time stepping With the additional assumptions introduced in this chapter, Equations [5-45] and [5-31] are the x and y direction momentum balance equations used in DAN3D. Numerical integration of these equations determines how the particles move in each time step. Again, the x direction is aligned with the local direction of motion, but when a particle is stationary this direction is undefined. In such a case, a predictor step, in which the local x-axis is arbitrarily oriented, is used to predict the direction of potential motion (based on the vector sum of local driving stresses) and the local axes are re-oriented accordingly. At the start of a time step, each particle occupies a new position determined in the previous time step. Due to the irregular 3D topography, the updated position of the particle may be temporary because it does not necessarily correspond with the new local sliding surface. Given its temporary global coordinates, the new position of the particle on the 3D sliding surface is approximated using a bicubic interpolation, which ensures smoothness of both elevation and slope values across the global reference grid (Press et al. 2002). For gradual bed curvature and a reasonably small time step, this correction is small. The particle velocity vector is then re-oriented so that it lies in the new local tangent plane. 115 This correction also provides information for approximating the local bed-normal radius of curvature of the path in the direction of motion: As. [5-55] R,= \u00E2\u0080\u0094 Pi where As,, is the distance travelled by the particle during the time step and /?; is the angle between the original and corrected velocity vectors. Following this correction, each particle occupies a position on the 3D sliding surface and is centred at a reference column. The particle smoothing length is calculated using Equation [5-23] (using flow depth values from the previous time step) and the positions of neighbouring particles are transformed into the local coordinate system (x,y,z)r The updated local flow depth is calculated using Equation [5-10] and the local depth gradient is calculated using Equations [5-13] and [5-14]. The local total bed-normal stress is calculated using Equation [4-29] and the local basal shear resistance is calculated using one of Equations [4-34] to [4-42], as governed by the specified local basal rheology and its associated parameters. The incremental tangential strain state is approximated using the simplified interpolation method described in Section 5.4 and the local stress coefficients are incremented using Equation [5-42]. If the updated stress state breaches the failure envelope, the stress coefficients are corrected accordingly using Equations [5-32] to [5-37]. The momentum balance Equations [5-45] and [5-31] are then used to solve for the instantaneous Dv Dv local acceleration, \u00E2\u0080\u0094 - and \u00E2\u0080\u0094 - , of the particle and its associated reference column. Dt Dt Assuming constant local acceleration during a short time step, At, the column/particle velocity is updated by a forward difference approximation: 116 [5-56] Dv. Dt DV Dt where velocities at the beginning of the time step are primed. The user-prescribed time step interval must be short enough to ensure stability, which is achieved empirically by running a few simulations and systematically decreasing At until no significant difference is observed in the model results. To minimize so-called \"particle penetration\" (Monaghan 1989), which tends to occur near strong shocks (cf, Chapter 6), especially in the absence of real or numerical viscosity, a velocity smoothing option is available in DAN3D. The following correction can be applied to the updated column/particle velocities (based on the \"XSPH variant\" proposed by Monaghan 1989): [5-57] v,=v,+CAv, where v(. is the velocity vector in the global coordinate system (and the vector on the right side is derived from Equation [5-56]), C is a user-specified velocity smoothing coefficient and Av;. is the following velocity correction: [5-58] A v , = \u00C2\u00A3 V, 7=1 h,+h; Equation [5-58] is simply a variation of Equation [5-8] that ensures conservation of linear momentum. Essentially, use of Equation [5-57] means that the particles move with a velocity that is closer to the average velocity in their neighbourhood (Monaghan 1992), and the magnitude of the velocity smoothing coefficient determines how much the velocities of 117 neighbouring particles influence the central particle. The sensitivity of the model to C is investigated in Chapter 6. Velocity smoothing introduces some numerical diffusion, which appears to smooth out strong shocks, increase stability and reduce the tendency for particles to line up in the downstream direction in channelized reaches of the path. But too much smoothing can impose unrealistically streamlined and coherent behaviour. Following velocity smoothing, the particle is advanced along the local tangent plane to its temporary position by a central difference approximation: The particle volume is updated using Equation [5-48] and the bed depth at neighbouring grid nodes is updated using Equation [5-49]. If the bed elevation change option is enabled, the bed elevation is correspondingly updated. When all of the particle positions are updated, the model proceeds to the next time step. 5.9 Discussion The numerical method described in this chapter includes a number of innovations that were driven by the objectives outlined in Chapter 2. Smoothed Particle Hydrodynamics, which was originally developed for the simulation of compressible flow in 3D, has been adapted for incompressible, depth-averaged flow across an irregular 3D surface. The method is fully-Lagrangian, which means that the relatively simple momentum balance equations derived in Chapter 4 can be solved efficiently. At the same time, mass balance is satisfied without the need for a mesh. A unique interpolation method based on strain gauge rosette theory permits strain estimation within this meshless framework, and because the method provides strain magnitude, rather than simply direction, the internal stresses can be incremented in proportion. This incremental approach has a physical basis but may also improve numerical 118 stability. The updated internal stresses are limited by a frictional yield criterion, which allows for non-hydrostatic, anisotropic conditions, and are coupled to the basal shear stress, which remains a generic function within the equivalent fluid framework. Bicubic interpolation of surface elevations and slopes ensures smoothness of these values across the global reference grid and permits the estimation of local path curvature and resulting centripetal acceleration, which can strongly influence the internal and basal stresses. In addition, a simple entrainment algorithm has been incorporated. The rates of mass and momentum transfer are controlled by a single parameter, the natural exponential growth rate, which has a physical significance and is related in a simple way to the established concept of erosion velocity. The availability of entrainable material is controlled by a spatially-variable bed depth, which permits the simulation of supply-limited conditions.- Rheology changes associated with entrainment can also be implemented. Finally, the SPH-based method can produce hazard intensity maps without requiring separate post-processing of data. This capability eliminates the need for supplementary software and improves the overall efficiency of the runout analysis process. A l l of these features have been designed to be modular and as simple as possible, which should facilitate future modifications and allow the model to evolve with increasing understanding of extremely rapid landslides. 5.10 Conclusion An original numerical method has been developed to solve the system of governing equations that were derived in Chapter 4. Together, the governing equations and the solution method satisfy all of the objectives that were outlined in Chapter 2, making DAN3D unique among existing landslide continuum dynamic models. With these capabilities, the model can be used to analyze a wide variety of landslide types at any scale, and the model output is suitable for direct use in landslide risk assessment. But before it can be used in practice, its accuracy must be validated using controlled tests. Model testing is the subject of the next chapter. 119 6 T E S T S 6.1 Introduction The DAN3D system of governing equations and corresponding numerical solution method were introduced in the previous two chapters. The program incorporates all of the essential features described in Chapter 2, but its ability to accurately solve the governing equations and its applicability to real cases must still be verified. Testing allows a somewhat clearer distinction to be made between model adaptability and model accuracy (e.g., Iverson 2003). This is particularly important within a calibration-based framework such as the equivalent fluid approach, in which arbitrary adjustment of resistance parameters can often produce good-looking results. i However, analytical solutions of the complete system of equations are not possible, and so only the most basic form of DAN3D (i.e., excluding the influences of irregular terrain, internal strength and material entrainment) can be tested this way. The full algorithm can only be evaluated (not strictly tested) by comparison with controlled experiments. At the same time, hypothetical experiments and parametric analyses can provide valuable insight into basic model behaviour and sensitivity. The purpose of this chapter is to validate and evaluate DAN3D by comparison with both analytical and experimental results. A series of parametric analyses are included to demonstrate the influences of the various control and rheological parameters. The ability of the model to provide accurate first-order runout predictions is also demonstrated. 6.2 ID dam-break 6.2.1 Description The ID dam-break problem provides a basic test for the computational algorithm presented in Chapter 5. It has been used previously by Hungr (1995), Wang and Shen (1999), 120 Mangeney et al. (2000) and Denlinger and Iverson (2004) to validate other dynamic models. Stoker (1957) presented an analytical solution to the classical I D dam-break problem, which involves zero bed slope, zero internal friction (i.e., hydrostatic internal stresses) and zero basal resistance. Mangeney et al. (2000) generalized the classical dam-break solution to account for bed slope and basal friction. In the general case, at time t after failure the bed-normal flow depth, h, at horizontal location X between the crest o f the flow at X = (2 -tjgh0 cos a H\u00E2\u0080\u0094(gsina-gcosatanb is the bulk basal friction angle (Equation [4-41]). The general solution converges on the classical solution when a and [K\u00C2\u00A3. For example, in Figure 6.1a the smoothing length at \u00C2\u00A3 = 30 s was approximately 5.3 m, and the flow depth at the lead particle, which had separated from the others, was therefore approximately 0.53 m; theoretically, it should have been lower. This issue was slightly exaggerated in these examples because the value of the smoothing length was influenced considerably by the large number of particles that remained closely spaced behind the crest, strongly skewing the smoothing length towards a lower value. The lOOx vertical exaggeration in the figures also exaggerates this discrepancy. However, it should be noted that the asymptotic profile predicted by the analytical solution is unrealistic, as a bore would be expected to develop at the flow front in a real fluid (Chow 1959) or a geological material. Second, the analytical solution predicted a sharp change in the profile at the crest (a singularity), which the model smoothed. Increasing the number of particles used in the simulation or employing a spatially-variable smoothing length, as in Wang and Shen (1999), would increase the resolution and improve the results. As mentioned in Chapter 5, these improvements are ongoing, but at present the number of particles is limited to 4000 and a satisfactorily robust spatially-variable smoothing method has not yet been developed. In spite of these limitations, the generally good correspondence between the simulated and analytical results demonstrates the ability of the basic algorithm to accurately solve the governing mass and momentum balance equations. No calibration was required; the input parameters for both the numerical and analytical methods were identical. 123 analytical \u00E2\u0080\u00A2 \u00E2\u0080\u00A2\u00E2\u0080\u00A2 numerical t = 30 s -500 i n r i i i i | r 500 1000 Horizontal Distance (m) a) a = 0 \u00C2\u00B0 , h =0\u00C2\u00B0 I I 1500 2000 \u00E2\u0080\u0094 analytical \u00E2\u0080\u00A2 \u00E2\u0080\u00A2\u00E2\u0080\u00A2 numerical CL Q -500 500 I I I | I 1000 1500 Horizontal Distance (m) b) a = 30\u00C2\u00B0, = 20\u00C2\u00B0 2000 analytical \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 numerical -500 \"i r n i | i i r n | i i i i | r 0 500 1000 1500 Horizontal Distance (m) c) a =30\u00C2\u00B0, = 10\u00C2\u00B0 2000 Figure 6.1. Analyses of the dam-break problem, (a) The classical problem with a = 0\u00C2\u00B0 and b = 10\u00C2\u00B0. The flow depth profiles are shown at 10 s intervals. 124 6.3 Slump test 6.3.1 Description The collapse of a cylindrical column of material, a 2D analogue of the dam-break problem, provides a simple test of grid-dependency (Denlinger and Iverson 2004). The results should be symmetrical regardless of the specified orientation of both the global and local reference systems. A similar test was used by Denlinger and Iverson (2004) to validate their dynamic model. 6.3.2 Methodology The full DAN3D algorithm was used in this test and in every other example that follows in this thesis. The initial radius and height of the cylindrical column were set to 0.354 m and 0.5 m, respectively, similar to the initial conditions used by Denlinger and Iverson (2004). The initial geometry is shown in Figure 6.2. To simulate the behaviour of a typical dry sand, a frictional basal rheology (Equation [4-41]) with a bulk basal friction angle of ., the basal rheological parameters (which depend on the selected basal rheology) and, i f applicable, the entrainment growth rate, Es. The comparative analyses presented in this section demonstrate the general behaviour of the model and its sensitivity to these parameters. 6.4.2 Methodology The sliding surface and source geometry were the same in each case in order to facilitate comparison. The simulations involved the motion of a hypothetical material across a simple 3D surface comprising a 45\u00C2\u00B0 incline and a horizontal plane, which were joined by a curved transition zone. The bed-normal source depths, h0 (in metres), were defined by: [6-5] h0 =0.15(l-(s 0 /0.3) 2 ) where (Ora. is the internal friction angle, d>b is the bulk basal friction angle (for frictional basal resistance), / is the basal friction coefficient (for Voellmy basal resistance), \u00C2\u00A3 is the turbulence parameter (for Voellmy basal resistance) and Es is the entrainment growth rate within a specified entrainment zone. 130 Figure 6.7. Initial configuration of the parametric analyses. The starting plane was inclined at 45\u00C2\u00B0 and connected to the horizontal plane by a 1 m radius smooth transition zone. The source depths (Equation [6-5]) are shown at 1 cm intervals. 6.4.3 Results and discussion: number of particles and particle smoothing As described in Chapter 5, continuum simulation is achieved through discretization of the governing equations, but a sufficiently large number of computational elements (particles) are required to capture the behaviour at every important location within the slide mass. At the lower limit, using only one particle, the model can only simulate the motion of the centre of mass, without accounting for internal deformation. Increasing the number of particles increases the resolution of the continuum method. The present upper limit is N = 4000, based on the speed and memory limitations of current personal computers. A l l of the analyses presented in this thesis were performed using either 4000 or 2000 particles (2000 was a previous upper limit). Interpolation plays a key role in this numerical method and the particle smoothing coefficient, B, which governs the interpolation distances and therefore the degree of smoothing, is expected to have a strong influence. As shown in Equation [5-23], its influence is in turn affected by the number of particles. This series of analyses demonstrates the related influences of control parameters N and B . The parameter values used in each 131 analysis are shown in Table 6.1. A frictional basal rheology (Equation [4-41]) was implemented and hydrostatic conditions were imposed (i.e., . = 0\u00C2\u00B0). The results are shown as depth contour maps in Figure 6.8 (parts (a) to (h) correspond with the respective cases listed in Table 6.1). With lower values of B, the interpolated depths were irregular (the contour lines appear wavy), and these irregularities appeared to propagate to subsequent time steps. As the particle smoothing coefficient increased, the smoothness of the depth interpolation at each time step also increased, resulting in a more stable and symmetrical simulation. But it appears to be possible to produce too much smoothing with a sufficiently high smoothing coefficient, an effect that was evident in run (h), which exhibited too much tapering at the margins. Decreased lateral spreading was also exhibited as the smoothing coefficient increased, a result of lower interpolated depth gradients. At the extreme, sufficiently high smoothing would produce negligible depth gradients and therefore negligible particle interactions, resulting in zero spreading. A high smoothing coefficient can also cause inefficiency, as the number of particles included in the necessary summations increases exponentially with the value of B . The key is to optimize B to produce sufficiently smooth depth interpolations while maximizing efficiency, but the results in Figure 6.8 suggest that the optimum value changes with the value of N . In this particular series of analyses, the optimum values appeared to be about 5 = 4 for N = 2000 (Figure 6.8d) and B = 6 for N = 4000 (Figure 6.8g). These two simulations produced very similar, reasonably smooth and stable results. The initial and final smoothing lengths (Equation [5-23]) for each simulation are listed in Table 6.2. The similarities between the smoothing lengths in simulations (d) and (g) are notable. It is possible that the ideal amount of smoothing is related to the smoothing length, irrespective of the number of particles. However, each case is unique and no general relationship can be drawn from these results. At present, it is therefore recommended that the particle smoothing coefficient, B, be adjusted by the user until the initial depth interpolation appears smooth, but still corresponds with the given geometry (probably in the range of B = 4 to 6). At the least, this approach limits instabilities that may compound as a result of irregularities introduced by the numerical method at startup. 132 Run N B P initial (tn) ^ final (m) (a) 4000 3 0.022 0.038 (b) 2000 3 0.031 0.052 (c) 4000 4 0.029 0.049 (d) 2000 4 0.041 0.067 (e) 4000 5 0.036 0.061 (f) 2000 5 0.052 0.081 (g) 4000 6 0.044 0.071 (h) 2000 6 0.062 0.093 Table 6.2. Initial and final smoothing lengths corresponding to analyses (a) to (h) presented in Section 6.4.3 and Figure 6.8. 6.4.4 Results and discussion: internal friction As described in Chapter 2, the internal strength of landslide material has an important mitigating influence on spreading; because of it, landslides do not spread out or contract as readily as fluids. This series of analyses demonstrates the influence of the internal friction angle, $ , which governs the limiting internal stress states in DAN3D (cf, Chapter 5). The parameter values used in each analysis are shown in Table 6.1. A frictional basal rheology (Equation [4-41]) was implemented. Note that run (g) in Section 6.4.3 and Figure 6.8 is considered part of this series, as the same set of parameter values were used, but with = 0\u00C2\u00B0 to simulate hydrostatic conditions. The results are shown as depth contour maps in Figure 6.9 (parts (a) to (d) correspond with the respective cases listed in Table 6.1). A higher value of $ resulted in slower initial spreading, as the internal stresses initially tended to the active state in both the longitudinal and lateral directions. However, longitudinal convergence during motion through the concave transition zone resulted in the development of passive stresses in the longitudinal direction, which increased the longitudinal spreading forces in this segment of the path. The result was a longer deposit (the final position of the flow front increased, while the final 133 position of the tail decreased). The final location of the centre of mass did not change very much, although it did increase slightly as . increased, due mainly to the degree of lateral spreading. In this case, the tendency to active lateral stress states along most of the path resulted in less lateral spreading with higher values of $ . Less lateral spreading in turn resulted in less energy dissipation, which allowed the centre of mass to travel farther (cf, Chapter 2). The prevailing stress states recorded at each time interval during run (c) are shown in Figure 6.10. The limiting values of the stress coefficients in this case are shown in Table 6.3. Simultaneous longitudinal contraction and lateral expansion through the concave transition zone resulted in anisotropic conditions (i.e., simultaneous longitudinal passive stresses and lateral active stresses). Note that passive longitudinal stress conditions near the flow front did not necessarily persist through the entire runout zone. Run (c) represents the \"base case\" for the rest of the analyses presented in this section, which all used $ = 35\u00C2\u00B0 . k x (mill/ max) k y {min/max) *, ( \u00E2\u0084\u00A2o-3.26 ^ ( m , \u00E2\u0080\u009E )=0 .91 *,(\u00C2\u00BB\u00E2\u0080\u00A2\u00C2\u00BB> =0.71 *\u00E2\u0080\u009E(-\u00C2\u00BB, =0.36 Table 6.3. Limiting values of the stress coefficients in run (c), in which b=25\u00C2\u00B0 and # = 35\u00C2\u00B0 (cf, Equations [5-32] to [5-34] in Chapter 5). 6.4.5 Results and discussion: internal stiffness As described in Chapter 5, the dimensionless stiffness coefficient, D, controls the strain-dependent rate of the transition between active and passive internal stress states. At one extreme, a value of D = 0 produces hydrostatic conditions, as the stresses cannot change from their initial values, which are assumed to be hydrostatic in DAN3D. At the other 134 extreme, an infinite value of D produces instantaneous stress response, with no intermediate states possible between the active and passive limits. In the original D A N model (Hungr 1995), the value of the stiffness coefficient was based on values measured in granular soil behind retaining walls (after Terzaghi and Peck 1967). Different values were used depending on whether the deformation was contractional or expansional: D = (kp -ka)/0.05 in contraction and D = [kp -ka)/0.025 in expansion, where kp and ka were user-specified values of the passive and active earth pressure coefficients, respectively. For the limiting stress coefficients that applied in this series of simulations (Table 6.3), the appropriate value of D should therefore fall within the range of about 25 to 100. A value of D = 10000, representing the approximate equivalent of infinite stiffness, commonly assumed by other workers (cf, Chapter 2), was also tested for comparison. The parameter values used in each analysis are shown in Table 6.1. A frictional basal rheology (Equation [4-41]) was implemented. Note that run (g) in Section 6.4.3 and Figure 6.8 and run (c) in Section 6.4.4 and Figure 6.9 are considered part of this series, as they represent comparable simulations using D = 0 and D - 200, respectively. The results are shown as depth contour maps in Figure 6.11 (parts (a) to (d) correspond with the respective cases listed in Table 6.1). A higher stiffness coefficient produced a longer deposit, due to the relatively faster transition between the active longitudinal stress state on the inclined plane and the passive state through the concave transition zone. The faster transition increased the duration of the passive state, which gave the flow front a longer \"push\" that caused it to travel slightly farther. The otherwise minor differences between these simulations suggest that the use of different stiffness values in contraction or expansion would not make a significant difference and is probably not necessary. Somewhat surprisingly, the use of a pseudo-infinite stiffness coefficient did not appear to introduce numerical instability in this case (cf, Chapter 5). On the contrary, run (d) appeared to produce even more orderly and symmetrical results than the other analyses. This may have been in part due to the very simple geometry, which caused only a couple of cycles 135 between active and passive states. The process of stress transition may play a more important role over more complex topography. But it is also possible that infinite stiffness actually improves stability by ensuring an immediate reduction in tangential stresses during expansion and an immediate increase in tangential stresses during contraction, which is an extremely responsive form of damping. On the other hand, infinite stiffness does not allow for the possibility of stress states intermediate between active and passive, which could possibly develop during pseudo-steady-state motion in long, uniform reaches of a landslide path. Unfortunately, a general conclusion on this issue cannot be drawn from this series of results alone. A default value of D = 200 is therefore recommended, pending further work. 6.4.6 Results and discussion: velocity smoothing As described in Chapter 5, velocity smoothing, which introduces some numerical diffusion, is an option that smoothes out strong shocks by reducing so-called particle penetration (Monaghan 1989). It may also increase numerical stability and improve the behaviour of the model in channelized reaches of the path by reducing the tendency for particles to line up in the downstream direction. This series of analyses demonstrates the influence of the velocity smoothing coefficient, C , which governs the degree of velocity smoothing in DAN3D (cf, Equation [5-57] in Chapter 5). The parameter values used in each analysis are shown in Table 6.1. A frictional basal rheology (Equation [4-41]) was implemented. Note that run (c) in Section 6.4.4 and Figure 6.9 is considered part of this series. The results are shown as depth contour maps in Figure 6.12 (parts (a) to (d) correspond with the respective cases listed in Table 6.1). An increase in the velocity smoothing coefficient tended to reduce the amount of initial spreading, as particles near the margins became increasingly influenced by the conditions further within the flow. This behaviour could potentially be exploited to simulate initial cohesion in real landslides, including rock slides in the early stages of fragmentation. Velocity smoothing also smoothed out the shock, or hydraulic jump (e.g., Chow 1959), that formed near the tail of the flow during deposition; the tail of the deposit was not as steep and did not travel as far with increased smoothing, since particles could not penetrate as easily into the material that had deposited ahead. 136 It appears that it is possible for this form of numerical damping to produce too much smoothing. For example, in run (d) the material may have exhibited too much cohesiveness, considering that the friction angles used in this case were appropriate for dry sand with no cohesion. Although a universal value of C cannot be recommended based on these results, a value up to about C = 0.01 is probably appropriate. It should be noted that the final location of the centre of mass was not affected significantly. A very slight decrease with increasing velocity smoothing (about 1 mm horizontally for each 0.005 increment of C ) can be attributed to the shock smoothing as well as a corresponding decrease in energy dissipation due to decreased lateral spreading. But otherwise, linear momentum was completely conserved. 6.4.7 Results and discussion: basal rheology Modifications to the basal rheology and its parameters can have a dominant influence on the results, a characteristic that is demonstrated in this series of analyses. In the first two analyses, the frictional rheology (Equation [4-41]) was used to simulate a range of granular flow behaviour between approximately saturated and dry conditions. In the last two analyses, the Voellmy rheology (Equation [4-42]) was used with friction coefficients and turbulence parameters within a range of common calibrated values (e.g., Hungr and Evans 1996). The parameter values used in each analysis are shown in Table 6.1. As described in Chapter 4, a number of basal rheologies have been implemented in DAN3D, but the frictional and Voellmy rheologies have been used exclusively in the examples presented in this thesis. The results are shown as depth contour maps in Figure 6.13 (parts (a) to (d) correspond with the respective cases listed in Table 6.1). With the frictional basal rheology, an increase in the bulk basal friction angle reduced both spreading and translation of the slide mass, resulting in a deeper deposit. With b = 20\u00C2\u00B0 it spread out relatively thinly and covered a larger area. In both 137 analyses, as in all of the previous frictional analyses presented in this section, the final location of the centre of mass was in the proximal part of the deposit. Significantly different behaviour resulted with use of the Voellmy basal rheology. While the frictional components of resistance in runs (c) and (d) were less than those in runs (a) and (b) (friction coefficients 0.1 and 0.2 are equivalent to bulk basal friction angles of about 5.7\u00C2\u00B0 and 11.3\u00C2\u00B0, respectively), the added velocity-dependent resistance mitigated spreading and downslope acceleration during the early stages of motion. This added resistance had a relatively higher influence near the flow margins, as the turbulence term in Equation [4-42] dominates when the flow depths are low (infinite resistance results when the flow depth is zero, which occurs directly at the margins). Gray and Tai (1998) recognized the existence of this singularity and proposed an alternative to the Voellmy model to permit motion of the margins. In DAN3D, the interpolated flow depths at particle locations are always non-zero, but this effect can still be significant at very low flow depths. In this case, it allowed the deeper part of the flow to attain relatively higher velocities on the source slope, which created a steep flow front followed by a tapered tail, similar to the classic profile of a debris flow surge, and allowed the centre of mass to travel farther into the distal part of the deposit, in contrast to the frictional examples. The relatively lower frictional resistance then allowed the material to spread out thinly over a large area. The dominant influence of the turbulence term is evident when comparing cases (c) and (d). Even though a smaller friction coefficient was used in case (c), the mass did not travel as far because the lower turbulence parameter significantly reduced its momentum. Note the inverse relationship between the turbulence parameter and resistance, and that the Voellmy model converges on the frictional model when the turbulence parameter is infinite. 6.4.8 Results and discussion: entrainment As described in Chapter 2, mass and momentum transfer during entrainment of path material can have an important influence on landslide dynamics. On one hand, the volume of the landslide increases, which can result in increased flow and deposit depths as well as 138 increased impact and deposit areas. On the other hand, momentum transfer contributes resistance, which can have an effect on the degree of spreading and translation. This series of simulations demonstrates the effects of mass and momentum transfer during entrainment. For comparison, two cases were considered: one involving no entrainment and the other involving entrainment that doubled the volume over a short distance of 0.3 m (Figure 6.7) using a growth rate of Es=2.22 m\"1. In this example, centripetal acceleration was neglected. It would have otherwise contributed a velocity-squared component to the internal and basal stresses in the concave transition zone, making it more difficult to highlight the differences between the two cases. The parameter values used in each analysis are shown in Table 6.1. A frictional basal rheology (Equation [4-41]) was implemented. The results are shown as depth contour maps in Figure 6.14. It is useful to compare these results with the idealized sliding block model described in Chapter 2 (both the single and multi-block cases). Recall that in the case of the single sliding block, the inclination of the energy line would equal the dynamic basal friction angle i f centripetal acceleration was neglected. Analogously, in run (a) the inclination of the line connecting the initial and final positions of the centre of mass should have been slightly less than I I I I 1 1 1 t = 0.5 s 1 1 1 I I I I I I I I I I 1 1 1 t= 1 .0s < I I I I ^ \u00E2\u0080\u00A2 I I I I I I I I I I I I t= 1.5 s I I I I i = 40\u00C2\u00B0 Figure 6.9 (continued). Parametric analyses varying the internal friction angle, (c) . = 35\u00C2\u00B0, (d) d). = 40\u00C2\u00B0. The results are discussed in Section 6.4.4. 146 Figure 6.10. Prevailing longitudinal (x direction) and lateral (y direction) stress states recorded at each time interval during run (c) shown in Figure 6.9. The results are discussed in Section 6.4.4. 147 t - 0 s sjm s lope 1 1 1 ulated \u00C2\u00A7 1 err transit I f low depths i intervals) on z o n e 0.5 m o I I I I 1 1 1 t = 0 .5 s 1 1 1 I I I I I I I I I l 1 1 1 t = 1.0 s < I I I I I I I I I I I I I I I I t = 1 . 5 s I I I I I I I I I \u00C2\u00AE I I I I I I I t = 2 .0 s I I I I I I I I I I I I I I I I t = 2 .5 s I I I I I t = 0 s I I I I I i i i I I I t = 0 . 5 s I I I I I I I I I i i i i I I I t = 1.0 s < I I I I ^ ' I I I I I ) I i i i I I I t = 1.5 s I I I I i i l I I I \u00C2\u00AE i i i i I I I t = 2 .0 s I I I I I I I I I i i i i I I I t = 2 . 5 s I I I I I s a) D = 2 5 b) D = 50 Figure 6 .11 . Parametric analyses varying the dimensionless stiffness coefficient, (a) D = 25 . (b) D = 5 0 . The results are discussed in Section 6 .4 .5 . 148 t - 0 s ejrn ^ / (d w slope 1 1 1 ulated \u00C2\u00A7 1 err ransiti I flow depths intervals) on zone 0.5 m O I I I I 1 1 1 t = 0.5 s 1 1 1 I I I I I I I I I I 1 1 1 t= 1.0 s I I I I I I I I t i l l I I I t= 1.5 s I I I i i I I I I s I I I I I I I t = 2.0 s i i I I I I I I I I t = 2.5 s i I I I I 0 i i i i i i i r c) D = 100 t = 0s I I I I I I I I I I I t = 0.5s I I I I I I I I I I I I I I I I t= 1.0 s < I I I I ^ I I I I I t i l l I I I t= 1.5s I i I I << I I I I I i i i i i I I I t = 2.0 s I I I I I i i i i 0 i i i i I I I t = 2.5 s I i i i i d) D = 10000 Figure 6.11 (continued). Parametric analyses varying the dimensionless stiffness coefficient, (c) D = 100. (d) D = 10000 (a pseudo-infinite value). The results are discussed in Section 6.4.5. 149 t ~ 0 S s j m Mr ((i slope 1 1 1 ulated \u00C2\u00A3 1 err xansiti I flow depths intervals) on zone 0.5 m o I I I I t = 0.5 s 1 1 1 I I I I I I I I 1 1 1 t= 1.0 s c I I I I : I I I I I I I I I t= 1.5 s I I I I { I I I I I I I I I I I I t = 2.0 s I I I I ( I I I I I \u00C2\u00A9 I I I I I I I t = 2.5 s I ( I I I I \u00C2\u00A9 t = 0s I I I I I I I I I I I t = 0.5 s I I I I I I I I I I I I I I I I t= 1.0s < I I I I I I I I I I I I I t= 1.5 s I I I I { I I I I I \u00C2\u00AE I I I I I I I t = 2.0 s I I I I ( I I I I I \u00C2\u00A9 I I I I I I I t = 2.5 s I ( I I I I \u00C2\u00A9 a) C = 0.005 b)C = 0.01 Figure 6.12. Parametric analyses varying the velocity smoothing coefficient, (a) C = 0.005 (b) C = 0.01. The results are discussed in Section 6.4.6. 150 t = 0 s sim 0tr (<* w s lope 1 1 1 ulated g 1 cm transit I f low depths i intervals) on z o n e 0.5 m O I I I I 1 t = 0 . 5 s 1 1 1 I I I I I I I I I I 1 1 1 t = 1.0 s I I I I I I I I I I I I I I I I t = 1.5 s I I I I I I I I I \u00C2\u00A9 I I I I I I I t = 2 . 0 s I I I I (j I I I I I \u00C2\u00A9 I I I I I I I t = 2 .5 s I (j I I I I \u00C2\u00A9 t = 0 s I I I | I I I I I I I t = 0 . 5 s I I I I I I I I I t i l l I I I t = 1.0 s < i i i I I I I I I I I I I i i t = 1.5 s I I I I ( I I I I I I I I I I I I t = 2 .0 s I I I I I I I I \u00C2\u00A9 I I I I I I I t = 2 . 5 s I. I I I I \u00C2\u00A9 c ) C = 0 . 0 1 5 d ) C = 0.02 Figure 6.12 (continued). Parametric analyses varying the velocity smoothing coefficient, (c) C = 0.015 . (d) C = 0.02. The results are discussed in Section 6.4.6. 151 t - 0 s sirr Mr ((* slope 1 1 1 ulated % 1 err xansiti I flow depths intervals) on zone 0.5 m O I I I I t = 0.5s 1 1 1 I I I I I I I I I 1 1 1 t= 1.0 s <( I I I | % I I I I I 1 1 1 1 I I I t= 1.5 s I I I I I 1 1 1 1 \u00C2\u00A9 1 1 1 1 I I I t = 2.0s I I I I I 1 1 1 1 \u00C2\u00A9 1 1 1 1 I I I t = 2.5 s I 1 1 1 1 \u00C2\u00A9 I L t = 0s t= 1.0s \u00E2\u0080\u0094 I \u00E2\u0080\u0094 I \u00E2\u0080\u0094 [ t = 1.5 s t = 2.0s t = 2.5 s J L a) frictional b= 20\u00C2\u00B0 i r b) frictional b= 30 i r Figure 6.13. Parametric analyses varying the basal rheology and its parameters. A frictional rheology was implemented in these two analyses, (a) (pb = 20\u00C2\u00B0. (b) I I I I t = 0 . 5 s 1 1 1 I ntre o I I I I I f m a s s I I I I 1 1 1 t = 1.0 s I I I I I I I I I J) I I I I I I I t = 1.5 s I I I I I I I I I @^ I I I I I I I t = 2 . 0 s I I I I I I I I I \u00C2\u00AE I I I I I I I t = 2 .5 s I I I I I \u00C2\u00AE a) E s = 0 rrf1 b) E s = 2 .22 rrf' Figure 6 .14. Parametric analyses neglecting centripetal acceleration and varying the entrainment growth rate, (a) Es = 0 . (b) Es - 2.22 m\" 1. The results are discussed in Section 6.4 .8 . 154 simulated erosion depths 1 c m intervals) \u00C2\u00AB ( - s mulated entrainment z o n e Figure 6.15. Erosion map corresponding to run (b) shown in Figure 6.14b. The results are discussed in Section 6.4.8. Preliminary simulations of real landslides using DAN3D showed that there is a tendency for particles to line up in the downstream direction in highly channelized reaches. In such cases, the model essentially reverts to 2D and is unable to properly account for lateral momentum balance (cross-stream gravity driving stresses are still accounted for but internal pressure gradients are not). This could be a problem, particularly in channel bends, where superelevation and potential avulsion may not be accurately simulated. Resolution is thought to be the main issue. As the particles stretch out in the downstream direction, gaps form between them. Since the particles are constantly drawn to the thalweg of the channel by gravity, when they find a gap they simply fall in line. Increasing the number of particles (i.e., increasing the resolution) decreases the gap spacing. Velocity smoothing may also help by forcing the particles to move together downstream, rather than individually towards the channel thalweg. These hypothetical experiments demonstrate the influence of the number of particles and velocity smoothing on the behaviour of the model in channels. 6.5.2 Methodology A hypothetical straight channel 350 m long by 50 m wide (in plan dimensions) was used in this test. The channel was inclined at 20\u00C2\u00B0 in the downstream (global X) direction and its parabolic cross-section was defined by: 6.5 Channelization 6.5.1 Description 155 [6-6] AZ = 0.0472 where AZ is the vertical elevation difference between the thalweg at Y = 0 and the sliding surface at cross-stream location Y. The centre of the source, 20 m long by 20 m wide (in plan dimensions), was located at global coordinates (X = 0, Y = 0) and its depth was specified to give zero initial depth gradients in both the downstream and cross-stream directions, producing initial conditions similar to a dam-break. The resulting source volume was approximately 1120 m 3 . The chute configuration and initial particle positions are shown in Figure 6.16. Figure 6.16. Initial configuration of the hypothetical particle channelization tests. The parabolic chute was inclined at 20\u00C2\u00B0. A frictional basal rheology (Equation [4-41]) with a bulk basal friction angle of . = 0\u00C2\u00B0 (to simulate hydrostatic conditions). The particle smoothing coefficient was fixed at 5 = 4. The number of particles, N, and the velocity smoothing coefficient, C, were varied and the following three cases were considered: a) N = 2000 and C = 0; b) N = 4000 and C = 0; and c) N = 4000 and C = 0.01. 156 6.5.3 Results and discussion The particle positions at time t = 30 s are shown in Figure 6.17. Particle alignment was evident in run (a), as the particles had formed a single file line by this time. In contrast, the increased resolution in run (b) allowed the stream to remain up 4 or 5 particles wide. The velocity smoothing in run (c) also appeared to help maintain streamlining and the general order of the flow. Of course, with a sufficiently long channel and continued downstream acceleration and spreading, the particles would eventually align regardless of resolution, but these improvements may be enough to prevent alignment in most practical situations, in which the length and uniformity of channelized reaches are limited. Similar improvements have been observed in simulations of full-scale events. Further sensitivity analyses, possibly using a supercomputer, may be beneficial. 25-oH -25-50 100 150 200 a) N = 2000, C = 0 250 300 25-0--25- ~1 1 1 1\u00E2\u0080\u0094 50 100 150 200 b) N = 4000, C = 0 250 300 25-0H -25-0 50 100 150 200 c) N = 4000, C = 0.01 250 300 Figure 6.17. Particle positions in channel at time / = 30 s. (a)yV = 2000 and C = 0. (b) /V = 4000 and C = 0. (c) # = 4000 and C = 0.01. The particle alignment problem improved with increased resolution and velocity smoothing. 157 6.6 Chute experiments 6.6.1 Description Gray et al. (1999) presented the results of small-scale laboratory experiments similar to the hypothetical experiments of Section 6.4. Dry quartz chips were released from rest on a 40\u00C2\u00B0 chute with parabolic cross-section, which was joined to a level runout surface by a smooth transition zone. The entire sliding surface was painted for an even finish. A series of companion experiments using different materials was presented by Wieland et al. (1999). These controlled experiments provide a simple test for the full DAN3D algorithm. 6.6.2 Methodology DAN3D was used to simulate the experiment involving dry quartz chips (Gray et al. 1999). Detailed descriptions of the experimental setup, presented by Wieland et al. (1999), were used to configure the model. The centre of the source, with a volume of approximately 21.6 litres, was located at global coordinates (X = 0,7 = 0). The initial conditions are shown in Figure 6.18. A frictional basal rheology (Equation [4-41]) was implemented and the bulk basal friction angle was calibrated so that the simulated final position of the front approximately matched that of the experiment. Experimental trimlines presented by Gray et al. (1999) were digitized to facilitate comparison. The internal friction angle was set to (j). = 40\u00C2\u00B0 to match the value measured by Gray et al. (1999) using independent tests. The following control parameters were used: N = 4000 particles, particle smoothing coefficient 5 = 6, velocity smoothing coefficient C = 0.01 and stiffness coefficient D = 200. 158 Figure 6.18. Initial configuration of the tests based on experiments presented by Gray et al. (1999) and Wieland et al. (1999). The parabolic chute was inclined at 40\u00C2\u00B0 and connected to the horizontal plane by a smooth transition zone. The source depths are shown at 1 cm intervals. 6.6.3 Results and discussion A comparison of the simulated and measured results is shown in Figure 6.19. A calibrated bulk basal friction angle of $,=28\u00C2\u00B0, close to the laboratory-measured value of ^ = 3 0 \u00C2\u00B0 reported by Gray et al. (1999), produced the best results (use of the laboratory-measured value produced slightly less runout). The difference between the calibrated and measured values is within the likely measurement error discussed by Wieland et al. (1999). However, it was necessary to offset the simulated results by + 0.1 s to synchronize them with the experiment. It is possible that the model, which assumes hydrostatic initial conditions, was inaccurate at startup, but experimental measurement error could have also contributed to this discrepancy. In the experiment, the material was initially contained under a perspex cap, which was raised \"rapidly\" (Gray et al. 1999) at time / = 0 s. But it is not clear exactly how rapidly the cap was raised, whether it might have remained in contact with the material for a split second after the timer had started, or whether it might have actually pushed the material downslope. The possibility of experimental error seems to be supported by the simulation results presented by Gray et al. (1999), which exhibited a similar time lag. Such time lags 159 are not important in hazard assessment, as long as the distribution of intensity can be accurately simulated. Figure 6.19. Simulation of the straight runout chute experiment presented by Gray et al. (1999). The dots represent the measured experimental trimline (digitized from Gray et al. 1999). Note that the simulation has been offset + 0.1 s for synchronicity with the experiment (i.e., / = 0.5 s as shown corresponds with / = 0.6 s in the simulation). Once synchronized, the model accurately simulated the general shape of the flow during motion. The predicted distribution of the final deposit was reasonably accurate, with the 160 model predicting slightly thicker proximal deposits starting closer to the toe of the chute and more radial spreading in the distal runout zone. The proximal discrepancies may be attributable to the shock (hydraulic jump) that developed in the experiment as the tail flowed into the main deposit, as discussed in Section 6.4. The velocity smoothing that was imposed in this case may have smoothed out the shock and reduced particle penetration, but the complex process of deposition near the tail, including intense transverse shearing and accompanying lateral spreading, was not accurately simulated. In the experiment, it is also possible that material deposited initially at the toe of the chute helped to smooth the slope change through the transition zone, reducing local centripetal acceleration and resistance and allowing the material to travel farther. Similar \"self ramping\" has been observed in granular runup experiments (e.g., Chu et al. 1995). The lower-than-measured calibrated bulk basal friction angle seems to support this notion, as a lower value may have been required to offset the influence of higher centripetal acceleration and resistance through the transition zone. Nevertheless, the generally close correspondence between the simulation and the experiment demonstrates the ability of the model to simulate real flows using tightly-constrained rheological parameters. Again, the bulk basal friction angle was the only parameter that was adjusted in this case. 6.7 Deflection experiments 6.7.1 Description Similar small-scale experiments, but involving flow deflection, were conducted at the University of British Columbia. These tests demonstrate the ability of the model to simulate strongly curved flow and to provide accurate first-order predictions following calibration. Dry, rounded, well-graded polystyrene beads were used in the experiments. Approximately 20 litres of material was released from a storage box onto a chute with variable slope, which allowed control of the approach width and velocity. The material flowed from the chute onto a 20\u00C2\u00B0 approach slope and was then deflected by an inclined plane oriented obliquely to the 161 approach direction. The deflection angle, 8 (the plan angle between the approach direction and the intersection of the deflection and approach planes), and the deflection plane true dip angle, atl, were adjusted by rotating the deflection plane about a fixed pivot point. Both planes were made of plywood sheets covered with smooth sheet metal, which was marked with a 10 cm square grid. A photograph of the experimental apparatus is shown in Figure 6.20. Four different experimental configurations were used, as shown in Table 6.4. Each configuration was repeated to verify precision. The experiments were photographed at 30 frames per second using a mounted high-speed video camera. The positions of the flow margins relative to the reference grid were recorded at each interval. Velocities were calculated from the positions of the flow front in successive frames. Figure 6.20. Photograph of the experimental apparatus. Experiment 8 (\u00C2\u00B0) o 1 and 2 60 33 3 and 4 45 33 5 and 6 45 45 7 and 8 60 45 Table 6.4. Configurations used in the deflected flow experiments. 8 is the deflection angle and ad is the deflection plane true dip angle. 162 6.7.2 Methodology Experiment 1 was arbitrarily chosen to be the prototype for calibration of the numerical model. The box used to contain and release the material at start-up could not be replicated by the digital sliding surface, due to its infinitely sloping sidewalls. As a result, a virtual release chute and initial distribution of material were used by the model (initial conditions similar to those described in Section 6.6). By trial-and-error adjustment of the virtual chute angle, the position, width and velocity of the simulated flow front were synchronized with the experiment at the start of the 20\u00C2\u00B0 approach slope. The observed approach velocity was approximately 4.5 m/s in all of the experiments. A frictional basal rheology (Equation [4-41]) was implemented and both the bulk basal and internal friction angles were systematically adjusted until the simulated maximum runup distance matched that of the prototype experiment. The calibrated model was then applied to prediction under the other three experimental configurations. 6.7.3 Results and discussion The calibrated simulation of Experiment 1 is shown in Figure 6.21. The model was calibrated to match the observed maximum runup distance of approximately 56 cm (within the measurement error of +/- 2 cm) using a bulk basal friction angle of b=\4\u00C2\u00B0. This low value (relative to that expected of dry broken rock) implies the presence of pore water pressure near the base of the flow, probably following the onset of entrainment downslope of the source area and almost certainly during impact on the floodplain deposits (in fact, liquefaction of the path material may have occurred at this point, resulting in complete loss of frictional resistance). As such, it represents an average value. The calibrated value of . =35\u00C2\u00B0 produces higher than hydrostatic tangential stresses in directions of compressional strain, corresponding with the passive failure condition (cf, Chapter 2). A trial comparative run using only hydrostatic stress conditions did not produce the same results, although similar runup could be achieved by lowering the basal resistance. Similar to Revellino et al.'s (2004) D A N analysis, the maximum simulated flow velocities recorded along the path also corresponded with independent velocity estimates. In the final frame of Figure 7.7, the simulated flow front reached the village and was still travelling at over 6 m/s. The simulation was stopped at this point because the D E M did not include details within the village (i.e., streets, buildings, etc.) that strongly influenced the real event. 183 Figure 7.7. Simulation of the Cervinara landslide using calibrated input parameters determined by Revellino et al. (2004). The flow/deposit depth contours are at 1 m intervals and the elevation contours are at 10 m intervals. The simulated flow front reached the village in the last frame, when it was still moving at over 6 m/s. 184 It must be stressed that this analysis included only a few changeable input parameters, and so the faithful simulation of the external dimensions of the flow did not result from arbitrary adjustments of multiple parameters. Furthermore, the parameters used were the same ones determined by Revellino et al. (2004) on the basis of D A N analyses of multiple events of similar descriptions, and the results again demonstrate good correspondence between the two models. The use of previously calibrated input parameters makes these results equivalent to a first-order prediction. 7.5 Quindici 7.5.1 Description A similar series of landslides occurred in the Campania region the previous year, in May 1998, causing a total of 161 deaths (Revellino et al. 2004). They were also triggered by prolonged rainfall and involved substantial entrainment of surficial pyroclastic deposits (cf, Figure 2.1 in Chapter 2 for a photograph of events in the area of Sarno, near Quindici). Many of these events were characterized by multiple surges caused by separate failures that were funnelled into the same creek channel. The subject of this analysis is a single surge from one of these events, which jumped out of its channel near the fan apex and impacted the south-eastern part of the village of Quindici. 7.5.2 Methodology Two cases were considered to demonstrate the sensitivity of the model to complex topography: one using input sliding surface data that was digitized from a contour map, and another with slight modifications to this data near the fan apex. As shown in the next section, these modifications were necessary to reproduce the observed channel avulsion. The initial failure was modelled as a 4200 m slide of approximately uniform 2 m depth. Similar to the preceding Cervinara analysis, the basal resistance in the source zone was modelled using a frictional rheology (Equation [4-41]) with a bulk basal friction angle of h = 30\u00C2\u00B0 is shown in Figure 7.15. As expected, the model predicted that a relatively deeper deposit, with an average depth of approximately 7.5 m and a maximum depth of approximately 17 m, would form at the foot of the source slope. The runout distance was underestimated by more than one kilometre. Clearly, volume and rheology change capabilities are essential to the successful simulation of this event. 194 Figure 7.15. Simulation of the Nomash River landslide neglecting entrainment and the associated change in resistance. Instead, a frictional basal rheology was used for the duration of motion. This simulation did not match the observed behaviour of the real event. The deposit depth contours are at 2 m intervals and the elevation contours are at 20 m intervals. 7.7 Zymoetz River The following analysis was undertaken with the help of Nichole Boultbee, Doug Stead and Jim Schwab. 7.7.1 Description At approximately 1:30 am PDT on June 8, 2002, a large debris flow severed the Pacific Northern Gas pipeline at the mouth of Glen Falls Creek in the Zymoetz River valley, 18 km east of the city of Terrace in northern British Columbia, Canada (Figure 7.16). The landslide originated more than 4 km from the pipeline, initiating as a slide in volcanic bedrock. The fragmenting rock mass overrode and entrained snow and saturated surficial material on the source slope and in the cirque basin at its toe, before continuing as a debris flow that travelled more than 3 km down a sinuous, confined valley. An aerial view of the landslide path is shown in Figure 7.17. On the source slope and in the upper part of the cirque basin the landslide deposited a thick, talus-like apron of fragmented rock, with material coming to rest on slopes exceeding 30\u00C2\u00B0. In the valley downslope of the cirque basin, the landslide exhibited much higher fluidity, with 195 sweeping superelevations and motion persisting on slopes lower than 6\u00C2\u00B0. Debris was deposited in lobes along Glen Falls Creek and in a large fan that extended across and dammed the Zymoetz River. Figure 7.17. Aerial view of the Zymoetz River rock slide - debris flow. (A) source; (B) cirque basin; (C) Glen Falls Creek; (D) Zymoetz River. (Aerial photograph taken June 12, 2002, courtesy of Jim Schwab, Northern Interior Forest Region, Ministry of Forests, Province of British Columbia. Copyright \u00C2\u00A9 2005 Province of British Columbia. A l l rights reserved. Reprinted with permission of the Province of British Columbia.) 196 The total economic losses associated with the landslide exceeded C D N $30 million, including the direct costs of repairing the pipeline as well as lost timber and forestry access time (Schwab et al. 2003). Most of the damage occurred downslope of the cirque basin and can be attributed to the high mobility of the event. With a total volume of approximately 1.4x10 m and a fahrboschung less than 17\u00C2\u00B0, the event plots well below the mean empirical relationships proposed for catastrophic landslides by Scheidegger (1973) and L i (1983), suggesting exceptional mobility for its size (Figure 7.18). However, the Zymoetz River rock slide - debris flow is only one of several recent events of this type and magnitude in British Columbia (e.g., Schwab et al. 2003; Hungr and Evans 2004a). Descriptions of the Zymoetz River landslide were presented by Schwab et al. (2003), Boultbee (2005) and Boultbee et al. (2006). The following description concentrates on aspects of the mobility of the landslide. 1 x 10 5 i i iiiiij n r m Volume (m3) ' 'l'x 109 Figure 7.18. Plot of landslide volume vs. the tangent of the fahrboschung (angle of the line connecting the crown of the source to the toe of the deposit, as measured along the curving path). The Zymoetz River event (volume = 1.4xl06 m 3 , tan(fahrboschung) = 0.31) plots below the empirical relationships proposed by (a) Scheidegger (1973) and (b) L i (1983). 197 An enlarged aerial view of the proximal part of the path is shown in Figure 7.19. The rock slide initiated in heavily jointed volcaniclastic rock mid-slope on the east wall of a cirque (location A on Figure 7.19). A prominent joint set dipping downslope at about 45\u00C2\u00B0 formed the main sliding surface, while a sub-vertical backscarp extended up to the crown at elevation 1390 m (Boultbee 2005). Regional precipitation and temperatures were above and below normal seasonal levels, respectively, in the months prior to failure, as reflected in a 200% above normal snowpack (Schwab et al. 2003), but no significant coinciding climatic trigger was recorded. Tectonic deformation, alteration, glacial over-steepening, frost action and pore pressure changes have all been suggested as possible factors contributing to the apparent long-term, progressive degradation of the rock mass, leading ultimately to failure (Schwab et al. 2003; Boultbee et al. 2006). There is evidence of at least one previous rock slide immediately adjacent to the present source scar (Schwab et al. 2003), which may have created an apron of coarse talus at the toe of the source slope. Figure 7.19. Enlarged aerial view of the proximal part of the path. (A) source; (B) exposed till; (C) thick deposits over snow; (D) thin deposits over snow in bend; (E) location of end moraine near elevation 880 m; (F) forested bedrock knob. (Aerial photograph taken June 12, 2002, courtesy of Jim Schwab, Northern Interior Forest Region, Ministry of Forests, Province of British Columbia. Copyright \u00C2\u00A9 2005 Province of British Columbia. A l l rights reserved. Reprinted with permission of the Province of British Columbia.) 198 Based on calculations using pre- and post-event digital elevation models (DEMs), which are described in more detail in the next section, approximately 720,000 m 3 of rock was displaced from the source area during the initial slide. Accounting for an estimated 25% bulking of the failure due to fragmentation (e.g., Sherard et al. 1963) gives a source volume of approximately 900,000 m 3 . The first phase of motion was over snow, colluvium (talus and old rock slide debris) and bare rock. The snowpack was up to 5 m deep in the cirque basin at the time of the event (Schwab et al. 2003). A snow-covered area of about 200,000 m 2 was traversed by the landslide within the cirque basin between the source zone and the location of the snowline, near a prominent end moraine located at approximately elevation 880 m. Some rock debris slid across the surface of the snow while some plowed into it and pushed it forward or became embedded. At the same time, snow was also entrained, as evidenced by exposed substrate at some locations. Figure 7.20 is a view up the source slope several months after the event. A ridge of till (a lateral moraine), which was covered with snow prior to the event, is exposed in the centre of the image (location B on Figure 7.19). Striations in the till that were carved by the passing rock suggest that it was stripped of surficial material, including snow and likely loose colluvium. On either side of the till ridge and further upslope are talus-like deposits of coarse rock avalanche debris (rockfall or secondary failures following the main event contributed to the amount of this material to an unknown degree). Approximately 100,000 m 3 of coarse rock debris was spread over the source slope, which averaged more than 30\u00C2\u00B0. Based on field measurements, another approximately 500,000 m of coarse rock debris, averaging about 3 m thick, came to rest in the upper part of the cirque basin immediately below the source area (Boultbee et al. 2006), some of it on top of snow (Figure 7.21 and location C on Figure 7.19). Plowed snow was deposited along the margins of this rock debris. 199 Figure 7.20. View from toe of source slope toward source area (at top of photograph). The till exposed in the centre (location B on Figure 7.19) is surrounded by coarse rock avalanche deposits. (Photograph courtesy of Nichole Boultbee) Figure 7.21. Rock avalanche deposits on top of snow in the upper part of the cirque basin (location C on Figure 7.19). (Photograph courtesy of Nichole Boultbee) 200 Contrasting with the coarse debris is a curving segment adjacent to its left margin on the outside of the first major curve of the path (location D on Figure 7.19). Here, the original snow cover remained in place, but was marked by strongly-defined curving striations, evidently caused by debris having travelled in an extremely rapid arching motion over the snow. The high point of the arch is almost 60 m above the valley thalweg. Using the Forced Vortex Equation for superelevation (e.g., Hungr et al. 1984), the flow velocity of the debris sliding over snow in this curve was estimated to be approximately 26 m/s (Boultbee et al. 2006). Little or no debris was deposited in this segment. The thicker rock avalanche deposits overlap its proximal and distal margins, suggesting that they arrived only when the curving motion over the snow had been completed. While evidently highly mobile, it is unlikely that this part of the landslide was responsible on its own for triggering the massive entrainment and motion of material into the distal segment of the path. The reason for this is that only a relatively small quantity of the landslide material was evidently involved, which slid over the snowpack but was unable to entrain much of it. The increase in mobility was probably caused by rapid, undrained loading of wet surficial material (e.g., Hutchinson and Bhandari 1971; Sassa 1985) by the leading edge of the more massive main flow. Such liquefaction of path material has long been invoked as a mechanism of long-runout (e.g., Buss and Heim 1881, Sassa 1988, Voight and Sousa 1994, Abele 1997). A very similar mechanism and sequence of transition between a rock slide and debris avalanche at Eagle Pass, British Columbia was described by Hungr and Evans (2004a). It is possible that an end moraine, which appears to demarcate the two main phases of motion (location E on Figure 7.19), affected local drainage of spring meltwater below the snowline, causing saturation of loose material accumulated in the path. At the same time, the thalweg upstream of the end moraine was probably filled with loose, fine-grained soil and organics. After crossing the end moraine, the landslide dropped out of the basin and continued downstream in a second, more fluid phase of motion (Figure 7.22). The flow divided around a forested bedrock knob near elevation 830 m (location F on Figure 7.19), with most of the material branching to the east. It rejoined near elevation 690 m and followed Glen Falls 201 Creek the remaining 2.6 km to the Zymoetz River. An enlarged aerial view of the distal part of the path is shown in Figure 7.23. Figure 7.22. Oblique aerial view of the cirque basin showing the change in deposit composition and morphology near the end moraine. (A) location of end moraine near elevation 880 m; (B) forested bedrock knob. (Photograph taken on June 8, 2002, courtesy of Jim Schwab, Ministry of Forests, Province of British Columbia) Figure 7.23. Enlarged aerial view of the distal part of the path. (A) channel widening due to entrainment at the lateral margins; (B) superelevation around bend; (C) upper end of gorge; (D) fan. (Aerial photograph taken June 12, 2002, courtesy of Jim Schwab, Northern Interior Forest Region, Ministry of Forests, Province of British Columbia. Copyright \u00C2\u00A9 2005 Province of British Columbia. A l l rights reserved. Reprinted with permission of the Province of British Columbia.) 202 The extremely rapid motion along Glen Falls Creek is evidenced by superelevations up to 40 m above the base of the creek channel. Velocities back-calculated from superelevation data in this segment ranged between 14 and 26 m/s (Boultbee et al. 2006). These high velocities continued far into the distal part of the path. For example, the velocity at location B on Figure 7.23 was estimated to be between 16 and 18 m/s (Boultbee et al. 2006). The high mobility of the debris flow was evidently aided by a significant amount of water, as mud splashes high above the creek channel were observed on trees next to the path (Boultbee 2005). The creek channel itself was widened substantially, by up to about 80 m in some locations, as colluvium, topsoil and vegetation were entrained along the margins of the flow (e.g., location A on Figure 7.23). Material was simultaneously deposited in discontinuous lobes up to 4 m thick (Boultbee 2005). Based on field measurements, approximately 200,000 m 3 of material was deposited along Glen Falls Creek (Boultbee et al. 2006), the majority of which stopped just upstream of a narrow bedrock gorge (location C on Figure 7.23). These deposits were finer than the material deposited in the cirque basin and included volcanic source material as well as additional volcanics, limestone, till and some organics derived from the basin, the channel banks and existing channel deposits (Boultbee 2005). The debris flow deposits can be seen in Figure 7.24, which is a view up the channel towards the cirque basin. Incorporation of finer material and additional water is evident in the character (gradation, lithology and water content) of deposits left in the channel downstream of the end moraine and in the distal fan deposit. It is likely that this change in the character of the flowing material increased the mobility of the landslide and allowed motion to continue beyond the cirque basin. 203 Figure 7.24. Oblique aerial view up the channel towards the source area (at top of photograph). The debris flow superelevated around numerous bends and left discontinuous deposits. (A) The flow overtopped a 20 m high ridge near the confluence with a tributary creek (Photograph taken on June 13, 2002, courtesy of Jim Schwab, Ministry of Forests, Province of British Columbia) The Pacific Northern Gas pipeline was severed near the mouth of the creek, igniting an intense fire. Approximately 600,000 m 3 of material reached the Zymoetz River and formed a large fan/dam that blocked the river for about 30 minutes (location D on Figure 7.23) (Schwab et al. 2003). A view of the fan is shown in Figure 7.25. The fine materials blanketing the fan were deposited during a second, smaller debris flow that was witnessed at 10:15 am PDT. This secondary event likely resulted from sidewall destabilization due to erosion by the main event. It appears to have originated from a debris slide that entered the channel below the cirque basin and may have remobilized fresh deposits in its path (Jim Schwab, personal communication). An estimated total volume of 1.4xl06 m 3 was mobilized during the event, about 500,000 m 3 of which was entrained along the path. 204 Figure 7.25. Oblique aerial view of the large fan created by the landslide, which temporary dammed the Zymoetz River. (Photograph taken on June 13, 2002, courtesy of Marten Geertsema, Ministry of Forests, Province of British Columbia) 7.7.2 Methodology Initial back-analyses using D A N (Hungr 1995) suggested that a more realistic simulation of the event could be achieved by accounting explicitly for the complex 3D geometry (Boultbee 2005). Pre- and post-event DEMs at 2 m resolution were provided. The source depths were approximated by subtracting the post- from the pre-event D E M and isolating the probable main failure zone. Data outside of this zone were filtered, leaving a displaced volume of approximately 720,000 m 3 . The isolated source depths were then subtracted from the pre-event D E M to approximate the initial sliding surface elevations. The failure was bulked approximately 25% at the source, for an initial volume of 900,000 m 3 . The data spacing was increased to 10 m for input into the model. 205 The two main phases of motion were simulated as one continuous event, starting as frictional motion of broken rock and transforming into debris flow at the point where significant entrainment of saturated material from the path began (estimated to be near location E on Figure 7.19). Above elevation 880 m, the approximate location of the snowline and the moraine, the basal resistance was governed by a frictional model (Equation [4-41]). The bulk basal friction angle was calibrated to produce a reasonable simulation of the distribution of coarse deposits in the proximal path. Entrainment was neglected during this phase. The first surge passing through the proximal area, which apparently slid over snow, could well have been analyzed as a separate event, even though it may have occurred only seconds before the main detachment. A n alternative method was used: a lower bulk basal friction angle was imposed on a relatively small volume of material at the flow front to simulate the influence of snow on the mobility of the leading edge of the rock slide and reproduce the superelevation against the opposite wall of the cirque basin (location D on Figure 7.19). A value of b = 18\u00C2\u00B0 (based on a separate DAN3D back-analysis using a smaller detachment volume, which was performed by Suzanne Chalindar) was assigned to the first 10% of the material (90,000 m3) to reach elevation 980 m, approximately where the real landslide began to spread out over the deep snowpack in the upper part of the cirque basin. This volume was chosen to represent the apparently small volume of material near the flow front, which ran over but did not entrain the snow. Trial runs showed that the model was not very sensitive to the selected flow front volume, within a reasonable range; a much smaller volume would not cover the observed impact area, while a much larger volume would likely produce large enough flow depths to cause entrainment of the snow. The relatively low value of b used to model this part of the event implies inherently lower dynamic friction between the rock and the snow, which may have also been mediated by pore fluid, possibly supplied by wet or rapidly melting snow incorporated into the flow. This methodology is consistent with the hypothesis that the relatively less mobile but more massive main body of the landslide was responsible for initiating the entrainment and change in behaviour in the cirque basin. Based on the estimated presence of loose, saturated soils in the path upstream of the end moraine, as described in the previous section, entrainment was enabled below elevation 880 206 m, and the Voellmy basal resistance model (Equation [4-42]) was implemented. The two Voellmy resistance parameters were adjusted by trial-and-error to achieve the best simulation of the second main phase of motion in terms of the estimated flow velocities, observed flow trimlines and observed distribution of deposits. To simulate entrainment during this phase, a volume growth rate of Es = 3.3 x 10\"4 m\"1 was specified below elevation 880 m (based on Equation [5-51] assuming Vf = 800,000 m 3 and VQ = 300,000 m 3 , based on the aforementioned volume balance estimates for the real event, and S = 3000 m). The internal friction angle was set to $ =35\u00C2\u00B0 and the following control parameters were used for the duration of motion: N = 4000, B = 4, C = 0.01 and D = 200. 7.7.3 R e s u l t s a n d d i s c u s s i o n The results are shown in Figure 7.26 (Figure 7.26a shows the first 60 s of motion at 10 s intervals, while Figure 7.26b shows 600 s of motion at 100 s intervals). The best simulation of the main proximal phase of motion, based mainly on the distribution of deposits near the toe of the source slope, was achieved using b = 30\u00C2\u00B0. Again, this value corresponds closely to Hsu's (1975) \"expected\" value for sliding rock, and it has been used in several of the preceding cases. In this case, it is also slightly lower than the average slope angle of the failure surface, and so its use approximates limiting equilibrium conditions at initiation. 214 The rheology change was implemented below elevation 1060 m, approximately where the more mobile lobe extended from the main deposit, and the Voellmy friction coefficient and turbulence parameter were set to / = 0.1 and 500 m/s2, respectively. These values were recommended by Hungr and Evans (1996) based on back-analyses of 23 rock avalanches using D A N . Similar calibrated values were determined by Hungr et al. (1998), Ayotte and Hungr (2000), Revellino et al. (2004) and Pirulli (2005). In this case, the friction coefficient approximates the average slope of the creek channel downstream of the main deposit, which allows the distal flow to spread out and deposit there. At the same time, the finite turbulence parameter, which is within the range of values used in the preceding cases, limits the flow velocities along the channel. Although entrainment did occur during the real event, volume change was neglected in this analysis, as the amount of surficial material available along the path was small relative to the source volume. As in previous cases, the internal friction angle was set to $ = 35\u00C2\u00B0 and the following control parameters were used: N = 4000, 5 = 4, C = 0.01 and D = 200. 7.8.3 Results and discussion The results are shown in Figure 7.30. A good simulation of the extent of the impact area and the distribution of deposits was achieved using the preliminary input parameter values. The relatively deep proximal and thin distal deposits were reproduced, with the distal lobe comprising approximately 1 M m 3 of material, as estimated in the real event. Similar to several of the preceding cases, comparative analyses have shown that the observed distribution of deposits cannot be simulated using constant rheology assumptions. The lack of significant runup and superelevation in the real event means that velocity estimates are not available for comparison. However, this feature was also reproduced by the model, suggesting that the simulated velocities were reasonable. The simulated velocity along the distal channel ranged up to about 25 m/s, which is within reason for a debris flow of this magnitude. A lower/higher simulated value could be achieved by reducing/increasing the value of the Voellmy turbulence parameter. 215 Figure 7.30. Simulation of the McAuley Creek landslide using \"off the shelf input parameters. The flow/deposit depth contours are at 5 m intervals and the elevation contours are at 20 m intervals. 216 These good preliminary results suggest that true prediction, involving objective input parameter selection from a database of back-analyses, may eventually be possible. In this case, the most subjective decision was where to implement the change from frictional to Voellmy behaviour. This decision would have been difficult to make without the post-event field evidence. In reality, there was probably no noticeable difference between the surficial deposits upstream and downstream of this location prior to the event. Any surficial deposits that may have contributed to the high mobility of the distal phase could also extend upstream beneath the present main deposit, which remained there perhaps because it was simply deep enough to bridge them. More experience analyzing this type of event may help reduce this subjectivity. 7.9 Kolka The following analysis was undertaken with the help of Olga Tutubalina, Sergey Chernomorets, Dmitry Petrakov, Valeriy Drobyshev and Steve Evans. 7.9.1 Description On September 20, 2002, a large ice/rock avalanche - debris flow occurred in the Caucasus Mountains in North Ossetia, Russia. The event involved the failure of Kolka Glacier at the head of the Genaldon River valley. It impacted a large area over a width of up to 500 m and a total distance of about 36 km, comprising an 18 km long proximal ice/rock avalanche phase, which stopped upstream of the Karmadon Gates (a narrow pass in the Skalistyy Range), and an 18 km long distal debris flow phase. Velocities up to 250 km/hr along the proximal path have been estimated based on flow cross-sections and seismic data (Drobyshev 2003b and 2005). The main deposit formed an ice and debris dam in the Karmadon depression that caused extensive flooding and damage to property upstream, as shown in Figure 1.2 in Chapter 1. Another view of the main deposit area is shown in Figure 7.31. An estimated 125 people were killed in the disaster. 217 Figure 7.31. View of the main deposit of the Kolka Glacier event, looking downstream from the right side of the Genaldon valley towards the Karmadon Gates. (Photograph taken on September 22, 2002, two days after the event, courtesy of Mr. Igor Galushkin, Ministry of Natural Resources, North Ossetia, Russia) The exact trigger for the event is still in dispute. Ice and rock falls onto Kolka Glacier from the northern face of Mt. Dzhimarai-khokh occurred in the two months leading up to the disaster, and a relatively large event of this type was originally thought to have triggered the catastrophic failure (Huggel et al. 2005). However, recently-discovered satellite images taken less than 9 hours prior to the event suggest that a large triggering slope failure did not occur, but that a surge of the rear part of Kolka Glacier was underway (Tutubalina et al. 2005). This suggests the possibility of a \"standing start\" (Olga Tutubalina, Sergey Chernomorets and Steve Evans, personal communication), with catastrophic failure and rapid acceleration perhaps caused by a combination of brittle failure of the ice and rapid loading of water and/or saturated material at its base or immediately downstream. Not disputed is the poor condition of Kolka Glacier at the time, and the possible involvement of a significant amount of meltwater. 218 The following back-analysis was undertaken to produce a simulation of the proximal movement phase in order to aid in the description and detailed analysis of the event. This work is ongoing. A simple two-parameter calibration method is demonstrated. 7.9.2 Methodology A D E M of the sliding surface at 50 m resolution was created using large-scale topographic maps. An outline of the Kolka Glacier source area was digitized from a detailed post-event geomorphic map (Figure 3b in Petrakov et al. 2004a). In accordance with the \"standing start\" scenario described above, a total bulked volume of 130 M m 3 , representing Kolka Glacier and the antecedent accumulation of ice and rock fall debris on its surface, was assumed to start from rest within this area (a uniform vertical source thickness of approximately 85 m was used). This volume is within the range of total volume estimates presented by Drobyshev (2003a), Petrakov et al. (2004b) and Huggel et al. (2005). The estimated volumes of material entrained and deposited along the path were approximately equal and were small relative to the source volume; entrainment was therefore neglected. A Voellmy resistance model (Equation [4-42]) with constant parameters was used for the following three reasons: 1) There was no evidence of a significant change in behaviour along the path (upstream of the Karmadon Gates) that would justify a user-imposed change in rheology or rheological parameters in the simulation. 2) The event initiated, travelled and deposited on very low slope angles, implying a persistently low frictional component of resistance. This behaviour can be simulated with the Voellmy model using a suitably low friction coefficient. 3) The turbulence component of the Voellmy model can be used to limit the simulated flow velocities, which (considering the above requirement for low frictional resistance) would otherwise be too high. 219 The Voellmy friction coefficient, / , and turbulence parameter, were systematically adjusted, as shown in the next section, to produce good simulations of the estimated velocities along the path, the estimated time of arrival of the flow front at the Karmadon Gates and the final distribution of deposits upstream (only a relatively small volume squeezed through the Gates to generate the distal debris flow). The internal friction angle was set to $ = 20\u00C2\u00B0, considered appropriate for a mass composed predominantly of flowing ice, and the following control parameters were used: /V = 4000, 5 = 4, C = 0.01 and D = 200. 7.9.3 Results and discussion The results of a model run using / = 0.05 and E, = 1000 m/s2 are shown in Figure 7.32. This combination of parameter values produced the best simulation of the estimated flow velocities along the Genaldon valley and the arrival time of the front at the Karmadon Gates (approximately 6 minutes). A 2D matrix (for various / and \u00C2\u00A3 ) of maximum simulated velocity plots, which was used to help make this judgment, is shown in Figure 7.33. A l l of these model runs reproduced the sinuous motion of the ice/rock avalanche along the Genaldon valley and the main deposit that dammed the valley in the Karmadon depression. As such, the independent velocity estimates provided a crucial secondary constraint. A comparison of estimated and simulated arrival times along the path is shown in Figure 7.34. A different reference starting point was used in each case, so the DAN3D results have been offset by + 45 s to filter out this discrepancy and facilitate comparison. The Drobyshev (2003b and 2005) curve is slightly steeper, reflecting slightly higher velocities. Similarly higher velocities were recorded by the model in the deeper part of the flow immediately behind the front. This is a characteristic of the Voellmy rheology; the turbulence term produces relatively higher resistance in shallower parts of the flow, including near the front. This depth-dependence results in surging behaviour, with the flow front constantly being overtaken and/or pushed forward by the deeper and faster trailing flow. When this behaviour is taken into account, the Drobyshev and DAN3D results are actually in good agreement. 220 Figure 7.32. Calibrated simulation of the Kolka event at 60 s intervals. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. The sinuous motion of the ice/rock avalanche along the Genaldon valley was reproduced. 221 Figure 7.32 (continued). Calibrated simulation of the Kolka event at 60 s intervals. The flow/deposit depth contours are at 10 m intervals and the elevation contours are at 100 m intervals. The ice and debris dam in the Karmadon depression was reproduced. 222 Figure 7.33. Matrix of maximum velocity plots for model runs with various combinations of / and \u00C2\u00A3 . The velocity contours are at 10 m/s intervals and the elevation contours are at 100 m intervals. The best simulation of the estimated velocities along the path was obtained using / = 0.05 and cf = 1000 m/s2 (middle left). 223 3200 \u00E2\u0080\u0094 i Drobyshev 2800 \u00E2\u0080\u0094 2400 \u00E2\u0080\u0094 E c o 2000 1600 \u00E2\u0080\u0094 1200 800 0 100 200 300 400 500 Time (s) Figure 7.34. Comparison of arrival times along the path estimated by Drobyshev (2003b and 2005, based on estimated maximum velocities) and DAN3D using / = 0.05 and \u00C2\u00A3 = 1000 m/s (based on simulated flow front arrival times). The DAN3D timing has been offset by + 45 s to facilitate comparison, as different reference starting points were used in each case. A comparative analysis using the same pair of input parameters was performed to examine the triggering scenario postulated by Huggel et al. (2005). An initial downslope velocity of 25 m/s was imposed on the 130 M m 3 source to account for momentum transfer from a large ice or rock avalanche from Mt. Dzhimarai-khokh. This had only a minor influence on the results. The flow accelerated slightly faster but soon reached similar maximum velocities in the Genaldon channel. The mass arrived at the Karmadon Gates sooner but formed an almost identical deposit. Note that this result may simply reflect the character of the Voellmy rheology and does not necessarily discount the possibility of a triggering slope failure. The Voellmy rheology is largely empirical, and so physical justification for parameter selection is often difficult. This case is no exception. The calibrated friction coefficient / = 0.05 is lower than the average slope gradient in the starting zone (approximately 0.1), which allowed the mass to accelerate rapidly from rest without the need to impose an initial velocity. At the same time, it is slightly higher than the average slope gradient in the Karmadon depression, which allowed the mass to deposit there. While justification for the 224 exact value of / used is still difficult, liquefaction of material at the base of the flow could account for such a low component of frictional resistance, and a continuous supply of water from ice melt during motion could explain its persistence. In contrast, the turbulence parameter is almost purely empirical, and can only be justified on the basis of estimated and simulated flow velocity comparisons. The matrix method demonstrated in Figure 7.33 is an effective way to perform these types of comparisons when two independent parameters must be calibrated. Similar matrices can be constructed to compare simulated trimlines and deposit distributions. Efficiency can be gained by stepping through the matrix systematically, rather than running every possible combination of parameters. This approach has been used in other back-analyses and sensitivity analyses. 7.10 Other cases A number of other cases have been analyzed for private consulting assignments and are detailed in various technical reports. A brief description of each case follows: 1) A calibration back-analysis of the 1946 Six des Eaux Froides rock avalanche in the Swiss Alps was performed. D A N (Hungr 1995) was used for preliminary calibration, with good correspondence to the subsequent DAN3D results. Various combinations of frictional and Voellmy rheologies were attempted (employing the matrix calibration method outlined in the previous section). The best simulation of the event, in terms of the observed extent of the impact area and distribution of deposits, was achieved using a constant Voellmy rheology with parameters / =0.13 and \u00C2\u00A3 = 450 m/s2, close to Hungr and Evans' (1996) recommended values of / =0.1 and E, = 500 m/s2. 2) A back-analysis of a small rock avalanche near Preonzo, Switzerland in 2002 and a forward-analysis of a potential similar failure nearby were performed using input parameters derived from an independent D A N analysis. The results were used to 225 cross-check the D A N results and to help assess the risk to facilities and the effectiveness of protective measures located near the base of the slope. A frictional rheology was implemented and (j)b was varied between 30\u00C2\u00B0 and 35\u00C2\u00B0. The results were very sensitive to the value of the bulk basal friction angle in a range close to the mean slope angle. 3) A calibration back-analysis of the 1915 Jane Camp landslide in British Columbia was performed and the results were used to guide a forward-analysis of a potential failure nearby. Similar to the preceding case, a parametric approach was adopted during the predictive phase. Various combinations of frictional and Voellmy rheologies were used with a range of parameter values, and a number of possible entrainment scenarios were examined. 4) A forward-analysis of a potential mine waste dump failure was performed to help assess the risk to expensive equipment located downslope. Different assumed source volumes and locations were analyzed. Simulated flow velocities and depths along a possible protective fence line were output for use in impact pressure and energy calculations. 5) A similar forward-analysis of a potential colluyium failure at another mine site was performed. In this case, a parametric approach was used to assess the sensitivity of the results to the assumed soil conditions. A subsequent analysis was performed to help assess the effectiveness of a berm in the distal path. 6) A calibration back-analysis of the January 2005 landslide in North Vancouver, British Columbia was performed and the results were used to guide a forward-analysis of a potential failure nearby. Forward-analyses using various assumed source volumes originating from within the 2005 scar were also analyzed to help assess the effectiveness of a berm/basin at the toe of the slope. The results were cross-checked using an established analytical method for predicting runup against adverse slopes 226 (e.g., Hungr et al. 1984). Hazard intensity maps of simulated flow velocities and depths were output for use in vulnerability estimation. As an important step in model development and distribution, DAN3D has also been used by other workers, and their input has led to several improvements. Alex Strouth (University of British Columbia) used the model to back-analyze the 2003 Afternoon Creek rock avalanche in Washington State. Suzanne Chalindar (Ecole Polytechnique Federate de Lausanne) used the model to back-analyze the Mystery Creek rock avalanche in British Columbia and the 1970 Huascaran debris avalanche in Peru, and continues to use it to analyze several cases in the Swiss Alps. Tyler Trudel (University of British Columbia) used the model to perform a detailed parametric analysis of the 1970 Huascaran event. Alexandra Campbell, Graham Janson and Heather Taylor (Queen's University, Kingston) used the model as part of an undergraduate design project to assess outburst flood risks at Ruapehu Volcano in New Zealand. 7.11 Discussion The back-analyses presented in this chapter demonstrate the applicability of DAN3D at full-scale to a wide variety of landslide types. The results are summarized in Table 7.1. The importance of each of the features outlined in Chapter 2 has been demonstrated in context. In particular, the analysis of rock slide - debris avalanche type events (Hungr and Evans 2004a) requires the full suite of model capabilities. In each case, the model has been calibrated in some way (either directly or by using D A N as a preliminary tool) to produce a satisfactory simulation of the real event, but it is important to note that this process is not arbitrary. In general, the frictional basal resistance model appears to produce good simulations of dry granular behaviour, while the Voellmy model produces good results when a significant amount of water is apparently involved. As such, the appropriateness of each model for a given case can usually be objectively determined. Although the associated input parameters are not considered real material properties, they are still (often tightly) constrained by the observed and/or estimated behaviour of the real event, 227 as defined by objective calibration criteria. At the same time, the parameters can often be related to both the initial and final limiting equilibrium conditions. There should never be a need to introduce any input that is not supported by field observations or in conflict with previous modelling experience. Case Study Landslide Type N B C D Y; (\u00C2\u00B0) Basal Rheology and Parameters (m 1 ) Before Entrainment After Entrainment Frank rock avalanche 2000 4 0 200 40 frictional Y b = l4\u00C2\u00B0-no change 0 V a l Pola rock avalanche 4000 6 0.01 200 35 frictional no change 0 Cervinara debris avalanche/ flow 2000 4 0 200 35 frictional Yb=30\u00C2\u00B0 Voellmy /=0.07 \u00C2\u00A3 = 2 0 0 m/s 2 0.01 Quindici debris avalanche/ flow 2000 4 0 200 35 frictional Yb=30\u00C2\u00B0 Voellmy /=0.07 \u00C2\u00A3 = 2 0 0 m/s 2 0.01 Nomash River rock slide - debris avalanche 2000 4 0 200 35 frictional Voellmy 7=0.05 \u00C2\u00A3 = 4 0 0 m/s 2 1.9 x 10~3 Zymoetz River rock slide - debris flow 4000 4 0.01 200 35 frictional ^ = 3 1 \u00C2\u00B0 Voellmy 7=0.1 \u00C2\u00A3 = 1 0 0 0 m / s 2 3.3 x io- 4 McAuley Creek rock avalanche 4000 4 0.01 200 35 frictional & = 3 0 \u00C2\u00B0 Voellmy 7=0.1 \u00C2\u00A3 = 5 0 0 m/s 2 0 Kolka ice/rock avalanche 4000 4 0.01 200 20 Voellmy 7=0.05 \u00C2\u00A3 = 1 0 0 0 m / s 2 no change 0 Table 7.1. Parameter values used in each case study. N is the number of particles, B is the particle smoothing coefficient, C is the velocity smoothing coefficient, D is the stiffness coefficient, . is the internal friction angle, b is the bulk basal friction angle (for frictional basal resistance), f is the basal friction coefficient (for Voellmy basal resistance), \u00C2\u00A3 is the turbulence parameter (for Voellmy basal resistance) and Es is the entrainment growth rate within a specified entrainment zone. 228 For this reason, the repeatedly good correspondence between the D A N and DAN3D results is encouraging, and suggests that the already extensive database o f DAN-calibrated cases w i l l continue to be a useful reference. A t the same time, significant efficiency can be gained by using the two models in a complementary manner. Sti l l , a high level of uncertainty remains when it comes to parameter selection, making accurate deterministic runout prediction difficult. A s a result, a parametric approach to forward-analysis has been adopted. Quantifying this uncertainty is the next step, but the database wi l l have to grow substantially before probabilistic parameter selection wi l l be possible. More experience with the model by a larger community of users may accelerate this process and lead to the development of a standardized analysis methodology that may also help reduce subjectivity. 7.12 Conclusion A variety of real landslides have been back-analyzed using DAN3D. The results have further validated the model and in turn provided insight that has aided in event description and analysis. They have also served to expand the existing database of calibrated cases. This experience has been translated into forward-analyses that have been used in practice to help assess landslide risks and provide guidance for the design o f protective measures, within the established framework outlined in Chapter 1. A move towards probabilistic runout prediction is desirable but w i l l require a considerable expansion of the database. A s such, continued back-analysis of real case studies remains a priority. 229 8 CONCLUSIONS 8.1 Introduction As described in Chapter 1, this thesis has been an effort to advance the state-of-the-art in dynamic modelling for the practical purposes of landslide risk assessment and the design of protective measures. To do so, specific objectives were defined and a new continuum dynamic model that satisfies these objectives, DAN3D, has been developed. The new model is based on a number of existing continuum dynamic models, most significantly D A N (Hungr 1995), but includes several innovations, incorporated into a unique system of governing equations and a unique numerical solution method, that make it substantially original. It has been validated using analytical and experimental methods and its applicability at full-scale has been demonstrated using back-analyses of real case studies. This work has been extended in practice to parametric forward-analyses of potential events. The purpose of this final chapter is to summarize these innovations and the main results of this research to date. Recommendations are also presented to guide future work with the model by others. 8.2 Summary of work completed Highlights of the work presented in Chapters 1 to 7 follow: 1) The need for effective methods of performing landslide runout analysis and the potential for continuum dynamic modelling to address this need were identified. Inherent limitations of the 2D model D A N , which has been used extensively, were also identified and a 3D extension of the model was therefore proposed. 2) The following five key objectives for model development were defined: a model intended for practical landslide runout analysis should be able to simulate motion across complex 3D terrain, simulate strain-dependent, non-hydrostatic, anisotropic 230 internal stresses, simulate mass and momentum transfer due to entrainment, simulate variations in rheology along the path and within the landslide and, at the same time, operate in an efficient manner. 3) A comprehensive review of existing continuum dynamic models in the context of these objectives was presented. 4) The depth-averaged governing equations of motion for an \"equivalent fluid\" (Hungr 1995) were formally derived. The resulting equations incorporate aspects of each of the five objectives for model development. 5) Smoothed Particle Hydrodynamics (SPH) was adapted for incompressible, depth-averaged flow across an irregular 3D surface. The method permits discretization of the governing equations within a Lagrangian framework and satisfies mass balance without the need for a computational mesh. 6) A strain interpolation method based on strain gauge rosette theory was developed. The method takes advantage of the particle-tracking nature of SPH and provides estimates of strain magnitude in addition to orientation. A simplified version of the original rigorous method is currently implemented in the model. 7) A corresponding stress redistribution method was developed, which increments the internal stresses in proportion to the estimated strain to account for finite material stiffness. The updated internal stresses are coupled to the basal shear stress (a generic function in the equivalent fluid framework) and limited by a frictional yield criterion, which permits non-hydrostatic, anisotropic stress states in accordance with classical Rankine earth pressure theory. 8) A simple entrainment algorithm based on the concept of natural exponential volume growth was developed. The rates of mass and momentum transfer are controlled by a user-specified growth rate and the availability of material is limited by a user-231 specified, spatially-variable bed-depth. Both parameters can be constrained by actual field evidence. Rheology changes associated with entrainment can also be implemented. 9) An intensity parameter interpolation method was developed to eliminate the need for separate post-processing of particle-centred data. This feature provides output in the form of hazard intensity maps that can be used directly in risk studies to estimate vulnerability, design protection and characterize potential secondary effects. 10) Established analytical solutions of the ID dam-break problem were used to test the accuracy of the basic numerical solution method, with satisfactory results. 11) A hypothetical experiment analogous to a 2D dam-break was used to evaluate the grid-dependency of the numerical solution method. The model did not introduce any grid-dependence beyond the grid-dependent initial particle distribution. 12) Additional hypothetical experiments were used to demonstrate the general behaviour of the model and evaluate its sensitivity to the various input parameters. Interim values for the input control parameters were recommended. 13) Controlled experiments were used to evaluate the ability of the model to produce accurate simulations of real granular flows. Satisfactory results were achieved with appropriate calibration of the input rheological parameters, although the calibrated values were very close to the measured material properties. 14) Additional controlled experiments involving deflected flow (performed at the University of British Columbia) were used to demonstrate the calibration-prediction methodology and the ability of the calibrated model to produce accurate first-order predictions, even when the geometry differs from that of the prototype. 232 15) A variety of full-scale landslide case studies were used to evaluate the full capabilities of the model and the importance of each of the five objectives for model development was demonstrated in context. 16) As in the case of the small-scale experiments, the calibrated model produced accurate simulations at full-scale and, in the process, proved to be a useful tool for event description and analysis. Various systematic and objective calibration methods were demonstrated. In general, independent observations and/or estimates of both the extent and velocity/duration of the real event are required to provide reasonable constraint. Both the initial and final limiting equilibrium conditions can provide further constraint. 17) In general, the frictional basal resistance model appears to be appropriate for simulating dry granular flows, while the Voellmy model appears to be more appropriate for simulating wet or saturated flows. The concept of modelling rock slide - debris avalanches/flows using a frictional proximal phase and a Voellmy distal phase (e.g., Hungr and Evans 2004a) was further demonstrated. 18) Good correspondence between D A N and DAN3D results was repeatedly observed, suggesting that D A N is a useful tool for efficient preliminary calibration of the new model and that the existing database of DAN-calibrated parameters is a useful reference. 19) Sensitivity to complex terrain was observed in several cases, highlighting the importance of accurate topographic input. 20) The seasonal influences of snow and associated meltwater on mobility were noted, which may be important to consider in risk studies. 233 8.3 Recommendations for future work Although the model is considered fit for interim use and has already been applied in practice, several limitations have been identified and further research is warranted. A number of recommendations follow: 1) Limited spatial resolution is an issue, as shown in the dam-break and channelization cases in Chapter 6. Resolution is sacrificed near flow margins and there is a tendency for particles to align in the downstream direction in highly channelized reaches of a path. An increase in the number of particles the model can use and the introduction of a velocity smoothing method has improved the situation, but more work can still be done. Optimization of the code and coinciding advances in computational capacity may permit a further increase in the maximum number of particles, which would be the simplest way to increase resolution. Tests using a supercomputer with a parallelized version of the code may be beneficial. Implementation of a spatially-variable smoothing length (e.g., Wang and Shen 1999) is another option, which has only been investigated briefly to date. In the meantime, these limitations need to be taken into account when determining the appropriateness of the model in a given case and when interpreting the results. 2) The model is not explicitly shock-capturing. Velocity smoothing smoothes out strong shocks but does not replicate actual energy losses. Explicit shock-capturing capabilities may be a good addition to the model, or at least provide a way to investigate the influence of neglecting such energy losses. 3) Validation testing should continue as more analytical methods and experimental results become available. 4) Comprehensive parametric analyses should be performed to expand on the preliminary study presented in Chapter 6. In particular, the influence of internal stiffness, as controlled by the dimensionless stiffness coefficient, D, requires further 234 investigation. Carefully controlled laboratory experiments and/or well-described full-scale case studies would make good prototypes for this type of work, which Suzanne Chalindar and Tyler Trudel have already initiated. Eventually, appropriate values for the control input parameters should be built into the model, reducing the number of user-adjustable parameters. 5) A more thorough look at the influence of different smoothing kernels may be beneficial. For example, the cubic kernel is very similar to the default Gaussian kernel (cf, Chapter 5), but may increase computational efficiency and does not require a cutoff depth to be specified at the flow margins. 6) The rigorous strain interpolation method described in Chapter 5 can be implemented once more realistic/justifiable stress redistribution assumptions can be made. 7) In general, more research into actual entrainment mechanisms and processes is needed. For example, laboratory experiments could be designed to isolate the main controlling factors. Eventually, more complex entrainment relationships based on this type of work and additional field observations can be implemented in the model, with the current assumption of natural exponential volume growth serving as a baseline. The possible influence of self-channelization due to erosion on the spreading behaviour of real events could also be investigated. 8) The influence of rheology changes caused by surficial material in the path of large rock avalanches and their relationship to the plowing mechanism of entrainment could be investigated in more detail. For example, the Frank and Val Pola analyses presented in Chapter 7, which used constant rheology assumptions, could be revisited. The Hope and Cheam rock avalanches near Vancouver, British Columbia are other potential case studies that exhibited this behaviour (e.g., Orwin et al. 2004). 9) The influence of flow heterogeneity requires more research. Rheology variations within the landslide can be implemented in DAN3D, but to date this feature has only 235 been used to a limited extent. The original SPH method can also handle spatial density variations, and this capability could be useful for investigating the effects of density variations caused by entrainment, dilation, comminution, fragmentation and general pore pressure and material property differences. 10) A method that accounts explicitly for airborne/freefall conditions and the associated momentum and energy changes on impact could be implemented. 11) A method that accounts for the initial cohesion of a landslide mass (e.g., prior to significant fragmentation of a rock slide) could be implemented. A progressively-decreasing velocity smoothing coefficient, C , could be used for this purpose, although a progressively-decreasing internal cohesion (additional to the assumed internal friction) may be a more elegant method. 12) A commercial Windows-based version of the model is currently at an advanced stage of development (by others), but will require considerable debugging and testing before it can be distributed. In the meantime, distribution of the original DOS-based version should continue, in order to expand the community of users. 13) The back-analyses presented in Chapter 7 form the basis for a collection of DAN3D-calibrated cases, but significantly more work is required. Continued back-analysis of real events in parallel with the other recommendations on this list is essential. 14) A probabilistic approach to input parameter selection should be investigated once the database of calibrated parameters has reached a sufficiently large size. In the meantime, a parametric approach to forward-analysis is recommended. 15) The application of DAN3D to the analysis of snow avalanche motion could be investigated. 236 8.4 Discussion Clearly, many more improvements can be made, but it is important to avoid getting mired in the process. With ongoing development in mountainous areas, the need for practical landslide risk assessment tools is immediate, and perpetual testing of continuum dynamic models does not serve this need very well. At the same time, it must be acknowledged that landslides are extremely complex phenomena, and computer models cannot be expected to account for all of this complexity. In the end, all models, no matter how complex, are subject to errors and uncertainty. Arguably, more headway can be gained by trying to quantify this uncertainty than by trying to perfect the governing equations or numerical solution methods. In keeping with the approach formalized by Hungr (1995), this thesis has been influenced to a large extent by two key concepts: 1) the idea that continuum dynamic models should be able to account for the most important characteristics of real landslides, while remaining as simple as possible; and 2) the idea that full-scale application is an essential component of the development process. Although there is significant room for improvement, which will inevitably result in increased complexity, I hope that future work with the model will remain faithful to this basic approach. 8.5 Conclusion The continuum dynamic model proposed in this thesis satisfies the overall objective of this work. A 3D extension of the existing 2D model D A N has been developed for the purposes of practical landslide runout analysis. The new model addresses the inherent limitations of the original model but retains all of its key features, which are essential for the simulation of real landslides. In the process, several innovations have been made. The model is based on a new system of governing equations, which are discretized and solved using a new and efficient numerical method. It has been validated at both laboratory and field-scale and its utility has been demonstrated in practice. With continued work, including ongoing back-analysis of real landslide case studies, the new model has the potential for true predictive capability. 237 REFERENCES Abele, G. 1997. Rockslide movement supported by the mobilization of groundwater-saturated valley floor sediments. Zeitschrift fur Geomorphologie, 41: 1-20. Amsden, A . A . and Harlow, F.H. 1970. The S M A C method: a numerical technique for calculating incompressible fluid flows. Los Alamos National Laboratory Report LA-4370. Ayotte, D. and Hungr, O. 2000. Calibration of a runout prediction model for debris-flows and avalanches. In Proceedings of the Second International Conference on Debris-Flow Hazards Mitigation, Taipei. Edited by G.F. Wieczorek and N.D. Naeser. A A . Balkema, Rotterdam, pp. 505-514. Bagnold, R A . 1954. Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear. Proceedings of the Royal Society of London, A , 225: 49-63. Benz, W. 1990. Smooth particle hydrodynamics: a review. In The numerical modelling of nonlinear stellar pulsations. Edited by J.R. Buchler. Kluwer Academic, Dordrecht, pp. 269-288. B G C Engineering Inc. 2000. Geotechnical Hazard Assessment - South Flank of Frank Slide, Hillcrest, Alberta. Report to Alberta Environment. Boultbee, N . 2005. Characterization of the Zymoetz River Rock Avalanche. M.Sc. Thesis, Simon Fraser University, Vancouver. Boultbee, N . , Stead, D., Schwab, J., Geertsema, M . 2006. The Zymoetz River rock avalanche, June 2002, British Columbia, Canada. Engineering Geology, 83(1): 76-93. Brufau, P., Garcia-Navarro, P., Ghilardi, P., Natale, L. and Savi, F. 2000. ID mathematical modelling of debris flow. Journal of Hydraulic Research, 38(6): 435-446. Bursik, M . , Martinez-Hackert, B., Delgado, H. and Gonzalez-Huesca, A . 2003. A smoothed-particle hydrodynamic automaton of landform degradation by overland flow. Geomorphology, 53: 25-44. Buss, E. and Heim, A. 1881. Der Bergsturz von Elm. Worster, Zurich, (in German) Calvetti, F., Crosta, G. and Tartarella, M . 2000. Numerical simulation of dry granular flows: from the reproduction of small-scale experiments to the prediction of rock avalanches. Rivista Italiana di Geotecnico, 34(2): 21-38. Canada Department of Mines. 1917. Frank, Alberta, Map 57 A. Canada Department of Mines, Ottawa. 238 Cannon, S.H. 1993. An empirical model for the volume-change behavior of debris flows. In Proceedings of Hydraulic Engineering '93, San Francisco. Edited by H.W. Shen, S.T. Su and F. Wen. American Society of Civil Engineers, New York, Vol. 2, pp. 1768-1773. Cannon, S.H. and Savage, W.Z. 1988. A mass-change model for the estimation of debris-flow runout. Journal of Geology, 96: 221-227. Cannon, S.H. and Savage, W.Z. 1990. A mass-change model for the estimation of debris-flow runout: a discussion. Journal of Geology, 99: 792-794. Carson, M . A . and Kirkby, M.J. 1972. Hillslope form and process. Cambridge University Press, London. Chen, H. and Lee, C F . 2000. Numerical simulation of debris flows. Canadian Geotechnical Journal, 37: 146-160. Chen, H. and Lee, C F . 2002. Runout analysis of slurry flows with Bingham model. Journal of Geotechnical and Geoenvironmental Engineering, December, pp. 1032-1042. Chen, H . and Lee, C F . 2003. A dynamic model for rainfall-induced landslides on natural slopes. Geomorphology, 51: 269-288. Chiou, M . C , Wang, Y . and Hutter, K. 2005. Influence of obstacles on rapid granular flows. Acta Mechanica, 175: 105-122. Chow, V.T. 1959. Open-channel Hydraulics. McGraw-Hill, New York. 680 pp. Chu, T., Hi l l , G., McClung, D.M. , Ngun, R. and Sherkat, R. 1995. Experiments on granular flows to predict avalanche runup. Canadian Geotechnical Journal, 32: 285-295. Cleary, P.W. and Prakash, M . 2004. Discrete-element modelling and smoothed particle hydrodynamics: potential in the environmental sciences. Philosophical Transactions of the Royal Society of London, A , 362: 2003-2030. Corominas, J. 1996. The angle of reach as a obility index for small and large landslides. Canadian Geotechnical Journal, 33: 260-271. Costa, J.E. 1991. Nature, mechanics, and mitigation of the Val Pola Landslide, Valtellina, Italy, 1987-1988. Zeitschrift fur Geomorphologie, 35(1): 15-38. Coussot, P. 1994. Steady, laminar flow of concentrated mud suspensions in open channel. Journal of Hydraulic Research, 32(4): 535-559. Crosta, G.B., Chen, H. and Lee, C F . 2004. Replay of the 1987 Val Pola Landslide, Italian Alps. Geomorphology, 60(11): 127-146. 239 Cruden, D . M . 1976. Major rock slides in the Rockies. Canadian Geotechnical Journal, 13:8-20. Cruden, D . M . and Krahn, J. 1973. A re-examination of the geology of the Frank Slide. Canadian Geotechnical Journal, 10: 581-591. Cruden, D . M . and Hungr, O. 1986. The debris of the Frank Slide and theories of rockslide-avalanche mobility. Canadian Journal of Earth Sciences, 23: 425-432. Cruden, D . M . and Varnes, D.J. 1996. Landslide types and processes. In Landslides: investigation and mitigation, Transportation Research Board Special Report 247, United States National Research Council. Edited by A . K . Turner and R.L. Schuster. National Academy Press, pp. 36-75. Cundall, P.A. and Strack, O.D.L. 1979. A discrete numerical method for granular assemblies. Geotechnique, 29(1): 47-65. Davies, T.R. 1982. Spreading of rockslide debris by mechanical fluidization. Rock Mechanics, 15: 9-24. Davies, T.R., McSaveney, M.J. and Hodgson, K . A . 1999. A fragmentation-spreading model for long-runout rock avalanches. Canadian Geotechnical Journal, 36: 1096-1110. Davies, T.R. and McSaveney, M.J. 2002. Dynamic simulation of the motion of fragmenting rock avalanches. Canadian Geotechnical Journal, 39: 789-798. Denlinger, R.P. and Iverson, R . M . 2001. Flow of variably fluidized granular masses across three-dimensional terrain, 2. Numerical predictions and experimental tests. Journal of Geophysical Research, 106(B1): 553-566. Denlinger, R.P. and Iverson, R . M . 2004. Granular avalanches across irregular three-dimensional terrain: 1. Theory and computation. Journal of Geophysical Research, 109: F01014. Dent, J.D. and Lang, T.E. 1980. Modeling of snow flow. Journal of Glaciology, 26(94): 131-140. Dent, J.D. and Lang, T.E. 1983. A biviscous modified Bingham model of snow avalanche motion. Annals of Glaciology, 4: 42-46. Drobyshev, V . N . 2003a. Glacial disaster in North Ossetia on 20 September 2002. Vladikavkaz, unpublished manuscript, (in Russian) Drobyshev, V . N . 2003b. Estimation of the velocity for the ice collapse movement along the Genaldon gorge. Vladikavkaz, unpublished manuscript, (in Russian) 240 Drobyshev, V . N . 2005. Analysis and interpretation of the process of Kolka glacier collapse in the Genaldon gorge on 20 September 2002 on the basis of records from Tsey and Fiagdon seismic stations. Vladikavkaz, unpublished manuscript, (in Russian) Egashira, S. and Ashida, K. 1997. Sediment transport in steep slope flumes. In Proceedings of the Roc Japan Joint Seminar on Water Resources. Egashira, S., Honda, N . and Itoh, T. 2001. Experimental study on the entrainment of bed material into debris flow. Physics and Chemistry of the Earth (C), 26(9): 645-650. Erismann, T.H. and Abele, G. 2001. Dynamics of Rockslides and Rockfalls. Springer-Verlag, Berlin. Erlichson, H. 1991. A mass-change model for the estimation of debris-flow runout: a second discussion: conditions for the application of the rocket equation. Journal of Geology, 99: 633-634. Evans, M.W. and Harlow, F.H. 1957. The Particle-in-Cell method for hydrodynamic calculations. Los Alamos National Laboratory Report LA-2139. 76 pp. Evans, S.G. 2001. Landslides. In A synthesis of geological hazards in Canada. Edited by G.R. Brooks. Geological Survey of Canada, Bulletin 548, pp. 43-79. Evans, S.G., Hungr, O. and Clague, J.J. 2001. Dynamics of the 1984 rock avalanche and associated distal debris flow on Mount Cayley, British Columbia, Canada; implications for landslide hazard assessment on dissected volcanoes. Engineering Geology, 61: 29-51. Evans, S.G., Couture, R., Turner, K. and Fuller, T. 2003. The 2002 rock avalanche at McAuley Creek, near Vernon, British Columbia; implications for regional landslide hazard assessment. In Proceedings of the 3 r d Canadian Conference on Geotechnique and Natural Hazards, Edmonton, p. 299. Fannin, R.J. and Wise, M.P. 2001. An empirical-statistical model for debris flow travel distance. Canadian Geotechnical Journal, 38: 982-994. Fell, R., Ho, K.K.S . , Lacasse, S. and Leroi, E. 2005. A framework for landslide risk assessment and management. In Proceedings of the International Conference on Landslide Risk Management, Vancouver. Edited by O. Hungr, R. Fell, R. Couture and E. Eberhardt. A . A . Balkema, London, pp. 3-25. Garcia, R., Lopez, J.L., Noya, M . , Bello, M.E. , Bello, M.T., Gonzalez, N . , Paredes, G., Vivas, M.I. and O'Brien, J.S. 2003. Hazard mapping for debris-flow events in the alluvial fans of northern Venezuela. In Proceedings of the 3 r d International Conference on Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment, Davos. Edited by D. Rickenmann and C L . Chen. Millpress, Rotterdam, pp. 589-599. 241 Gauer, P. and Issler, D. 2004. Possible erosion mechanisms in snow avalanches. Annals of Glaciology, 38: 384-392. Ghilardi, P., Natale, L. and Savi, F. 2001. Modeling of debris flow propagation and deposition. Physics and Chemistry of the Earth, 26(9): 651-656. Gingold, R A . and Monaghan, J.J. 1977. Smoothed Particle Hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 181: 375-389. Gonzalez, E., Herreros, M.I., Pastor, M . , Quecedo, M . and Fernandez Merodo, J A . 2003. Discrete and continuum approaches for fast landslide modelling. In Proceedings of the 1s t International PFC Symposium, Gelsenkirchen. Edited by H . Konietsky. Swets and Zeitlinger, Lisse. pp. 307-313. Govi, M . 1989. The 1987 Landslide on Mount Zandila in the Valtellina, Northern Italy. Landslide News, 3: 1-3. Gray, J.M.N.T. and Tai, Y . C . 1998. On the inclusion of a velocity-dependent basal drag in avalanche models. Annals of Glaciology, 26: 277-280. Gray, J.M.N.T., Wieland, M . and Hutter, K. 1999. Gravity-driven free surface flow of granular avalanches over complex basal topography. Proceedings of the Royal Society of London, A , 455: 1841-1874. Gray, J.M.N.T., Tai, Y . C . and Noelle, S. 2003. Shock waves, dead zones and particle-free regions in rapid granular free-surface flows. Journal of Fluid Mechanics, 491: 161-181. Greve, R. and Hutter, K. 1993. Motion of a granular avalanche in a convex and concave curved chute: experiments and theoretical predictions. Philosophical Transactions of the Royal Society of London, Physical Sciences and Engineering, 342(1666): 573-600. Greve, R., Koch, T. and Hutter, K. 1994. Unconfined flow of granular avalanches along a partly curved surface, Part 1: Theory. Proceedings of the Royal Society of London, A , 445: 399-413. Guthrie, R.H., Evans, S.G. and Hungr, O. 2003. A rock slide - debris avalanche of May 1999, Nomash River, Vancouver Island. In Proceedings of the 3 r d Canadian Conference on Geotechnique and Natural Hazards, Edmonton, p. 289. Harlow, F.H. 1957. Hydrodynamic problems involving large fluid distortions. Journal of the Association for Computing Machinery, 4(2): 137-142. Harlow, F.H. and Welch, J.E. 1965. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, 8(12): 2182-2189. 242 Heim, A . 1932. Bergsturz und Menschenleben (Landslides and Human Lives). Translated by N . Skermer. Bitech Press, Vancouver, 1989. Heinrich, Ph., Boudon, G., Komorowski, J.C., Sparks, R.S.J., Herd, R. and Voight, B. 2001. Numerical simulation of the December 1997 debris avalanche in Montserrat, Lesser Antilles. Geophysical Research Letters, 28(13): 2529-2532. Hsu, K.J . 1975. Catastrophic debris streams (sturzstroms) generated by rockfalls. Geological Society of America Bulletin, 86: 129-140. Huber, A . 1980. Schwallwellen in Seen als Folge von Felsstiirzen. Mitteilung No. 47 der Versuchsanstalt fur Wasserbau, Hydrologie und Glaziologie an der ETH Zurich, pp. 1-122. (in German) Hiibl, J. and Steinwendtner, H. 2001. Two-dimensional simulation of two viscous debris flows in Austria. Physics and Chemistry of the Earth, 26(9): 639-644. Huggel, C , Zgraggen-Oswald, S., Haeberli, W., Kaab, A. , Polkvoj, A. , Galushkin, I. and Evans, S.G. 2005. The 2002 rock/ice avalanche at Kolka/Karmadon, Russian Caucasus: assessment of extraordinary avalanche formation and mobility, and application of QuickBird satellite imagery. Natural Hazards and Earth System Sciences, 5: 173-187. Hungr, O. 1981. Dynamics of rock avalanches and other types of slope movements. Ph.D. thesis, The University of Alberta, 506 pp. Hungr, O. 1990. A mass-change model for the estimation of debris-flow runout: a reply. Journal of Geology, 98:791 -792. Hungr, O. 1995. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Canadian Geotechnical Journal, 32: 610-623. Hungr, O. 1997. Some methods of landslide hazard intensity mapping. In Proceedings of the Landslide Risk Workshop, Honolulu. Edited by D . M . Cruden and R. Fell. A . A . Balkema, Rotterdam, pp. 215-226. Hungr, O., Sun, H.W. and Ho, K.K.S . 1998. Mobility of selected landslides in Hong Kong -Pilot back-analysis using a numerical model. Report of the Geotechnical Engineering Office, Hong Kong SAR Government. Hungr, O. 2000. Analysis of debris flow surges using the theory of uniformly progressive flow. Earth Surface Processes and Landforms, 25: 1-13. Hungr, O. 2002. Analytical models for slides and flows. In Proceedings of the UNESCO/IGCP International Symposium on Landslide Mitigation and Protection of Cultural and Natural Heritage, Kyoto. Edited by K . Sassa. Disaster Prevention Research Institute, Kyoto University, pp. 559-586. 243 Hungr, O. 2004. Landslide hazards in BC: achieving balance in risk assessment. Innovation, Journal of the Association of Professional Engineers and Geoscientists of BC, 8(3): 12-15. Hungr, O., Morgan, G.C. and Kellerhals, R. 1984. Quantitative analysis of debris torrent hazards for design of remedial measures. Canadian Geotechnical Journal, 21: 663-677. Hungr, O. and McClung, D . M . 1987. An equation for calculating snow avalanche run-up against barriers. In Proceedings of the Symposium on Avalanche Formation, Movement and Effects, Davos. International Association of Hydrological Sciences, Publication 162, pp. 605-612. Hungr, O. and Evans, S.G. 1996. Rock avalanche runout prediction using a dynamic model. In Proceedings of the 7 t h International Symposium on Landslides, Trondheim. Edited by K . Senneset. A A . Balkema, Rotterdam, pp. 233-238. Hungr, O., Evans, S.G., Bovis, M . and Hutchinson, J.N. 2001. A review of the classification of landslides of the flow type. Environmental and Engineering Geoscience, 7(3): 1-18. Hungr, O., Dawson, R., Kent, A. , Campbell, D. and Morgenstern, N.R. 2002. Rapid flow slides of coal mine waste in British Columbia, Canada. In Catastrophic landslides: effects, occurrence and mechanisms. Edited by S.G. Evans and J.V. DeGraff. Geological Society of America, Reviews in Engineering Geology 15. pp. 191-208. Hungr, O. and Evans, S.G. 2004a. 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. and Evans, S.G. 2004b. The occurrence and classification of massive rock slope failure. Felsbau, Rock and Soil Engineering, Journal for Engineering Geology, Geomechanics and Tunnelling, 2: 16-23. Hurlimann, M . , Corominas, J., Moya, J. and Copons, R. 2003. Debris-flow events in the eastern Pyrenees: preliminary study on initiation and propagation. In Proceedings of the 3 r d International Conference on Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment, Davos. Edited by D. Rickenmann and C L . Chen. Millpress, Rotterdam, pp. 115-126. Hutchinson, J.N. and Bhandari, R.K. 1971. Undrained loading, a fundamental mechanism of mudflow and other mass movements. Geotechnique, 21: 353-358. Hutter, K. and Savage, S.B. 1988. Avalanche dynamics: The motion of a finite mass of gravel down a mountain side. In Proceedings of the 5 t h International Symposium on Landslides, Lausanne. Edited byC. Bonnard. A A . Balkema, Rotterdam, pp. 691-697. Hutter, K. and Koch, T. 1991. Motion of a granular avalanche in an exponentially curved chute: experiments and theoretical predictions. Philosophical Transactions of the Royal Society of London A , 334:93-138. 244 Hutter, K., Siegel, M . , Savage, S.B. and Nohguchi, Y . 1993. Two-dimensional spreading of a granular avalanche down an inclined plane, Part 1: Theory. Acta Mechanica, 100: 37-68. Iverson, R . M . 1997. The physics of debris flows. Reviews of Geophysics, 35(3): 245-296. Iverson, R . M . 2003. How should mathematical models of geomorphic processes be judged? In Prediction in Geomorphology, Geophysical Monograph Series, Vol . 135. Edited by P.R. Wilcock and R . M . Iverson. American Geophysical Union, Washington, pp. 83-94. Iverson, R .M. , Schilling, S.P. and Vallance, J.W. 1998. Objective delineation of lahar inundation zones. Geological Society of America Bulletin, 110: 972-984. Iverson, R . M . and Denlinger, R.P. 2001. Flow of variably fluidized granular masses across three-dimensional terrain, 1. Coulomb mixture theory. Journal of Geophysical Research, 106(B1): 537-552. Iverson, R .M. , Logan, M . and Denlinger, R.P. 2004. Granular avalanches across irregular three-dimensional terrain: 2. Experimental tests. Journal of Geophysical Research, 109: F01015. Jaboyedoff, M . , Couture, R. and Locat, P. 2004. Structural analysis of Turtle Mountain (Alberta) using digital elevation model (DEM). In Proceedings of the 57 t h Canadian Geotechnical Conference, Quebec City. Section 1C, pp. 7-13. Jakob, M . 2005. Debris-flow hazard analysis. In Debris-flow Hazards and Related Phenomena. Edited by M . Jakob and O. Hungr. Praxis-Springer Publishers, Heidelberg, pp. 411-443. Jakob, M . , Anderson, D., Fuller, T., Hungr, O. and Ayotte, D. 2000. An unusually large debris flow at Hummingbird Creek, Mara Lake, British Columbia. Canadian Geotechnical Journal, 37: 1109-1125. Jaky, J. 1944. The coefficient of earth pressure at rest. Journal of the Society of Hungarian Architects and Engineers, October, pp. 355-358. Jeyapalan, J.K. 1981. Analyses of flow failures of mine tailings impoundments. Ph.D. thesis, University of California, Berkeley. Jeyapalan, J.K., Duncan, J .M. and Seed, H.B. 1983a. Analyses of flow failures of mine tailings dams. Journal of Geotechnical Engineering, 109(2): 150-171. Jeyapalan, J.K., Duncan, J .M. and Seed, H.B. 1983b. Investigation of flow failures of tailings dams. Journal of Geotechnical Engineering, 109(2): 172-189. Johnson, A . M . 1970. Physical Processes in Geology. Freeman, Cooper and Company, San Francisco. 577 pp. 245 Johnson, N.L. 1996. The legacy and future of CFD at Los Alamos. In Proceedings of the 1996 Canadian CFD Conference, Ottawa. King, J. 1996. Tsing Shan debris flow. Geotechnical Engineering Office, Hong Kong Government, Special Project Report SPR 6/96. Kowloon, Hong Kong. 100 pp. Kobayashi, Y . and Kagawa, T. 1987. The prediction of hazards from debris avalanches and rockfalls with the aid of computer simulations. In Proceedings of the International Symposium on Engineering Geological Environment in Mountainous Areas, Beijing, Vol . 1, pp. 567-572. Koch, T., Greve, R. and Hutter, K. 1994. Unconfined flow of granular avalanches along a partly curved surface, Part 2: Experiments and numerical computations. Proceedings of the Royal Society of London, A , 445: 415-435. Korner, H.J. 1976. Reichweite und Geschwindigkeit von Bergsturzen und FlieBschneelawinen. Rock Mechanics, 8: 225-256. Laigle, D. and Coussot, P. 1997. Numerical modeling of mudflows. Journal of Hydraulic Engineering, 123(7): 617-623. Lang, T.E., Dawson, K . L . and Martinelli Jr., M . 1979. Application of numerical transient fluid dynamics to snow avalanche flow, Part 1: Development of computer program A V A L N C H . Journal of Glaciology, 22(86): 107-115. Lang, T.E. and Martinelli Jr., M . 1979. Application of numerical transient fluid dynamics to snow avalanche flow, Part 2: Avalanche modeling and parameter error evaluation. Journal of Glaciology, 22(86): 117-126. L i , J. and Yuan, J. 1983. The main features of the mudflow in Jiang-Jia Ravine. Zeitschrift fur Geomorphologie, 27: 325-341. L i , T. 1983. A mathematical model for predicting the extent of a major rockfall. Zeitschrift fur Geomorphologie, 27: 473-482. Lied, K. and Bakkeh0i, S. 1980. Empirical calculations of snow-avalanche runout distance based on topographic parameters. Journal of Glaciology, 26(94): 165-177. Liu, K.F. and Mei, C C . 1989. Slow spreading of a sheet of Bingham fluid on an inclined plane. Journal of Fluid Mechanics, 207: 505-529. Lucy, L.B. 1977. A numerical approach to the testing of the fission hypothesis. Astrophysical Journal, 82(12): 1013-1024. 246 MacArthur, R.C. and Schamber, D.R. 1986. Numerical methods for simulating mudflows. In Proceedings of the 3 r d International Symposium on River Sedimentation, Mississippi, pp. 1615-1623. Major, J J . and Iverson, R . M . 1999. Debris-flow deposition: Effects of pore-fluid pressure and friction concentrated at flow margins. Geological Society of America Bulletin, 111:1424-1434. Mangeney, A. , Heirich, P. and Roche, R. 2000. Analytical solution for testing debris avalanche numerical models. Pure and Applied Geophysics, 157: 1081-1096. Mangeney-Castelnau, A. , Vilotte, J.P., Bristeau, O., Perthame, B., Bouchut, F., Simeoni, C. and Yemeni, S. 2003. Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme. Journal of Geophysical Research, 108(B11): 2527. McClung, D .M. 2000. Extreme avalanche runout in space and time. Canadian Geotechnical Journal, 37: 161-170. McClung, D . M . 2001. Superelevation of flowing avalanches around curved channel bends. Journal of Geophysical Research, 106(B8): 16,489-16,498. McClung, D . M . and Mears, A.I. 1991. Extreme value prediction of snow avalanche runout. Cold Regions Science and Technology, 19: 163-175. McClung, D .M. and Mears, A.I. 1995. Dry-flowing avalanche run-up and runout. Journal of Glaciology, 41: 359-372. McConnell, R.G. and Brock, R.W. 1904. Report on the great landslide at Frank, Alberta, 1903. In Part VIII, Annual Report 1903, Department of the Interior, Dominion of Canada. McEwen, A.S. and Malin, M.C. 1989. Dynamics of Mount St. Helens' 1980 pyroclastic flows, rockslide-avalanche, lahars, and blast. Journal of Volcanology and Geothermal Research, 37: 205-231. McLellan, P.J. and Kaiser, P.K. 1984. Application of a two-parameter model to rock avalanches of the Mackenzie Mountains. In Proceedings of the 4 l International Symposium on Landslides, Toronto. Vol . 1, pp. 135-140. Mears, A.I. 1981. Design criteria for avalanche control structures in the runout zone. United States Department of Agriculture, Forest Service, Rocky Mountain Forest and Range Experimental Station, General Technical Report RM-84, 28 pp. Monaghan, J.J. 1989. On the problem of penetration in particle methods. Journal of Computational Physics, 82: 1-15. 247 Monaghan, J.J. 1992. Smoothed Particle Hydrodynamics. Annual Reviews in Astronomy and Astrophysics, 30: 543-574. Morgenstern, N.R. and Sangrey 1978. Methods of slope stability analysis. In Landslides, Analysis and Control. Transportation Research Board Special Report 176, United States National Academy of Sciences. Edited by R.J. Schuster and R.J. Krizek. National Academy Press, pp. 155-171. Naaim, M . , Faug, T. and Naaim-Bouvet, F. 2003. Dry granular flow modelling including erosion and deposition. Surveys in Geophysics, 24: 569-585. Nagasawa, M . and Kuwahara, K. 1993. Smoothed particle simulations of the pyroclastic flow. International Journal of Modern Physics, B, 7(9/10): 1979-1995. Nicoletti, P.G. and Sorriso-Valvo, M . 1991. Geomorphic controls of the shape and mobility of rock avalanches. Geological Society of America Bulletin, 103: 1365-1373. O'Brien, J.S. 2003. Reasonable assumptions in routing a dam break mudflow. In Proceedings of the 3 r d International Conference on Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment, Davos. Edited by D. Rickenmann and C L . Chen. Millpress, Rotterdam, pp. 683-693. O'Brien, J.S., Julien, P.Y. and Fullerton, W.T. 1993. Two-dimensional water flood and mudflow simulation. Journal of Hydraulic Engineering, 119(2): 244-261. Orwin, J.F., Clague, J.J. and Gerath, R.F. 2004. The Cheam rock avalanche, Fraser Valley, British Columbia, Canada. Landslides, 1: 289-298. Pariseau, W.G. 1980. A simple mechanical model for rockslides and avalanches. Engineering Geology, 16: 111-123. Patra, A . K . , Bauer, A . C , Nichita, C C , Pitman, E.B., Sheridan, M.F. , Bursik, M . , Rupp, B., Webber, A. , Stinton, A.J . , Namikawa, L . M . and Renschler, C S . 2005. Parallel adaptive numerical simulation of dry avalanches over natural terrain. Journal of Volcanology and Geothermal Research, 139: 1-21. Perla, R., Cheng, T.T. and McClung, D . M . 1980. A two-parameter model of snow-avalanche motion. Journal of Glaciology, 26(94): 197-207. Petley, D.N., Dunning, S.A. and Rosser, N.J. 2005. The analysis of global landslide risk through the creation of a database of worldwide landslide fatalities. In Proceedings of the International Conference on Landslide Risk Management, Vancouver. Edited by O. Hungr, R. Fell, R. Couture and E. Eberhardt. A . A . Balkema, London, pp. 367-373. 248 Petrakov, D.A., Tutubalina, O.V. and Chemomorets, S.S. 2004a. Changes of glacial features and terrain after the 2002 Genaldon Disaster: assessment and forecast. In Proceedings of the Conference on High Mountain Hazard Prevention, Vladikavkaz and Moscow. Petrakov, D.A., Tutubalina, O.V. and Chemomorets, S.S. 2004b. The 2002 Genaldon glacial catastrophe: one year later. Kriosfera Zemli, VIII(l): 29-39. (in Russian) Pierson, T.C. 1986. Flow behaviour of channelized debris flows, Mount St. Helens, Washington. In Hillslope Processes. Edited by A .D. Abrahams. Allen and Unwin, Boston, pp. 269-296. Pirulli, M . 2005. Numerical modelling of landslide runout: A continuum mechanics approach. Ph.D. thesis, Politechnico di Torino, Italy. 204 pp. Pirulli, M . , Preh, A. , Roth, W., Scavia, C. and Poisel, R. 2003. Rock avalance run out prediction: combined application of two numerical methods. In Proceedings of the International Symposium on Rock Mechanics, South African Institute of Mining and Metallurgy, Johannesburg, pp. 903-908. Pirulli, M . , Scavia, C. and Hungr, O. 2004. Determination of rock avalanche run-out parameters through back analyses. In Proceedings of the 9 t h International Symposium on Landslides, Rio de Janeiro. Edited by W A . Lacerda, M . Ehrlich, S.A.B. Fontoura and A.S.F. Sayao. A . A . Balkema, London, pp. 1361-1366. Pitman, E.B, Nichita, C C , Patra, A. , Bauer, A. , Sheridan, M . and Bursik, M . 2003a. Computing granular avalanches and landslides. Physics of Fluids, 15(12): 3638-3646. Pitman, E.B., Nichita, C C , Patra, A . K . , Bauer, A . C , Bursik, M . and Webb, A . 2003b. A model of granular flows over an erodible surface. Discrete and Continuous Dynamical Systems, Series B, 3(4): 589-599. Pouliquen, O. 1999. Scaling laws in granular flows down rough inclined planes. Physics of Fluids, 11(3): 542-548. Press, W.H., Teukolshy, S.A., Vetterling, W.T. and Flannery, B.P. 2002. Numerical recipes in C++, the art of scientific computing. Cambridge University Press, Cambridge. Pudasaini, S.P. and Hutter, K. 2003. Rapid shear flows of dry granular masses down curved and twisted channels. Journal of Fluid Mechanics, 495: 193-208. Pudasaini, S.P., Wang, Y . and Hutter, K. 2005. Rapid motions of free-surface avalanches down curved and twisted channels and their numerical simulation. Philosophical Transactions of the Royal Society of London, A , 3 63: 1551-1571. Rankine, W.J.M. 1857. On the stability of loose earth. Philosophical Transactions of the Royal Society of London, 147. 249 Read, R.S., Langenburg, W., Cruden, D .M. , Field, M . , Stewart, R., Bland, FL, Chen, Z., Froese, C.R., Cavers, D.S., Bidwell, A . K . , Murray, C , Anderson, W.S., Jones, A. , Chen, J., Mclntyre, D., Kenway, D., Bingham, D.K., Weir-Jones, I., Seraphim, J., Freeman, J., Spratt, D., Lamb, M . , Herd, E., Martin, D., McLellan, P. and Pana, D. 2005. Frank Slide a century later: the Turtle Mountain monitoring project. In Proceedings of the International Conference on Landslide Risk Management, Vancouver. Edited by O. Hungr, R. Fell, R. Couture and E. Eberhardt. A . A . Balkema, London, pp. 713-723. Revellino, P., Hungr, O., Guadagno, F .M. , Evans, S.G. 2004. Velocity and runout simulation of destructive debris flows and debris avalanches in pyroclastic deposits, Campania region, Italy. Environmental Geology, 45: 295-311. Rickenmann, D. 1999. Empirical relationships for debris flows. Natural Hazards 19(1): 47-77. Rickenmann, D. and Koch, T. 1997. Comparison of debris flow modelling approaches. In Proceedings of the 1s t International Conference on Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment. Edited by C L . Chen. ASCE, New York. pp. 576-585. Rickenmann, D., Weber, D. and Stepanov, B. 2003. Erosion by debris flows in field and laboratory experiments. In Proceedings of the 3 r d International Conference on Debris-flow Hazards Mitigation: Mechanics, Prediction and Assessment, Davos. Edited by D. Rickenmann and C L . Chen. Millpress, Rotterdam, pp. 883-893. Sailer, R., Rammer, L. and Sampl, P. 2002. Recalculation of an artificially released avalanche with SAMOS and validation with measurements from a pulsed Doppler radar. Natural Hazards and Earth System Sciences, 2: 211-216. Salm, B. 1966. Contribution to avalanche dynamics. In Proceedings of the Symposium on Scientific Aspects of Snow and Ice Avalanches, Davos. International Association of Scientific Hydrology Publication 69. pp. 19-29. Salm, B., Burkard, A . and Gubler, H . 1990. Berechnung von Fliesslawinen eine Anleitung fur Praktiker mit Beispielen. Mitteilungen des Eidgenossischen Institut fur Schnee- und Lawinenforschung, No. 47, 37 pp. Sampl, P. 1993. Current status of the A V L avalanche simulation model - Numerical simulation of dry snow avalanches. In Proceedings of the Pierre Beghin International Workshop on Rapid Gravitational Mass Movements, Grenoble. Edited by L. Buisson and G. Brugnot. pp. 269-276. th Sassa, K. 1985. The mechanism of debris flows. In Proceedings of the 11 International Conference on Soil Mechanics and Foundation Engineering, San Francisco. Vol . 1, pp. 1173-1176. 250 Sassa, K . 1988. Geotechnical model for the motion of landslides. In Proceedings of the 5 International Symposium on Landslides, Lausanne. Edited by C. Bonnard. A . A . Balkema, Rotterdam, pp. 37-56. Sassa, K . 2002. Mechanism of rapid and long travelling flow phenomena in granular soils. In Proceedings of the UNESCO/IGCP International Symposium on Landslide Mitigation and Protection of Cultural and Natural Heritage, Kyoto. Edited by K . Sassa. Disaster Prevention Research Institute, Kyoto University, pp. 11-30. Savage S.B. and Hutter, K. 1989. The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics, 199: 177-215. Savage S.B. and Hutter, K . 1991. The dynamics of avalanches of granular materials from initiation to runout, Part 1: Analysis. Acta Mechanica, 86: 201-223. Schamber, D.R. and MacArthur, R.C. 1985. One-dimensional model for mud flows. In Proceedings of the A S C E Specialty Conference on Hydraulics and Hydrology in the Small Computer Age. American Society of Civil Engineers, New York. Vol . 2, pp. 1334-1339. Scheidegger, A .E . 1973. On the prediction of the reach and velocity of catastrophic landslides. Rock Mechanics, 5: 231-236. Schwab, J.W., Geertsema, M . and Evans, S.G. 2003. Catastrophic rock avalanches, west-central British Columbia. In Proceedings of the 3 r d Canadian Conference on Geo technique and Natural Hazards, Edmonton, pp. 291-299. Selby, M.J. 1985. Earth's changing surface. Oxford University Press, Oxford. Sherard, J.L., Woodward, R.J., Gizienski, S.F. and Clevenger, W A . 1963. Earth and earth-rock dams. John Wiley and Sons, New York, 722 pp. Sheridan, M.F. and Patra, A . K . 2005. Modeling and simulation of geophysical mass flows. Journal of Volcanology and Geothermal Research, 139: vii. Sheridan, M.F. , Stinton, A.J . , Patra, A. , Pitman, E.B., Bauer, A . and Nichita, C C . 2005. Evaluating Titan2D mass-flow model using the 1963 Little Tahoma Peak avalanches, Mount Rainier, Washington. Journal of Volcanology and Geothermal Research, 139: 89-102. Sousa, J. and Voight, B. 1991. Continuum simulation of flow failures. Geotechnique, 41(4): 515-538. Sovilla, B. and Bartelt, P. 2002. Observations and modelling of snow avalanche entrainment. Natural Hazards and Earth System Sciences, 2: 169-179. Stadler, R. 1986. Stationares, schnelles FlieBen von dicht gepackten trockenen und feuchten Schiittgutern. Dr. Ing. Dissertation, University of Karlsruhe. 251 Stoker, J J . 1957. Water Waves. Interscience, New York. Tai, Y . C . and Gray, J.M.N.T. 1998. Limiting stress states in granular avalanches. Annals of Glaciology, 26: 272-276. Takahashi, T. 1978. Mechanical characteristics of debris flow. Journal of the Hydraulics Division of the American Society of Civi l Engineers, 104(HY8): 1153-1169. Takahashi, T. 1991. Debris flow. International Association for Hydraulic Research monograph. A . A . Balkema, Rotterdam. Takahashi, T. and Yoshida, H. 1979. Study on the deposition of debris flows Part 1: Deposition due to abrupt change of bed slope. Annals of the Disaster Prevention Research Institute, Kyoto University, 22B-2. (in Japanese) Takahashi, T., Nakagawa, H., Harada, T. and Yamashiki, Y . 1992. Routing debris flows with particle segregation. Journal of Hydraulic Engineering, 118(11): 1490-1507. Terzaghi, K . 1920. Old earth pressure theories and new test results. Engineering News Record, 85: 632. Terzaghi, K . and Peck, R.B. 1967. Soil mechanics in engineering practice. John Wiley and Sons, New York. Thurber Engineering Ltd. and Golder Associates Ltd. 1993. Cheekye River terrain hazard study. Unpublished report to the British Columbia Ministry of Environment, Lands and Parks, Victoria. Trunk, F.J., Dent, J.D. and Lang, T.E. 1986. Computer modeling of large rock slides. Journal of Geotechnical Engineering, 112: 348-360. Turnbull, B. and Bartelt, P. 2003. Mass and momentum balance model of a mixed flowing/powder snow avalanche. Surveys in Geophysics, 24: 465-477. Tutubalina, O.V., Chernomorets, S.S. and Petrakov, D.A. 2005. Kolka Glacier before the 2002 collapse: new data. Kriosfera Zemli, IX(4): 62-71. (in Russian) VanDine, D.F. 1985. Debris flows and debris torrents in the southern Canadian Cordillera. Canadian Geotechnical Journal, 22: 44-68. Voellmy, A . 1955. Uber die Zerstorungskraft von Lawinen. Schweizerische Bauzeitung, 73: 212-285. Voight, B. and Pariseau, W.G. 1978. Rockslides and Avalanches: A n Introduction. In Rockslides and Avalanches. Edited by B. Voight. Elsevier, New York. Vol . 1, pp. 1-67. 252 Voight, B. and Sousa, J. 1994. Lessons from Ontake-san: a comparative analysis of debris avalanche dynamics. Engineering Geology, 38: 261-297. Voight, B., Komorowski, J.C., Norton, G.E., Belousov, A .B . , Belousova, M . , Boudon, G., Francis, P.W., Franz, W., Heinrich, P., Sparks, R.S.J, and Young, S.R. 2002. The 26 December (Boxing Day) 1997 sector collapse and debris avalanche at Soufriere Hills Volcano, Montserrat. In The Eruption of Soufriere Hills Volcano, Montserrat from 1995 to 1999. Edited by T.H. Druitt and B.P. Kokelaar. Geological Society of London Memoir, 21: 363-407. Wang, Z. and Shen, H.T. 1999. Lagrangian simulation of one-dimensional dam-break flow. Journal of Hydraulic Engineering, 125(11): 1217-1220. Wang, Y . , Hutter, K. and Pudasaini, S.P. 2004. The Savage-Hutter theory: A system of partial differential equations for avalanche flows of snow, debris, and mud. Zeitschrift fur Angewandte Mathematik und Mechanik, 84(8): 507-527. Wieland, M . , Gray, J.M.N.T. and Hutter, K. 1999. Channelized free-surface flow of cohesionless granular avalanches in a chute with shallow lateral curvature. Journal of Fluid Mechanics, 392: 73-100. Wise, M.P., Moore, G. and VanDine, D. 2004. Land Management Handbook 56: Landslide Risk Case Studies in Forest Development and Planning Operations. British Columbia Ministry of Forests, Research Branch, Victoria. 124 pp. 253 "@en . "Thesis/Dissertation"@en . "10.14288/1.0052928"@en . "eng"@en . "Geological Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en . "Graduate"@en . "A new continuum dynamic model for the analysis of extremely rapid landslide motion across complex 3D terrain"@en . "Text"@en . "http://hdl.handle.net/2429/18500"@en .