Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Mitchell, Andrew David
Master of Applied Science - MASc
University of British Columbia
Rock falls are a common hazard affecting people and infrastructure below natural cliffs and in man-made environments such as roads, railways, quarries and open pit mines. Computer models have become standard in assessing the hazards posed by rock falls, and there are a wide variety of models currently available. The main difference between these models is the way in which they represent the impact process. New models utilizing rigorous solutions for rigid body impacts have become more common, however, the complexity of these models makes them difficult to calibrate and difficult to apply in practice. The PIERRE 2 model, presented here, returns to a simplified lumped-mass model assuming collinear impact conditions. The natural variability of rock falls is represented using the following features: • The slope angle is varied by a stochastic roughness angle to approximate the effects of non-collinear impacts. • The conservation of momentum during an impact is determined using hyperbolic relationships so that high energy impacts will have greater losses than low energy impacts. • A stochastic shape factor is used to vary the sphere dimensions at impact to approximate the rotational effects of real, non-spherical boulders. Impact mechanics theory was used as the basis for these features; however their validity is demonstrated through 2D and 3D model calibration. The model was calibrated using data from experiments with detailed kinematic information on two natural slopes, Vaujany, France and Ehime, Japan, and excavated slopes at three Austrian quarries. The results from the calibrations were then applied to a natural slope, Tornado Mountain, and an excavated slope at Nicolum quarry, both in British Columbia. The results of this calibration and validation program showed the model can reproduce runout and kinematic behaviours of rocks based on relatively limited site characterization data. Select results were also compared to design guidelines. Using reliability engineering principles, the probability of the design parameters derived from the model being exceeded was found to be acceptably low. The model is also capable of reproducing observed relationships between impact incidence angle and the normal restitution coefficient reported by others.
COMPUTER MODELING OF ROCK FALL: IMPROVEMENT AND CALIBRATION OF THE PIERRE 2 STOCHASTIC ROCK FALL SIMULATION PROGRAM by Andrew David Mitchell B.Sc., University of Alberta, 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Geological Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) May 2015 © Andrew David Mitchell, 2015 ii Abstract Rock falls are a common hazard affecting people and infrastructure below natural cliffs and in man-made environments such as roads, railways, quarries and open pit mines. Computer models have become standard in assessing the hazards posed by rock falls, and there are a wide variety of models currently available. The main difference between these models is the way in which they represent the impact process. New models utilizing rigorous solutions for rigid body impacts have become more common, however, the complexity of these models makes them difficult to calibrate and difficult to apply in practice. The PIERRE 2 model, presented here, returns to a simplified lumped-mass model assuming collinear impact conditions. The natural variability of rock falls is represented using the following features:  The slope angle is varied by a stochastic roughness angle to approximate the effects of non-collinear impacts.  The conservation of momentum during an impact is determined using hyperbolic relationships so that high energy impacts will have greater losses than low energy impacts.  A stochastic shape factor is used to vary the sphere dimensions at impact to approximate the rotational effects of real, non-spherical boulders. Impact mechanics theory was used as the basis for these features; however their validity is demonstrated through 2D and 3D model calibration. The model was calibrated using data from experiments with detailed kinematic information on two natural slopes, Vaujany, France and Ehime, Japan, and excavated slopes at three Austrian quarries. The results from the calibrations were then applied to a natural slope, iii Tornado Mountain, and an excavated slope at Nicolum quarry, both in British Columbia. The results of this calibration and validation program showed the model can reproduce runout and kinematic behaviours of rocks based on relatively limited site characterization data. Select results were also compared to design guidelines. Using reliability engineering principles, the probability of the design parameters derived from the model being exceeded was found to be acceptably low. The model is also capable of reproducing observed relationships between impact incidence angle and the normal restitution coefficient reported by others. iv Preface The work presented in this thesis is the result of a computer model calibration and validation. Previously published field experimental results from Canada, France, and Japan as well as previously unpublished experiments from Canada and Austria were used in the calibration process. The primary goal of this work is to demonstrate the applicability of the PIERRE 2 rockfall dynamics model as a viable engineering tool. The secondary goal is to present model inputs for first order predictions of rockfall dynamics from limited slope characterization data. During much of the work presented in this thesis, I worked as a member of a team of researchers. The paragraphs below describe my personal contribution to the variety of tasks. Chapters 3 and 4 represent an unpublished extended draft of a journal manuscript (at the time of submission of this thesis), written by the author. I have carried out the calibration of the model and improved the model algorithm by adding the normal roughness angle distribution, roughness reduction function and shape factors. The figures presented are all original to the author, and drafted for the manuscript. The model in its original form was developed by O. Hungr, who also assisted in editing the manuscript. A version of Chapter 5 has been published. Gischig, V., Hungr, O., Mitchell, A., Bourrier, F. 2015. Pierre3D – a 3D stochastic rock fall simulator based on random ground roughness and hyperbolic restitution coefficients. Canadian Geotechnical Journal, in press. I assisted in supplying a preliminary 2D calibration to set the ranges of inputs for the Monte Carlo simulations. I also assisted in determining the best objective function for comparing the results of the simulations. I also assisted in revising the text of the manuscript. Implementing the 3D version of the model in MATLAB was carried out by V. Gischig, as well as completing the v various calibrations in the 3D model and drafting the manuscript. O. Hungr supplied the conceptual model for the 3D version based on the original 2D model, as well as editing the manuscript. F. Bourrier assisted with editing the text and supplied insights into the testing conducted in 2003. A version of Chapter 6 has been submitted for publication. Preh, A., Mitchell, A., Hungr, O., Kolenprat, B. 2015, in review. I completed the initial parametric study referenced in the paper. I also completed the study of the optimal stopping criteria. I did all the work on applying the model and the optimal parameters to the Nicolum quarry simulation and drafted the text for that section of the paper. Figures 7, 13, 15 and 16 are original to the author. A. Preh performed the final calibrations for the Austrian quarry sites, drafted the main text, and produced the remaining figures. O. Hungr supplied the original version of the model, as well as assisted in editing the text. B. Kolenprat assisted in conducting the field testing. Chapter 7 consists of original research conducted by the author. Apart from work on the model, I also visited the IRSTEA rock fall testing site at Vaujany, France, during a 3-day period of testing in September, 2014. I assisted the IRSTEA team in conducting the experiments. I also helped Dr. D. Wyllie to carry out additional field tests at the Nicolum Quarry site in the winter of 2015. While data from those experiments are not presented here, they have been invaluable for giving me additional insights into the complexities of obtaining information on falling rocks in the field. vi Table of Contents Abstract .......................................................................................................................................... ii Preface ........................................................................................................................................... iv Table of Contents ......................................................................................................................... vi List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................. xiv List of Symbols ......................................................................................................................... xxiv Acknowledgements .................................................................................................................. xxvi Chapter 1: Introduction ................................................................................................................1 1.1 Background ....................................................................................................................... 1 1.2 Statement of Purpose and Objectives ............................................................................... 5 1.3 Organization of the Thesis ................................................................................................ 6 Chapter 2: Literature review ........................................................................................................7 2.1 Impact Mechanics ............................................................................................................. 7 2.1.1 Stereomechanical impacts ....................................................................................... 7 2.1.2 Contact laws .......................................................................................................... 11 2.1.3 Impacts on granular soil ........................................................................................ 13 2.1.4 Fragmentation ........................................................................................................ 15 2.2 Coefficients of Restitution .............................................................................................. 16 2.3 Computer Modeling for Rock Fall ................................................................................. 23 2.3.1 Spatial representation ............................................................................................ 23 2.3.2 Lumped mass models ............................................................................................ 24 2.3.3 Rigid body models ................................................................................................ 25 vii 2.3.4 Stochastic models .................................................................................................. 25 2.3.5 Colorado Rockfall Simulation Program ................................................................ 26 2.3.6 Rocfall ................................................................................................................... 28 Chapter 3: PIERRE 2 model description ..................................................................................30 3.1 Slope ............................................................................................................................... 30 3.2 Trajectory ........................................................................................................................ 31 3.2.1 Flight ..................................................................................................................... 31 3.2.2 Sliding ................................................................................................................... 32 3.3 Collision .......................................................................................................................... 33 3.4 Slope Properties .............................................................................................................. 34 3.4.1 Roughness ............................................................................................................. 34 3.4.2 Interface friction angle .......................................................................................... 39 3.4.3 Reference normal deformation energy .................................................................. 39 3.4.4 Reference tangential deformation energy ............................................................. 43 3.5 Boulders .......................................................................................................................... 44 3.6 Impact Model .................................................................................................................. 47 3.7 Stopping Criterion .......................................................................................................... 49 3.8 Apparent Restitution Values ........................................................................................... 51 3.9 Program Outputs ............................................................................................................. 53 Chapter 4: Calibration and validation of the PIERRE 2 stochastic rock fall dynamics simulation program .....................................................................................................................54 4.1 Introduction ..................................................................................................................... 54 4.2 Data ................................................................................................................................. 55 viii 4.2.1 Calibration ............................................................................................................. 55 4.2.2 Validation .............................................................................................................. 62 4.3 Methodology ................................................................................................................... 66 4.3.1 Parametric studies ................................................................................................. 67 4.3.2 Velocity versus elevation data ............................................................................... 70 4.4 Initial Model Calibration ................................................................................................ 72 4.4.1 Vaujany – open talus slope .................................................................................... 72 4.4.2 Preg – strong rock quarry ...................................................................................... 81 4.5 Detailed Model Calibration ............................................................................................ 92 4.5.1 Ehime – open slope with weak bedrock and talus ................................................ 92 4.5.2 Vaujany ................................................................................................................. 99 4.5.3 Preg quarry .......................................................................................................... 102 4.6 Model Validation .......................................................................................................... 103 4.6.1 Tornado Mountain – firm talus slope .................................................................. 103 4.6.2 Nicolum Quarry – exposed hard rock ................................................................. 107 4.7 Use of Model Result as Design Inputs ......................................................................... 109 4.8 Discussion and Conclusions ......................................................................................... 112 Chapter 5: PIERRE3D - a 3D stochastic rock fall simulator based on random ground roughness and hyperbolic restitution coefficients ...................................................................114 5.1 Introduction ................................................................................................................... 114 5.2 Method .......................................................................................................................... 123 5.2.1 Flight trajectory ................................................................................................... 123 5.2.2 Collision .............................................................................................................. 123 ix 5.2.3 Applying stochastic roughness ............................................................................ 125 5.2.4 Goldsmith model ................................................................................................. 127 5.2.4.1 Hyperbolic restitution factors ................................................................ 128 5.2.4.2 Evaluation of the friction limit angle .................................................... 129 5.2.4.3 Contained slip impacts .......................................................................... 130 5.2.4.4 Uncontained slip impacts ...................................................................... 131 5.2.5 Model implementation ........................................................................................ 132 5.3 Data ............................................................................................................................... 132 5.4 Model Calibration ......................................................................................................... 136 5.5 Results........................................................................................................................... 138 5.5.1 First calibration attempt, single homogenous ground material ........................... 138 5.5.2 Second calibration attempt, three ground materials ............................................ 139 5.6 Discussion ..................................................................................................................... 142 5.6.1 Model and data uncertainties ............................................................................... 142 5.6.2 Probabilistic rock fall hazard maps ..................................................................... 145 5.7 Conclusions................................................................................................................... 148 Chapter 6: Stochastic analysis of rock fall dynamics on quarry slopes ................................150 6.1 Introduction: Rock Fall Hazards in Quarries and Open Pit Mines ............................... 150 6.2 Full Scale Testing at Quarry Sites ................................................................................ 152 6.2.1 Test procedure ..................................................................................................... 152 6.2.2 Summary of field test results ............................................................................... 154 6.3 PIERRE 2D and 3D, Lumped-mass Stochastic Model of Rock Fall Ballistics ............ 156 6.4 Calibration of the Model ............................................................................................... 161 x 6.4.1 General ................................................................................................................ 161 6.4.2 Calibration results, Austrian sites ........................................................................ 163 6.4.2.1 Influence of the degree of detail of the modeled slope geometries on the calculated hazard distances ................................................................... 171 6.4.3 Calibration results, Canadian quarry ................................................................... 172 6.5 Discussion and Conclusion ........................................................................................... 176 Chapter 7: Integrating size and shape effects in three dimensions .......................................179 7.1 Model inputs ................................................................................................................. 179 7.2 Calibration .................................................................................................................... 180 7.2.1 Vaujany ............................................................................................................... 180 7.2.2 Nicolum ............................................................................................................... 188 7.3 Sensitivity to Shape Factor ........................................................................................... 190 7.4 Discussion ..................................................................................................................... 193 Chapter 8: Conclusion ...............................................................................................................197 8.1 Applicability of the PIERRE 2 Model .......................................................................... 198 8.2 Model Limitations ........................................................................................................ 201 8.3 Future Research ............................................................................................................ 202 Bibliography ...............................................................................................................................205 Appendices ..................................................................................................................................213 Appendix A Roughness Angle Distribution Calibration ....................................................... 214 A.1 Minimum Number of Simulations to Reduce Variability ................................... 214 A.2 Input Mean and Standard Deviation for Truncated Normal Distribution ........... 217 A.3 Size Dependent Roughness Reduction ................................................................ 220 xi Appendix B Stopping Criterion ............................................................................................. 223 Appendix C Ehime Additional Calibration Plots ................................................................... 229 Appendix D Statistical Functions for Reliability Calculations .............................................. 235 xii List of Tables Table 1.1: Factors contributing to an area being a rock fall source (adapted from Volkwein et al., 2011). .............................................................................................................................................. 3 Table 2.1: Summary of normal and tangential coefficients of restitution (en and et, respectively), related to slope material type (adapted from Heidenreich, 2004). ................................................ 18 Table 4.1: Calculation table to determine the SEL and MEL for Screen 1 and Screen 2 locations..................................................................................................................................................... 111 Table 4.2: Probability of failure from exceeding the net's impact capacity or interception height...................................................................................................................................................... 112 Table 4.3: Summary of best fit input parameters from calibration process. ............................... 112 Table 5.1: Best fit parameter sets used to reproduce observed rock fall data from the Vaujany test site. .............................................................................................................................................. 141 Table 6.1: Input parameters required by the model and their values, determined by calibration, used in all the analyses ................................................................................................................ 162 Table 6.2: Comparison of calculated and observed Runout Distance statistics, obtained with the optimal set of input parameters, shown in Table 1. (Note: in the last line, the distance RD is divided by the slope height, H, see Figure 6.6) .......................................................................... 167 Table 6.3: Comparison of calculated and observed First Impact distance statistics, obtained with the optimal set of input parameters, shown in Table 1. (Note: in the last line, the distance ID is divided by the slope height, H, see Figure 6.6) .......................................................................... 168 Table 6.4: Comparison of runout distance statistics using a fine and coarse (smoothed) slope geometry ..................................................................................................................................... 172 Table 7.1: Model input parameters from Vaujany 3D recalibration. .......................................... 181 xiii Table 7.2: Comparison of the design velocity values from the 3D and 2D simulations. ........... 188 Table 8.1: Best-fit model parameters for both the 2D and 3D models. ...................................... 201 xiv List of Figures Figure 2.1: Ratio of rebound over incident kinetic energy as a function of the interface friction minus the friction limit angle. ....................................................................................................... 11 Figure 2.2: Graphical comparison of the relations described by Equation 2.6 and Equation 2.7. 22 Figure 3.1: Modification of the local slope angle by the roughness angle. .................................. 35 Figure 3.2: Comparison of roughness angle distributions using a uniform distribution for the longitudinal and tangential roughness, a and b, and using a normal distribution, c (truncated) and d..................................................................................................................................................... 38 Figure 3.3: Conceptual diagram showing the force-displacement diagrams during a) the start of the incident phase, b) the point of maximum force (Fmax) and penetration (αmax), and c) the end of the restitution phase. The shaded area represents the total plastic work for the impact. .............. 41 Figure 3.4: Hyperbolic relationship between deformation energy and normal restitution, with the scaling point highlighted. .............................................................................................................. 43 Figure 3.5: Hypothetical impact of an irregularly shaped natural rock, a) contrasting equivalent spherical radius to the actual distance to the contact point, and b) components contributing to the moment about the contact point. ................................................................................................... 45 Figure 3.6: Representation of non-collinear impact using an effective radius. ............................ 46 Figure 3.7: Relationship between the incident and rebound kinetic energy as a function of vnin (m/s) and log mass (mass in kg). .................................................................................................. 49 Figure 4.1: Distribution of boulder sizes used in the 2003 Vaujany experiments a) on the basis of mass and b) diameter. ................................................................................................................... 56 Figure 4.2: Vaujany test site slope profile and the locations of Screen 1, Screen 2 and the forest road. .............................................................................................................................................. 57 xv Figure 4.3: Vaujany test site, looking downslope from the release location. ............................... 57 Figure 4.4: Preg Quarry test site profiles, a) Section 1, b) Section 2, c) Section 3, and d) Section 4..................................................................................................................................................... 58 Figure 4.5: Preg Quarry test site (photo A. Preh). ........................................................................ 59 Figure 4.6: Distribution of boulder sizes used in the Preg experiments a) Section 1 mass, b) Section 1 diameter, c) Section 2 mass, d) Section 2 diameter, 3) Section 3 mass, f) Section 3 diameter, g) Section 4 mass, h) Section 4 diameter. ..................................................................... 60 Figure 4.7: Ehime test site profile. ................................................................................................ 61 Figure 4.8: Distribution of natural rock sizes used in the Ehime tests a) mass, and b) diameter. 62 Figure 4.9: Location of tree impacts along the profile used to estimate jump heights and velocities for boulders A and B from Tornado Mountain. ............................................................ 63 Figure 4.10: Nicolum quarry test site profile. ............................................................................... 64 Figure 4.11: Distribution of boulder sizes used in the Nicolum tests on the basis of a) mass, and b) diameter. ................................................................................................................................... 65 Figure 4.12: Nicolum quarry test site with experimental fence located at the toe of the slope (Photo O. Hungr). ......................................................................................................................... 66 Figure 4.13: Initial parametric surface plots for runout distance. a) RMSE of the runout distribution and b) the absolute value of the residual of the maximum runout. ........................... 68 Figure 4.14: Initial parametric surface plots for runout distance. a) KS test statistic for the runout distribution and b) p-value for the KS test. ................................................................................... 69 Figure 4.15: Parametric surface plots for the runout distribution RMSE, a) for the initial wide range of inputs, and b) for the refined input parameters. .............................................................. 73 xvi Figure 4.16: Parametric surface plots for runout distance. a) RMSE of the runout distribution and b) the absolute value of the residual of the maximum runout. ...................................................... 74 Figure 4.17: Parametric surface plots of Screen 1 jump height. a) RMSE of the jump height distribution and b) residual of the maximum jump height. ........................................................... 74 Figure 4.18: Parametric surface plots of Screen 1 velocity. a) RMSE of the velocity distribution and b) residual of the maximum velocity. .................................................................................... 75 Figure 4.19: Parametric surface plots of Screen 2 jump height. a) RMSE of the jump height distribution and b) residual of the maximum jump height. ........................................................... 75 Figure 4.20: Parametric surface plots of Screen 2 velocity. a) RMSE of the velocity distribution and b) residual of the maximum velocity. .................................................................................... 76 Figure 4.21: Comparison of observed and simulated runout distances from Vaujany experiments using the truncated normally distributed roughness angle. ........................................................... 79 Figure 4.22: Screen 1 Jump Height and Velocity Distributions using μ = σ = 0.5 and θscale = 0.72........................................................................................................................................................ 80 Figure 4.23: Screen 2 Jump Height and Velocity Distributions using μ = σ = 0.5 and θscale = 0.72........................................................................................................................................................ 80 Figure 4.24: Preg quarry parametric surface plots of RMSE for varying quarry face θscale, a) first impact distribution for Section 1, and b) first impact distribution for Section 2. ......................... 82 Figure 4.25: Preg quarry parametric surface plots for varying quarry face θscale, a) residual maximum first impact location for Section 1, and b) residual maximum first impact location for Section 2........................................................................................................................................ 82 Figure 4.26: Preg quarry parametric surface plots for quarry face θscale = 0.5, a) runout distribution for Section 1, and b) runout distribution for Section 2. ............................................. 83 xvii Figure 4.27: Preg quarry parametric surface plots for quarry face θscale = 0.7, a) runout distribution for Section 1, and b) runout distribution for Section 2. ............................................. 83 Figure 4.28: Cumulative distributions comparing the observed distributions to the simulations using the optimum parameters from the parametric study for a) Section 1 runout, b) Section 2 runout, c) Section 1 first impact, and d) Section 2 first impact. ................................................... 85 Figure 4.29: Results of optimum parameters identified in the parametric study applied to a) Section 3 runout, b) Section 4 runout, c) Section 3 first impact, and d) Section 4 first impact. .. 86 Figure 4.30: Comparisons of the observed runout versus boulder mass for Preg sections 1 to 3. 88 Figure 4.31: Comparison of runout versus mass without roughness reduction and with roughness reduction for Preg Section 1, solid line indicating a linear regression fit to the simulated runout data. ............................................................................................................................................... 90 Figure 4.32: Runout versus mass for Preg Section 2, solid line indicating a linear regression fit to the simulated runout data. ............................................................................................................. 90 Figure 4.33: Runout versus mass for Preg Section 3, solid line indicating a linear regression fit to the simulated runout data. ............................................................................................................. 91 Figure 4.34: Runout distance versus duration with Hstop = 0.25R for Preg Section 3. ................. 92 Figure 4.35: a) Plot of observed and simulated incident velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. The solid and dashed lines on b) and c) indicate the mean and ± 2 standard deviations on the residual distributions. .................................................................................................................... 94 xviii Figure 4.36: a) Plot of observed and simulated rebound velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data.95 Figure 4.37: a) Plot of observed and simulated rebound angular velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. ............................................................................................................................... 96 Figure 4.38: Comparison of rebound translational kinetic energy to angular kinetic energy from Ehime cube dataset, dashed lines are from Ushiro et al. (2006). a) Without the inclusion of a shape factor, and b) with a shape factor uniformly distributed between 0.85 and 1.3. ................ 97 Figure 4.39: a) Plot of observed and simulated angular velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data.99 Figure 4.40: Comparison of observed and simulated runout distances from Vaujany experiments using the optimum input parameter set. ...................................................................................... 100 Figure 4.41: Comparison of observed and simulated jump height, translational velocity and kinetic energy at the Vaujany Screen 1 location. ........................................................................ 101 Figure 4.42: Comparison of observed and simulated jump height, translational velocity and kinetic energy at the Vaujany Screen 2 location. ........................................................................ 102 Figure 4.43: Runout distributions for boulders A and B from Tornado Mountain using the input parameters found in the Vaujany and Ehime calibrations, with the observed runout shown by the dashed line. ................................................................................................................................. 103 xix Figure 4.44: Runout distributions for boulders A and B from Tornado Mountain with the roughness scale reduced to 0.6, with the observed runout shown by the dashed line. ............... 104 Figure 4.45: Comparison of simulated and observed jump heights versus the runout distance for Boulders A and B from Tornado Mountain. ............................................................................... 104 Figure 4.46: Comparison of estimated and simulated velocities versus the runout distance for boulder A (left column) and boulder B (right column) from Tornado Mountain. ...................... 105 Figure 4.47: Comparison of Pierre 2 calculated normal restitution factors relative to the impact incidence angle for the Tornado Mountain simulation compared to incidence angle and normal restitution values presented by Wyllie (2014). Also shown are similar data from four sites in Australia collected by Spadari et al. (2011). ............................................................................... 107 Figure 4.48: a) Plot of observed and simulated incident velocities vs elevation for the Nicolum quarry dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data...................................................................................................................................................... 108 Figure 4.49: a) Plot of observed and simulated rebound velocities vs elevation for the Nicolum quarry dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data...................................................................................................................................................... 109 Figure 5.1: A conceptual diagram of contact force, F and displacement during impact (Falcetta, 1985) ........................................................................................................................................... 119 Figure 5.2: Shapes of three relationships between the restitution factor and approach momentum, discussed in the text. Red: Equation (2) (Falcetta, 1985). Green: Equation (3) (Pfeiffer and Bowen, 1989). Black: Equation (4), Hyperbolic (Bourrier and Hungr, 2011). ........................ 120 xx Figure 5.3: Schematic overview of the 3D collision algorithm. ................................................. 124 Figure 5.4: a) Orthophoto of the Vaujany rock fall test site in France. Also shown are the end locations of 71 blocks out of a total of 100 blocks released from the upper forest road. b) Map showing end locations, position of observation screens, definition of the scatter angle as used in the calibration as well as the release point. c) Distribution of block weights. d) Run-out distance distribution. e) Distribution of scatter angle. f) and g) Velocity and jump height distribution observed at screen 1. h) and i) Velocity and jump height distribution observed at screen 2. ..... 134 Figure 5.5: Vaujany test site photos. a) Shallow gully in the mid-reaches of the path, upslope of the observation screens. b) A talus segment outside the gully. c) Talus forming the lowest segment of the path, down-slope from the road. d) Forest road crossing the lower third of the path (Photos O. Hungr). .............................................................................................................. 135 Figure 5.6: Cumulative distribution function (CDF) and probability density functions (PDF) of the observed quantities shown in Figure 5.4d-i (black line). The red curves are the results of the best-fit model using only one parameters set for the entire slope. For the orange curve the parameters of the lower forest road was changed such that fewer blocks stop on it (see Figure 5.7a). The blue curve represents the best-fit model for which a third material was introduced for the talus cone below the forest road (see Figure 5.7b). .............................................................. 137 Figure 5.7: Maps of end positions of the three best-fit model versions shown in Figure 5.5. Also indicated are the observed end positions in magenta. ................................................................. 139 Figure 5.8: Maps showing the spatial distribution of ground material (Figure 5.5). The parameters for the corresponding material numbers are given in Table 5.1. a) Different properties were assigned to the forest road to inhibit blocks from stopping on the lower road. b) A third material was assigned for the slope below the lower forest road. .................................. 140 xxi Figure 5.9: Model results for the case with three ground materials. The variability of the results related to low sampling numbers is derived by randomly drawing 100 realizations from 10 000 model realization. This is repeated 1000 times. The blue shading represents the 95% percentile of the 1000 random CDFs. a) to f) as in Figure 5.6. ....................................................................... 143 Figure 5.10: Comparison of model results based on two different calibration strategies; calibration giving weight to all parameters (as described in Section 4), and giving weight to screen velocities only. a) to f) as in Figure 6. ............................................................................. 144 Figure 5.11: Probabilistic hazard maps. a) Impact rate derived from 100 model realizations (block releases). b) Impact rate derived from 10 000 model realizations (case with three materials). b) Maximum jump height per cell. c) Maximum velocity. d) Maximum kinetic energy. e) Rate of exceeding a velocity of 20m/s. f) Rate of exceeding jump heights of 5 m. For b – f 10 000 model realizations were used. .................................................................................... 147 Figure 6.1: Photo of a rock fall accident at a quarry in Austria (Kolenprat, 2012). ................... 150 Figure 6.2: Danger areas associated with a quarry face (Kolenprat, 2012). ............................... 151 Figure 6.3: Rock fall drop tests in progress at the Klöch Quarry in Austria (Photo A. Preh). ... 152 Figure 6.4: Nicolum quarry test site, British Columbia, Canada (Photo O. Hungr). .................. 154 Figure 6.5: Summary of all drop tests in nine Austrian quarries. Histograms and probability density functions of (a) the first impact distance (412 measurements) and (b) runout distance (589 measurements), measured from the slope toe and normalized by slope height. ................ 155 Figure 6.6: Cross-sections used for the back analyses of Austrian quarry tests. ........................ 164 Figure 6.7: Mass distributions of boulders used in the field tests and in the calculations. a) Klöch, b) Pauliberg, c) Preg 1, d) Preg 2, e) Preg 3, f) Preg 4, and g) Nicolum. ....................... 164 xxii Figure 6.8: Example locus of calculated trajectories from the Pauliberg site (80 Boulders, each simulated 10 times). .................................................................................................................... 165 Figure 6.9: Example linear and rotational kinetic energy lines, Pauliberg site. ......................... 166 Figure 6.10: Comparison of calculated and observed hazard distances, Klöch site (four combined sections, 97 throws observed, 970 throws calculated) ................................................................ 169 Figure 6.11: Comparison of calculated and observed hazard distances, Pauliberg site (80 throws observed, 800 throws calculated) ................................................................................................ 169 Figure 6.12: Comparison of calculated and observed hazard distances, Preg sites (25 throws observed, 250 throws calculated per site). .................................................................................. 170 Figure 6.13: Comparison of simulated and observed 50% of exceedance runout durations. ..... 171 Figure 6.14: Example locus of calculated trajectories from the Nicolum site (34 Boulders). .... 174 Figure 6.15: Comparison of observed and simulated maximum bounce heights at different elevations along the Nicolum quarry profile. ............................................................................. 174 Figure 6.16: Frequency of impacts per square meter obtained from 250 model runs within the 154 m to 102 m elevation range (note the logarithmic scale). Red points correspond to observed impact locations. ......................................................................................................................... 176 Figure 7.1: Vaujany runout distributions using the 3D simulation with the truncated normal distribution and size dependent roughness reduction. ................................................................ 182 Figure 7.2: Jump height and linear velocity comparisons for Vaujany Screen 1. ...................... 183 Figure 7.3: Jump height and linear velocity comparisons for Vaujany Screen 2. ...................... 184 Figure 7.4: Vaujany runout distributions using the 3D simulation a) raw runout distribution, and b) filtered to remove the blocks stopping within the first 10 m of the release point. ................. 185 xxiii Figure 7.5: Jump height and linear velocity comparisons for Vaujany Screen 1 using the shape factor. .......................................................................................................................................... 186 Figure 7.6: Jump height and linear velocity comparisons for Vaujany Screen 2 using the shape factor. .......................................................................................................................................... 187 Figure 7.7: a) Plot of observed and simulated incident velocities vs elevation for the Nicolum quarry dataset using the 3D model. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data......................................................................................................................... 189 Figure 7.8: a) Plot of observed and simulated rebound velocities vs elevation for the Nicolum quarry dataset using the 3D model. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data......................................................................................................................... 190 Figure 7.9: Effect of the shape factor on the simulated jump height and velocity with constant roughness, column a) is Screen 1 and column b) is screen 2. ..................................................... 191 Figure 7.10: Comparison of scatter angles for a) Shape factor = 1.5, b) No shape factor, c) optimum simulation without the shape factor and d) cumulative distributions calculated for a), b), and c). .................................................................................................................................... 192 Figure 7.11: Impact rate map for the Vaujany site using the Shape Factor, aggregate of 10 000 trajectories. .................................................................................................................................. 193 Figure 7.12: Spatial distribution of trajectories with modeled linear kinetic energy greater than a) 500 kJ, b) 1000 kJ, and c) 1500 kJ. Maps produced using an aggregate of 500 model trajectories...................................................................................................................................................... 195 xxiv List of Symbols Notes: Bold values are used to designate vectors and matrices. Subscript values are used to designate directions Superscript values are used to designate incident or rebound conditions e “Classical” coefficient of restitution, ratio of rebound over incident normal velocity for an impact D Equivalent spherical diameter of a rock E Young’s modulus EK Linear kinetic energy Eω Rotational kinetic energy G Modified Goldsmith equations for the collinear impact of a sphere on a plane surface in Impact incidence, the moment just before an object makes contact with a surface k Restitution factor, alternative to coefficient of restitution used in the PIERRE 2 model n Direction normal to an impact surface n′ Direction normal to the impact surface rotated by the surface roughness angle R Equivalent spherical radius for a rock Rθ Rotation matrix to transform incident vectors from the n - t plane to the n′- t′ plane re Impact rebound, the moment just after the restitution phase of an impact S Random number sampled from a uniform distribution between 0 and 1. t Direction tangential to an impact surface t′ Direction tangential to an impact surface rotated by the surface roughness angle v Linear velocity vector component xxv v Velocity vector = (vn, vt, ω) α Impact angle β Local slope angle γ Friction limit angle, inclination of the incident impulse vector relative the vector tangential to the impact point. θ Surface roughness angle ν Poisson’s ratio ρ Density ω Angular velocity . xxvi Acknowledgements I would like to offer my thanks and gratitude to my supervisor, Dr. Oldrich Hungr. It has been a great privilege to learn from you and work with you. I would also like to thank all those who have offered advice and suggestions during committee meetings and in the review of this thesis, Dr. Erik Eberhardt, Dr. Davide Elmo, Dr. Roger Beckie, and Dr. Eldad Haber. Thanks to my colleagues within the Geological Engineering group for everything, academic and social. And a special thanks to my office-mates for the impromptu opinions and brain storming sessions. I give sincere thanks to my parents, David and Maxine, and my siblings, Evan and Tara, for providing support and encouragement throughout my life. I would also like to thank Mila Mello for being such a great companion to me and making our home such a wonderful place to come back to during the writing of this thesis. I would also like to show my gratitude to my co-authors in Chapter 5. It was a joy to work with Dr. Valentin Gischig, I learned a lot and I greatly enjoyed our time together. Thanks to Dr. Franck Bourrier for his helpful suggestions during the editing of the paper and for facilitating my trip to France. The authors are grateful to the research team at IRSTEA, Grenoble, led by Dr. Fred Berger, for providing the Vaujany test data. Professor Alexander Preh of the Vienna University of Technology participated in the development of the hyperbolic impact model. The study was partly funded by the Swiss National Science Foundation, Project P300P2_151291, as well as by a National Science and Engineering Research Council of Canada (NSERC) Accelerator Grant to O. Hungr. An anonymous reviewer provided constructive suggestions. Many thanks to Dr. Alexander Preh for taking the lead for the paper included as Chapter 6. Your insights and enthusiasm for the project was a great source of motivation. The research xxvii was supported by an NSERC operating grant to O.Hungr. For their support of the research project the following are gratefully acknowledged: Asamer Kies- und Betonwerke GmbH, Lafarge Perlmooser Gmbh, Schotter- und Betonwerk Karl Schwarzl Betriebsgesellschaft mbH, Pornat Steinbruch Preg GmbH, Wopfinger Baustoffindustrie GmbH, Dolomit Eberstein Neuper GmbH, Klöcher Basaltwerke GmbH & Co KG (Asamer Kies- und Betonwerke GmbH), Basaltwerk Pauliberg GmbH & Co KG, VA Erzberg GmbH, Cemex Austria AG, and the British Columbia Ministry of Transportation and Highways, Canada. 1 Chapter 1: Introduction The following Chapter provides an introduction to what rock fall is and why it is important to understand the hazards and risks associated with rock fall. The individual components that define rock fall hazard are reviewed to provide the larger context for the research program detailed in the remainder of the thesis. The use of computer models to assess the reach and intensity of rock fall events is introduced. The specific objectives of this research program are stated, and the overall organization of the thesis is provided. 1.1 Background Rock fall is the process of detachment, falling, bouncing and rolling of rock fragments, either singly or in clusters (Hungr et al., 2014). A key characteristic of this process that distinguishes it from other landslide types is that the motion of the rock fragments is determined by the dynamic interaction of the fragments with the slope substrate, with interactions between individual fragments having negligible effect (Hungr et al., 2014). Rock falls are a common hazard in mountainous areas, as well as man-made environments such as road cuts, open pit mines and quarries. The pervasiveness of this hazard has made rock fall an area of extensive research, and a continuing challenge for design engineers seeking to minimize the risks posed by rock fall. Risk is defined as the product of hazard and vulnerability (Lee and Jones, 2004). Hazard is the possibility of an event occurring, resulting in an interaction with an element at risk (Lee and Jones, 2004). The hazards posed by rock falls to a location can be broken down into three basic components, as described by Volkwein et al. (2011). The first is the onset, which is the probability that an unstable block of a given size detaches from a slope at a given source location. The next is the reach, which is the probability that the rock fall reaches the location at 2 risk. The final element is the event intensity, which is a measure of the bounce height and kinetic energy of the block when it reaches the location at risk (Volkwein et al., 2011). Equation 1.1 gives the mathematical description for the hazard, H, during a time interval T, at location i, a distance of L from the source, from an event with magnitude j and intensity k. Equation 1.1: 𝐻𝑖𝑗𝑘 = 𝑃(𝐿)𝑗 × 𝑃(𝑇|𝐿)𝑖𝑗𝑘 The vulnerability of an element at risk describes the outcome of an interaction (Lee and Jones, 2004). Assessing the vulnerability requires assigning a probability of some potential harm to a given intensity of event, which is often difficult to do (Hungr et al, 1999). Delving into the complexities of a quantitative assessment of vulnerability (c.f. Massey et al, 2012) is beyond the scope of this thesis. However, it is important to recognize that having a good understanding of the intensity of a rock fall event is critical in understanding the hazard the event poses and developing an appropriate hazard reduction strategy. The following paragraphs present a more detailed examination of the methods to determine the terms of the hazard equation. Assessing the onset probability first requires the identification of source areas. The basic requirements for rock fall are a slope steep enough to allow the blocks to move by gravity and a trigger to detach the unstable block from the source area (Volkwein et al, 2011). Typical site and environmental factors that lead to rock falls are summarized in Table 1.1. 3 Table 1.1: Factors contributing to an area being a rock fall source (adapted from Volkwein et al., 2011). Slope Characteristics Environmental Factors Morphology slope angle, height, exposure Weather & climate rainfall and snowmelt, freeze thaw cycles Geology rock types, structure (bedding, folding, fracture orientation, length and persistence), previous weathering Active tectonics & seismicity uplift exposing more rock to weathering, seismic ground motions Material properties cohesion, friction, strength, durability Weathering & erosion combination of previous factors results in varying weathering and erosion rates Hydrogeology matrix and joint permeability, infiltration rates, static groundwater level Human activities quarrying, road cuts Traditionally, rock fall onset has been assessed retroactively by determining the frequency and size of observed rock falls in known hazardous areas and creating a “rock fall inventory”, which establishes qualitative hazard levels (Pierson et al., 1990). This inventory method can also be extended to create magnitude-frequency relationships, following the observation that small volume rock falls happen frequently and large rock falls happen less frequently. This behavior can be approximated by a power law relationship (Hungr et al., 1999). Rock fall onset can be determined proactively using analytical methods such as limit equilibrium analysis or key block analysis, where jointing and other physical characteristics of the rock mass are used to calculate the factor of safety (FS) for movement (Bruckno et al, 2010). This method can provide a very accurate assessment of stability and the expected size of unstable 4 blocks, but requires a large amount of data collection and modeling by competent professionals, which means it is usually only applied to small, high consequence areas (Bruckno et al, 2010). To examine the potential source areas within a large region, GIS based methods are typically used. These will use general parameters, such as the slope angle compared to some threshold angle, and may include a variety of other parameters, including presence of faults or other structural geology information, proximity of talus slopes, and cliff heights (Volkwein et al, 2011). Incorporating detailed structural geology information, such as the orientation, spacing and trace length of joint sets into large scale GIS models has been proposed to identify areas of instability (Jaboyedoff et al., 2004), however this has yet to be widely used in practice. The work detailed in this thesis has used experimental or case study data with known rock sizes. Determining the expected distribution of rock sizes for a site as an additional forward analysis is beyond the scope of this thesis. Rock fall reach has traditionally been assessed using empirical methods. The simplest method uses the concept of the fahrböschung angle, which is the angle formed by the elevation difference between the highest point on the source area over the maximum horizontal distance traveled by a rock (Evans and Hungr, 1993). This concept was originally proposed by Albert Heim for predicting the runout of rock avalanches, and has been applied to the study of fragmental rock falls, but the angles found at different sites tend to vary widely (Turner and Duffy, 2012). An alternative approach is to look at the angle formed by the change in elevation from the apex of the talus deposit below the source area to the farthest run out boulder and the horizontal distance from the talus apex to the boulder, referred to as the minimum shadow angle (Evans and Hungr, 1993). This concept takes into consideration the observation that large initial fall heights tend to result in major energy loss when the boulder first contacts the talus deposit, 5 meaning that the kinetic energy available at the start of the runout is largely independent of the initial fall height. This helps to reduce the observed variability obtained when using a fahrböschung angle. A minimum shadow angle of 27.5° is recommended as a first estimate for runout distances (Evans and Hungr, 1993). Determining the intensity of a rock fall event is not a trivial matter. As the motion of a rock fragment is governed by its interaction with the slope, the processes happening during an impact with the slope are a primary concern. In the context of rock fall, these processes are exceedingly complex and achieving exact, deterministic solutions is not feasible. The current state of practice is to use computer models for the determination of rock fall reach and intensity (Bourrier et al., 2012). Due to the difficulties in accurately measuring the kinematics of the particles, information such as linear and angular velocities and bounce heights have generally not been incorporated into model calibrations (Bourrier et al, 2012). 1.2 Statement of Purpose and Objectives The purpose of this thesis is to aid in the development of an improved computer model for the prediction of rock fall trajectories and dynamics. The objectives of this work are:  To develop an impact model that is consistent with the established principles of impact mechanics and observed behavior of rock falls, with a focus on simultaneously being able to represent runout and kinematic behaviors of the rocks in 2D and 3D.  To make the model as robust as possible, using a minimum number of model parameters, so that reliable results can be achieved with an amount of site information that corresponds to normal site investigation practice, such as topography and general material descriptions. 6  To calibrate the model using a variety of case studies, focusing on sites with detailed information on the particle kinematics, in order to provide model input parameters based on general slope descriptions at a range of locations for forward analysis. 1.3 Organization of the Thesis Chapter 2 provides a review of literature relevant to this subject, starting with a review of basic impact mechanics, moving to the application of these principles to rock fall impacts, and finally to computer modeling of rock falls. Chapter 3 outlines the major features of the PIERRE 2 model and the theoretical basis of those features. Chapter 4 presents the calibration of the PIERRE 2 stochastic rock fall simulation program. Chapter 5 outlines the development and calibration of the 3D version of the model, PIERRE 3D, in the form of a journal manuscript. Chapter 6 presents the results of a model calibration using full scale experiments from three Austrian quarries and one Canadian quarry, also in the form of a journal manuscript. Chapter 7 provides the results of work with the 3D model that has incorporated features outlined in Chapters 3 and 4, but not included in Chapter 5. In Chapter 8, the main discussion points resulting from the research are reviewed, and areas for future research are presented. 7 Chapter 2: Literature review 2.1 Impact Mechanics Rock fall is characterized by the dynamic interactions of a falling rock and the slope, thus understanding rock fall requires an understanding of impact mechanics concepts. The study of physical impact phenomena began simultaneously with the study of mechanics, with pioneering studies on rigid body impacts being performed by Galileo and then Newton in the 17th Century (Goldsmith, 1960). Since then, the study of impacts has incorporated other concepts from the study of mechanics, including vibrational effects, stress-strain behavior, and elastic, inelastic or plastic material responses. Despite these advances, the representation of impacts occurring between rigid bodies with correction factors applied for energy losses remains in common usage to the present day (Goldsmith, 1960). The following section introduces the concepts of stereomechanical impacts and coefficients of restitution, which have changed little since the time of Isaac Newton. A brief examination of the stresses and strains occurring during impacts, and the effect of material properties and impact geometry is then presented. Results from experimental studies involving projectiles impacting granular soil are discussed to aid in understanding the theoretical constructs of the contact laws in the context of rock fall dynamics. Finally, issues relating to fragmentation during impact are introduced. 2.1.1 Stereomechanical impacts The fundamentals of impact mechanics are described by stereomechanics, which is based on the impulse-momentum law for rigid bodies (Goldsmith, 1960). Stereomechanical theory allows the prediction of rebound velocities and applied linear and rotational impulses for an impact from an object with a given incident velocity. For an ideal, perfectly elastic collision the 8 law of conservation of momentum can be applied directly to find the final velocities. Likewise an ideal, perfectly plastic collision will have final velocities of zero in all directions. For an inelastic collision, a coefficient of restitution, designated “e”, is introduced. Goldsmith defines the coefficient of restitution for an object impacting a massive, rigid surface as the ratio of the rebound normal velocity over the incident normal velocity. Assuming the object loses no mass during the impact, this is exactly equal to the ratio of rebound and incident impulses in the plane normal to the impact plane, summarized in Equation 2.1. The coefficient of restitution is used to represent the degree of plasticity in the collision. The theory has a number of limitations, some of which will be discussed further, however, it performs well for the collision of two spheres or a sphere with a large rigid mass (Goldsmith, 1960). Equation 2.1: 𝑒 = 𝑣𝑛𝑟𝑒𝑣𝑛𝑖𝑛= 𝑃𝑛𝑟𝑒𝑃𝑛𝑖𝑛 A special case of a stereomechanical impact is presented for the scenario where a rotating sphere is moving in a 2D plane and impacts a rigid wall (Goldsmith, 1960). In this case, the mass of the wall can be considered infinite with an incident velocity of zero, resulting in a rebound velocity of zero for the wall. A local coordinate system can be defined, spanned by vectors normal, n, and tangential, t, to the point of contact. If the vector normal to the contact point intersects the centre of mass of the sphere, the impact is collinear (Stronge, 2000). Collinearity allows for the decoupling of the linear and rotational momentum, thus the equations for linear and rotational momentum conservation can be solved separately (Stronge, 2000). This simplifies the mathematics of finding the rebound velocities greatly and allows for greater computational efficiency. 9 Relative to this local n-t coordinate system, the sphere has an incident velocity vector vin = (vnin, vtin, Rωin), and the contact point of the sphere and wall has a coefficient of restitution, e, and a contact friction coefficient, φ. The inclination of the incident impulse vector, referred to as the friction limit angle γ, can be calculated using impulse-momentum theory, as shown in Equation 2.2 (Goldsmith, 1960). Equation 2.2: tan 𝛾 = 2(𝑣𝑡𝑖𝑛 + 𝑅𝜔𝑖𝑛)7𝑣𝑛𝑖𝑛(1 + 𝑒) Tangential and rotational motion relative to the contact point will result in the sphere sliding at the contact point, referred to as slip. For a given set of impact conditions, there are two possible situations, the first being the case in which φ > γ, termed a “contained slip impact”, and the second in which φ < γ, termed an “uncontained slip impact”. For the contained condition, the slip at the contact point is contained within the time of impact so that the tangential and angular velocities become perfectly synchronized. In this case the coefficient of restitution alone is used to calculate the rebound velocity vector, shown in Equation 2.3 (Goldsmith, 1960). Equation 2.3: |𝑣𝑡𝑟𝑒𝑣𝑛𝑟𝑒𝑅𝜔𝑟𝑒| = ||570 −270 −𝑒 0−57027|| × |𝑣𝑡𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝜔𝑖𝑛| For the uncontained slip condition, the slip has not ended at the end of the contact period, thus the tangential and angular velocities have not synchronized. The coefficient of restitution and contact friction coefficient are used to calculate the rebound velocity vector, and Equation 2.4 can be used to predict the final velocity vector of the sphere. 10 Equation 2.4: |𝑣𝑡𝑟𝑒𝑣𝑛𝑟𝑒𝑅𝜔𝑟𝑒| = |1 − tan𝜑 (1 + 𝑒) 00 −𝑒 00 −52tan𝜑 (1 + 𝑒) 1| × |𝑣𝑡𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝜔𝑖𝑛| During an impact, the angular velocity will either increase or decrease in order to get closer to, or achieve, synchronization with the tangential velocity, i.e. vtre = Rωre. If the motion of the particle is synchronized already during incident motion, no frictional energy loss due to sliding at the contact point will occur. If it is not synchronized, frictional energy will be expended to achieve the synchronization by slip between the surface of the sphere and the substrate. In the contained slip condition, which is the most common condition for this impact model, as the rotational velocity increases with all other terms in the equation held constant, the friction limit angle will increase, and the amount of energy required to synchronize the rotation with the tangential velocity will decrease. This leads to more of the linear kinetic energy being conserved in the impact. A point will be reached where φ < γ, and the impact will transition to an uncontained slip condition. In this case, no additional energy transfer from the linear motion is needed to synchronize the rotation, and as can be seen in Equation 2.4, the angular velocity no longer influences the calculation of the rebound tangential velocity. This leads to a constant conservation of translational kinetic energy when the other terms in the equations are held constant. This is shown graphically in Figure 2.1, which was obtained by plotting the total kinetic energy corresponding to the incident and rebound motion vectors in Equation 2.3 and Equation 2.4. Assuming the particle is spherical, the linear and angular kinetic energy, EK and Eω, respectively, are found using Equation 2.5. . 11 Equation 2.5: 𝐸𝐾 = 12𝑚𝑣2 𝐸𝜔 = 12𝐼𝜔2 𝑤ℎ𝑒𝑟𝑒 𝐼 = 25𝑚𝑅2 Figure 2.1: Ratio of rebound over incident kinetic energy as a function of the interface friction minus the friction limit angle. The equations developed for 2D stereomechanical impacts can also be extended into three dimensions, as demonstrated by Goldsmith (1960). The underlying concepts remain the same, and were exploited in the formulation of the 3D model described in Chapter 5. 2.1.2 Contact laws One of the limitations of the stereomechanical approach is that it gives no indication of the stresses and deformations at the impact point. This is of practical significance for applications such as determining the force applied to a protective structure. The Hertz Law of 12 Contact attempts to address some of the limitations of assuming instantaneous impacts for a stereomechanical impact. The Hertz Law describes the force – deformation relationship for static compression of two isotropic, elastic bodies with smooth surfaces (Goldsmith, 1960). This theory has been used for impact situations where permanent deformations were produced, so long as plastic yielding extends only within the vicinity of the contact point. The compressive stress distribution must be in equilibrium, thus:  Total applied force must equal total resisting force supplied by the vertical pressure within the contact area;  Displacements in the interior of the body are governed by general elastic stress – strain relationships, with the stress decreasing rapidly as distance from the contact area increases;  Components of displacement vanish at infinity, so that they can practically be neglected at distances from the contact area that are large relative to the dimensions of the contact area; and,  Along the contact area the shear stresses must vanish, normal stresses acting on the two bodies must balance, and outside the contact area the normal stresses must vanish. From these constraints, the contact area geometry and the distribution of stresses over the contact area can be calculated for a given applied force. Considering the case of a sphere impacting an extensive plane surface the contact area geometry and stress distribution can be calculated from the Young’s modulus, E, Poisson’s ratio, ν, and density, ρ, of the materials, and the radius of the sphere (Goldsmith, 1960). The complete derivation is presented by Goldsmith, 13 and will not be reproduced here. If we consider impacts of the same material, i.e. E, ν and ρ are constant, the contact geometry and stresses are a function of the radius. When these contact geometry and stress relations are combined with Newton’s second law of motion, a force-time relationship can be approximated for an idealized elastic impact (Goldsmith, 1960). To do this the initial velocity of the sphere that results in the applied force is the only additional parameter required. The Hertz Law of Contact provides a useful conceptual model of the stress distribution over a contact area; however it has some significant limitations. The most obvious is that it neglects plastic yielding around the contact area. The Meyer Law addresses some of these limitations, as described by Goldsmith (1960). The Meyer Law describes three phases of contact:  Initial elastic compression (as described by Hertz law) until a critical stress, p0, is attained.  Plastic zone moving outward from the centre of the contact area with constant pressure p0, with an elastic annulus, until the relative velocity of the bodies is zero.  Restitution phase where the elastic strain energy from the deformed area is recovered. The Meyer law requires the same input parameters as the Hertz law, with the addition of a new material property, the critical stress for plastic deformations to occur, p0.This is an improvement on the Hertz contact law because it can predict permanent deformations. However, this is still a significant simplification of the actual processes occurring during an impact (Goldsmith, 1960). 2.1.3 Impacts on granular soil The Hertz and Meyer contact laws were both derived based on the results of experiments involving metal objects impacting one another. These ideas were applied to soils in experiments 14 conducted with projectiles impacting a granular soil surface (Forrestal and Luk, 1992). The authors formulated equations to describe the penetration of projectiles with regular geometry into a granular soil, assuming the projectile penetrates by forming a spherically symmetric cavity in the soil. For the projectile to penetrate deeper into the granular surface the cavity must expand. The expansion of the cavity for a given amount of stress on the contact area can be determined by the frictional properties of the soil. Using the Mohr-Coloumb and Tresca criteria, the authors described the force and stress generated by the projectile at impact, and the formation of elastic and plastic zones in the soil. They then tested this theoretical framework with experiments and found a good agreement between the theory and experimental results (Forrestal and Luk, 1992). Full scale experiments dropping granitic boulders onto gravel pads were conducted by Pichler et al. (2005). The objective was to test the relations between three key physically quantities used in the design of cushioning layers for rock fall protective structures, the rock penetration depth, impact duration and impact force, resulting from differences in the boulder mass, fall height, and indentation resistance of the gravel substrate. Dimensional analysis and Pi-theorem were used to establish dimensionless relationships for the penetration depth, impact duration and impact force. The experiments were set up using test boulders of the same material with an approximate cubic shape dropped so that they would impact on a point, and a standard gravel material was used for the substrate. This allowed the authors to make the assumptions that the test boulders were geometrically similar over the impact area, and the boulders and gravel have constant density. Combined with the dimensional analysis, this allows for description of the impact process with the mass, incident velocity and characteristic diameter of the impact area of the boulder, and the indentation resistance of the gravel, which is a strength-like parameter (Pichler et al., 2005). 15 From this brief examination of contact laws and impacts on granular soils, it can be seen that the various methods for determining the stresses over a contact area have some similar characteristics. All depend on the incident velocity and size of the impacting projectile as well as some material properties to describe the strength and deformation properties of the projectile and/or substrate. These concepts will be utilized in the upcoming Chapter on the model development. 2.1.4 Fragmentation The previous sections of this Chapter deal with particles that may deform, but do not break as the result of an impact. Neglecting fragmentation of boulders during impacts is commonly done in design practice, despite fragmentation during rock falls being a widely observed phenomenon (Giacomini et al., 2009). The perils of this were highlighted by comparing the kinematic properties of different sizes of rocks from full scale testing, reported by Agliardi and Crosta (2003). Using their kinematic model, the authors were able to correctly predict maximum velocities for rocks with volumes ranging from 0.05 m3 to 2 m3. Despite that success, it was found that an area on the test slope where fragmentation occurred frequently had a number of relatively small fragments which were ejected from impacts at velocities greater than the modeled maximum (Agliardi and Crosta, 2003). Situations such as this lead to the so called “bullet effect” (Buzzi et al., 2011), which means that smaller fragments impacting a barrier or fence at high velocity can cause localized failures. Although the fragments have a lower kinetic energy due to their smaller mass, the small impact area results in very high stress concentrations (Buzzi et al., 2011). Typically large particles are the ones considered in a rock fall analysis, however, the possibility of fragmentation should be considered as it may be the more hazardous scenario. 16 Much of the research on rock fragmentation has focused on excavations in rock. An example of this was a study conducted by Zhang et al. (2000) which examined the effect of loading rates on rock fracture formation to optimize things like the impact speed for percussive drilling or the rotation speed of a boring machine. Zhang et al. separated the energy absorbed by a rock specimen into two components, the fracture and damage energy, EFD, and kinetic energy of flying fragments, EK. As one would expect, increasing the impact velocity or loading rate, thus increasing the initial kinetic energy, will result in both EFD and EK increasing (Zhang et al., 2000). In an attempt to apply these ideas to rock fall modeling the concept of a fragmentation energy threshold has been put forward, discussed by Giacomini et al. (2009). When the incident energy of an impact exceeds the fragmentation energy threshold, the block is randomly divided into smaller fragments with the rebound kinetic energy distributed amongst them. Full scale tests of schistose rocks dropped in quarries were conducted by Giacomini et al. to test the effect of different impact angles and kinetic energy on fragmentation for rocks with anisotropic strength. The results of the experiment showed that the observed fragmentation is very sensitive to the impact angle, which suggested approaches using threshold fragmentation energy were not adequately representing the actual fragmentation process for rocks with anisotropic strength (Giacomini et al. 2009). In the present development of the model discussed here, particle fragmentation has not been considered. 2.2 Coefficients of Restitution The concept of a coefficient of restitution as a constant material property dates back to the work of Isaac Newton, and has been retained by many practitioners to this day. A number of authors have suggested values for coefficients of restitution to be used in rock fall analyses. In 17 the following discussion it will be demonstrated that the momentum and energy conservation processes during a rock fall event are too complex to be accurately represented by a constant restitution coefficient. One of the recurring challenges through the history of the study of rock fall dynamics is accurately representing the coefficients of restitution for impacts. It should be noted here that coefficient of restitution may refer to the ratio of rebound over incident velocity, kinetic energy or impulse (Turner and Duffy, 2012). Many authors have attempted to find representative coefficients of restitution for different slope material types. A summary of some of these proposed coefficients of restitution (all based on the ratio of rebound over incident velocity) are given in Table 2.1 (Heidenreich, 2004). 18 Table 2.1: Summary of normal and tangential coefficients of restitution (en and et, respectively), related to slope material type (adapted from Heidenreich, 2004). Material en et Reference Smooth hard rock/ pavement 0.8 – 0.9 0.65 – 0.75 Piteau & Clayton, 1977 0.95 0.9 Heierli, 1985 0.4 – 0.53 0.9 – 0.99 Hoek, 1987 0.05 – 0.35 0.5 – 1 Urciuoli, 1988 0.37 – 0.42 0.87 – 0.92 Pfeiffer & Bowen, 1989 0.5 0.95 Giani, 1992, Barbierei et al, 1988 0.9 0.9 Evans & Hungr, 1993 0.17 – 0.43 0.45 – 0.88 Gerber, 1995 0.1 – 0.35 Kamijo, 2000 0.6 – 1.0 0.9 – 1.0 Jones et al., 2000 0.2 0.53 Budetta & Santo, 1994 Boulders/ Boulders covering bedrock 0.5 – 0.8 0.45 – 0.65 Piteau & Clayton, 1977 0.35 0.85 Hoek, 1985 0.33 – 0.37 0.83 – 0.87 Pfeiffer & Bowen, 1989 0.35 0.85 Giani, 1992, Barbierei et al, 1988 0.15 – 0.30 0.75 – 0.95 Jones et al., 2000 Talus 0.4 – 0.5 0.35 – 0.45 Piteau & Clayton, 1977 0.45 0.2 Heierli, 1985 0.32 0.8 – 0.82 Hoek, 1985 ~ 0 0.24 Urciuoli, 1988 0.30 – 0.33 0.80 – 0.87 Pfeiffer & Bowen, 1989 0.3 0.7 Giani, 1992, Barbierei et al, 1988 0.7 Evans & Hungr, 1993 0.12 – 0.2 0.65 – 0.95 Jones et al., 2000 Vegetated/ Soft Soil 0.2 – 0.4 0.2 – 0.3 Piteau & Clayton, 1977 0.3 0.8 Hoek, 1987 0.28 – 0.32 0.78 – 0.83 Pfeiffer & Bowen, 1989 0.25 0.55 Giani, 1992, Barbierei et al, 1988 0.1 – 0.2 0.5 – 0.8 Jones et al., 2000 19 It can be seen from Table 2.1 that the values for the coefficients of normal and tangential restitution vary widely within a material type, different studies found contradicting values for the coefficients, and several materials have overlapping ranges. Some of the studies summarized in Table 2.1 looked at factors that can affect the observed coefficient of restitution, which can be broadly grouped into block properties, slope properties and the kinematics of the collision (Heidenreich, 2004). To summarize, increasing block mass tends to decrease the normal restitution. The studies looking at the effect of the rock’s angularity had conflicting results. Increasing the slope’s Young’s modulus and dry density both resulted in increased coefficients of restitution. Increasing the impact angle resulted in a decrease in the normal restitution for every study that looked at that factor; however it had little effect on the tangential restitution values. It was also found that increasing the impact velocity would decrease the restitution coefficients (Heidenreich, 2004). Recent work has looked more deeply into these factors that affect the observed coefficients of restitution. Labiouse and Heidenreich (2009) conducted half scale laboratory experiments of impacts on sandy slopes to examine the effects of the ground properties, block properties and the kinematics on the normal, tangential and energetic coefficients of restitution. In this study, blocks with masses ranging from 100 to 1000 kg were tested on a sandy test slope with an inclination from 0 to 30° with drop heights up to 10 m. This work confirmed, in a controlled setting, the relationships between block mass, impact angle, and impact velocity discussed in the previous paragraph. However, it was also found that the relationship between increasing mass and decreasing normal coefficient of restitution started to level out to an almost constant value for masses larger than 500 kg. These tests used both spheres and cylindrical blocks with a hemispherical bottom. Comparison of the coefficients of restitution for spheres and a cylindrical 20 block with similar masses showed the cylindrical block had higher coefficients of restitution in both the normal and tangential directions. This was attributed to the larger radius of the base of the cylindrical block which reduced the penetration of the block into the surface allowing for easier translational movement and a greater normal coefficient of restitution. It was also observed that the cylinders tended to rotate more than the spheres, especially at high incident angles. The authors concluded that the motion of the test projectiles were governed by mechanisms of penetration, sliding and rotation, which to some extent are interdependent of one another. This demonstrates some of the problems in the representation of impact processes simply with a normal and tangential restitution coefficient. (Labiouse and Heidenreich, 2009). Asteriou et al. (2012) carried out lab testing with 5 different rock types with different material hardness with free fall drops and parabolic arcs with different impact angles and initial velocities. Limited field testing was also performed for verification purposes. Normal and kinematic coefficients of restitution were found to increase with material hardness and impact angle. No sensitivity was detected for the tangential coefficient of restitution. The field testing conducted generally agreed with the results of the lab testing (Asteriou et al., 2012). Ferrari et al. (2013) analysed video records from full scale experiments to determine the coefficients of restitution for impacts in a natural setting. The authors used angular boulders with a volume of approximately 0.8 m3, obtained from the talus deposit of the cliff where they conducted the tests. In these tests they observed the effect of the boulder size on its motion. If the boulder was larger than the substrate over which it was travelling it would tend to gain momentum until it encountered a flatter section of slope, or material of approximately its own size. If similarly sized material was encountered in the trajectory, the boulder would start to decelerate until it became trapped in a void within the talus. The video records of the impacts 21 were also analysed to determine the coefficients of restitution. This analysis resulted in normal restitution values greater than 1 in many cases. This phenomenon was found to be strongly affected by the block shape, with elongated blocks tending to impact in such a way that the longer moment arm at the impact point tended to propel the blocks away from the slope. Normal restitution, en, values greater than 1 were also more common for impacts with a shallow incidence angle (Ferrari et al., 2013). The phenomenon of impacts with en > 1 was also observed by Spadari et al. (2011). The results of video analysis from four sites in Australia indicated that the average en was slightly greater than 1 for all four sites (Spadari et al., 2011). However, when the tangential restitution was considered with the normal, there was always a net loss of kinetic energy (Spadari et al., 2013). The value of en was related to the impact angle, αin, using a linear relationship, shown in Equation 2.6. Equation 2.6: 𝑒𝑛 = −0.04𝛼𝑖𝑛 + 1.8 + ∆,𝑤𝑖𝑡ℎ ∆ ∈ [−0.5, 0.5] The correlation between the incidence angle and the normal restitution has also been described using a power law, shown in Equation 2.7 (Wyllie, 2014). Equation 2.7: 𝑒𝑛 = 19.5(𝛼𝑖𝑛)−1.03 Both Spadari and Wyllie found that the proposed relationships are accompanied by very strong variability, especially at low values of incident angle, where the range of the en values is nearly equal to the average. It is interesting to note that the two trend lines plot very close to one 22 another for 10° < αin < 30°, shown in Figure 2.2, indicating the two trends are in fact quite similar. Figure 2.2: Graphical comparison of the relations described by Equation 2.6 and Equation 2.7. Chau et al. (2002) performed laboratory testing utilizing impacts of formed plaster spheres and a plaster slope. In these experiments, it was also found that decreasing the incident angle had the effect of increasing the normal coefficient of restitution. The authors examined the ratio of rotational kinetic (Eω) energy over linear kinetic energy (Ek) and experimentally found it to have an upper bound of 0.4, which is in agreement with a theoretical impact involving a sphere. The upper bound corresponds to the case when the rebound angle is negligibly small, i.e. all of the rebound velocity is in the tangential direction, and the rotation perfectly synchronizes with the translational velocity during the impact, thus ωre = vtre/R. By substituting this upper bound relationship into the equation for the angular kinetic energy of a sphere, then dividing that 23 by the translational kinetic energy, the resulting ratio is Ek/Eω = 2/5, or 0.4. The amount of rotational kinetic energy induced at each impact was found to be a function of the impact angle as well, with the rotational energy increasing with αin up to 40°, and then decreasing to a negligible value with αin approximately equal to 70° (Chau et al., 2002). In summary, detailed studies of restitution “coefficients” indicate strong variability, and defy attempts to make specific classifications on the basis of material type. The conservation of momentum during an impact shows a strong dependence on the incident impact angle. Under certain impact conditions, the apparent normal restitution coefficient can be greater than 1, violating another commonly held assumption. In light of this evidence, the concept of a constant restitution coefficient is not tenable. 2.3 Computer Modeling for Rock Fall The use of computer simulations in determining rock fall trajectories dates to the 1970’s (Piteau and Clayton, 1977). Since then a multitude of models have been developed. The following is by no means a comprehensive review of all the models that have been developed, but it does provide an overview of some of the significant developments in computer modeling. A brief overview of the general approaches to computer modeling is presented. Some specific examples of common rock fall simulation programs are given, with a discussion of the model approaches they adopt. 2.3.1 Spatial representation Three broad categories are used to define the spatial representation for rock fall models: 2D, 2.5D or 3D (Volkwein et al., 2011). Some 2D models examine the planar trajectory, i.e. in the x and y directions, by representing the boulder as a sliding block moving down the path of steepest descent from the source area. This representation cannot provide any information on the 24 bounce heights associated with the trajectory. The more common 2D representation is to calculate the trajectory between bounces along a profile section, in the x and z directions. A 2.5D model uses a concatenation of the previous 2 representations, where a fall path is defined by the path of steepest descent from the elevation model, which is then used as a 2D profile for the trajectory calculation in the x and z directions (Volkwein et al., 2011). In 3D models, the slope is represented by a 3D surface and the trajectory is allowed to move in all 3 directions (Volkwein et al., 2011). The primary advantage of this is the ability to explicitly represent the effects of converging or diverging topographic features (Volkwein et al., 2011) . 2.3.2 Lumped mass models The simplest representation of the rock fall process is using a lumped mass model, where the momentum of the boulder is represented as being concentrated at a single, dimensionless point (Volkwein et al., 2011). For the most basic model definition the rock is represented only by its mass and translational velocity vector (Leine et al., 2014). If the dimensions of the particle are considered, it is only for the determination of contacts with the slope or obstacles. An extension of the lumped mass model is the sphere model, where the boulder is represented as a solid, homogenous sphere which is represented by the mass and moment of inertia, as well as the translational and angular velocity for the sphere (Leine et al., 2014). Many recent developments with lumped mass models have focused on incorporating them in GIS. Examples include STONE (Guzzetti et al., 2002) and Rockfall Analyst (Lan et al., 2007) which use simple impact models in order to produce rock fall hazard maps for large areas. 25 2.3.3 Rigid body models The rigid body approach represents the boulder as a geometric shape, typically a sphere, cube, cylinder or ellipsoid (Volkwein et al., 2011). For a rigid body model the boulder is described by its geometry, its mass, the location of the centre of mass within the block geometry, the inertial tensor, as well as the orientation of the block, and the translational and angular velocity vectors (Leine et al., 2014). 2.3.4 Stochastic models Many models use stochastic methods to represent the natural variability in rock falls (Turner and Duffy, 2012). Parameters describing any or all of the model inputs, such as the block shape, block mass, location of source areas and mechanical properties of the slope, can be varied from simulation to simulation or impact to impact by randomly sampling from predefined probability distributions for the stochastic parameters (Turner and Duffy, 2012). In contrast to the approach of using stochastic variables which are subsequently fed into a deterministic model, a fully stochastic method was developed by Bourrier et al. (2009a). In this method the incident velocity vector is modified by a coefficient matrix, A, in which each element, a1 to a9, is a normal random variable, shown in Equation 2.8. Correlations between the elements in A are also considered in this model (Bourrier et al., 2009a). Equation 2.8: |𝑣𝑡𝑟𝑒𝑣𝑛𝑟𝑒𝑅𝜔𝑟𝑒| = |𝑎1 𝑎2 𝑎3𝑎4 𝑎5 𝑎6𝑎7 𝑎8 𝑎9| × |𝑣𝑡𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝜔𝑖𝑛| The calculation of A requires a vector for the mean of each coefficient, a vector of standard deviations, as well as a 3 x 3 covariance matrix (Bourrier, 2008). The model was calibrated using the results of a numerical model which was in turn calibrated to small scale lab 26 tests (Bourrier et al., 2009b). The fundamental assumption of this calibration is that the scale effects are not sufficiently large as to preclude the use of this method for full scale modeling. When tested against the data from tests conducted at the Vaujany test site in France (also used in the calibration of the present model, see Chapter 4 and Chapter 5), the model was able to reproduce the results of the field tests well. One of the primary challenges of using this model is the stochastic variables are obtained from numerical analyses which are calibrated using lab experiments. These lab experiments must be applicable to the field situation which the user would like to model (Bourrier et al., 2009a). The calibration of the stochastic model is complicated and unique and would need to be repeated for each type of slope surface. This represents a practical disadvantage of the method. 2.3.5 Colorado Rockfall Simulation Program The first version of the Colorado Rockfall Simulation Program, CRSP, developed by Pfeiffer at the Colorado School of Mines, was released in 1989, with a total of 4 subsequent versions released (Andrew et al., 2012). The original version of CRSP used a 2D slope representation and a sphere model (Jones et al., 2000). Slopes were represented by a series of line segments, and at an impact the local slope angle is modified by a random surface roughness angle controlled by a prescribed probability distribution. The impact properties are functions of the friction, normal coefficient of restitution and tangential coefficient of restitution (Jones et al., 2000). In this model, the surface roughness is defined as the angle formed by the tangent of the perpendicular deviation from the average plunge line within a distance equal to the radius of the rock, divided by the radius of the rock (Jones et al., 2000). The average plunge line minus the 27 roughness angle is used to calculate the normal and tangential components of the incident velocity vector. To account for plastic energy losses in heavier impacts, an empirical relationship was developed that described the normal rebound velocity as a function of a constant coefficient of restitution scaled by the incident normal velocity, the metric version shown in Equation 2.9 (Jones et al., 2000). The reference normal velocity of 9.14 m/s is a conversion of the original value of 30 ft/s. Equation 2.9: 𝑣𝑛𝑟𝑒 = 𝑣𝑛𝑖𝑛𝑒𝑛1 + (𝑣𝑛𝑖𝑛9.14)2 The rebound tangential velocity is described as a function of the incident tangential velocity, the mass, radius, and moment of inertia of the block, an empirical friction function, f(F), and an empirical scaling factor, SF. This function is based on the principle of the conservation of energy, shown in Equation 2.10 (Jones et al., 2000). Equation 2.10: 𝑣𝑡𝑟𝑒 = √𝑅2 (𝐼(𝜔𝑖𝑛)2 +𝑀(𝑣𝑡𝑖𝑛)2) 𝑓(𝐹)𝑆𝐹𝐼 + 𝑀𝑅2 The rebound angular velocity is solved using the assumption that the block’s angular and tangential velocities perfectly synchronize during the impact (contained slip condition is assumed), resulting in the relationship shown in Equation 2.11 (Jones et al., 2000). Equation 2.11: 𝜔𝑟𝑒 = 𝑣𝑡𝑟𝑒𝑅 28 The most recent version of CRSP utilizes a 3D rigid body model (Andrew et al., 2012). The model utilizes the Discrete Element Method, DEM, assuming a linear elastic contact force-displacement relationship between the block and the slope. The blocks themselves are discretized into a system of three dimensional spheres, with forces applied at each contact. The model explicitly considers the forces resulting from the penetration of the rock into the slope at the contact point. The normal force is a function of the relative normal velocity at the contact point, the normal contact stiffness, the viscosity of the material, and the relative displacement of the contact points, i.e. the penetration into the slope. The tangential force is a function of the relative tangential velocity at the contact point, shear contact stiffness, friction coefficient, relative displacement, and the previously calculated normal force. The model solves the equations of motion using a numerical integration technique. At each time step, t, the boundary forces for each element ,i, in the model are solved as a function of the position, x, velocity, v, and acceleration, a, vectors, the mass, M, and damping, C, matrices, and the resultant contact force, P, vectors for each element, shown in Equation 2.12: (Andrew et al., 2012). This does not require separate considerations of flight and impact phases, as the contact forces simply drop out of the equations and the acceleration due to gravity remains. Equation 2.12: 𝒇𝑖,𝑡 = 𝑴𝒂𝑖,𝑡 + 𝑪𝒗𝑖,𝑡 + 𝑷𝒊,𝒕𝒙𝒊,𝒕 2.3.6 Rocfall Rocscience Inc. has released a total of 5 versions of the program Rocfall to date, all of which have a 2D, lumped mass boulder representation based on the work of Stevens (1998). In his thesis, Stevens identified the slope geometry, material properties and the rock fall source as being both inherently variable, and important for the accurate representation of rock fall events. 29 The inherent variability was addressed by performing probabilistic analyses. The slope geometry has a probability distribution associated with each vertex, which leads to a random variance in the slope profile from run to run. The variability of the material properties was addressed by assigning a probability distribution to the coefficients of restitution, and the mass of the boulders and the source location could also be varied using independent probability distributions (Stevens, 1998). Later versions of Rocfall have included a consideration of the rotation of the blocks, however the available documentation does not indicate how the rotation is calculated (Rocscience). Rocscience has also adopted the normal restitution scaling factor originally developed in CRSP, shown in Equation 2.9. They have also introduced an empirical mass scaling factor, which is equivalent to Equation 2.9, but with the relationship scaled to the mass divided by a constant as opposed to the normal velocity divided by a constant. A statistical variability in the slope angle, the surface roughness can also be defined (Rocscience). The most recent version of Rocfall (Version 5) has incorporated a 2D rigid body model, based on the work of Ashayer (2007). The impact models proposed by Descouedres et al. (1987), Kobayashi (1990), and Azzoni et al. (1995) have been implemented into this model (Ashayer, 2007). To summarize from Azzoni et al., the impacts are modeled as inelastic, using an assumed, idealized block shape, and the contact area is assumed to be infinitesimal, thus it is represented as a point. These assumptions, combined with energetic coefficients of restitution are used to calculate the rebound velocities (Azzoni et al., 1995). A roughness angle is also used to simulate the effects of slope variations, similar to the method used in the 2D model (Rocscience). 30 Chapter 3: PIERRE 2 model description The new model, called PIERRE 2, utilizes the principles of the lump mass sphere solution, similar to the original CRSP model. However, several important modifications have been introduced. The following Chapter provides definitions and details for the key features of the PIERRE 2 model. The descriptions are given for the 2D model, with specific changes necessary for the 3D version discussed later. The model uses simple representations of the slope, boulders, trajectory and impact geometry combined with the equations describing the collinear impact of a sphere with a planar surface neglecting fragmentation, shown in Section 2.1.1. The model allows movement by flight or sliding. Although rolling motion is not explicitly represented in the model, it can be approximated by a number of small jumps, discussed in Section 3.4.4. The model has been calibrated primarily with datasets where the rocks have not had a significant amount of fragmentation, thus in this context, the assumption of constant mass is valid. Only impacts with the slope surface are considered. The calibration has focused on sites where tree impacts did not occur, or had a negligible effect on the trajectory. 3.1 Slope Slopes are defined by planar distance and elevation (x and z) coordinates for the 2D model. The x and z coordinates are connected using linear line segments. Slope material properties, discussed further in Section 3.4, must be assigned to each line segment (with the default of all segments being “Material 1”). A local slope angle, designated as β, is calculated for each line segment. In the 3D version of the model, slopes are defined by a regular grid (x, y) of elevation points (z), which is used to create a surface using bilinear interpolation. Material properties are 31 also defined using a regular grid of points which must have the exact dimensions as the elevation grid, where each grid centre corresponds to the index of a material to be applied to any impact within that cell. 3.2 Trajectory The model begins with the particle either in flight or sliding from the user defined source point. An initial velocity (vx, vz, ω) may also be specified at this point. The model calculates the position and velocity of the particle at each time step, which has been set at 0.02 s (real time). This time step size was determined in the calibration of a previous version of the 2D model, where the time step was made progressively smaller until the effect of further reducing the time step was not observable. This process was repeated for the 3D model as well, refer to Chapter 5. The model continues this process of updating the position and velocity until it detects a collision with the slope surface, see Section 3.3, at which point the particle velocity vector is passed to the functions for calculating the impact processes. At the end of an impact, the rebound velocity vector is used as the new initial velocity for the next trajectory. The trajectory algorithms describe the vast majority of the simulated particle’s movement. 3.2.1 Flight If the particle source is set above the ground surface, the trajectory will begin in flight. At the end of each impact the particle will return to flight with an initial planar velocity of (vxre, vzre). The flight is described using a classic ballistic parabola neglecting air resistance, so that gravity, g, acting in the z-direction, is the only external force acting on the particle. The position of the particle, (x, z) at the time from the start of the trajectory phase, t, can be found using Equation 3.1. 32 Equation 3.1: (𝑥, 𝑧)(𝑡) = (𝑣𝑥𝑟𝑒𝑡, 𝑣𝑧𝑟𝑒𝑡 − 12𝑔𝑡2) The angular velocity remains constant throughout each flight phase of the trajectory as well, but it has no effect on the position of the particle with the assumed conditions. The three-dimensional trajectory model also has gravity acting in the z-direction as the only external force. Thus the spatial coordinates can be described by the rebound velocity in the x, y, and z directions multiplied by t, as well as the gravitational acceleration in the z direction, as shown in Equation 3.2. Equation 3.2: (𝑥, 𝑦, 𝑧)(𝑡) = (𝑣𝑥𝑟𝑒𝑡, 𝑣𝑦𝑟𝑒𝑡, 𝑣𝑧𝑟𝑒𝑡 − 12𝑔𝑡2) 3.2.2 Sliding If the particle begins its trajectory at the ground surface it will start moving by sliding. If the source point is set by the user below the ground surface the program will automatically place the particle at the ground surface and begin in sliding as well. The driving force is the downslope (tangential) component of the particle’s weight, which is the cosine of 90 minus the local slope angle, β. The resisting force is the normal component of the particle’s weight, sin(90 – β), multiplied by the tangent of the interface friction angle, φ. From Newton’s second law of motion, the position of particle can be determined using the initial velocity of the block and the acceleration, a, due to the driving force minus the resisting force divided by the mass of the particle, mp. 33 Equation 3.3: (𝑥, 𝑧)(𝑡) = (𝑣𝑥𝑟𝑒𝑡 +12(sin(90 − 𝛽) 𝑎)𝑡2, 𝑣𝑧𝑟𝑒𝑡 − 12(cos(90 − 𝛽) 𝑎)𝑡2) 𝑎 = g(cos(90 − 𝛽) − sin(90 − 𝛽) tan𝜑) The particle will move by sliding until the normal force becomes zero due to a convex slope shape, the “ski jump” effect, or the velocity exceeds the arbitrary maximum sliding velocity of 5 m/s (at which point the program forces the particle into a flight trajectory), or the particle stops. In the 3D model, sliding is handled in an equivalent manner. To determine the direction of the vector normal to the local slope, the 3 nearest points to the particle are found, a triangular surface is defined, and the vector normal to that surface is calculated. From the normal vector the downslope and normal components of the weight can be calculated. 3.3 Collision Collisions are detected by checking when the calculated particle elevation is below the elevation of the ground surface. When this occurs the position of the particle is reverted back to its position at the previous time step so that the particle is just above the ground surface. The position of the particle at impact is compared to the program inputs to determine the slope angle and the material properties to be passed to the collision function. At this point the translational velocity vector is transformed from the global x-z coordinate system into a local coordinate system that is normal and tangential to the local slope, referred to here as the n-t coordinate system. At the end of the collision the rebound total energy head is compared to the stopping criterion to determine if the particle will return to flight or not, detailed in Section 3.7. 34 3.4 Slope Properties In the PIERRE 2 model, slopes are described using a surface roughness scale factor, the boulder-substrate interface friction angle, and the scaling point for the reference unit deformation energies for determining the normal and tangential restitution factors. These input parameters are discussed in detail in the following sections. 3.4.1 Roughness Irregularities in the slope and the surface of the particle are represented by a single roughness angle, designated θ. One of the key assumptions of this model is that the effects of non-collinear impacts of irregular particles can be represented within the model using an idealized collinear impact of a sphere on a plane, with the inclination of the plane being perturbed by a random angle. This is similar to the roughness used in the CRSP model, however the roughness used in PIERRE is not meant to represent the variability in the slope angle alone. The concept developed here is to use a single parameter to represent the net effect of the variation between the actual slope angle relative to the local slope input in the model, as well as the effect of the non-spherical shape and non-collinear impacts. A mathematical demonstration of how this simple angle can allow the prediction of more complex interactions is given in Section 3.8. However, the validity of this assumption is demonstrated through calibration of the model. At the point of contact with the slope a roughness angle θ is subtracted from the local slope angle, i.e. the slope is made less steep, as shown in Figure 3.1. The roughness angle is treated as a purely stochastic variable calculated for each impact. The local n-t coordinate system is rotated again by the roughness angle, designated n′- t′. The components of the incident 35 velocity vector (vn′, vt′, ω) relative to this twice rotated coordinate system are used for the subsequent calculations. Figure 3.1: Modification of the local slope angle by the roughness angle. In theory any distribution could be used to describe the roughness angle. A uniform distribution ranging from zero to some user input maximum angle was tested in the initial version of the model and provided reasonable initial results (Bourrier and Hungr, 2012). The equation for the uniformly distributed stochastic roughness is simply a sample from a uniform distribution in the range of 0 to 1, S, multiplied by the user input scaling value, θscale, shown in Equation 3.4. Equation 3.4: tan 𝜃 = 𝜃𝑠𝑐𝑎𝑙𝑒 ∗ 𝑆 Although there was initial success using the uniform distribution, detailed calibration efforts indicated the need for a more complex distribution in order to better represent the 36 extremes of bounce height and velocity while not over-predicting the runout. The Box-Mueller approximation for a normal distribution was used to create a truncated normal distribution of θ. This requires a mean and standard deviation, μ and σ respectively, to be specified along with the scaling value. The implementation of this requires two successive samples, Si and Sj, to be drawn from a uniform distribution in the range of 0 to 1, shown in Equation 3.5 (Scott, 2011). For numerical stability, values of θ that are negative are discarded and new random samples Si and Sj are drawn until a positive value for θ is obtained. Not doing so can result in a case where the perturbation to the slope can result in the trajectory being below the ground surface. Equation 3.5: tan 𝜃 = 𝜃𝑠𝑐𝑎𝑙𝑒 [𝜇 + 𝜎(cos 2𝜋𝑆𝑗)√−2(ln(𝑆𝑖))] Equation 3.5 is recommended for the calculation of the roughness, as detailed in Chapter 4. Observations of natural talus deposits show that larger boulders will tend to travel farther downslope than smaller ones, because larger boulders are able to override smaller perturbations on the slope surface. Thus, the surface roughness value is scaled to the boulder size. A linear relationship was assumed to describe the effective roughness angle, θeff, from a theoretical maximum roughness, b, corresponding to the roughness multiplier as the boulder diameter approaches 0, and decreasing at a rate of m multiplied by the boulder equivalent diameter D, shown by Equation 3.6. Equation 3.6: 𝜃𝑒𝑓𝑓 = 𝜃(𝑏 −𝑚𝐷) With a large enough boulder, this equation could give a negative roughness value. To deal with this, a minimum value for θeff was set to 0.3. Note: the simple linear Equation 3.6 is 37 probably valid only within a relatively narrow range of particle sizes. As a result, one should be careful about extending any calibration results to a very wide range of sizes. The calibration efforts completed during the present work show that the model is sufficiently accurate within a range of particle sizes used in the test programs, as given below In the 3D model two roughness angles are considered, one acting in the direction of the azimuth of the trajectory, designated θ, as before, the other acting perpendicular to the first, the transverse roughness, designated ψ. The roughness reduction is applied to the longitudinal roughness, θ, only. The value of the transverse roughness is a random number uniformly distributed between - ψscale and ψscale, shown in Chapter 5, or using a normal distribution (Equation 3.5) with a mean of 0, shown in Chapter 7. A comparison of the uniform and normal roughness distributions, generated from simulations using the same values for θscale and ψscale, is shown in Figure 3.2. The longitudinal roughness distribution is the same as the 2D roughness distribution. 38 Figure 3.2: Comparison of roughness angle distributions using a uniform distribution for the longitudinal and tangential roughness, a and b, and using a normal distribution, c (truncated) and d. The roughness is the primary stochastic element in the model. Representing much of the variability in rock fall trajectories using this single stochastic parameter was a deliberate choice to constrain the problem. More stochastic elements would involve a greater amount of calibration effort in order to define the necessary distribution properties, and increase the overall uncertainty in the model with the increase in the number of parameters. It is for this reason that the interface friction angle has been set as a constant, and the reference deformation energies are calculated in a deterministic manner, described in the following sections. 39 3.4.2 Interface friction angle The dynamic coefficient of friction between the slope surface and the particle, referred to as the interface friction angle, φ, is used in the calculation of sliding motion and to determine if an impact has contained or uncontained slip 3.4.3 Reference normal deformation energy Instead of the commonly used “coefficient of restitution”, the term “restitution factor” is adopted here to reflect the fact this is not a constant material property (refer to Section 2.2). The definition of the ratio of rebound over incident velocity is used for the restitution factors. As the model does not consider fragmentation of the boulder during impact, this is equivalent to the ratio of rebound over incident momentum. Using the principles outlined in Section 2.1, a novel method for calculating restitution factors has been implemented in this model. The impact mechanics principles outlined in Section 2.1 provide the theoretical basis for the development of the deformation energy relationship. Consider a rigid spherical boulder with mass m, diameter D, and incident velocity vector vin, having a normal impact with a surface with a constant indentation resistance IR (strength like parameter described by Pichler et al., 2005). The constant indentation resistance proposed by Pichler et al. is conceptually similar to the constant plastic pressure of the Meyer Contact Law. Fragmentation does not occur in the impact, and the force of gravity applied on the boulder over the impact period is relatively small and can be neglected. A conceptual force-time and force-penetration history can be derived from these inputs. The surface applies a force to the boulder that is equal to the impact resistance multiplied by the contact area. The time integral of the resisting force and the contact duration gives the resisting impulse. The boulder’s initial momentum is reduced by the resisting impulse. The maximum penetration occurs when the resisting impulse is equal to the initial momentum, i.e. 40 the boulder’s velocity is instantaneously zero. At this point, elastic strains stored in the surface are recovered, and the surface exerts a rebound impulse on the boulder that is a function of the elastic strain recovered and the duration of the rebound phase. The force-penetration plot for this conceptual impact model is shown in Figure 3.3. The integral of the force-penetration curve allows for the identification of two regions. The area bounded by the upper curve and the rebound curve represents the plastic work (permanent deformations) that occurred. The area bounded by the rebound curve and the vertical line from the maximum force point represents the elastic work done during the impact. 41 Figure 3.3: Conceptual diagram showing the force-displacement diagrams during a) the start of the incident phase, b) the point of maximum force (Fmax) and penetration (αmax), and c) the end of the restitution phase. The shaded area represents the total plastic work for the impact. 42 The concave shape of the force-penetration diagram is a reflection of the assumed constant yield strength of the impacted surface, so that the total force is a product of the contact area, which for a sphere, scales in proportion to D2. The difference between the incident and rebound momentum is a product of the plastic work done. Plastic work is equivalent to the change in kinetic energy, which is 1/2mΔv2. The mass is proportional to D3. Thus, if we assume all of the particles have approximately equal densities, the change in momentum during an impact is a function of the incident kinetic energy, which is proportional to v2D3. A larger-diameter particle will have a greater contact area for a given penetration depth relative to a smaller diameter particle. This concept is similar to the “bullet effect” discussed in Section 2.1.4 (Buzzi et al., 2011). To represent this, the incident kinetic energy is distributed over the contact area, which is proportional to D2. The term “unit deformation energy”, Ed, is used to represent the incident normal kinetic energy distributed over the contact area, which has characteristic dimensions of vn2D. A hyperbolic relationship is proposed as a simplification to these theoretical constructs. The hyperbola is scaled using the unit deformation energy, Ed. In an infinitesimally light impact the collision can be considered to be fully elastic, so that Ed ≈ 0 and kn ≈ 1. Likewise, an infinitely heavy impact results in a fully plastic condition, where Ed approaches infinity and kn ≈ 0. The validity of the hyperbolic restitution hypothesis is difficult to determine theoretically, however the calibrations shown in Chapters 4, 5, and 6 demonstrate its validity. The program calculates the restitution factor normal to the impact plane rotated by the roughness angle, kn′, for each impact from a user input reference normal deformation energy, En, and the corresponding restitution factor for the input En, shown in Equation 3.7. The work 43 detailed in the following chapters uses reference unit deformation energy values corresponding to kn′ = 0.5, which is designated En,50, shown in Figure 3.4. Equation 3.7: 𝑘𝑛′ =𝐸𝑛,50𝐸𝑑 + 𝐸𝑛,50 Figure 3.4: Hyperbolic relationship between deformation energy and normal restitution, with the scaling point highlighted. 3.4.4 Reference tangential deformation energy The tangential restitution factor is determined in an equivalent manner to the normal restitution factor. The assumption here is that the non-frictional tangential momentum losses, characterized by kt′, are also scaled by the value of the unit deformation energy. These losses are meant to account for cratering and displacement of material from the impact site. Note that the reference unit deformation energies are not the same for normal and tangential velocities, i.e. 44 En,50 ≠ Et,50. However, for both kn′ and kt′ the impact unit deformation energy value, Ed, is the same. Equation 3.8: 𝑘𝑡′ =𝐸𝑡,50𝐸𝑑 + 𝐸𝑡,50 The 3D model applies the tangential deformation energy to the rebound velocity vector for both the longitudinal and transverse directions. The hyperbolic scaling of the restitution factors allows the model to approximate rolling motion. When the modeled boulder has a shallow incident angle on a low roughness surface it will tend to move in a series of short, low altitude bounces with the angular velocity synchronized with the tangential velocity. The relatively low normal and tangential restitution factors in this scenario take the place of rolling friction. An advantage of this is that the roughness angle is still applied at each modeled impact, which allows the modeled particle to resume a higher trajectory, similar to a real rolling boulder encountering a small change in topography that results in it returning to a bouncing motion. 3.5 Boulders Boulders are described by their mass and a shape factor. The mass of the boulder is used to calculate the equivalent diameter for a sphere with a constant density of 2.65 Mg/m3. This is used as the characteristic diameter for detecting impacts with the slope surface and the calculation of restitution factors. The program assumes the particle is a sphere, but natural rocks have a wide variety of shapes. One consequence of this is that the distance from the impact point to the centre of mass will likely not be the equivalent spherical radius, R. A sketch of this is shown in Figure 3.5 a), 45 where the distance between the centroid and the contact point, d, is greater than R. The scenario shown in this sketch is not a collinear impact as the centre of mass does not lie on the line normal to the contact point. The sum of the moments about the contact point are a function of all of the incident vector components as well as the moment arms Δn and Δt, illustrated in Figure 3.5 b. The magnitude of Δn and Δt will be a function of the orientation of the rock. Furthermore, the linear and angular momentum balances cannot be decoupled. This can be contrasted to the case of a collinear impact where the sum of the moments is simply a function of the incident tangential or normal velocity (depending on if it is a contained or uncontained slip impact), the angular velocity and the radius. A collinear impact is much simpler to solve, however, it was found in the calibration of the model this representation has a tendency to over-predict the angular velocity. Figure 3.5: Hypothetical impact of an irregularly shaped natural rock, a) contrasting equivalent spherical radius to the actual distance to the contact point, and b) components contributing to the moment about the contact point. 46 Maintaining the mathematical simplicity of the lump mass model assuming collinear impacts is desirable. To address the over-prediction of angular velocity, a shape factor is hypothesized based on the differences in distance from the centroid to the contact point, utilizing the Goldsmith equations for collinear impacts. It is proposed that varying the radius used in the calculations of rebound tangential and angular velocity from one impact to the next can be used to represent the effects of non-collinearity on the angular velocity. Figure 3.6 b) shows the equivalent spherical radius calculated from the mass of the irregular boulder shown in a), as well as the effective radius used in the impact calculation to represent the rotational effects of the non-collinear impact. Figure 3.6: Representation of non-collinear impact using an effective radius. This is implemented in the model by the user specifying a maximum radius Rmax, called the shape factor, in the boulder table. This is used to create a range of possible radius modifications. At each impact, a uniformly distributed random number, S, is multiplied by the range of the radius modification, shown in Equation 3.9. 47 Equation 3.9: 𝑅𝑒𝑓𝑓 = 𝑅𝑚𝑖𝑛 + (𝑅𝑚𝑎𝑥 − 𝑅𝑚𝑖𝑛)𝑆, 𝑤𝑖𝑡ℎ 𝑆 ∈ [0, 1] The assumption is that the non-collinear impact in Figure 3.6 a) can be represented by a collinear impact. The effect of the moments applied by the incident normal and tangential velocity acting through the moment arms Δn and Δt is approximated by an effective radius. A theoretical proof to show the effective radius is a good approximation of a real non-collinear impact is difficult to demonstrate. Instead, the validity of this simplification will be demonstrated through the calibration process described in Chapter 4. 3.6 Impact Model The impact model used is a modified version of the stereomechanical impact model described by Equation 2.2, Equation 2.3, and Equation 2.4. The process described in Section 2.1.1 considers a single coefficient of restitution acting in the normal direction. The modification described here uses both the normal and tangential restitution factors determined from the hyperbolic relationships. The inclusion of a tangential restitution factor is meant to represent the non-frictional energy losses due to plastic deformations and displacement of material at an impact site. Following the procedure outlined by Goldsmith (1960), the friction limit angle, γ, is calculated as the inclination of the incident impulse vector. Taking the tangential restitution factor into Equation 2.2, considering all calculations are done relative to the n′- t′ plane, results in Equation 3.10. Equation 3.10: tan 𝛾 = 2(𝑘𝑡′𝑣𝑡′𝑖𝑛 + 𝑅𝑒𝑓𝑓𝑘𝑡′𝜔𝑖𝑛)7𝑣𝑛′𝑖𝑛(1 + 𝑘𝑛′) 48 As with the original Goldsmith equations, the friction limit angle, γ, is compared to the interface friction angle, φ, with φ > γ being a contained slip impact. The system of equations describing a contained slip impact using the modified equations is shown in Equation 3.11. Equation 3.11: |𝑣𝑡′𝑟𝑒𝑣𝑛′𝑟𝑒𝑅𝜔𝑟𝑒| = ||57𝑘𝑡′ 0 −27𝑘𝑡′0 −𝑘𝑛′ 0−57𝑘𝑡′ 027𝑘𝑡′|| × |𝑣𝑡′𝑖𝑛𝑣𝑛′𝑖𝑛𝑅𝑒𝑓𝑓𝜔𝑖𝑛| If ϕ < γ the impact has uncontained slip and the interface friction angle must be considered in the calculation of the rebound velocities, shown in Equation 3.12. Equation 3.12: |𝑣𝑡′𝑟𝑒𝑣𝑛′𝑟𝑒𝑅𝜔𝑟𝑒| = |𝑘𝑡′ − tan𝜑 (1 + 𝑘𝑛′) 00 −𝑘𝑛′ 00 −52𝑡𝑎𝑛𝜑(1 + 𝑘𝑛′) 𝑘𝑡′| × |𝑣𝑡′𝑖𝑛𝑣𝑛′𝑖𝑛𝑅𝑒𝑓𝑓𝜔𝑖𝑛| The difference between contained and uncontained slip impacts for a stereomechanical impact was described in Section 2.1.1. The hyperbolic restitution coefficients introduced here add another level of complexity to the simple relationship presented graphically in Figure 2.1. The normal and tangential restitution factors depend on the incident normal velocity and the diameter of the impacting particle. Using the mass as a representation of the diameter, the rebound linear kinetic energy can be plotted as a function of the incident normal velocity and boulder mass with the interface friction angle, tangential velocity, and the ratio of radius multiplied by rotational velocity divided by the tangential velocity held constant. Figure 3.7 shows the interaction between vnin and boulder mass. The discontinuity in the surface shows the transition from uncontained to contained slip conditions, with the uncontained slip conditions 49 having a higher energetic restitution. This discontinuity in the surface is a result of the resolution of the plot and the sensitivity of the normal restitution to the mass and velocity. If this was plotted with a fine enough resolution the transition would become smooth. This simple set of equations utilizing the assumption that the impacts are collinear is able to represent a wide range of possible energy restitution situations when the hyperbolic restitution factors are introduced. Figure 3.7: Relationship between the incident and rebound kinetic energy as a function of vnin (m/s) and log mass (mass in kg). 3.7 Stopping Criterion Models generally use an arbitrary minimum velocity to determine when the particle has stopped (c.f. Stevens, 1998, Lan et al., 2007). Using the energy and work principles outlined in 50 the development of the deformation energy, it can be shown that larger rocks require a greater amount of rebound energy to keep them moving following an impact. The rock must accumulate sufficient elastic work in the slope in the normal direction to overcome the normal component of the weight of the rock if it is to become airborne again at the end of the restitution phase. If a rock fails to become airborne, it may still move by rolling or sliding. Plastic work in the tangential direction resulting from friction and other losses during the impact will reduce the tangential and angular kinetic energy. If the net incident tangential and angular kinetic energy is greater than the net plastic work done over the contact interval, the rock will continue to move. As the model treats impacts as instantaneous, it is not possible to solve these force balances directly, as all of these processes would have to be integrated over time. However, from this conceptual framework, the resisting forces acting over the contact area can all be related to the size of the block, the weight for the normal direction, and the contact area in the tangential direction. The stopping criteria used in the model compares the total energy head, linear plus angular kinetic energy divided by the block weight, to a fraction of the block dimension, Hstop (measured in m). The rebound velocities calculated for the end of the restitution phase of an impact are used to find the total energy head, which is equivalent to the vertical distance the block would be above the impact point if all the rebound kinetic energy were converted to potential energy. Hstop is calculated from the particle radius multiplied by a dimensionless constant. If the rebound total energy head is less than Hstop the model will stop the trajectory, shown in Equation 3.13. 51 Equation 3.13: 𝒗𝑟𝑒 = (𝑣𝑡𝑟𝑒𝑣𝑛𝑟𝑒𝜔𝑟𝑒) 𝑓𝑜𝑟 𝐸𝐾𝑟𝑒 + 𝐸𝜔𝑟𝑒𝑚𝑔 ≥ 𝐻𝑠𝑡𝑜𝑝 𝒗𝑟𝑒 = (000) 𝑓𝑜𝑟 𝐸𝐾𝑟𝑒 + 𝐸𝜔𝑟𝑒𝑚𝑔 < 𝐻𝑠𝑡𝑜𝑝 If the rebound energy head does not exceed some minimum value, the boulder is unlikely to be able to bounce, roll or slide out of its impact crater. The value of the constant used in the calculation of Hstop was determined empirically by calibration, detailed in Chapter 4. 3.8 Apparent Restitution Values The restitution factors calculated using Equation 3.7 and Equation 3.8 are always greater than zero and less than one. As discussed in Section 2.2, the observed normal restitution values are dependent on the incidence angle of the impact, and can be greater than 1. This is represented in the model with the rotation of the incident and rebound vectors using the roughness angle. As the values of kn′ and kt′ are calculated in the n′- t′ plane, the values for kn and kt, which are relative to n - t plane, will be different. The values of kn and kt are referred to as the apparent restitution factors, while kn′ and kt′ are referred to as the true restitution factors. The apparent restitution factors are used when comparing the model restitution factors to those from experiments, such as those detailed in Section 2.2. The normal and tangential components of the apparent incident vectors in the true incident vectors are shown in Equation 3.14. 52 Equation 3.14: 𝑣𝑡′𝑖𝑛 = 𝑣𝑡𝑖𝑛 cos 𝜃 + 𝑣𝑛𝑖𝑛 sin 𝜃 𝑣𝑛′𝑖𝑛 = −𝑣𝑡𝑖𝑛 sin 𝜃 + 𝑣𝑛𝑖𝑛 cos 𝜃 The rotation of the local coordinate system results in the incident normal velocity having an apparent effect on the tangential rebound velocity, and the incident tangential velocity having an apparent effect on the normal rebound velocity. The rotation of the incident and rebound vectors can be shown in the rotation matrix Rθ, Equation 3.15. Equation 3.15: 𝑹𝜽 = |cos 𝜃 sin 𝜃 0− sin 𝜃 cos 𝜃 00 0 1| Since the output from PIERRE is given relative n-t plane, but the calculations are performed in the n′- t′ plane, both sides of the modified Goldsmith system of equations, designated G, must be multiplied by the rotation matrix Rθ, shown in Equation 3.16. Equation 3.16: 𝑹𝜽𝒗𝒓𝒆 = 𝑮𝑹𝜽𝒗𝒊𝒏 Finally, to isolate the rebound velocities vector, both sides must be multiplied by the inverse of Rθ. Carrying through this calculation using the matrix for a contained slip impact, where the 3 x 3 matrix is the result of Rθ -1G Rθ, the matrix shown in Equation 3.17 is found. Equation 3.17: |𝑣𝑡𝑟𝑒𝑣𝑛𝑟𝑒𝑅𝜔𝑟𝑒| = ||57𝑘𝑡 cos2 𝜃 − 𝑘𝑛 sin2 𝜃57𝑘𝑡 cos 𝜃 sin 𝜃 +𝑘𝑛 sin 𝜃 cos 𝜃 −27𝑘𝑡 cos 𝜃57𝑘𝑡 sin 𝜃 cos 𝜃 + 𝑘𝑛 sin 𝜃 cos 𝜃57𝑘𝑡 sin2 𝜃 − 𝑘𝑛 cos2 𝜃 −27𝑘𝑡 sin 𝜃−57𝑘𝑡 cos 𝜃 −57𝑘𝑡sin 𝜃27𝑘𝑡|| × |𝑣𝑡𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝑒𝑓𝑓𝜔𝑖𝑛| 53 This shows how the roughness angle changes the simple Goldsmith equations, thus allowing more complex behaviors to be modeled. This can be compared to the fully stochastic matrices presented by Bourrier (2008). The coefficient matrices of the two methods are not directly comparable, as the matrix shown in Equation 3.17 does not consider covariance. However, the inclusion of the stochastic variables θ and Reff result in all 9 elements of the coefficient matrix having a stochastic element, which is consistent with Bourrier’s observations. 3.9 Program Outputs The program records the local slope angle, roughness angle, incident and rebound velocity components, coefficients of restitution in both the “apparent” n-t directions and the “true” n′- t′ directions and the maximum jump height for each impact. The maximum jump height is taken as the maximum vertical distance from the center of mass of the particle to the ground surface. These results are stored in the TR.dat datafile. The trajectory is stopped when the block’s total energy head drops below the stopping criterion discussed previously, or when the block passes the model boundary. At the end of each trajectory the start and end coordinates, as well as the total duration of the simulated trajectory, are recorded. This information is stored in the DIS.dat datafile. The user may also specify an observation point at a point on the x axis. At the observation point the jump height and velocity will be recorded for each boulder that passes the designated point on the x axis and recorded in the OBS.dat datafile. 54 Chapter 4: Calibration and validation of the PIERRE 2 stochastic rock fall dynamics simulation program 4.1 Introduction Rock fall model calibration often focuses on matching runout distances (c.f. Guzzetti et al., 2002, Lan et al., 2007), which is important, but it only covers the reach term of the hazard equation (Equation 1.1). Being able to accurately represent the particle kinematics is needed to properly evaluate the event intensity, thus it is critical to determine rock fall hazard for an area. These kinematic properties are also required for the design of protection methods, such as catchment areas or barriers. Catchments and barriers must be adequately sized so that incoming rocks are intercepted, and barriers, which either deflect or contain rock fall, must be able to withstand or attenuate the impact energy (Turner and Schuster, 2012). A realistic prediction of the kinematics of the design boulder is necessary to accurately estimate the impact energy at a barrier. The rotational motion of a rock is also frequently overlooked in model calibration and design practice; however, a significant portion of the total energy of a falling rock is due to the rotation (c.f. Chau et al., 2002). The challenge in modeling rockfall events is not in the flight portion of the trajectory, which can be accurately represented by a ballistic trajectory, but in the representation of the impacts. A great deal of research has been done in the field of impact mechanics, and the physical processes during an impact have been quite well described (Goldsmith, 1960, Stronge, 2000). The material properties for a falling rock and the impact site, such as the Young’s modulus or strength properties, are variable, and cannot be accurately characterized in a practical application. This limits the usefulness of the established impact mechanics relationships, 55 described in Section 2.1, for application to rock fall modeling. However, this research does provide a conceptual framework for the modeling of rockfall impacts, which was developed in Chapter 3. This chapter describes the calibration process of the PIERRE 2 model described in Chapter 3. The calibration has focused on sites which have detailed information on the rock kinematics, including one site with detailed measurements of the rock angular velocity. Details of the calibration process to find model inputs for both firm talus and exposed rock slopes are presented, along with general model inputs for these broad slope types. The talus slope input parameters are applied to another site for validation purposes. The results are also compared to design standards to demonstrate the applicability of the PIERRE 2 program as a design tool. 4.2 Data 4.2.1 Calibration The initial calibration of the 2D model was done using data from detailed observations of multiple rock fall parameters during a real-size rock fall experiment performed by the IRSTEA Institute (formerly CEMAGREF) of Grenoble, near Vaujany in the Savoy Alps. Details of the experiments were described by Dorren et al. (2005). The initial 2D calibration was done concurrently with the 3D calibration presented in Chapter 5. A total of 100 blocks were released from the same spot. Figure 4.1 a) shows the distribution of block masses. They range from 1.5 Mg to 3.5 Mg, which corresponds to equivalent block diameters of about 1 to 1.4 m at a rock density of 2650 kg/m3, shown in Figure 4.1 b). 56 Figure 4.1: Distribution of boulder sizes used in the 2003 Vaujany experiments a) on the basis of mass and b) diameter. The total runout distances, and the velocity and jump height of each block, as recorded by cameras at Screens 1 and 2, and the translational kinetic energy obtained from the measured block mass and velocity, were used in the calibration (CEMAGREF, 2003). The slope profile and locations of Screens 1 and 2 are shown in Figure 4.2. The slope surface had an average angle of 32° and consisted of talus of gravel to cobble-sized particles, with patches of small to medium boulders, all of which is covered by a thin layer of forest litter. A view of the test slope looking down from the boulder release location is shown in Figure 4.3. 57 Figure 4.2: Vaujany test site slope profile and the locations of Screen 1, Screen 2 and the forest road. Figure 4.3: Vaujany test site, looking downslope from the release location. 58 A second calibration was done using data from the Preg Quarry site in Austria. Four slope sections along one bench in the quarry were used in this testing. The sections ranged from slightly convex to slightly concave, shown in Figure 4.4. The surfaces involved were rough, steep faces created by uncontrolled blasting of strong rock and flat, horizontal quarry floors of highly compacted rock debris, see Figure 4.5. The experimental protocol and additional details on the site were presented by Preh et al. (2015), included as Chapter 6 in this thesis. Figure 4.4: Preg Quarry test site profiles, a) Section 1, b) Section 2, c) Section 3, and d) Section 4. 59 Figure 4.5: Preg Quarry test site (photo A. Preh). Twenty-five rocks were released from the top of each of the 4 slopes, for a total of 100 rocks. First, the volume of the boulder to be thrown was estimated from measurements of the three dominant boulder axes, and the volume was used to calculate the mass. The mass of the boulders ranged from 0.14 Mg up to 8.13 Mg, shown in Figure 4.6. An excavator was used to move the boulders gently over the edge of the slope. At the end of each experiment, the stopping point of the boulders was measured. The impact distances were estimated by video analyses. After 5 throws, the boulders on the bench were removed by a wheel-loader in order to prevent boulder interaction during the following throws. Rock fall trajectories were filmed using two digital cameras, installed at the top and the toe of the slope. Additional details of this test site are available in Preh et al. (2015), included as Chapter 6 in this thesis. 60 Figure 4.6: Distribution of boulder sizes used in the Preg experiments a) Section 1 mass, b) Section 1 diameter, c) Section 2 mass, d) Section 2 diameter, 3) Section 3 mass, f) Section 3 diameter, g) Section 4 mass, h) Section 4 diameter. The third calibration dataset is from the Ehime prefecture in Japan. The experimental procedure and test site for their 2003 rockfall study were described by Ushiro et al. (2006). In 61 the testing both cast concrete blocks and natural rocks were launched down a relatively open slope using a hydraulic excavator. The trajectories of the blocks were tracked using high speed cameras, and the blocks were instrumented with triaxial accelerometers and angular rate sensors. This allowed for tracking the position, and determining the linear and angular velocity of all the test blocks. The test site is located on a 42 m high rock and talus slope covered by thin, exposed soil with sparse grass and shrubs in the test track. The upper 26 m of the slope is thinly mantled sandstone and mudstone with an average face angle of 44°. Below this is a 16 m high firm talus slope at angle of 35°, shown in Figure 4.7 (Ushiro et al., 2006). Figure 4.7: Ehime test site profile. 62 A total of 23 trajectories were recorded for the Ehime tests. A total of 4 cast concrete spheres with masses of 194 kg and diameters of 0.52 m, 7 concrete cubes with masses of 520 kg and diameters of 0.72 m, and 12 natural rocks, with masses ranging from 120 to 2060 kg, shown in Figure 4.8, were launched. Figure 4.8: Distribution of natural rock sizes used in the Ehime tests a) mass, and b) diameter. 4.2.2 Validation The fourth data set was from Tornado Mountain in British Columbia, which was described in detail by Wyllie (2014). In 2004, two separate rockfalls originating on the rock face travelled distances of approximately 340 m vertically and 610 m horizontally and came to a stop on a horizontal bench excavated for the railway. The source area was identified as a 50 m high rock face in very strong, blocky limestone. The runout path was along a uniform talus slope of gravel to small boulder size fragments, similar in description to the Vaujany site. The slope ranges from 35° on the upper slope to 22° on the lower slope, with no significant irregularities. 63 The slope is vegetated with sparse conifers with diameters up to 500 mm. There was no evidence of recent rockfalls on the talus slope prior to the 2004 event. The impact points on the slope were identified from a field survey, and the co-ordinates were recorded with a GPS unit and a laser rangefinder. Broken trees were noted at certain locations along the trajectories, and it was assumed that the trees broke at the elevation that they were impacted. Those elevations were used to reconstruct the trajectory and velocity of the rocks between the impact points before and after, by fitting a ballistic parabola to the three points. The key assumption here is that breaking the tree did not significantly slow or deflect the rock, or otherwise affect the trajectory. The analysis presented by Wyllie (2014) used this assumption, with the justification that the trees were quite small relative to the boulders. That assumption is carried over in this work. The slope profile and the locations of the tree impacts are shown in Figure 4.9. Figure 4.9: Location of tree impacts along the profile used to estimate jump heights and velocities for boulders A and B from Tornado Mountain. 64 Fragments of both boulders were observed along their runout paths. The final masses for the boulders were about 3,750 kg (maximum dimension 1.6 m), and 5,600 kg (maximum dimension 2.5 m) respectively (Wyllie, 2014). The final stopping positions of the two blocks, and jump heights and velocities, where available, were used for the validation. The final validation dataset was from experiments conducted by the British Columbia Ministry of Transportation and Highways at Nicolum Quarry, 200 km east of Vancouver (49°22’49”N, 121°22’31”W). The quarry is developed in quartz diorite belonging to the Coast Plutonic Complex. The rock is extremely strong and the rock mass is of excellent quality. The slope used for the testing is approximately 80 m high with an average slope of 52°, shown in Figure 4.10. Figure 4.10: Nicolum quarry test site profile. 65 The falling rocks were filmed by two analog video cameras. A total of 34 trajectories were recorded for the Nicolum quarry tests. The blocks used ranged in mass from 0.03 Mg to 3.04 Mg (equivalent diameters of 0.28 m to 1.3 m), shown in Figure 4.11. Figure 4.11: Distribution of boulder sizes used in the Nicolum tests on the basis of a) mass, and b) diameter. The trajectory of each boulder was traced on the videos with reference to survey marks on the slope, so that the geometric position and timing of each impact could be recorded. This allows estimation of incident and rebound velocities using ballistic calculations. Runout of the boulders was not measured at Nicolum because there was an experimental rock fence at the foot of the slope which stopped most of the rocks. The site at the time of this testing is shown in Figure 4.12. The fence impact parameters are not reported in this thesis. 66 Figure 4.12: Nicolum quarry test site with experimental fence located at the toe of the slope (Photo O. Hungr). 4.3 Methodology The five sites used in this study all have different data available, which meant that the calibration or validation methodology had to be adjusted for the data available to each site. The general procedure in this process was to group slopes with similar properties together. Vaujany, Ehime and Tornado Mountain were used as the “firm talus/weak rock” group, and the Preg and Nicolum quarries made up the “hard, exposed rock” group. Parametric studies using select data sets, one from each group, were conducted to establish the model response to a wide range of inputs, and to establish narrower zones of inputs for more detailed calibrations. The detailed calibrations were performed using a trial and error method, with a goal of providing a reasonable 67 prediction of the complete observed data sets for each site within each group. Details of the specific calibration and validation methods are included in the following sections. 4.3.1 Parametric studies The Vaujany and Preg Quarry sites were selected for the initial parametric studies for an open talus slope and exposed hard rock slope, respectively. Preliminary testing of the model found En,50 and θscale were the most sensitive, thus they were used as the variables, with constant Et,50 and φ, in the parametric study. Using two variables allows for the objective function to be plotted as a surface, with the minimum of the surface corresponding to the best fit. The following objective functions were examined for the comparison of the simulations to the experimental data:  Mean;  Root mean squared error;  Mean bias error;  Two sample Kolmogorov-Smirnov test; and  Maximum values. The mean was judged to be too coarse of a comparison, as many different distributions can end up giving the same mean value. The root mean squared error, RMSE, gives an indication of the overall fit of the simulated distribution compared to the experimental distribution, and will always result in a value greater than zero, with lower values indicating a better fit. The mean bias error, MBE, is also a measure of the overall fit, and can be positive or negative, with positive values indicating the model is under-predicting the result on average, and vice-versa. The Kolmogorov-Smirnov (KS) test gives a test statistic, which is a measure of the probability that the two samples are from the same population, with a value near zero indicating a high 68 probability, and a value near 1 indicating a low probability. The KS test also gives a p-value between 0 and 1, which is the significance level of the test statistic (GraphPad, 2014). The maximum values were also considered, as these are important to consider for design. Initially the Vaujany model was run with a wide range of inputs to find the model response to a very broad range of parameters, and to select the appropriate means of comparison. The runout was used for this initial phase. The RMSE and MBE are similar tests, so it was decided that one or the other would be used. An overall good fit was desired, so the ablity to differentiate over from under-prediction using the MBE provided no specific advantage. The RMSE provided easier interpretation of the parametric plots as it has a minimum value of zero, thus it was examined further. To select the objective functions for the following parametric studies, plots of the RMSE, observed minus simulated maximum values and a two-sample KS test were compared, shown in Figure 4.13 and Figure 4.14, respectively. Figure 4.13: Initial parametric surface plots for runout distance. a) RMSE of the runout distribution and b) the absolute value of the residual of the maximum runout. 69 Figure 4.14: Initial parametric surface plots for runout distance. a) KS test statistic for the runout distribution and b) p-value for the KS test. The KS test statistic agrees with the RMSE generally; however the p-value gives low significance for the tests (all tests had less than 2% significance). More gradation in the results can be observed using the RMSE plot as well, whereas with the KS test, the test statistic increases rapidly moving away from the best fit area and stays near 1 for the majority of the parameter space. The RMSE and residual maximum value were selected as the objective functions for the parametric studies. The distributions for the simulated results were calculated using the same break points as the ones provided in the data reports for the experiments (CEMAGREF, 2003, Preh et al., 2015). The comparisons between the observed and simulated datasets were made on the basis of the cumulative distributions. The reasoning for this was to focus the comparison on optimizing the fit for the entire distribution. The root mean squared error for the predicted distribution, Ŷ, and the observed distribution, Y, divided into n bins is calculated using Equation 4.1. 70 Equation 4.1: 𝑅𝑀𝑆𝐸 = √1𝑛∑(𝑌?̂? − 𝑌𝑖)2𝑛𝑖=1 The maximum prediction for the runout, jump height, velocities or first impact point were also compared, depending on the simulation being tested, as these are the values that will likely govern the design of rock fall protection system, and will represent the highest consequence events in a risk assessment. To simplify the following graphical representations, the absolute value of the residual of the maximum prediction, Ŷmax, compared to the maximum observation, Ymax, was calculated using Equation 4.2. Equation 4.2: 𝑌𝑟𝑒𝑠𝑖𝑑 = 𝑎𝑏𝑠(?̂?𝑚𝑎𝑥 − 𝑌𝑚𝑎𝑥) 4.3.2 Velocity versus elevation data The velocity data available for the Ehime and Nicolum tests were obtained at specific elevations on the slope corresponding to impacts. Due to the random nature of the process, the observed points were distributed irregularly on the profile, thus attempting to directly compare the observation data to the simulation was not feasible. However, examination of the data shows that there appears to be some relationship between the elevation (i.e. distance from the release point), and the velocity. The relationship between the elevation and velocity could not be defined with a simple relationship, such as a linear or polynomial functions fit to the data. To address this, a non-parametric local regression method was implemented using the “loess” regression function in R (R Core Team, 2013). The primary advantage of this method is a relationship can be fit to the data without having to make broad assumptions about the form of the relationship. 71 The theory behind the method is described in detail by Cleveland et al. (1992), and is briefly summarized here for the case with one numeric predictor variable. The curve is fit to a certain point using the observed points in its neighbourhood, weighted by the distance of the observed point to the point being fit. The proportion of the total points considered in the fitting is set by the span, which for the general case is a number greater than 0 and less than or equal to 1, with 1 being all the points in the observation dataset considered. The span may be set to greater than 1, however that case was not considered in this calibration. As the span is increased, the fit curve will become smoother. The weight that an observation point within the span is given to the fit curve at any location, T(u; t), is inversely proportional to the distance of the observation point to the point being fit, u, over the length of the span, t, shown in Equation 4.3 (Cleveland et al., 1992). Equation 4.3: 𝑇(𝑢; 𝑡) = {(1 − (𝑢𝑡)3)3, 𝑓𝑜𝑟 0 ≤ 𝑢 < 𝑡0 , 𝑓𝑜𝑟 𝑢 ≥ 𝑡 Essentially, it was assumed that the velocity is related to the distance from the release point. While the form of that relationship is unknown, the observation points nearer to a certain location on the slope have a greater influence on the relationship. This was used as an objective method to compare the simulation results to the observed data for Ehime and Nicolum. First a loess regression curve was fit to the observed data, and the residuals of the observed points compared to the loess curve were found. Next the residuals of the simulated points compared to the loess curve fit to the observed data were found. If the residuals from the simulated points have a mean near zero, that is an indication that the simulation is unbiased. The two residual distributions can also be compared to see if they have a 72 similar range. Finally, a loess curve is fit to the simulation data to see if the two curves fall near one another for the entire elevation profile, or if there is a tendency to over or under-predict velocities for certain elevation bands. Further specific details are provided in the Ehime and Nicolum calibrations. The calibration and validation process has been performed iteratively. For example, the size dependent roughness scaling was calibrated using the Preg quarry data, and the shape factor was calibrated using the Ehime data, and the parameters identified were then applied to the Vaujany simulation. This iterative method was required to achieve the goal of general parameter sets applicable to broadly similar slopes. 4.4 Initial Model Calibration 4.4.1 Vaujany – open talus slope As described in the previous section, a parametric study was conducted by varying the reference normal deformation energy and roughness angle. The losses due to the tangential coefficient of restitution are, in general, much smaller than those due to the normal coefficient of restitution. The observed values also generally fall in a smaller range (Turner and Schuster, 2012), which is why variability in the reference tangential unit impact energy was not considered in the initial parametric study. A tangential deformation energy value for 95 % energy conservation, Et,95 = 50 m3/s2 was used for all trials. Variation in the contact friction angle was also not considered, and a constant value of 30°, reported by Bourrier et al. (2012), was used. An initial range of normal deformation energy values for 50 % energy conservation, En,50 = 1 to 50 m3/s2 and a roughness scale, θscale = 0.1 to 1.0 was used to establish the response of the model over a wide range of inputs. For simplicity, a uniform distribution was used to represent the roughness angle distribution, Equation 3.4, so that only θscale is needed to describe the 73 distribution. The initial run looked at finding the minimum RMSE for the runout distance distribution to establish an optimum zone to focus on for more detailed analysis. The best results for the runout distribution were found at the lower ranges of the normal reference deformation energy, and high roughness. A comparison of the parametric runout plots for the initial range of inputs and then the refined search are shown in Figure 4.15. Figure 4.15: Parametric surface plots for the runout distribution RMSE, a) for the initial wide range of inputs, and b) for the refined input parameters. A more refined parametric search was conducted for En,50 from 5 to 20 m3/s2 and θscale from 0.7 to 1.2. In this range, the RMSE for the distributions of total runout distance, jump height and translational velocity at Screens 1 and 2 were calculated for all parameter sets tested. Since design of catchment or protective structures will be governed by the highest predicted valyes, the difference between the simulated maximum jump height and maximum velocity and the observed maxima were calculated. The parametric plots from the refined search area are shown in Figure 4.16 to Figure 4.20. 74 Figure 4.16: Parametric surface plots for runout distance. a) RMSE of the runout distribution and b) the absolute value of the residual of the maximum runout. Figure 4.17: Parametric surface plots of Screen 1 jump height. a) RMSE of the jump height distribution and b) residual of the maximum jump height. 75 Figure 4.18: Parametric surface plots of Screen 1 velocity. a) RMSE of the velocity distribution and b) residual of the maximum velocity. Figure 4.19: Parametric surface plots of Screen 2 jump height. a) RMSE of the jump height distribution and b) residual of the maximum jump height. 76 Figure 4.20: Parametric surface plots of Screen 2 velocity. a) RMSE of the velocity distribution and b) residual of the maximum velocity. From a visual inspection of the parametric plots, some general observations can be made:  The plots show a banded behavior, indicating there is not a unique best fit for any output using a two parameter set;  Increasing roughness tends to decrease the runout distance, jump height and velocity;  Increasing normal reference deformation energy tends to increase the runout distance and maximum jump height, with a negligible effect on the velocity; and  The model is most sensitive to the roughness. There is not one clear area where the model performs best for all outputs, so a compromise was made to try to capture the general behavior of the system while also capturing the maximum values. An En,50 of 15 m3/s2 with θscale = 0.9 on the forest slope section and θscale = 0.5 on the road were selected. Using these inputs, the model would give acceptable predictions of the jump height and velocity distributions. The simulated maxima were close to the observed, 77 although the distribution did not show the same frequency of the higher velocities or jump heights. The simulations would also tend to have a slight over-prediction for the runout. Because of the stochastic roughness element in the model there will be some variation from run to run. The model was run a number of times with the same input parameters. The variation between simulated runout distributions was minimal, but the variation between the simulated Screen results was more pronounced. It was found that using the aggregate of 5 runs (each of the 100 blocks thrown 5 times in the simulation) provided repeatability of the mean and constrained the range of predicted maximum Screen values compared in this study. Details of this process have been included in Appendix A.1. All of the subsequent results use 5 combined model runs per simulation. Following to the observation that the model is most sensitive to roughness, further investigation was made into the roughness angle distribution as a means to improve the overall model performance. The effect of using the truncated normal distribution described previously (Equation 3.5) was tested. By using a truncated normal distribution the objective was to have more impacts with a low roughness to result in higher translational velocities, with the occasional high roughness value to improve the fit of the runout distribution and to better simulate extreme values. To test the truncated normal distribution, a number of different mean (μ) and standard deviation (σ) values were used with all other model values held constant. Means ranging from 0.3 to 0.5 were used with standard deviations ranging from 0.15 to 0.45. Increasing the input mean will tend to decrease the runout distance, average jump height, and average velocity. It was found that increasing the standard deviation tends to reduce the runout, and increase the variation in the jump height and velocity distributions, which also results in a slight reduction in 78 the mean values. It was determined that the agreement between runout, jump height and velocity were best where the mean was approximately equal to the standard deviation. For convenience, a mean and standard deviation of 0.5 were chosen, which could then be scaled by the roughness scaling factor, θscale, input by the user. Additional plots from the truncated normal distribution testing are included in Appendix A.2 Using the truncated normal roughness distribution, the optimal values from the Vaujany simulation were recalculated as:  Normal reference deformation energy, En,50 = 15;  Tangential reference deformation energy, Et,95 = 50; and,  Talus slope roughness, θscale = 0.72; and,  Road roughness, θscale = 0.4. The results from the Vaujany calibration are summarized in the following figures. The runout distribution, Figure 4.21, agrees much better with the observed distribution using the truncated normal roughness distribution. The Screen 1 and 2 jump height and velocity results, Figure 4.22 and Figure 4.23 respectively, are very similar to those obtained with the uniformly distributed roughness. 79 Figure 4.21: Comparison of observed and simulated runout distances from Vaujany experiments using the truncated normally distributed roughness angle. 80 Figure 4.22: Screen 1 Jump Height and Velocity Distributions using μ = σ = 0.5 and θscale = 0.72. Figure 4.23: Screen 2 Jump Height and Velocity Distributions using μ = σ = 0.5 and θscale = 0.72. 81 4.4.2 Preg – strong rock quarry The process of conducting an initial parametric study was repeated for two slope sections in the Preg quarry in Austria with the goal of establishing parameters for exposed hard rock faces. The parameters identified from the initial study using the first two slopes were then applied to the other two slopes to see if the parameters were valid for those slopes. The truncated normal distribution was used for the roughness angle calculation. An initial range of En,50 from 5 to 25 m3/s2 and θscale from 0.1 to 1.0 was used to establish the response of the model over a range of inputs similar to the optimal range for the Vaujany simulation. The distribution of the first impact locations was used to calibrate the quarry face parameters, with the distribution of runout distances used to calibrate the quarry floor parameters. The RMSE for both distributions was used as the objective parameter, similar to the Vaujany calibration. The residual maximum value for the first impact location was also calculated. Sections 1 and 2 were selected for this study. Since one of the goals of this study is to have a model that performs well with relatively few input parameters, the normal deformation energy was set using the same reference values for both the quarry face and quarry floor, meaning the roughness scale value is the only difference between the two parameter sets. The parametric plots for the first impact location for Section 1 and 2 are shown in Figure 4.24 and Figure 4.25, respectively. The parametric surface plots for the runout distribution found by varying the quarry floor roughness for Sections 1 and 2 with a constant θscale = 0.5 on the face are shown in Figure 4.26. Equivalent plots with θscale = 0.7 on the face are shown in Figure 4.27. 82 Figure 4.24: Preg quarry parametric surface plots of RMSE for varying quarry face θscale, a) first impact distribution for Section 1, and b) first impact distribution for Section 2. Figure 4.25: Preg quarry parametric surface plots for varying quarry face θscale, a) residual maximum first impact location for Section 1, and b) residual maximum first impact location for Section 2. 83 Figure 4.26: Preg quarry parametric surface plots for quarry face θscale = 0.5, a) runout distribution for Section 1, and b) runout distribution for Section 2. Figure 4.27: Preg quarry parametric surface plots for quarry face θscale = 0.7, a) runout distribution for Section 1, and b) runout distribution for Section 2. The locations of the minima do not match for the two slopes, however the inputs En,50 = 10 to 20 and quarry face θscale = 0.4 to 0.8 give results with similarly low errors. The model 84 performs well with a quarry floor roughness between 0.3 and 0.5 over the range of normal unit deformation energies tested. In order to further constrain the range for the quarry face θscale, the difference between the maximum observed and simulated first impact locations were also examined. This constrained the quarry face roughness to a range of 0.5 to 0.8. From this analysis, an optimum set of input parameters were identified as:  En,50 = 20;  Face θscale = 0.7; and  Bench θscale = 0.3. The cumulative distributions for the runout and first impact distances for Sections 1 and 2 are shown in Figure 4.28. 85 Figure 4.28: Cumulative distributions comparing the observed distributions to the simulations using the optimum parameters from the parametric study for a) Section 1 runout, b) Section 2 runout, c) Section 1 first impact, and d) Section 2 first impact. After the optimal parameters for Sections 1 and 2 were identified from the parametric study, these parameters were applied to Sections 3 and 4. The cumulative distributions for the runout distances and first impact locations are shown in Figure 4.29. 86 Figure 4.29: Results of optimum parameters identified in the parametric study applied to a) Section 3 runout, b) Section 4 runout, c) Section 3 first impact, and d) Section 4 first impact. The magnitude of the errors is similar for Sections 3 and 4 compared to Sections 1 and 2. The simulations of Sections 1 and Section 3 tended to over-predict the number of blocks stopping in the first 10 m from the toe of the quarry face; however, all the simulations gave a reasonable estimate of the maximum runout. Matching the first impact location was more 87 difficult, with the model tending to under-predict the number of blocks having a first impact in the first 2 m from the toe, but still under-predicting the maximum first impact distance. As with the Vaujany calibration, a number of runs were necessary to get repeatable results. It was found that using the aggregate of each of the 25 blocks thrown 10 times in the simulation provided repeatability of the values compared in this study, also detailed in Appendix A.2. All of the subsequent results are the averages of 10 combined model runs. With no size scaling of the roughness, the smaller blocks tend to have a further runout. This is inconsistent with the observations from the field tests, with the results for Sections 1 through 3 shown following, where a linear model fit to the data results in a nearly horizontal trend with very high scatter (low R2 values). From this result, the negative trend in the simulated runout versus mass is not supported. This is also contrary to observations of natural slopes where larger boulders tend to be found toward the toes of talus slopes. 88 Figure 4.30: Comparisons of the observed runout versus boulder mass for Preg sections 1 to 3. To address this issue, the linear roughness scaling (Equation 3.6) was introduced. To calibrate the relationship the simulated and observed runout versus boulder mass were compared. As the previous calibration of the Vaujany model was achieved without the roughness scaling, the relationship was set so that it would average to 1 (no net change on average) for the Vaujany boulders. 89 The Vaujany test used blocks with diameters ranging from 1.0 m to 1.4 m, with a mean of 1.2 m, thus the effective roughness is set to equal 1 at a diameter of 1.2 m. Using Equation 3.6, with the constraint that θeff = 1.0 at D = 1.2, the slope of the line, m, can be related to the intercept, b, reducing the problem to a single parameter for calibration. The b parameter corresponds to the maximum roughness multiplier at a diameter approaching zero. Equation 4.4: 𝑏 = 1 + 1.2 × 𝑚 The calibration process looked at three cases, (m, b) = (0.1, 1.12), (0.2, 1.24), and (0.3, 1.36). It should be noted that this relationship can yield a negative value for θeff where D > b/m. This issue was not encountered with the range of boulder sizes, slope and intercept parameters tested. However, to ensure this does not happen, the program was set so that any reduction calculated less than 0.3 would be set at a constant 0.3. This was not rigorously calibrated as there is a lack of records for blocks of this extremely large size (61 tonne minimum boulder mass to result in θeff = 0.3 using m = 0.3 and b = 1.36). These values for slope and intercept for the roughness relationship were tested against the Preg Sections 1 to 3 models. Plots for the best match using Equation 4.5 are shown following. Equation 4.5: 𝜃𝑒𝑓𝑓 = 𝜃(1.36 − 0.3𝐷) 90 Figure 4.31: Comparison of runout versus mass without roughness reduction and with roughness reduction for Preg Section 1, solid line indicating a linear regression fit to the simulated runout data. Figure 4.32: Runout versus mass for Preg Section 2, solid line indicating a linear regression fit to the simulated runout data. 91 Figure 4.33: Runout versus mass for Preg Section 3, solid line indicating a linear regression fit to the simulated runout data. Without the roughness reduction, the model tended to over-predict the runout for small boulders and under-predict the runout for larger boulders. As desired, introducing the roughness reduction reduced the average runout for the smaller boulders and increased the average runout for the larger boulders. The Vaujany simulations were re-run with Equation 4.5, and, as desired, the effect on the runout, jump height and velocity distributions are negligible. Example plots for Vaujany are shown in Appendix A.3. The duration of each test was compared to the simulated duration to calibrate the stopping criterion. Values of Hstop = 0.1R, 0.25R, 0.5R and 0.75R were tested. The experimental values for the runout duration were plotted against the runout distance. A linear trend line was fit to the data, and the 95 % prediction interval was calculated. The simulated results were plotted and the criterion that resulted in the best correspondence to the 95 % prediction interval was selected. 92 The effect of the stopping criteria on the total runout distance was found to be negligible for the range of values tested. The optimum criterion was selected with Hstop = 0.25R, with an example comparison plot is shown in Figure 4.34. Complete plots are included in Appendix B . Figure 4.34: Runout distance versus duration with Hstop = 0.25R for Preg Section 3. The results of this initial analysis were used as a starting point for the more detailed calibration of these 4 slopes and 2 others in Austria, detailed in Chapter 6. 4.5 Detailed Model Calibration 4.5.1 Ehime – open slope with weak bedrock and talus The Ehime data set provided incident and rebound linear and rotational velocities for each impact point along the test slope. The optimum parameters for Vaujany combined with the size dependent roughness scaling were used for the initial model. The masses used in the 93 simulations were exactly those provided for the spheres and cubes. As the natural rock masses were given as ranges, the mass for all the blocks within a size range was taken as the midpoint of the range. As with Vaujany, each simulation used 5 throws per every real throw. The loess regression technique described in Section 4.3.2 was used in the analysis of the data. Different values for the span used in the regression were tested, ranging from 0.5 to 1. It was found that the range of the residuals was not strongly influenced by the span chosen, but the higher values for the span resulted in a smoother curve. As a result, a span of 1 was chosen, and is used in all the following loess regressions. To compare the distributions of the residuals, ± 2 standard deviation lines were plotted over the residual distributions, approximating the 95% prediction interval. The residuals are non-normally distributed; however the distributions are fairly symmetrical, so the comparison based on standard deviations is considered valid. The initial model run using the optimum Vaujany inputs resulted in a good fit for the translational velocity profiles, shown in Figure 4.35 and Figure 4.36. 94 Figure 4.35: a) Plot of observed and simulated incident velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. The solid and dashed lines on b) and c) indicate the mean and ± 2 standard deviations on the residual distributions. 95 Figure 4.36: a) Plot of observed and simulated rebound velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. The means for the residuals are close to zero, meaning that the simulation is unbiased overall. It can also be seen that the shape of the residual distributions are similar for the incident and rebound velocities, indicating that the model is representing the scatter away from the mean velocity well. The width of the 95 % prediction interval for the simulations is very similar to that of the observations. All of the data plots below the free fall velocity line, indicating the assumption of neglible initial velocity is valid. These results gave a significant over prediction of the angular velocity, and over predicted the range of angular velocities, shown in Figure 4.37. 96 Figure 4.37: a) Plot of observed and simulated rebound angular velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. The shape factor described in Section 3.5 was introduced to attempt to reduce the bias in the angular velocity prediction. The results presented by Ushiro et al. (2006) gave plots of the linear kinetic energy, EK, versus the angular kinetic energy Eω, which were used in the calibration of the shape factor. The ratio of Eω / EK plots as a linear relationship if the only random element is the roughness angle, but it will become more scattered if the radius is varied at impact. The geometry of the particles can be used to obtain a first estimate for the shape factors. Theoretically, the spheres do not require a shape factor. The shortest and longest distances from the centroid to the surface of a cube divided by the radius of a sphere with an equivalent volume are 0.81 and 1.4, respectively. The moment of inertia of a cube is within 10% of the moment of 97 inertia of a sphere with an equivalent volume, thus the effect of representing a cube as a sphere is considered negligible. Due to the variable geometry of the natural rocks it was not attempted to theoretically estimate shape factors, they were determined by model calibration only. An example of results without and with the shape factor is shown for the cube dataset in Figure 4.38. The plots for the spheres and the natural rocks are shown in Appendix C Figure 4.38: Comparison of rebound translational kinetic energy to angular kinetic energy from Ehime cube dataset, dashed lines are from Ushiro et al. (2006). a) Without the inclusion of a shape factor, and b) with a shape factor uniformly distributed between 0.85 and 1.3. The calibration process focused on matching the angular velocity versus elevation profiles, as well as the ratio of Eω / Ek. Calibration showed that it was necessary to include a relatively small shape factor to replicate the experimental results for the sphere dataset. It is believed that this is required due to the modification of the spheres during the trajectory due to fragmentation as well as other, small magnitude processes occurring at the impact location not 98 explicitly dealt with by the model. The shape factor for the cubes was slightly less than the theoretical range, with the effects of small magnitude processes also being the probable cause of this. The natural rocks required the highest shape factor to replicate their observed behavior. After matching the observed angular velocities, the incident and rebound linear velocity profiles were recalculated to ensure the improved fit for the angular velocity wasn’t at the expense of the prediction of the linear velocity. The optimal results were found using a minimum value of 1.0 and a maximum shape factor of 1.2 for the spheres, a minimum of 0.85 and a maximum of 1.3 for the cubes, and a minimum of 1.0 and a maximum of 1.5 for the natural rocks. A plot of the angular velocity versus elevation for the natural rock dataset is shown in Figure 4.39. The loess curves fit to both data sets plot much closer to one another than they did without the shape factor, and the mean residual for the simulation comparison is close to zero. The range of simulated angular velocities is still much wider than the observed data set. The fit for the linear velocity profiles was similar to that of the original simulations. This is typical for the spheres and cubes as well, with the complete set of plots shown in Appendix C . This indicates that the shape factors described improve the model fit, but do not completely address the issues resulting from the representation of the boulders as spheres. 99 Figure 4.39: a) Plot of observed and simulated angular velocities vs elevation for the Ehime natural rock dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. 4.5.2 Vaujany All of the model features tested previously have been applied to the Vaujany model and the model has been re-run using the optimum parameters from the Ehime detailed calibration. The following figures show the results from this simulation. Figure 4.40 compares the modeled to the observed runout distribution. The model generally produces the overall runout behaviour, as can be seen in the cumulative distribution. The observed data shows many blocks becoming trapped in relatively small, discrete sections, where as the model shows more of a smooth distribution. This result is unsurprising considering the runout distribution was generated using a relatively coarse topography and generalized, stochastic material properties. 100 Figure 4.40: Comparison of observed and simulated runout distances from Vaujany experiments using the optimum input parameter set. Figure 4.41 and Figure 4.42 compare the jump height, velocity and linear kinetic energy distributions for Screen 1 and Screen 2, respectively. The general shapes of the distributions are similar for each output, which is confirmed by examining the cumulative distributions, which all 101 plot near one another. Weibull distributions were also fit to each output dataset, which are used in the reliability analysis of the model results as design inputs, detailed in Section 4.7. Details of why the Weibull distribution was chosen and how it was fit are given in Appendix D. Figure 4.41: Comparison of observed and simulated jump height, translational velocity and kinetic energy at the Vaujany Screen 1 location. 102 Figure 4.42: Comparison of observed and simulated jump height, translational velocity and kinetic energy at the Vaujany Screen 2 location. These parameters were then used for the validation simulation using the Tornado mountain model. 4.5.3 Preg quarry The detailed calibration for the Preg quarry is presented as part of Chapter 6. From this detailed calibration, the following optimum parameter set was identified for exposed, hard rock, which has been used with the Nicolum quarry validation data set: 103  Normal reference deformation energy, En,50 = 5;  Tangential reference deformation energy, Et,95 = 50;  Face roughness, θscale = 0.65; and,  Bench roughness, θscale = 0.35. 4.6 Model Validation 4.6.1 Tornado Mountain – firm talus slope Since there were only two trajectories available for comparison at Tornado Mountain, this was used as a validation dataset. To attempt to get a statistically significant model result, 50 simulated trajectories were used for each boulder. Using the Vaujany/Ehime parameters directly results in some of the simulated blocks having a runout equal to or slightly greater than the observed runout, shown in Figure 4.43. Decreasing the roughness from 0.7 to 0.6 increases the probability of the simulated blocks matching the observed runout, as shown in Figure 4.44. Figure 4.43: Runout distributions for boulders A and B from Tornado Mountain using the input parameters found in the Vaujany and Ehime calibrations, with the observed runout shown by the dashed line. 104 Figure 4.44: Runout distributions for boulders A and B from Tornado Mountain with the roughness scale reduced to 0.6, with the observed runout shown by the dashed line. The reduced roughness also gives a good match to the observed jump heights and velocities, shown in Figure 4.45 and Figure 4.46, respectively. In this case, the data is compared visually, as the observed data set was too small for fitting a loess curve with confidence. Figure 4.45: Comparison of simulated and observed jump heights versus the runout distance for Boulders A and B from Tornado Mountain. 105 Figure 4.46: Comparison of estimated and simulated velocities versus the runout distance for boulder A (left column) and boulder B (right column) from Tornado Mountain. The masses used in the simulation are the final mass. The initial mass of the boulders was greater, as there was evidence of fragmentation along the trajectory. If the initial mass is used, the roughness will be reduced due to the larger diameter, which could result in a greater runout 106 distance. Significant amounts of fragmentation were not noted for the Vaujany or Ehime tests, so this may explain the necessity to reduce the roughness to achieve the same results. The relationship between the apparent normal restitution and the incidence angle was also examined. The Tornado simulation was used for this comparison as 20 values for the incidence angle and normal restitution were available from the trajectories with tree impacts (Wyllie, 2014), which allows for a direct comparison to the simulation. The values from Spadari et al. (2011) obtained from a four slopes in Australia were also included in the comparison. The raw data points were used instead of the trend lines fit to the data. Our model does not explicitly take αin into the restitution calculation; however, the roughness value applied to the slope combined with the restitution factors can produce the same behavior. The hyperbolic restitution calculations result in values of kn’ and kt’ between 0 and 1, which are applied in the relative to the roughness angle. When these values are rotated back into the plane normal and tangential to the input slope, the apparent normal restitution factor, kn, can be greater than 1, while ensuring there is total energy lost during impact because both kn’ and kt’ are less than 1. This relationship is shown in Figure 4.47. It can be seen there is good agreement between the two observed datasets and the simulation results. 107 Figure 4.47: Comparison of Pierre 2 calculated normal restitution factors relative to the impact incidence angle for the Tornado Mountain simulation compared to incidence angle and normal restitution values presented by Wyllie (2014). Also shown are similar data from four sites in Australia collected by Spadari et al. (2011). 4.6.2 Nicolum Quarry – exposed hard rock The optimal parameters identified for the Austrian sites were used for the Nicolum quarry model to test if these same parameters would accurately reproduce the observed velocity profiles. From the impact coordinates and impact times obtained from the video records, the incident and rebound velocities for each impact were calculated. To do this, a ballistic parabola was fit between the impact coordinates, assuming that air resistance was negligible during flight. The loess curve fitting process was used to compare the simulated and observed data sets. To get an indication of the potential error in the data the theoretical free fall velocity along the profile was calculated. Since the rocks were launched with a negligible initial velocity, 108 none of the calculated velocities should plot outside the free fall line. A significant proportion of the points for both the incident and rebound velocity profiles plotted outside the free fall line above 154 m elevation. This area was the furthest from the camera and at a greater angle, decreasing the accuracy of the impact coordinates, thus the calculated velocities. The impacts above 154 m and below 102 m elevation (approximate catchment fence elevation) were excluded from the analysis, so that only the points in the open slope that the camera had a clear line of sight to were analysed. A set of 5 simulations were used for the comparison of the incident and rebound velocity profiles, shown in Figure 4.48 and Figure 4.49, respectively. Figure 4.48: a) Plot of observed and simulated incident velocities vs elevation for the Nicolum quarry dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. 109 Figure 4.49: a) Plot of observed and simulated rebound velocities vs elevation for the Nicolum quarry dataset. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. The means of the residuals are close to zero, meaning that the simulation is unbiased on average. It can also be seen that the shape of the residual distributions are similar for the incident and rebound velocities, indicating that the model is representing the scatter away from the mean velocity well. The width of the 95 % prediction interval for the simulations is narrower than for the observations. However, even within the filtered observation data some of the points do fall outside of the free fall velocity line, indicating some uncertainty in the observation data. 4.7 Use of Model Result as Design Inputs One of the motivations of making the model as simple as possible is to develop an engineering design tool that is will provide reliable results when limited site information is 110 available. A hypothetical design for placing flexible barriers at the Vaujany Screen 1 and Screen 2 locations was conducted following the ETAG 027 design guideline. The design kinetic energy, determined from the velocity and block mass, and the barrier height were determined from the PIERRE 2 analyses following the procedure outlined in Peila and Ronco (2009). The inputs to determine the design velocity, vp, are the 95th percentile of the computed velocities, vt, and a safety factor, γF. The design block mass, mp, is determined from the rock density, ρ, the design block volume, volb, and a safety factor, γM. The design interception height, hi, is determined from the 95th percentile of the forecast bounce height, hp, plus an additional clearance, f, of not less than ½ of the average block diameter. The computed velocity was calculated as the 95th percentile value of a Weibull distribution fit to the velocity distributions. A safety factor of γF = 1.16 was used, corresponding to the case with the velocity determined by computer simulation and the slope characterized by a precise topographic survey. The maximum block volume, 1.26 m3, and rock density, 2800 kg/m3, were given by CEMAGREF (2003), which were used in the calculation of the block mass, ρ volb. A safety factor of γM = 1.05, corresponding to a high level of confidence in the rock size, was selected. Table 4.1 presents the results of the calculations to determine the maximum energy levels (MEL) and serviceability energy levels (SEL). The 95th percentile of the bounce height was also calculated from Weibull distributions fit to the simulation data. The mean block diameter was 1.2 m, so the minimum clearance height is 0.6m. 111 Table 4.1: Calculation table to determine the SEL and MEL for Screen 1 and Screen 2 locations Screen vp = vt γF (m/s) mp = ρ volb γM (kg) SEL = 0.5mpvp2 (kJ) MEL = 1.3(0.5mpvp2) (kJ) hi ≥ hp + f (m) 1 24.0 3700 1070 1390 4.02 2 24.9 3700 1150 1490 4.16 These design values were then compared to the distributions of the experimental data. The goal was to look at the reliability of the predictions, i.e., what is the probability that the actual factor of safety is less than 1.0 (Duncan, 2000). Both the SEL and MEL exceed the maximum observed kinetic energies at Screen 1 and Screen 2, which were 786 kJ and 958 kJ, respectively. Log-normal and Weibull distributions were fit to the observed kinetic energy distributions for Screens 1 and 2. The Weibull distribution had a better agreement with the observed values at the higher values in the kinetic energy distribution; as a result, the Weibull distribution was used for the reliability calculations. This distribution fitting process was carried out so that the probability of exceedance for a value outside of the range of observations could be calculated. Using the SEL and MEL values from Table 4.1, and the Weibull distribution fit to the observations, the probability of the design values being exceeded was calculated, given in Table 4.2. The minimum interception height was less than the observed maximum height for both Screen locations. However, the maximum simulated bounce height did exceed the observed value for both Screen locations. The probability of exceeding the calculated interception height, hi, and the simulated maximum, hmax, are shown in the following table. More details on this process are included in Appendix D. 112 Table 4.2: Probability of failure from exceeding the net's impact capacity or interception height. Screen P(EK > SEL) P(EK > MEL) P(h > hi) P(h > hmax) 1 0.072 % 0.004 % 2.4 % 0.006 % 2 0.099 % 0.006 % 6.2 % 0.39 % The results obtained for the prediction of the SEL or MEL show using the factors of safety given in the ETAG 027 guideline will result in a low probability of failure. Using the minimum clearance specified by the guideline results in an appreciable probability of exceedance. However, a more appropriate level of safety can be used if the clearance height, f, is set so hi is equal to the maximum height for the simulation. 4.8 Discussion and Conclusions The results of the calibration and validation demonstrate that this simple model can reproduce the kinematic behavior of rockfalls on both talus and exposed rock slopes. The best fit parameter sets for firm talus slopes and exposed rock faces are given in Table 4.3. Table 4.3: Summary of best fit input parameters from calibration process. Parameter Talus Rock Face Quarry Floor Normal deformation energy, E0.5,n 15 5 5 Tangential deformation energy, E0.95,t 50 50 50 Roughness scale value, θscale 0.6 – 0.7 0.65 0.35 Shape factor for natural rocks 1.0 – 1.5 NC NC Often exposed rock faces are considered to have higher normal restitution factors than talus, see Table 2.1, so it is counter-intuitive that the normal deformation energy should be 113 higher for the talus than the rock face. One explanation for this could be that there is a greater amount of energy loss due to fragmentation for the rock on rock contact. The exposed rock faces are also steeper, which leads to more of the impacts having a high incidence angle, thus reducing the effect of the normal deformation energy. The initial calibration of the model using the roughness angle as the only stochastic element within the lumped mass model demonstrated that it could be used to accurately predict ranges of runout distances, velocities and jump heights simultaneously. The combination of a stochastic contact roughness and particle shape factor in the model can be used to produce a more realistic range of angular velocities, as well as achieving an improved fit for the distributions of runout distances, velocities and jump heights. The ability to model impacts with normal restitution factors greater than 1 is also made possible with this stochastic representation of the impact. These concepts can also be effectively applied to a 3D model space, as will be demonstrated in the following Chapters. A hypothetical design calculation using the ETAG 027 guidelines for the Vaujany test site was conducted, and compared to the experimental data using reliability engineering principles. From this comparison it was found the model’s results for translational velocities would provide design kinetic energy values with a probability of exceedance of less than 1 in 1000 for all cases considered. To achieve a similar level of safety for the design interception height, the maximum values predicted for bounce heights should be used instead of the 95th percentile plus the minimum clearance value specified in ETAG 027. 114 Chapter 5: PIERRE3D - a 3D stochastic rock fall simulator based on random ground roughness and hyperbolic restitution coefficients1 5.1 Introduction Rock fall is a type of landslide consisting of detachment, fall, rolling, and bouncing of rock fragments (Hungr et al., 2014). It may occur as single events or in clusters, but each of the moving fragments interacts dynamically only with the slope substrate (path). Rock fall occurs on natural cliffs and mountain sides, on eroded soil slopes or on artificially excavated slopes. Of special practical importance are rock fall hazards threatening buildings, roads and railways. Although a variety of techniques to prevent detachment of rock fragments from slopes exist, such mitigation measures are technically difficult and cannot be applied everywhere, due to the profusion and inaccessibility of potential rock fall source areas. An alternative approach towards mitigation of rock fall hazards is to identify hazard zones that are subject to potential impact and keep human activities within areas where the risk is acceptable. Alternatively, the hazard zones can be eliminated or trimmed back by intercepting potential rock fall using earthworks or rigid or flexible barriers. Either of these possible mitigation measures requires that the rock fall hazard intensity parameters, including movement velocity, kinetic energy, trajectory path and height of bounces be predicted for potential events. Ideally, a rock fall hazard intensity map should evaluate these parameters in a probabilistic manner within the boundaries of the hazard zone (Jaboyedoff et al., 2005; Volkwein et al., 2011). When intersected with vulnerability data, such a map can be used 1 A version of this Chapter has been published: Gischig, V., Hungr, O., Mitchell, A., Bourrier, F. 2015. Pierre3D – a 3D stochastic rock fall simulator based on random ground roughness and hyperbolic restitution coefficients. Canadian Geotechnical Journal, in press. Refer to the Preface for details of the author’s contribution. 115 to prepare a quantitative risk map required for risk assessment. At the same time, a map or diagram of bounce heights and energies can be used to optimize the location of barriers or earthworks and provide parameters necessary for their design. Since detailed quantitative field observations of rock fall are rare, the sole existing means for evaluating intensity parameters is through dynamic modelling of rock fall trajectories. This type of computer-based modelling had begun to develop in the nineteen-seventies (e.g. Piteau and Clayton, 1977, Azzimi and Desvarreux, 1977), and has now become routine as part of geotechnical hazard mitigation practice. Nevertheless, the state-of-art of rock fall dynamic modelling is still far from well-established, and many models, including some of those used in practice, are insufficiently well validated and calibrated. All models use the same method of calculating movement of rock fragments through the air, based on well-known Newtonian ballistic equations (Equation 5.5). Air resistance is routinely neglected. The main scientific challenge and the difference between various models lie in the treatment of periodic contacts (impacts) of the particle and the slope substrate (Bourrier and Hungr, 2012). The earliest models were two-dimensional and used the simplified "lumped mass" approach, neglecting shape of the particle and its rotation. Later models, such as the widely-used CRSP model of Pfeiffer and Bowen (1989), retained the lumped mass approach, but accounted for the angular momentum of the particle. In European practice, the most popular model based on similar assumptions is that of Spang and Rautenstrauch (1988), while in North America it is the Rocscience (TM) program, based originally on the work of Stevens (1998). All of these programs are being constantly updated and improved. Three-dimensional extensions of the lumped mass model, allowing the user to trace the rock fall trajectory over 3D terrain, were developed by Guzzetti et al. (2002), Agliardi and Crosta (2003), Lan et al. (2007), Dorren (2012) 116 and others. A comprehensive review of all rock fall trajectory models has recently been compiled by Volkwein et al., (2011). A consideration of rotational energy in the framework of lumped mass analysis is only possible with the assumption of a collinear impact. In such an impact, the normal to the contact surface passes through the centre of gravity of the particle. Thus, no moment is transmitted by the normal contact force and this permits separate balancing of momentum normal and parallel to the impact surface (Stronge, 2000). In practical terms, the assumption of collinearity means that the particle must be of spherical shape and the impact equations take the form of the Goldsmith (1960) equations, described below. With any other shape, most impacts are non-collinear and the normal, tangential and rotational momentum equations must be solved simultaneously for each impact, considering the precise position of the centre of gravity of the body at the moment of impact, respective to the contact point. This is the basis of "rigorous" methods of analysis, which follow the motion and impact of non-spherical particles. Surprisingly, several early three-dimensional models dealt with non-spherical particles (Falcetta, 1985, Bozzolo and Pamini, 1986, Descoeudres and Zimmermann, 1987). At the present time, with the unprecedented increase in computing power, many sophisticated new two and three-dimensional rigorous models are appearing (Glover et al., 2012, Vijayakumar et al., 2012, Andrew et al., 2012). The rigorous solution of the impact equations for a non-spherical particle is much more mathematically complicated than the simplified lumped mass solution. This, however, is not a crucial disadvantage of rigorous methods, given the high power of contemporary computers. More important concern is the availability of relevant data. Every impact of a rock particle is a random process (Bourrier et al., 2009a). The first source of randomness is the irregularity of the slope surface, presence of obstacles such as boulders or trees and the variability of the surface 117 strength and stiffness. The second source is the irregularity of the shape of the particle. For a non-spherical particle, the impact conditions change dramatically with every slight shift of particle orientation at the point of contact. An irregular particle encounters similar random impact conditions when interacting with a featureless surface, as a sphere would experience on a rough surface. In general, both random effects occur at the same time. In our opinion, randomness is the main characteristic of rock fall impacts. Therefore, it appears counterproductive to model each impact rigorously, at considerable expense in terms of data requirements and computing effort. The authors propose instead to return to the premise of the lumped mass approach, and to devote the majority of their effort to the refinement of the stochastic impact model and its calibration. Thus, the approach detailed in this paper does not differ substantially from some previous solutions. However, we introduce several new concepts which we consider important for correct simulation of rock fall behaviour: The first concept concerns the definition of "roughness". Our hypothesis is that stochastic perturbation of the slope surface topography, combined with a solution of impacts of a spherical particle using the Goldsmith Equations, can simulate both sources of randomness described in the previous paragraph. Previous solutions such as CRSP (Pfeiffer and Bowen, 1989) also perturbed the surface angle through geometric and randomized functions. In our case, we do not consider the perturbation merely in a geometric sense. We assume that the transformations required to implement roughness approximate the effects of non-collinearity of irregular particles impacting an irregular surface. We rely on calibration results to provide a proof of this hypothesis. The present version of the model makes the assumption that each impact is completely random and independent of other impacts. Irregular particle shapes may introduce dependence: 118 successive impacts of the same particle, although random, may be related by the particle shape. In the future, this may be accounted for by applying two perturbations with different distributions - one to account for particle shape and the other for surface roughness - to handle more irregular block shapes where the effects of impact eccentricity are more pronounced. The two roughness distributions will need to be calibrated by analysing experiments using particles of various shapes (Volkwein and Klette, 2014). For now, however, both effects are included in the same stochastic model, assuming that particle shapes used in field testing are reasonably similar. A very important consideration regards the treatment of non-elastic energy losses during impact. Clearly, light impacts will occur under relatively elastic conditions. Heavy impacts, brought about by greater mass of the particle and higher approach velocity, will consume energy through plastic deformation and fracturing. This subject has been addressed by a number of authors, but not in a conclusive way. The original concept of constant "restitution coefficient" can be explained by plotting a contact-force-displacement diagram as shown in Figure 5.1 (Falcetta, 1985). 119 Figure 5.1: A conceptual diagram of contact force, F and displacement during impact (Falcetta, 1985) The simple model assumes that, during light impacts, the contact force, F, increases linearly with displacement (i.e. the change of distance from the projectile center of gravity to the contact point). The slope of the linear relationship is C1. On rebound, the contact is assumed to recover elastically, with a steeper "rebound stiffness" C3. Considering the overall loss of energy in a light impact, the normal "elastic" momentum restitution coefficient, kn,el is: Equation 5.1: 𝑘𝑛,𝑒𝑙 = √𝐶1𝐶3 This is the standard, constant "coefficient of restitution" used in the earliest rock fall dynamics models. Falcetta (1985) assumed that, after the contact force exceeds a certain maximum (yield) limit F0, the stiffness of the contact will decrease to a yield slope C2. The current, variable normal restitution factor can now be calculated as a function of the maximum impact force F: 120 Equation 5.2: 𝑘𝑛 = 𝑘𝑛,𝑒𝑙√𝐶2[𝐶1 + (𝐹0𝐹 )2(𝐶2 − 𝐶1)] Note that, because kn is not a constant, it should not be called a "coefficient" to avoid confusion. Restitution “factor” may be a better term, and will be used in the following. The result of Equation 5.2 is plotted schematically by the red line in Figure 5.2. The plot is presented on a scale of incident momentum, which is proportional to the maximum impact force. The bi-linear relationship was used in a rigorous model by Falcetta (1985) and in a lumped mass model by Hungr and Evans (1989). Figure 5.2: Shapes of three relationships between the restitution factor and approach momentum, discussed in the text. Red: Equation (2) (Falcetta, 1985). Green: Equation (3) (Pfeiffer and Bowen, 1989). Black: Equation (4), Hyperbolic (Bourrier and Hungr, 2011). An alternative function to account for plastic energy losses was derived by purely empirical means by Pfeiffer and Bowen (1989) and used in the CRSP program, as well as a 121 number of other codes (e.g. Dorren, 2012, Rocscience, 2010). Here, the restitution factor is scaled against the normal approach velocity, vn: Equation 5.3: 𝑘𝑛 = 𝑘𝑛,𝑒𝑙11 + (𝑣𝑛9.15)2 The empirical factor 9.15 is simply the equivalent of the original empirical constant of 30 ft/s, expressed in SI units (m/s). If we assume that the reduction scales against momentum, rather than velocity, the term vn can be considered as momentum of a particle with a mass of unity. With that assumption, the relationship of Equation 5.3 can be plotted by the green line in Figure 5.2. Note that the specific numerical parameters used in the plots are arbitrary; the intent is to compare the shape implied by the functions. Evidently, there is some similarity between the two relationships, although one is developed from simple theoretical consideration, the other empirical. Both relationships discussed so far act as reduction functions related to a constant restitution "coefficient" kn,el. Bourrier and Hungr (2012) suggested a hyperbolic relationship, which tends to unity when the impact momentum approaches zero. (Of course, by implication, the contact stiffness is non-linear). As shown by the black line in Figure 5.2, the hyperbola can be fully defined by a single parameter M0.5, being the user-specified value of incident momentum at which the value of the restitution factor equals 0.5: Equation 5.4: 𝑘𝑛 = 𝑀0.5𝑀 +𝑀0.5 The hyperbolic formulation has several positive qualities: 1) It allows very light impacts to occur without expending significant amount of energy - a phenomenon often observed in 122 nature when rock fragments proceed in a series of low bounces over a steep slope, maintaining steady velocity over long distances. The low impact momentum in these cases occurs because the impact angle (the inclination of the approach vector with respect of the slope) is low (cf. Bourrier et al., 2009a). 2) It correctly predicts large energy expenditure in heavy impacts. 3) It is controlled by a single parameter which facilitates calibration. In the present work, the hyperbolic function is scaled against impact energy, as outlined in Section 5.2.4. The last novel concept used in the present research concerns the changes in direction over 3D topography. All 3D models use the concept of mirror reflection, to determine the steering of a particle trajectory by impacts. Most also include a random direction perturbation (e.g. Guzzetti et al., 2002, Dorren, 2012, Lan et al., 2007). In the present model, the randomness of the direction change is provided again by stochastic perturbation of the slope angle transverse to the plane of the approach trajectory. In effect, a longitudinal perturbation provides randomness of rebound velocity, while a similar transverse perturbation influences the rebound direction. The momentum equilibrium is resolved separately in both vectorial directions, using the Goldsmith Equations. Of course, the above concepts rely on many assumptions, and would be difficult to verify theoretically. Our approach is that, if these concepts are valid, they should result in a successful calibration of the resulting model. A complete calibration should be able to show good results in several dimensions, including the direction and length of trajectories, lengths and heights of bounces, linear and rotational velocities. For this reason, we use data from one of the most detailed rock fall observation programs ever carried out, conducted at Vaujany in the French Alps by the IRSTEA Institute (former CEMAGREF, Berger and Dorren, 2006). Ideally, calibrated parameters should produce good results on more than one observation site. An 123 extensive program of calibration using a variety of field data is presently in progress, using both a 2D and 3D version of the new model. 5.2 Method 5.2.1 Flight trajectory The trajectory between collisions is computed using the ballistic flight equations (neglecting air resistance). The incremental change in position is: Equation 5.5: ∆𝒙 = 12𝒂(∆𝑡)2 + 𝒗∆𝑡 where a = (0, 0, -g) is the gravitational velocity vector. Because the only acceleration is gravity, the trajectory lies in a vertical plane containing the rebound velocity vector of the last collision. The chosen calculation time steps of Δt = 0.02 s are sufficiently small to ensure repeatability; it was found that results did not change notably at time steps smaller than 0.05s. At each time step, the height h above ground is computed by subtracting the elevation of the vertical projection of the particle onto the DEM. As soon as h is negative, i.e. collision is detected, the time step is iteratively changed in order to find a more precise impact location. This iterative collision detection is terminated when the current location is less than 1 cm below the surface (i.e. -0.01 m < h < 0). 5.2.2 Collision At collision of the particle with the ground, it is necessary to define a local coordinate system (Figure 5.3a and b) and compute the two tangential components (i.e. parallel to the ground), and the normal component of the incident flight velocity and angular velocity. 124 Figure 5.3: Schematic overview of the 3D collision algorithm. Incident linear and angular velocity vectors in the original xyz-coordinate system are: Equation 5.6: 𝒗 = (𝑣𝑥𝑣𝑦𝑣𝑧) ,𝝎 = (𝜔𝑥𝜔𝑦𝜔𝑧) We span a local coordinate system at the collision site with a basis vector en corresponding to the unit vector normal to the local slope segment. The other two basis vectors are a unit vector eL parallel to the local slope and pointing in the direction of the incident vector, and a unit vector eT parallel with the local slope and perpendicular to the en and eL vectors (Figure 5.3b). Hence, the plane described by the vectors en and eL is perpendicular to the local slope segment. The slope-parallel vectors eL and eT are derived as follows: The velocity component normal to the slope is computed by the vector product of the velocity vector v and the basis vector en. Equation 5.7: 𝒗𝑛 = 𝒗 ∙ 𝒆𝑛 The velocity vector parallel to the slope in line with the incident vector is: vL = v - vn∙en, and the basis vector parallel to this velocity vector is: 125 Equation 5.8: 𝒆𝐿 = 𝒗𝐿|𝒗𝐿| referred to as longitudinal basis vector. Finally, the transversal basis vector eT can be computed as the cross product of the other two basis vectors Equation 5.9: 𝒆𝑇 = 𝒆𝐿 × 𝒆𝑛 Thus, eT, eL and en define a right-handed system. 5.2.3 Applying stochastic roughness Both irregularities of the slope surface as well as the irregular shape of the falling boulder result in an apparent random behaviour at impact. Such random behaviour is considered in our model in a stochastic manner by randomly changing the slope angle, which mimics roughness of both the surface and the boulder. Roughness is applied in that way in two directions: 1) in the longitudinal direction (parallel with the incident vector) and 2) in the transverse direction. In case of the former, the local slope angle is rotated by subtracting a random angle θ, such that the dip angle of the slope becomes flatter. In the present version of the code, the angle θ is obtained by drawing a random value for tan(θ) from a uniform distribution between 0 and θscale , here termed the longitudinal roughness coefficient. In the latter case, the slope is rotated in the transverse direction by an angle ψ. Again the angle ψ is obtained by drawing a random value for tan(ψ), however, this time from a uniform distribution between -ψscale and ψscale. Eventually, the collision takes place in a twice rotated coordinate system [eT″, eL″, en″] instead of the original local coordinate system [eT, eL, en] (Figure 5.3 c and d). For the first rotation A, the local coordinate system [eT, eL, en] is rotated around an axis perpendicular to the vertical plane that contains the particle trajectory. The underlying idea is 126 that one part of roughness arises from the rotation of the irregular body, which predominantly revolves around the horizontal axis. This vector is computed as follows: Equation 5.10: 𝒓𝐴 = 𝒗 × 𝒆𝑧|𝒗 × 𝒆𝑧| where ez = (0, 0, 1). The second rotation B, is done around a horizontal vector parallel within the vertical plane of the particle trajectory. Equation 5.11: 𝒓𝐵 = 𝒗 − 𝒆𝑧𝑣𝑧|𝒗 − 𝒆𝑧𝑣𝑧| The rotation of each basis vector ej, j = (T, L, n) around the rotation vectors rA, rB is performed by the following sequence of two matrix multiplications: Equation 5.12: (𝒆𝑇′, 𝒆𝐿′, 𝒆𝑛′) = 𝑹𝐴 ∙ (𝒆𝑇 , 𝒆𝐿 , 𝒆𝑛) (𝒆𝑇′′, 𝒆𝐿′′, 𝒆𝑛′′) = 𝑹𝐵 ∙ (𝒆𝑇′, 𝒆𝐿′, 𝒆𝑛′) Here, the rotation matrices RA and RB are Equation 5.13: 𝑹𝑋 = [𝑟12(1 − cos 𝛼) + cos 𝛼 𝑟1𝑟2(1 − cos 𝛼) − 𝑟3 sin 𝛼 𝑟1𝑟3(1 − cos𝛼) + 𝑟2 sin 𝛼𝑟1𝑟2(1 − cos 𝛼) + 𝑟3 sin 𝛼 𝑟22(1 − cos 𝛼) + cos 𝛼 𝑟2𝑟3(1 − cos 𝛼) − 𝑟1 sin 𝛼𝑟1𝑟3(1 − cos 𝛼) − 𝑟2 sin 𝛼 𝑟2𝑟3(1 − cos𝛼) + 𝑟1 sin 𝛼 𝑟32(1 − cos 𝛼) + cos𝛼] in which X refers to A or B, r1, r2, and r3 correspond to components of the rotation axes rA and rB and α to the rotation angles θ and ψ, respectively. The order or the two rotations (i.e. first around a transversal vector rA by an angle θ, then around a longitudinal vector rB by an 127 angle ψ) is chosen arbitrarily. A reversed order would result in a different rotated coordinate system. We obtain the normal vector and two tangential vectors of each of the linear incident velocity vector v and the angular velocity vector ω in the new coordinate frame [eT″, eL″, en″], by computing the scalar product of v and ω with each basis vector of the rotated coordinate frame, respectively. The vectors and their components are termed: Equation 5.14: 𝒗′′ = (𝑣𝑇𝑣𝐿𝑣𝑛) ,𝝎′′ = (𝜔𝑇𝜔𝐿𝜔𝑛) In this initial version of the model, both longitudinal and transverse roughness values are determined by uniform distributions between limits described above. Somewhat surprisingly, this extremely simple distribution, being the sole stochastic element in the model, yields good results when compared with the statistical records of particle run-out, as well as flight velocity and bounce heights (see Section 5.5 below). However, experiments with other distributions are being carried out and will be reported on in a future article. 5.2.4 Goldsmith model At collision with the ground the model by Goldsmith, (1960) is adopted. It was derived for the impact and rebound of a spherical body with a planar surface, e.g. a billiard ball on a smooth table. It considers not only normal and tangential impact velocity, but also the rotational components are included in the treatment of the impact. In the Goldsmith’s model, the normal impulse of the boulder is modified during the impact by a normal restitution factor kn. In our model, we assume the tangential impulse is also reduced by a tangential restitution factor kt. This purpose of this second restitution factor is to account for non-frictional 128 momentum losses in the longitudinal direction. The frictional losses are already included in the Goldsmith Equations, as explained below. 5.2.4.1 Hyperbolic restitution factors Unlike Bourrier and Hungr (2012, Equation 4), we here use restitution factors that are dependent on a measure of the impact energy distributed over the contact area, referred to as the plastic deformation energy. The underlying idea is that the greater the energy applied to a given area, the higher the resultant stress at the contact area, which will result in more local plastic yielding. Loss of momentum during impact is here assumed to correlate hyperbolically with the deformation work per unit contact area that is largely expended in local yielding of the soil or fracturing of the rock. Thus, the loss of momentum at impact depends on the impact deformation energy, which is proportional to the particle mass and the square of incident normal velocity, ∝ ρR3vn2, distributed over the impact area measured by the square of the particle characteristic dimension (radius) ∝ R2. Hence, assuming that density is a constant, we find that Ed ∝ Rvn2. Similar to equation 4, the functional relationship between the restitution factor and the deformation energy follows a hyperbolic function such that at high impact velocities the momentum loss is higher than at very low impact velocities: Equation 5.15: 𝑘𝑛 = 𝐸𝑛,0.5𝐸𝑑𝑖𝑛 + 𝐸𝑛,0.5 Here, En,0.5 [m3/s2] is a reference deformation energy value, i.e. the value at which kn equals to 0.5. For a given impact velocity vn, kn is computed using Edin = Rvn2. Note that the restitution factor becomes 1.0 when the deformation energy approaches zero and tends to zero at infinite values of energy. The same formulation is also applied to kt: 129 Equation 5.16: 𝑘𝑡 = 𝐸𝑡,0.5𝐸𝑑𝑖𝑛 + 𝐸𝑡,0.5 The assumption here is that the non-frictional tangential momentum losses, characterized by kt, are also scaled by the value of the normal impact energy. These losses are meant to account for cratering and displacement of material from the impact site. Note that the reference deformation energies are not the same for normal and tangential velocities, i.e. En,0.5 ≠ Et,0.5. However, for both kn and kt the impact deformation energy value is the same and depends only on particle radius and vn. 5.2.4.2 Evaluation of the friction limit angle The first step of the Goldsmith model is to compute an angle γ that represent the inclination of the incidental impulse vector at the surface of the particle. Equation 5.17: tan 𝛾 =2𝑘𝑡𝑣𝑡𝑖𝑛7𝑣𝑛𝑖𝑛(1 + 𝑘𝑛) where R is the particle radius. vtin is the tangential component of the velocity vector on the surface of the boulder, i.e. it consists of the tangential velocities vL and vT, as well as on the rotational velocity components on the boulder surface vLrot and vTrot. These can be computed from the equation vrot = ω″x R, where R = (0, 0, -R) in the coordinate frame [eT″, eL″, en″]: Equation 5.18: 𝒗𝑟𝑜𝑡 = (𝑣𝑇𝑟𝑜𝑡𝑣𝐿𝑟𝑜𝑡𝑣𝑛𝑟𝑜𝑡) = (−𝜔𝐿𝑅𝜔𝑇𝑅0) The tangential velocity at the surface of the boulder can then be computed: 130 Equation 5.19: 𝑣𝑡𝑖𝑛 = √(𝑣𝑇 + 𝑣𝑇𝑟𝑜𝑡)2 + (𝑣𝐿 + 𝑣𝐿𝑟𝑜𝑡)2 5.2.4.3 Contained slip impacts If the friction angle φ between the boulder and the ground is greater than the angle γ (i.e. φ > γ), the slip at the contact if fully contained during the collision and the linear and rotational velocities during rebound become perfectly synchronized. This condition is referred to by Goldsmith (1960) as the "non- slip" [contained slip] condition and the impact momentum loss is fully accounted for by the two restitution factors, which are independent of the contact friction angle. In the present model, the Goldsmith equations are applied in the longitudinal and transversal directions separately. Hence, the impact is described by the following three equations: In transversal direction: Equation 5.20: (𝑣𝑇𝑟𝑒𝑅𝜔𝐿𝑟𝑒) = (57𝑘𝑡−27𝑘𝑡−57𝑘𝑡27𝑘𝑡) ∙ (𝑣𝑇𝑖𝑛𝑅𝜔𝐿𝑖𝑛) In longitudinal direction: Equation 5.21: (𝑣𝐿𝑟𝑒𝑅𝜔𝑇𝑟𝑒) = (57𝑘𝑡−27𝑘𝑡−57𝑘𝑡27𝑘𝑡) ∙ (𝑣𝐿𝑖𝑛𝑅𝜔𝑇𝑖𝑛) In the normal direction Equation 5.22: 𝑣𝑛𝑟𝑒 = −𝑘𝑛𝑣𝑛𝑖𝑛 131 It can be shown that, unless the incident rotational velocity of the particle is perfectly synchronized with the translation velocity, Equation 5.20 and Equation 5.21 predict tangential momentum loss, even if kt were taken as 1.0. This loss is clearly due to frictional work, expended in synchronizing the rotational and linear velocities. Because the frictional impulse is fully contained in these impacts, the friction angle does not figure in the equations. As stated earlier, kt accounts for additional, non-frictional momentum losses. 5.2.4.4 Uncontained slip impacts If the angle γ exceeds the friction angle φ (i.e. φ < γ), slip persists throughout the period of contact and the energy that is lost depends on the magnitude of the contact friction coefficient [uncontained slip condition]. The following three equations apply in this situation: In transversal direction: Equation 5.23: (𝑣𝑇𝑟𝑒𝑅𝜔𝐿𝑟𝑒) = (𝑘𝑡 − tan𝜑 (1 + 𝑘𝑛) 00 −52tan𝜑 (1 + 𝑘𝑛) 𝑘𝑡) ∙ (𝑣𝑇𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝜔𝐿𝑖𝑛) In longitudinal direction: Equation 5.24: (𝑣𝐿𝑟𝑒𝑅𝜔𝑇𝑟𝑒) = (𝑘𝑡 − tan𝜑 (1 + 𝑘𝑛) 00 −52tan𝜑 (1 + 𝑘𝑛) 𝑘𝑡) ∙ (𝑣𝐿𝑖𝑛𝑣𝑛𝑖𝑛𝑅𝜔𝑇𝑖𝑛) In the normal direction Equation 5.25: 𝑣𝑛𝑟𝑒 = −𝑘𝑛𝑣𝑛𝑖𝑛 In this case, the tangential and angular rebound velocities depend on both φ and kt. 132 5.2.5 Model implementation The current version of the model is implemented in MATLAB (MATLAB, 2011). As input it requires a digital elevation model (DEM) defined on a regular grid, a ground material map defined on the same grid, and a table specifying initial coordinates of the release blocks. Collision detection relies on 2D interpolation, for which standard interpolation schemes implemented in MATLAB are used. Similarly, graphics tools built-in in MATLAB are used for visualization. Computation time strongly depends on the size of the DEM and the chosen input parameters. For the calibration case study presented in the following, a DEM of 560x400 m size with 35 840 cells and a cell size of 2.5 m was used. With the time step of 0.02s the average computation time for one realization (i.e. one block release) was 2.5 s. However, due to the stochastic nature of the method, computation time strongly varies among realizations (in our case from 0.2 – 4.8s) depending on if the block stops in the first meters of the slope or runs out to the slope toe with strong deviation to the slope maximum gradient. 5.3 Data The first version of the model was validated against detailed observations of multiple rock fall parameters during a real-size rock fall experiment performed by the IRSTEA Institute (former CEMAGREF) of Grenoble near Vaujany in the Savoy Alps (lat 45°12’, long 6°3’, Figure 5.4 a). Details of the experiments were described by Dorren et al. (2005). The data available to us were taken from a report by Berger (2003), which was intended for a comparative test study of different rock fall trajectory models (see also Berger and Dorren, 2006). In the experiment, granitic blocks were released down a steep slope, relatively free of trees (a snow avalanche path) using a hydraulic excavator. For each block a number of parameters were recorded, including the run-out final location. In addition, two instrumented screens were located 133 near the mid-point of the path, in which the flight velocity and height were recorded by high speed cameras. The site was located on a hill slope mantled by inactive talus deposits. Most of the surface is underlain by small boulders, embedded in a sandy matrix and covered by a thin forest soil (Figure 5.4 a and Figure 5.5 a). The middle part of the path contains discontinuous patches of somewhat coarser talus, with boulders up to approximately 30 cm diameter and less matrix (Figure 5.5 b). The lowest one-third of the path contains predominantly coarse talus, with typical grain sizes of 30 to 70 cm and no matrix (Figure 5.5 c). The slope is twice intersected by a forest road (Figure 5.5 d). Blocks were released from the upper forest road. Between this starting point and the lower forest road the path has the morphology of a shallow channel, approximately 10 m wide and up to 1 m deep (Figure 5.5 a). The distance between the release point and the slope toe is about 470 m (horizontal distance 400 m). The average slope dip angle along this distance is 32°. 134 Figure 5.4: a) Orthophoto of the Vaujany rock fall test site in France. Also shown are the end locations of 71 blocks out of a total of 100 blocks released from the upper forest road. b) Map showing end locations, position of observation screens, definition of the scatter angle as used in the calibration as well as the release point. c) Distribution of block weights. d) Run-out distance distribution. e) Distribution of scatter angle. f) and g) Velocity and jump height distribution observed at screen 1. h) and i) Velocity and jump height distribution observed at screen 2. 135 Figure 5.5: Vaujany test site photos. a) Shallow gully in the mid-reaches of the path, upslope of the observation screens. b) A talus segment outside the gully. c) Talus forming the lowest segment of the path, down-slope from the road. d) Forest road crossing the lower third of the path (Photos O. Hungr). A total of 100 blocks were released from the same spot. Figure 5.4 c shows the distribution of block masses. They range from 1500 kg to 3500 kg, which corresponds to equivalent block diameters of about 1 to 1.4 m at a rock density of 2600 kg/m3. For 71 blocks, the coordinates of the final landing position could be recovered from the maps given by Berger (2003). The remaining 29 blocks could not be located at the end of their trajectories. We used the end position along profile distance as well as the scatter angle with respect to the start position (Figure 5.4d and e). Velocity and jump height of each block, as recorded by cameras at Screens 1 and 2 are shown in Figure 5.4b (Berger, 2003). Histograms of data recoded at these screens are 136 shown in Figure 5.4f - i. The digital elevation model used for our study consists of regularly gridded elevation data with grid spacing of 2.5 m (Figure 5.4 b). 5.4 Model Calibration To calibrate the model against observations at the Vaujany test site, we performed a Monte-Carlo search. We randomly varied the four parameters En,0.5, Et,0.5 (normal and tangential reference deformation energies), θscale , and ψscale (i.e. the longitudinal and transversal roughness coefficients) to find a parameters set that best fits the observations. As a criterion to assess fit quality, we compare the cumulative distribution functions (CDF; Figure 5.6), and computed the root mean square error (RMSE) of the difference between the modeled and the observed CDFs. In a first Monte-Carlo step, we varied the parameters 2000 times using a uniform distribution in a generous range: 3 < En,0.5 [m3/s2]< 25; 20 < Et,0.5 [m3/s2] < 100; 0.3 < θscale < 2; 0.3 < ψscale < 2. From these 2000 model realizations, we chose the best-fit parameters set and initiated a second Monte-Carlo step, in which we drew parameters from a normal distribution with the mean being the best-fit parameters from the previous step, and a standard deviation of 10% around the mean. Again, we ran the model with 2000 random parameter sets. From these realizations we chose the best-fitting parameter set. 137 Figure 5.6: Cumulative distribution function (CDF) and probability density functions (PDF) of the observed quantities shown in Figure 5.4d-i (black line). The red curves are the results of the best-fit model using only one parameters set for the entire slope. For the orange curve the parameters of the lower forest road was changed such that fewer blocks stop on it (see Figure 5.7a). The blue curve represents the best-fit model for which a third material was introduced for the talus cone below the forest road (see Figure 5.7b). The single parameter set corresponding to only one, homogeneous ground material over the entire extent of the slope was able to fit the screen velocities and jump heights, but not the entire run-out distribution at the same time. In the first calibration attempt, we focussed on the screen observations and the run-out distribution in the first 200 m, corresponding to the channeled slope topography. As a first priority, we chose models that fit the screen observations well. Hence, in the aforementioned two-step calibration procedure, we chose the 10% of the 138 2000 model realizations (i.e. 20 model realizations) that best fit screen observations, and then, out of these, chose the ones that best fit the run-out distribution in the first 200 m of the slope. Hence, the ultimate best-fit is a compromise between fitting screen observations and run-out distance. A second calibration attempt was conducted with three parameter sets, fitted to approximate distribution of materials over the site. This led to a further improvement of the results. 5.5 Results 5.5.1 First calibration attempt, single homogenous ground material In a first attempt to fit the observations, we assume that the entire slope can be represented by a single ground material, i.e. we used only one parameter set. The run-out distribution as well as the screen observations can be fit reasonably well in the first 200 m of the slope (red curves in Figure 5.6). Further downslope, the run-out distribution does not fit well as many blocks (20-30%) stop on the forest road in the model (Figure 5.6a). Those blocks that pass the road mostly run down until they reach the bottom of the valley. This pattern can also be observed in the map of run-out distribution in Figure 5.7a. The model was able to reproduce the two run-out lobes that were observed below the forest road (Figure 5.6b and Figure 5.7a). However, the model over-predicts the scatter angle of the southern lobe (negative values in Figure 5.6b). Also, in the north (positive values), the model predicts blocks that run-out further north than observed. However, in this latter case, the model may not be inconsistent with the experiment. During the experiment, the end position was determined for only 71 out of 100 released blocks. Some of the blocks jumped out of the channel and into the forest at the top of the slope, and thus were not tracked to their final position (10 blocks). For some other blocks the final position was not determined, simply because the blocks could not be found (19 blocks). 139 Note that also about 10% of all blocks are predicted to travel out of the channel in Figure 5.6b, which thus is in agreement with the experiment. Figure 5.7: Maps of end positions of the three best-fit model versions shown in Figure 5.5. Also indicated are the observed end positions in magenta. The modeled block velocity distributions for Screen 1 and 2 are in close agreement within the lower range (i.e. for the 60-70% lowest velocities), but under-predict the observed maximum velocities (Figure 5.6c and e). In contrast, the maximum jump heights are predicted reasonably well by the model, but the frequency of the lowest jump heights is under-predicted (Figure 5.6d and f). Recall that the chosen best-fit model attempts to simultaneously fit screen velocities, jump height, and the run-out distribution in the first 200 m, and thus is a compromise. It was not possible to predict larger velocities, and, at the same time, match the run-out distribution in the first 200 m as will be discussed in Section 5.6.1. 5.5.2 Second calibration attempt, three ground materials In the second calibration attempt, three different sets of parameters were assigned to various zones of the slope, where the air-photo and field observations (Figure 5.4 and Figure 5.5) show ground conditions strongly contrasting with those of the upper part of the path. To account for the smooth, compacted road surface (Figure 5.5d), we assigned high reference velocities and low roughness values to the road (Figure 5.8a). These values were not derived by rigorous calibration, but only chosen extreme enough to allow most blocks to pass the road (i.e. low 140 roughness coefficients, high reference deformation energies). The resulting run-out distribution in the lower portion of the slope is still not in agreement with the observations (orange curve in Figure 5.6a). Most blocks that pass the road run out until they reach the valley bottom. Note that model result for the screen observations (velocity and jump height) have not changed, because the same ground material was used for the area around the screens. Since it is not possible to match the entire run-out distribution with two ground materials, we introduced a third material to apply to the coarse talus surface below the forest road. Figure 5.8: Maps showing the spatial distribution of ground material (Figure 5.5). The parameters for the corresponding material numbers are given in Table 5.1. a) Different properties were assigned to the forest road to inhibit blocks from stopping on the lower road. b) A third material was assigned for the slope below the lower forest road. Adding a third set of material properties is in agreement with the change to a substrate that contains blockier and less matrix-supported talus material (Figure 5.5c). The resulting run-out distribution is in good agreement with the observed run-out distribution (blue curve in Figure 5.6a). The best-fit parameters for the material below the forest road correspond to a material with higher impact energy loss, i.e. lower restitution factors due lower reference velocities, and higher roughness angles. Note that the discontinuous patches of talus, present on the upper part of the 141 slope (Figure 5.5b), were not mapped in detail. Thus, the primary generalized material of the slope applies to the entire area upslope of the logging road, including a range of surfaces such as shown in both Figure 5.5a and 5b. Only the surface in Figure 5.5c and the surface of the road were differentiated in the three-material model. Table 5.1: Best fit parameter sets used to reproduce observed rock fall data from the Vaujany test site. One ground material: Material 1 Material 2* Material 3 En,50 [m3/s2] 15 - - Et,50 [m3/s2] 68 - - θscale 1.16 - - ψscale 0.93 - - Two ground materials: Material 1 Material 2* Material 3 En,50 [m3/s2] 15 45 - Et,50 [m3/s2] 68 206 - θscale 1.16 0.39 - ψscale 0.93 0.31 - Three ground materials: Material 1 Material 2* Material 3 En,50 [m3/s2] 15 50 5 Et,50 [m3/s2] 68 206 34 θscale 1.16 0.55 1.58 ψscale 0.93 0.23 1.03 *Arbitrarily set to 3 time higher/lower values than material 1 to avoid block stopping on the forest road 142 5.6 Discussion 5.6.1 Model and data uncertainties For the best-fit model using three materials, we estimated the uncertainty in the cumulative distribution function due to the stochastic nature of the model. Note that for the modeled CDF curves we used 1000 released blocks, i.e. each of the 100 blocks was released 10 times. Such oversampling compared to the real experiment minimizes random fluctuations in the CDF curves that result from the limited number of block releases. However, it also allows assessing the uncertainty associated with a potentially statistically insufficient number of observations. To do so, we performed a model run with 10 000 realizations (each block released 100 times) and randomly drew 100 of these 10 000 realizations, and computed the CDF for this subset of blocks. The procedure was repeated 1000 times resulting in a bundle of 1000 CDFs for each observation parameter (Figure 5.9). The shading shown in Figure 5.9 corresponds to the 95% percentile of the CDF, i.e. the shading envelopes 95% of the bundle of 1000 CDFs. We extract for each modeled median and maximum value an uncertainty estimate representing the variability in the stochastic results. Also indicated in Figure 5.9 are the median and maximum values for the screen velocities and jump heights both from the observations and the model results. In brackets, the 95% percentile for each value is indicated. The observed median screen velocities are well within the modeled 95% percentile. However, the maximum velocities are higher than indicated by the 95% percentile. Hence, the discrepancy between the model and experiment cannot be explained by an insufficient number for blocks released in both the model and experiment. Nevertheless, the modeled uncertainty range indicates that the range of possible maximum velocities can vary between -20% and +20% at repetition of the experiment with 100 blocks. In contrast to the maximum velocities, the maximum jump heights are well contained 143 within the uncertainty range, but the median values are over-predicted by the model. Again the range of possible maximum jump heights may vary considerably between about -30% and +30%. Figure 5.9: Model results for the case with three ground materials. The variability of the results related to low sampling numbers is derived by randomly drawing 100 realizations from 10 000 model realization. This is repeated 1000 times. The blue shading represents the 95% percentile of the 1000 random CDFs. a) to f) as in Figure 5.6. Although there may be a potential under-sampling bias in the experimental data, the 95% percentile of maximum velocities derived from the model indicates that there still remains a discrepancy (although not a severe one) between the modeled and observed extreme velocities. In an attempt to reproduce higher maximum velocities, we re-calibrated our model, but this time by only giving weight to the screen velocities. Figure 5.10 shows a comparison between the calibrated model weighting observed parameters as described in Section 5.4, and the one weighting only velocities. For the latter, the entire CDF shifts towards higher velocities by about 144 2 m/s (Figure 5.10c and e). However, the maximum velocities can still not be reproduced. Note that the run-out distribution also changes such that fewer blocks arrest along the first 200 m of the slope (Figure 5.10a). Higher velocities come at the expense of higher discrepancies in run-out distribution. Note that during the experiment it was observed that some of the blocks that stopped in the first part of the slope tended to be more irregular in shape. These impacted the ground such that rotational and translational momentum strongly declined – more than a model assuming collinear impacts can account for. Figure 5.10: Comparison of model results based on two different calibration strategies; calibration giving weight to all parameters (as described in Section 4), and giving weight to screen velocities only. a) to f) as in Figure 6. A model being able to account for such non-collinear impacts – either in a rigorous manner, or preferably as an additional simple stochastic element - would allow to produce higher velocities and still have more stops in the uppermost part of the slope. To address the 145 problem of limited model variability, especially regarding extreme velocities, three model modifications will be considered in future research: 1) as mentioned before, assuming a uniform distribution for roughness angles θ and ψ may be too simple. Other roughness distributions, such as normal distribution, could yield improved results as more extreme roughness angles are allowed. 2) The block radius has been assumed constant so far, and was derived from the assumption that blocks are spheres in the first approximation. Effects of block irregularity on impact (like the aforementioned observations during the experiment) require accounting for more complex, rigorous impact physics. However, such effects may also be approximated with an additional stochastic element. For instance, randomly varying block radius at impact could help to introduce more model variability, hence more extreme velocities (although this approximation does not intend to explicitly simulate non-collinearity). 3) The assumption of a uniform set of parameters over the large area of the upper path is a major generalization. More spatial variability in ground properties may yield more variability and may result in a better representation of extremes. These issues are subject of current ongoing research, with very promising initial results. 5.6.2 Probabilistic rock fall hazard maps Considering ground and block roughness at collision as a random change of slope orientation reproduces inherent randomness in impact physics. With this stochastic element, the model is well suited for probabilistic rock fall hazard assessment. In Figure 5.11, we present a set of maps that can be produced from the best-fit model presented in Figure 5.6 and Figure 5.7c. A key property of stochastic modeling approaches in the context of hazard assessment is the capability of considering also very rare events that constitute residual hazard and risk. As Figure 5.11a shows, releasing only 100 blocks in the model as in the experiment produces an impact 146 pattern similar to the one observed: a generally narrow run-out path, a double-lobed path at the bottom as a result of the local topography, and a few impact traces outside the channel at the top of the slope (that are also in agreement with the observations; see Section 5.5.1). However, if we exceed the number of block releases by far, i.e. by using 10 000 blocks, the impact rate of very rare events can be estimated as shown in Figure 5.11b (note the logarithmic scale). The area affected by block impacts is now larger compared to the observations. This area includes an estimate of residual hazard – the component of total hazard that can contribute heavily in risk calculations, because elements at risk are often situated in marginal areas. Hence, for the following maps 10 000 blocks are always used. Figure 5.11c and d show maximum jump height per cell and maximum velocity, respectively. We can derive for each cell a rate of exceeding a certain velocity (Figure 5.11e) or jump height (Figure 5.11f). If a rate of rock fall events (events per year) for a given slope is known, the maps can be translated into hazard maps providing an estimate of the annual probability of impact, velocity or jump height. Intersection of the maps with infrastructure exposure, vulnerability and associated (monetary or human) loss estimates would yield a quantitative risk estimate for the site. 147 Figure 5.11: Probabilistic hazard maps. a) Impact rate derived from 100 model realizations (block releases). b) Impact rate derived from 10 000 model realizations (case with three materials). b) Maximum jump height per cell. c) Maximum velocity. d) Maximum kinetic energy. e) Rate of exceeding a velocity of 20m/s. f) Rate of exceeding jump heights of 5 m. For b – f 10 000 model realizations were used. 148 5.7 Conclusions We presented a novel 3D rock fall modeling approach with the following key characteristics: 1) it is based on the lumped-mass assumption, 2) rotational energy transfer at impact is considered using Goldsmith’s equations, 3) restitution factors depend on impact deformation energy (including particle mass) via a hyperbolic relationship, 4) ground roughness and particle irregularities are accounted for by a random perturbation of the slope angle at the impact location, and 5) lateral run-out scatter is enabled by randomly perturbing transverse slope orientation at impact. The model is controlled by only four parameters (for each material type), which makes it a pragmatic and flexible tool for rigorous calibration and forward run-out simulation. Through its stochastic nature, the model is ideal for probabilistic rock fall hazard studies. We calibrated the model against a comprehensive dataset of a rockfall experiment performed at the Vaujany rock fall test site in France (Berger and Dorren, 2006), which includes run-out end locations as well as velocities and jump heights at two observation screens for up to 100 blocks released from the same location. It was possible to reasonably reproduce run-out distribution, and velocity and jump heights distributions in the uppermost 200 m of the slope with a single material and single set of parameters. Better simulation of the overall runout distribution was obtained by using three surface materials - still very strongly generalized over the extensive area of the experimental site. The model has elements that tend to be conservative, since run-out distance, lateral scatter and jump heights are more extreme than in the experiment. Velocity distributions at the control screens are also well reproduced, except that the predicted velocity extremes are somewhat lower than measured. Further developments of the model intended to improve the maximum velocity predictions are suggested in Section 5.6.1. Given the 149 simplicity of the model and the small number of input parameters utilized, the performance of the model is remarkably good in all possible dimensions available for comparison. Future work will include testing other possible roughness distributions. For example, a low- kurtosis normal distribution gives promising results, particularly with respect to modelling extreme bounces. Further important improvement will include two roughness perturbations, one connected with the surface roughness and the other with the particle shape. It is possible that comparisons with “rigorous” models, which account for particle shape directly, may be used to calibrate the shape perturbation distribution. Hence, the work presented here is part of a suite of ongoing studies attempting to find well-calibrated parameters sets for different site conditions and to fine-tune details of the physics and the stochastic processes the model is based on. This work is currently continuing and its results will be reported in future articles. 150 Chapter 6: Stochastic analysis of rock fall dynamics on quarry slopes2 6.1 Introduction: Rock Fall Hazards in Quarries and Open Pit Mines Falling blocks and rock mass falls are a common cause of accidents in quarries and open pit mines (Figure 6.1). In Austria, a new regulation concerning health and safety at work plants for quarries came into force in 2011 (Tagbauarbeitenverordnung, Kolenprat, 2012). This regulation includes the requirement of specific risk assessment of geologic hazards with the purpose to determine areas in which danger of rock fall exists. These dangers are often assessed incorrectly and the areas at the toes of the slopes affected by rock falls are estimated wrongly. Figure 6.1: Photo of a rock fall accident at a quarry in Austria (Kolenprat, 2012). The dimensions required for assessment of hazards to personnel and equipment are defined in Figure 6.2. They include the distance from the slope toe to the location of the first 2 A version of this Chapter has been submitted for publication: Preh, A., Mitchell, A., Hungr, O., Kolenprat, B. 2015, in review. Refer to the Preface for details of the author’s contribution. 151 impact on the quarry floor (“First Impact distance”, ID) and the “Runout Distance” (RD), measured from the toe of the face to the point of deposition of the boulder. The size and form of the blocks, the height of the fall, the slope geometry, the surface roughness and the damping parameters (factors of restitution) of the ground control the mechanism of the fall (sliding, rolling, bouncing and falling) and cause different impact and runout distances. Similar safety regulations exist in Canada, although to date, quantitative risk analysis is not required there. Figure 6.2: Danger areas associated with a quarry face (Kolenprat, 2012). The purpose of this paper is to present results of full-scale tests of rock falls at several Austrian and Canadian quarries and to evaluate predictions of the basic hazard dimensions ID and RD using a new computer model of rock fall dynamics. 152 6.2 Full Scale Testing at Quarry Sites 6.2.1 Test procedure A total of ten different sites with several rock types (calcareous rocks, granite, basaltic rock and iron ore, etc.) were investigated, 9 of which are in Austria and one in British Columbia, Canada. A total of 589 boulders were released individually on 14 different slopes with slope heights between 9 and 75m (Figure 6.3), of which 311 tests on 7 different slopes were selected for detailed back-analysis during this research. The mean slope angles of all the quarry faces ranged from 53 to 71º. A metric 3D-image of the cut slope was used to generate representative cross sections for the analysis. The 3D-image was created using the ShapeMetriX3D system (3G Software & Measurement GmbH, Austria) which is based on the principles of classical photogrammetry. Figure 6.3: Rock fall drop tests in progress at the Klöch Quarry in Austria (Photo A. Preh). 153 During all the Austrian experiments the protocol was identical. First, the volume of the boulder to be thrown was estimated by measuring the three dominant boulder axes and recorded. The mass of the boulders varied from 0.002 Mg up to 44 Mg (16 m³). An excavator was used to move the boulders gently over the edge of the slope. At the end of each experiment, the stopping point of the boulders was measured. The impact distances were estimated by video analyses. After every 5 throws, the boulders deposited on the quarry floor were removed by a wheel-loader, in order to prevent boulder interaction during the following throws. Rock fall trajectories were filmed using two digital cameras, installed at the top and near the toe of the slope. In the Canadian quarry (Nicolum), the falling fragments were filmed by two analog video cameras (Figure 6.4). A total of 34 trajectories were recorded for the Nicolum quarry tests. The blocks used ranged in mass from 0.03 Mg to 3.04 Mg (equivalent diameters of 0.28 m to 1.3 m). The trajectory of each boulder was traced on the videos with reference to survey marks on the slope, so that the geometric position and timing of each impact could be recorded. This allows estimation of flight heights using ballistic calculations. Runout of the boulders was not measured at Nicolum, because there was an experimental rock fence at the foot of the slope which stopped most of the rocks. The fence impact parameters are not reported in this paper. 154 Figure 6.4: Nicolum quarry test site, British Columbia, Canada (Photo O. Hungr). 6.2.2 Summary of field test results One objective of the Austrian component of the research project was to document the runout distance as a function of the slope height. A summary histogram of all measured runout distances, RD, normalized by slope height, is shown in Figure 6.5 a) and fitted by the Weibull and log-normal distributions. The log-normal distribution shows the better fit, with a Kolmogorov–Smirnov test statistic of 0.037 and p–coefficient of 0.372, indicating a good fit and 155 a high confidence in the fit. Similarly, the summary of height-normalized impact distances, ID, can be approximated by the log-normal function, as shown in Figure 6.5 b). These summary plots and statistical distributions can be used for preliminary empirical estimates of the hazard dimensions on typical quarry slopes, corresponding to desired probability levels. Figure 6.5: Summary of all drop tests in nine Austrian quarries. Histograms and probability density functions of (a) the first impact distance (412 measurements) and (b) runout distance (589 measurements), measured from the slope toe and normalized by slope height. The following general observations may be drawn from the field test project in both countries:  Higher slopes typically produce larger average runout and impact.  Higher slopes produce impact distances that are more variable  Steep and/or overhanging slopes (convex shaped cross sections) typically produce larger impact distances and lesser average runout distances  More gentle or concave- shaped cross sections produce impacts near the toe, long runout distances and outliers 156  The surface roughness has a strong influence on bouncing and therefore on the impact distance  High bouncing could generally not be observed  Heavy blocks show lower rebound heights compared to light blocks  Light blocks show larger impact distances compared to heavy blocks  The type and strength of the rock does not strongly influence the interaction, except if fragmentation takes place  A strong influence of the block shape on the runout process could not be observed (possibly because few strongly elongated or tabular blocks were tested).  Slab shaped blocks (non-spherical) show slightly higher bouncing and more variation in the runout direction compared to equidimensional blocks  The dominant motion is bouncing; rolling occurs only at the beginning and the end of the process. 6.3 PIERRE 2D and 3D, Lumped-mass Stochastic Model of Rock Fall Ballistics Computer-based dynamic models of rock fall dynamics have been in existence since mid-1970’s (e.g. Azimi and Desvarreux (1977) and Piteau and Clayton (1977)) and a large number of models are now available, both in two and three dimensions (for a recent review, see Volkwein et al., 2011). Following some preliminary testing of several existing codes, the authors focused on a stochastic model developed recently at the University of British Columbia, Canada. While the recent trend in rock fall ballistics research uses complex rigorous modelling methods (e.g. Glover et al., 2012, Vijayakumar et al., 2012, Andrew et al., 2012), the new model, PIERRE2D and 3D returns to the simpler lumped-mass concept similar to the widely used CRSP model (Pfeiffer and Bowen, 1989), with some significant modifications. 157 Like all other models, PIERRE uses classical ballistics equations to describe the flight of the rock fall fragments, neglecting air resistance. The impacts of fragments on the slope surface are described by a combined mechanical/stochastic model. The mechanical part is based on the exact solution of a collinear impact of a sphere on a planar surface, derived by Goldsmith (1960). Two types of impact are recognized, differentiated by a friction-limiting angle γ, which represents the inclination of the incidental impulse vector: Equation 6.1: )1(7)(2tanninnintinttkvRkvk  Here, R is the particle radius, vt, vn and ω are the incident (superscript in) tangential, normal and rotational velocities respectively, while kt and kn are tangential and normal restitution factors, discussed further below. The type of impact depends on whether the friction angle between the substrate and the surface of the particle, φ, is greater of smaller than γ. In every impact, there must be sliding between the particle surface and the substrate, except in case of perfect synchronization, where Rω = vt. If φ > γ, however, the sliding process ends before the termination of contact. Such a “limited sliding” impact is described by: Equation 6.2: ininnintttnttrerenretRvvkkkkkRvv720750072075 (Here, the superscript re indicates rebound quantities). It is of interest that Equation 6.2 does not imply the absence of energy losses due to friction. Comparing the kinetic energy represented by the rebound tensor in Equation 6.2 with 158 the incident tensor, it is clear that energy is spent in proportion to the difference between the incidental rotational velocity Rωin and tangential velocity vtin. However, because the slip is contained within the contact period, the frictional impulse is independent of the friction angle and the energy loss is fully accounted for by the two restitution coefficients: This finding is strictly valid only for collinear impacts. In the case where φ < γ, the slip is not fully contained and the rebound conditions depend on φ: Equation 6.3: ininninttnnntrerenretRvvkkkkkRvv)1(tan250000)1(tan The normal restitution factor represents the fractional retention of normal momentum and is identical to the “restitution coefficient” used by Goldsmith (1960), Pfeiffer and Bowen (1989) and many others. The tangential restitution factor, kt, is added in this project to account for non-frictional momentum losses that may result from cratering, breakage and similar contact phenomena, not considered in the simple theory. During impacts, energy is spent by deformation, fracturing, friction and mass displacement (cratering), all of which intensify with increased incident kinetic energy. To account for this, the restitution factors in the PIERRE model are scaled by a hyperbolic function, proposed originally by Bourrier and Hungr (2012). While the original model by Bourrier and Hungr (2012) scaled the factors in proportion to the incident momentum, in the present code it is scaled to the unit incident normal kinetic energy, Enin, calculated as Enin = Rvn2. The reasoning here is that contact deformation is proportional to the normal kinetic energy, with characteristic 159 dimensions vn2R3, distributed over the contact area, with a dimension R2. The scaled restitution factors are calculated as: Equation 6.4: ninnnn EEEk5.0,5.0' tinntt EEEk,5.0,5.0' E0.5,n and E0.5,t are reference unit energy constants (input parameters), defined as the unit energy at which each respective factor evaluates to 0.5. Note that kn and kt are both scaled against the normal unit energy, as this is assumed to predominantly control the impact deformation. The hyperbolic Equation 6.4 tend to 1.0 in very light impacts and to 0 at extremely strong impacts. Each is determined by a single input parameter. The stochastic part of the theory is based on the definition of “roughness”, being a perturbation of the local slope angle at the impact site, designated as θ. In contrast to previous models, which ascribed a geometric meaning to roughness (e.g. a sinusoidal function in CRSP, Pfeiffer and Bowen, 1989), the present algorithm considers the roughness perturbation strictly as a stochastic parameter. A hypothesis is made, that random perturbation of the contact angle can mimic the effects of randomness introduced in non-collinear impacts between an irregular particle and a rough surface. While it is difficult to evaluate this hypothesis theoretically, the authors rely on calibration results to confirm its validity. While it is acknowledged that other aspects of rock fall impacts than the contact inclination are random, in the present model, the roughness is the only random element in the theory. Again, this assumption must be validated by calibration. 160 The model is very sensitive to the roughness angle, with an increase in the mean roughness angle resulting in a decrease in the runout distance, a decrease in the distance from the toe for the first impact point, and a decrease in the mean velocity. However, higher roughness angles result in more extreme values for the location of the first impact point and velocity, and create a more realistic distribution of stopping distances. Initial calibrations using a uniform distribution for the roughness had limited success in reproducing the maximum values for the outputs tested. To address this limitation, a skewed distribution with a lower mean, but a small number of high roughness values was tested. After some experimentation, the Box-Mueller approximation of a normal distribution, Equation 6.5, was selected (Scott, 2011). The tangent of the roughness value is calculated at each impact from a constant mean and standard deviation, μ and σ respectively, and two random numbers, normally distributed between 0 and 1, Si and Sj. The resulting distribution is linearly scaled by the user input θscale, which is a number greater than zero and typically 1 or less. Equation 6.5: tan 𝜃 = 𝜃𝑠𝑐𝑎𝑙𝑒 [𝜇 + 𝜎(cos 2𝜋𝑆𝑗)√−2(ln(𝑆𝑖))] For numerical stability, values of θ that are negative are discarded and new random samples Si and Sj are drawn until a positive value for θ is obtained. Not doing so can result in a case where the perturbation to the slope produces a trajectory that lies below the ground surface. The roughness angle distribution was tested on several sites, and a mean of 0.5 and standard deviation of 0.5 consistently performed well. After initial calibration work, it was found that the performance of the model could be improved by also taking the size of the particle into consideration when calculating roughness, as larger particles will be able to override smaller perturbations on the slope surface. After testing 161 runout versus mass for a number of simulations, the linear empirical Equation 6.6 with the fixed values of b = 1.36 and m = 0.3 was selected to simulate this relationship. Equation 6.6: tan 𝜃𝑒𝑓𝑓 = 𝑡𝑎𝑛𝜃(𝑏 − 𝑚𝐷) Note that all of the parameters used in Equation 6.5 and Equation 6.6 are fixed empirical constants in the model and not a part of the input parameter set to be dimensioned by calibration. The stopping criterion for the model was set when the total energy head (translational plus rotational kinetic energy divided by mass of the block) became less than 0.25 of the block diameter, following an impact. While the above-described relationships were largely derived empirically, albeit with due consideration of observed physical trends, they are validated by the calibration results reported below. A three-dimensional version of PIERRE, based on the same collision algorithm, has been developed simultaneously and is described in Gischig et al. (2015) [Chapter 5]. 6.4 Calibration of the Model 6.4.1 General Quarry faces and floors are fairly similar from one location to another in terms of geometry (roughness and angle) and material character. The goal of the calibration effort in this project is to produce a tool capable of first order approximate prediction of rock fall behavior in all dimensions. Therefore, it is desirable to minimize the number of required input parameters. The surfaces under consideration were of only two types: quarry face (rock face produced by uncontrolled blasting) and quarry floor (relatively smooth surface formed of heavily compacted crushed rock). Both surfaces are visible on Figure 6.3 and Figure 6.4. The input parameters for the two surfaces are independent of particle shapes and it is assumed that variation in particle 162 shape is encompassed in the stochastic roughness function. This assumption is partly justified by the fact that most test particles were selected so as to avoid extreme shape variability. In addition to the mass of the particle, the simple lumped mass algorithm described above requires the input parameters listed in Table 6.1, together with values adopted at the conclusion of the calibration process. These parameter sets were used to calculate all the results reported below, for all the sites in Austria and Canada. Table 6.1: Input parameters required by the model and their values, determined by calibration, used in all the analyses Parameter Quarry face Quarry floor Unit energy for normal restitution, E0.5,n 5 5 Unit energy for tangential restitution, E0.95,t 50 50 Roughness scale value, θscale 0.65 0.35 The present method of calibration uses a trial-and-error approach assuming the variables are independent. However, some guidance was obtained in the initial stages of the search from a more systematic calibration approach, based on plotting the Root Mean Square Error, RMSE, of the ID and RD distributions as functions of the input parameters (refer to Section 4.3.1). The objective here was to find a set of parameters that would minimize the RMSE for both the first impact distance (as an index related to bounce height) and the runout distance. This process will be described in a subsequent article. The final selection of the optimal parameters reported in Table 6.1 was made by trial and error. 163 6.4.2 Calibration results, Austrian sites The quarries Preg, Pauliberg and Klöch were selected from the nine investigated sites for the analysis. The three sites represent typical prototypes of cut slopes in Austrian quarries in terms of the slope height, the slope angle and the surface roughness. The Preg site is a serpentinite quarry (serpentinised dunites and peridotites) located in the Eastern Alps (Styria) at 47°17'03.3"N 14°56'51.8"E, developed in the Kraubath ultramafic massif. The rock is extremely strong (Approximate UCS of 170 MPa) and the rock mass forming the quarry face is of good quality. The sites Pauliberg and Klöch belong to the Transdanubian Volcanic Region (historical term Styrian Volcanic Arc) and both are basalt quarries. Pauliberg is located at the end of the volcanic arc in the transition zone between Eastern Alps and Pannonian Basin (47°35'05.0"N 16°20'21.1"E), the Klöch quarry is located at 46°46'17.0"N 15°57'53.7"E in the northern part of the formation. The hard durable basaltic rock at both sites exhibited a high UCS of between 200 and 250 MPa. The cross-sections used to back-analyse the Austrian cases are shown in Figure 6.6. In each case, the computer analysis consisted of throwing boulders with the same masses as recorded during the field experiments, as shown in Figure 6.7. However, to increase the statistical significance of the stochastic analyses, each boulder was “thrown” 10 times (from the same initial location, determined by the position of the excavator at the crest of the slope). 164 Figure 6.6: Cross-sections used for the back analyses of Austrian quarry tests. Figure 6.7: Mass distributions of boulders used in the field tests and in the calculations. a) Klöch, b) Pauliberg, c) Preg 1, d) Preg 2, e) Preg 3, f) Preg 4, and g) Nicolum. 165 Four cross-sections were used to carry out tests at both Preg and Klöch quarries. At the latter site, as shown in Figure 6.6 a, the four cross-sections were nearly identical. In the analysis, therefore, they were replaced by a single cross-section marked as b) in the Figure. An example locus of trajectories analysed at one of the sites (Pauliberg) is shown in Figure 6.8. Figure 6.9 shows the linear and rotational kinetic energy lines from the same site, corresponding to the trajectory yielding the maximum Impact Distance, ID and the maximum Runout Distance, RD. Figure 6.8: Example locus of calculated trajectories from the Pauliberg site (80 Boulders, each simulated 10 times). 166 Figure 6.9: Example linear and rotational kinetic energy lines, Pauliberg site. The normal unit deformation energy was set as being the same for both the quarry face and quarry floor, meaning that the roughness scale value is the only difference between the two parameter sets. Initial testing examined the following ranges of input parameters:  En,50 = 2.5 to 25  θscale = 0.1 to 1.0. After an initial examination of RMSE distribution plots, the next stage of calibrations focused on the probability of exceedance values for the first impact and runout distances, ID and RD. The probabilities were calculated directly from the test and data sets. For each site, the median value of the two parameters was calculated and the value representing the 95% percentile (5% probability of exceedance). These values were then used to compare the calculated and observed values. The input variables listed in Table 6.1 were obtained by adjusting the input parameter set used in all the analyses, so as to obtain the optimal agreement between calculations 167 and measurements in all the tests (including the Canadian results discussed below). The comparisons are reported in Table 6.2 and Table 6.3. The “goodness of fit” was judged subjectively, primarily by comparing the average values shown in the last lines of Table 6.2 and Table 6.3, but also considering similarity of other comparisons, as discussed below. Equal weight was given to both the impact distance, which is a measure of velocity and height of trajectory bounces on the face, and the runout distance, which depends on all the conditions of a given solution. Table 6.2: Comparison of calculated and observed Runout Distance statistics, obtained with the optimal set of input parameters, shown in Table 1. (Note: in the last line, the distance RD is divided by the slope height, H, see Figure 6.6) Slope (number of throws) RD 50% (m) RD 95% (m) RD maximum (m) Calc. Obs. Calc. Obs. Calc. Obs. Klöch, combined (97) 5.4 5.6 15.6 15.1 31.4 20.5 Pauliberg, comb. (80) 5.1 5.0 12.5 10.5 24.3 20.8 Preg 1 (25) 7.9 10.5 16.7 15.8 20.1 17.0 Preg 2 (25) 8.0 6.0 16.7 21.9 23.7 24.0 Preg 3 (25) 5.5 6.5 15.4 14.0 21.6 18.0 Preg 4 (25) 7.4 5.0 16.3 13.6 25.0 16.5 Average (normalized by H) 0.39 0.38 0.92 0.90 1.45 1.15 168 Table 6.3: Comparison of calculated and observed First Impact distance statistics, obtained with the optimal set of input parameters, shown in Table 1. (Note: in the last line, the distance ID is divided by the slope height, H, see Figure 6.6) Slope (number of throws) ID 50% (m) ID 95% (m) ID maximum (m) Calc. Obs. Calc. Obs. Calc. Obs. Klöch, combined (97) 1.8 0.8 3.4 2.6 5.0 5.2 Pauliberg, comb. (80) 3.8 3.0 6.4 6.0 8.1 7.5 Preg 1 (25) 1.1 1.0 2.2 3.8 3.3 4.3 Preg 2 (25) 1.6 1.0 3.0 2.4 3.9 3.0 Preg 3 (25) 1.0 0.8 2.4 1.5 3.3 4.0 Preg 4 (25) 2.2 1.0 4.2 2.4 5.4 5.0 Average (normalized by H) 0.11 0.07 0.21 0.18 0.28 0.28 Detailed results of all the back-analyses from the Austrian sites are presented in the form of histograms in Figure 6.10 through Figure 6.12. In each case, the final analyses were conducted with the input parameter set from Table 6.1 and using the actual tables of boulder masses, per Figure 6.7. In the analyses, 10 throws were simulated for each boulder. The runout distance histograms show distributions that can be fit reasonably well with the lognormal distribution, similar to the fit of the summary data shown in Figure 6.5. The First Impact distance data are more difficult to describe statistically, probably because this quantity is more variable and the number of throws at each site is limited. The mean fall velocities in the field experiments were estimated by recording the overall duration of the trajectory, from release of the particle, to its deposition. The median observed 169 trajectory durations for all the tests are compared with those calculated by the model in Figure 6.13. Figure 6.10: Comparison of calculated and observed hazard distances, Klöch site (four combined sections, 97 throws observed, 970 throws calculated) Figure 6.11: Comparison of calculated and observed hazard distances, Pauliberg site (80 throws observed, 800 throws calculated) 170 Figure 6.12: Comparison of calculated and observed hazard distances, Preg sites (25 throws observed, 250 throws calculated per site). 171 Figure 6.13: Comparison of simulated and observed 50% of exceedance runout durations. 6.4.2.1 Influence of the degree of detail of the modeled slope geometries on the calculated hazard distances A key question, regarding the presented stochastic approach, relates to how many uncertainties can be captured using the stochastic roughness function. In particular, how important is the degree of detail of the input slope geometries in terms of the surface roughness. To answer this question the sites Klöch, Pauliberg and Preg 1 have been re-analyzed, using a coarse and a fine slope geometry. The fine slope geometries were derived from the 3D-images (created using the ShapeMetriX3D system) having a resolution of about 3 cm and are shown in 172 Figure 6.6. Coarse geometries were manually smoothed by three to five straight lines. The results of the comparison are presented in Table 6.4. Table 6.4: Comparison of runout distance statistics using a fine and coarse (smoothed) slope geometry Slope (number of throws) RD 50% (m) RD 95% (m) RD maximum (m) Fine Coarse Fine Coarse Fine Coarse Klöch, combined (97) 5.4 5.9 15.6 15.7 31.4 23.5 Pauliberg, comb. (80) 5.1 5.5 12.5 13.1 24.3 19.5 Preg 1 (25) 7.9 8.1 16.7 17.4 20.1 25.7 Median values, values representing the 95% percentile as well as the maximum runout distances show a fairly good agreement. This preliminary analysis leads us to conclude that the stochastic roughness function is capable to reasonably approximate the influence of the surface roughness on the calculated hazard distances independently, without it being necessary to take in account the precise surface roughness in the input slope geometry. However, further study regarding the optimal approach to representing the slope geometry is required. 6.4.3 Calibration results, Canadian quarry The Canadian experiments were conducted by the British Columbia Ministry of Transportation and Highways at Nicolum Quarry, 200 km east of Vancouver (49°22’49”N, 121°22’31”W). The quarry is developed in quartz diorite belonging to the Coast Plutonic Complex. The rock is extremely strong and the rock mass is of excellent quality, although controlled blasting was not used (Figure 6.4). As stated earlier, runout or impact distances were not recorded in the Canadian field tests, because there was a fence located at the toe of the slope. 173 Instead, the lengths of jumps of each projectile along the length of the slope were measured by processing the impact points on the slope, as recorded by the videos. The same optimal input parameters identified for the Austrian sites (Table 6.1) were used also for the Nicolum quarry model. The comparison between observations and analysis were made in terms of the maximum bounce heights, hmax, plotted as a function of elevation along the steep slope profile. The impact coordinates and impact times were determined from the video records as described above. In the processing stage, a ballistic parabola was fit between each contiguous pair of impact coordinates, assuming that air resistance was negligible during flight. The maximum bounce height above the slope surface was then determined from the geometry of the parabola. The same quantity was calculated for each bounce during the computer analysis. Figure 6.14 shows the results of the analysis, with each boulder thrown once (note that the trajectories shown in this figure are only a random sampling – like the field tests themselves). The data used for comparison purposes was the maximum bounce height vs. elevation (Figure 6.15). The area above el. 154 was the furthest from the camera and at a greater angle, decreasing the accuracy of the impact coordinates, thus the bounce heights. The impacts above el. 154 and below el. 102 (near the catchment fence) were excluded from the analysis, so that only the points in the open slope that the camera had a clear line of sight were analysed. 174 Figure 6.14: Example locus of calculated trajectories from the Nicolum site (34 Boulders). Figure 6.15: Comparison of observed and simulated maximum bounce heights at different elevations along the Nicolum quarry profile. 175 Again, the analysis used the real masses of the tested boulders. Each test was simulated five times. Figure 6.15 shows the comparison of the simulated and observed hmax profiles. Both the observed and simulated bounce height data were fitted by linear regression, to indicate the average trends. The observed and calculated maximum bounce heights show similar average trends and similar variability, despite the fact that a three-dimensionally complex slope had to be represented by a single cross-section. The analysis of the Nicolum experiments was repeated using the 3D version of the PIERRE model, described by Gischig et al. (2015). In 3D, it is necessary to define “roughness” both in longitudinal (downslope) and transverse (cross slope) directions. The longitudinal and transverse roughness are designated θscale and ψscale, respectively. The calibration performed by Gischig et al. (2015) on a site involving a talus slope found that the ratio of ψscale to θscale ranged from 0.65 to 0.8. The 3D model was modified from the version published, to include the truncated normal distribution for the roughness angle. The calibrated parameters from the 2D model were applied directly to the 3D model, and then the calibration focused on the longitudinal and transverse roughness scale. The calibrated 3D parameters are as follows:  Reference normal deformation energy, En,0.5 = 5 m3/s2  Reference tangential deformation energy, Et,0.95 = 50 m3/s2  Longitudinal roughness scale, θscale = 0.45  Transverse roughness scale, ψscale = 0.3 The values for the longitudinal and transverse roughness result in a norm of 0.54, comparing favourably with the roughness value of 0.65 used in the 2D analyses (see the first column of Table 6.1). As shown in Figure 6.16, there is generally good agreement between the 176 observed impact locations and the areas with a higher simulated impact density. There was not a significant variation in the distributions of scatter angles for the simulations, so this was not used to select the best fit. The agreement between analysis and observation is not very good near the upper limit of the trajectories, above elevation 154, where the slope was not well oriented relative to the sight line of the camera. Figure 6.16: Frequency of impacts per square meter obtained from 250 model runs within the 154 m to 102 m elevation range (note the logarithmic scale). Red points correspond to observed impact locations. 6.5 Discussion and Conclusion The research reported in this article represents an attempt to calibrate a simple dynamic model, with a single set of input parameters, to measured results at a number of quarry faces. As can be seen in the summary Table 6.2 and Table 6.3 and in the comparison plots of Figure 6.10 177 through Figure 6.16, the agreement between observation and analysis is never perfect, but largely consistent. Natural rock fall is a strongly random process and a “perfect” agreement between observed and calculated trajectories can never be achieved. But the model presented here provides approximate simulation of every dimension of the phenomenon, including both average behaviour (see last lines of Table 6.2 and Table 6.3), distributions of output variables (Figure 6.10 through Figure 6.16) and, importantly, the values of extremes. Both distances and velocities (as reflected in the trajectory durations and the heights of bounces) are represented. A three-dimensional analysis of one of the sites also shows good approximation of the lateral and longitudinal distribution of impact points on the rock slope. To the authors’ knowledge, such a comprehensive comparison of analysis and measurement for multiple rock fall sites has never been done. Our success points to the possibility of making practical first order stochastic predictions within reasonable limits of approximation at similar sites, using a simple model with a minimal number of input parameters. The randomness of the rock fall process has many sources. As shown by the above comparisons, our stochastic model is capable of simulating both the average behavior and the variability of the process, in spite of simplifications such as representing a faceted, highly irregular three-dimensional rock face by a single idealized two-dimensional cross-section. It was also found that the results are not too sensitive to the precision of the geometrical input, presumably because the stochastic “roughness” model has a dominant influence. The results underline the importance of accounting for certain aspects of rock fall: a) Particle mass has strong influence on the energy conservation during impacts. The hyperbolic model, scaling the restitution factors against unit incidental kinetic energy (Equation 6.4) is shown to perform well. b) The stochastic model of “roughness”, using the Box-Muller normal 178 distribution approximation, represented by Equation 6.5, is essential and forms the sole stochastic element in the theory. c) The results are improved by including a simple, empirical size correction, using Equation 6.6. d) The hypothesis that the assumed “roughness” model can reasonably approximate the effects of both particle irregularity and non-collinear impact conditions is confirmed. Of course, particle shape is a parameter that could systematically affect the stochastic model and violate the implicit assumption of independence. It is possible to incorporate a certain index of particle shape into the “roughness” model, to simulate this systematic effect. In this modification, each group of particles of a similar shape “class” would have a separate roughness distribution. This would be a worthwhile subject for further research, but it will require a large amount of field data, parsed into certain shape categories (e.g. Volkwein and Klette, 2014). Further research regarding this aspect of rock fall is in progress. 179 Chapter 7: Integrating size and shape effects in three dimensions In general, this model has been first developed and tested using a 2D spatial representation, then extended to three dimensions. The normally distributed roughness, size dependent roughness scaling and shape factor described in Chapter 4 were still in development when the paper included as Chapter 5 was submitted for publication. The analysis of Nicolum quarry presented in Chapter 6 focused on the areal distribution of the impact locations and the simulated bounce heights. The ballistic parabola used to determine the bounce heights can also be used to determine incident and rebound velocities. The total energy head stopping criterion had been implemented in the 3D model for Nicolum, however it had not been tested in an open talus slope model. These model features have been incorporated into the 3D version for completeness. This Chapter details the implementation of the size and shape dependent model inputs in the 3D model, the 3D analyses of Vaujany using the updated model, as well as the velocity data from the Nicolum quarry 3D simulations. This demonstrates how the model features developed and tested in 2D can then be applied to the 3D model. The values obtained from the hypothetical design shown in Section 4.7 have been compared with equivalent results from the 3D model to ensure the two are comparable. The application of the 3D model as a design aid is also discussed. 7.1 Model inputs The longitudinal roughness has also been modeled using a truncated normal distribution with a mean of 0.5 and standard deviation of 0.5, multiplied by θscale, exactly like the 2D model. Since the transverse roughness can perturb the slope in either direction the program uses a mean of 0 and a standard deviation of 0.5, multiplied by ψscale. This implicitly makes zero transverse roughness the highest probability output. 180 The roughness reduction function, Equation 3.6, has been applied only in the longitudinal direction. The calibrations shown in the following section produced acceptable matches to the observed data, thus it was determined it was not necessary to apply the roughness in the transverse direction. The shape factor is applied to the combined translational velocity, i.e. in both the longitudinal and transverse directions. Thus, the incident angular velocity terms in Equation 5.20 and Equation 5.21 are multiplied by Reff as opposed to R, with the calculation of vT and ωL shown in Equation 7.1 as an example. The calculation for vL and ωT is equivalent. Equation 7.1: (𝑣𝑇𝑟𝑒𝑅𝜔𝐿𝑟𝑒) = (57𝑘𝑡−27𝑘𝑡−57𝑘𝑡27𝑘𝑡) ∙ (𝑣𝑇𝑖𝑛𝑅𝑒𝑓𝑓𝜔𝐿𝑖𝑛) 7.2 Calibration 7.2.1 Vaujany The optimal inputs identified in Chapter 5 for three slope materials were used as a starting point for this recalibration (refer to Table 5.1). The model has been recalibrated with the normally distributed roughness and the size dependent roughness reduction, without the inclusion of the shape factor so that the effects of the changes to the roughness representation can be isolated. The model reference points for the restitution factors were made equal so that this analysis would be more similar to the 2D analysis, and the roughness was adjusted to the new optimum, summarized in Table 7.1. 181 Table 7.1: Model input parameters from Vaujany 3D recalibration. Material 1 – Upper Slope Material 2 - Road Material 3 – Lower Slope No Shape Factor En,50 [m3/s2] 15 15 15 Et,50 [m3/s2] 50 50 50 θscale 0.7 0.4 1.5 ψscale 0.6 0.3 1.0 With Shape Factor 1.0 – 1.5 En,50 [m3/s2] 15 15 15 Et,50 [m3/s2] 50 50 50 θscale 1.0 0.4 1.5 ψscale 0.9 0.3 1.0 Each of the real blocks was simulated 5 times in the program to reduce the random variation. The runout distribution calculated for the new optimum without the inclusion of the shape factor is shown in Figure 7.1. As with Chapter 5, the runout distribution is determined by finding the elevation of each simulated boulder end point and projecting that elevation onto the “average” 2D profile utilized in Chapter 4. 182 Figure 7.1: Vaujany runout distributions using the 3D simulation with the truncated normal distribution and size dependent roughness reduction. The jump height and velocity distributions obtained using these parameters are shown in Figure 7.2 and Figure 7.3 for Screens 1 and 2, respectively. The distributions were determined by isolating the point where each trajectory crossed the ground elevations of the Screens. 183 Figure 7.2: Jump height and linear velocity comparisons for Vaujany Screen 1. 184 Figure 7.3: Jump height and linear velocity comparisons for Vaujany Screen 2. A shape factor ranging from 1.0 to 1.5 was applied to the boulders, equivalent to that of the 2D model. There was an issue with blocks stopping shortly after the release point, Figure 7.4a, but if the blocks stopping in the first 10 m are filtered out, there is an excellent agreement with the runout distribution, Figure 7.4b. 185 Figure 7.4: Vaujany runout distributions using the 3D simulation a) raw runout distribution, and b) filtered to remove the blocks stopping within the first 10 m of the release point. The real boulders would have some initial velocity due to the way they are released by the excavator in the experiments which is not represented in the simulation. This is a potential physical explanation for the blocks that stop very shortly after release. The observed and 186 simulated jump heights and velocities at Screens 1 and 2 show very good agreement, shown in Figure 7.5 and Figure 7.6, respectively. Figure 7.5: Jump height and linear velocity comparisons for Vaujany Screen 1 using the shape factor. 187 Figure 7.6: Jump height and linear velocity comparisons for Vaujany Screen 2 using the shape factor. Compared to Figure 5.6, or Figure 7.2 and Figure 7.3, the fit of the model has been improved by combining the truncated normally distributed roughness and shape factor. The SEL and MEL calculations presented in Section 4.7 use the P95 of the distributions for the velocity. The P95 values for the velocity in the 2D model could be used to provide design energies with a 188 low probability of exceedance. The P95 values for the velocity distributions at each screen in the 3D model compared to the ones obtained from the 2D simulation are summarized in Table 7.2. Table 7.2: Comparison of the design velocity values from the 3D and 2D simulations. Screen 1 Screen 2 3D P95 Velocity 20.2 m/s 23.0 m/s 2D P95 Velocity 20.8 m/s 21.5 m/s Velocity difference -0.6 m/s / -2.9% 1.5 m/s / 7.0% SEL difference - 50 kJ / -4.7% 170 / 15% MEL difference -70 / -5.0% 220 / 15% The prediction of the jump height was better at both Screen locations using the 3D model. The P95 jump height increased from 3.42 m to 4.16 m at Screen 1 and from 3.56 m to 4.31 m at Screen 2. 7.2.2 Nicolum The velocity profiles from the 3D model were compared to the observed data using the loess regression curves, as was done for the 2D model. The incident and rebound velocities for each impact were plotted against the elevation for the impact, and the data from the simulation was compared to the regression curve fit to the observed data. Values for θscale ranging from 0.45 to 0.6 were combined with values of ψscale ranging from 0.3 to 0.5 such that the norm of the two vector components ranged of 0.54 to 0.78, approximately the same as the roughness used in the 2D analysis. The ratio of ψscale / θscale ranged from 0.67 to 0.83, consistent with the results from Chapter 5. The quality of fit for the model was determined on the basis of the incident and rebound velocity versus elevation 189 profiles. There was not a significant variation in the distributions of scatter angles for the simulations, so this was not used to select the best fit. The best fit for the 3D analysis was found using the following parameters:  Reference normal deformation energy, En0.5 = 5 m3/s2;  Reference tangential deformation energy, Et0.95 = 50 m3/s2;  Longitudinal roughness scale, θscale = 0.45; and,  Transverse roughness scale, ψscale = 0.3. Plots of the incident and rebound velocity profiles are shown in Figure 7.7 and Figure 7.8, respectively. Figure 7.7: a) Plot of observed and simulated incident velocities vs elevation for the Nicolum quarry dataset using the 3D model. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. 190 Figure 7.8: a) Plot of observed and simulated rebound velocities vs elevation for the Nicolum quarry dataset using the 3D model. b) Distribution of residuals about the loess curve fit to the observed data. c) Distribution of residuals from the simulation relative to the loess curve fit to the observed data. As with the 2D simulation, the parameters resulted in an unbiased estimate of the mean velocity, with the simulated data tending to be clustered closer to the mean. The rebound velocity profile also has the same behavior of first under-predicting, then over-predicting the rebound velocity. The model was also tested using a shape factor of 1.0 to 1.5. This was found to have almost no effect on the results, thus it was not tested further. 7.3 Sensitivity to Shape Factor The sensitivity of the model to the shape factor was tested using the Vaujany simulation. Using the constant roughness parameters, the model was tested with shape factors of 1.25 and 191 1.0. It was found that with constant roughness, increasing the shape factor tends to increase the jump height and velocity, shown in Figure 7.9, and an increase in the runout distance. Figure 7.9: Effect of the shape factor on the simulated jump height and velocity with constant roughness, column a) is Screen 1 and column b) is screen 2. The Shape = 1.5 and No Shape cases as well as the optimum simulation without the shape factor were used to see if the shape factor has an effect on the scatter angle (defined in 192 Chapter 5). Plots for the scatter angles calculated for each trajectory are shown for the optimum fit with the shape factor, as well as no shape factor with the same roughness, and no shape factor with reduced roughness to match the runout better are shown in Figure 7.10. Figure 7.10: Comparison of scatter angles for a) Shape factor = 1.5, b) No shape factor, c) optimum simulation without the shape factor and d) cumulative distributions calculated for a), b), and c). 193 There is not a large difference between the results with and without shape factor; however the shape factor does result in slightly more lateral scatter in the trajectories. The two scatter angle distributions without the shape factor plot very close to one another, indicating the reduced lateral dispersion observed for the second case was not simply a result of the simulation tending to have shorter runout distances. While there is a slight effect, looking at the distribution of impact points, the trajectories are still concentrated in the channel on the upper slope, with the trajectories then breaking into three “fingers” in the lower half of the slope, shown in Figure 7.11. Figure 7.11: Impact rate map for the Vaujany site using the Shape Factor, aggregate of 10 000 trajectories. 7.4 Discussion The re-calibration of the Vaujany simulation with the changes to the roughness calculations and the inclusion of the shape factor have demonstrated the validity of these 194 concepts applied in 3D following their development in 2D. The results in Chapter 5 indicated it was difficult to match the bulk of the distributions while capturing the higher jump heights and velocities as well. The inclusion of the shape factor has helped reduce the error in the upper end of the distributions. To be consistent with the 2D calibration, the input reference normal and tangential deformation energies were made constant for the three materials. The input values for the two roughness scales, θscale and ψscale, had to be reduced from the calibrated values from Chapter 5 to match the runout, velocity and jump height distributions. The velocity and jump height distributions are similar for the 2D and 3D models. The difference between the SEL and MEL determined using the 2D or 3D models was within 15 %, which was also the error on the kinetic energy estimates given by CEMAGREF (2003) for the experimental data. The prediction of the jump heights is significantly better for the 3D model relative to the 2D model. To expand on the hypothetical design presented in Section 4.7, the results of the 3D modeling can be used to determine the required width of a barrier or catchment system. A simple way to do this is to determine an acceptable residual hazard and use the impact rate map, Figure 7.11, to isolate the areas with an impact probability greater than the residual hazard. This would give the required width for the barrier. The energy can also be mapped to determine the spatial limits of areas with trajectories that exceed a given energy level, as shown in Figure 7.12. This figure shows the frequency of simulations having a kinetic energy greater than a certain threshold taken from 500 throws. These results could then be combined with a onset probability to establish return periods for the 195 exceedance of a certain energy at any point within the model area to create hazard curves for a fully quantitative hazard assessment, as recommended by Abbruzzese and Labiouse (2014). Figure 7.12: Spatial distribution of trajectories with modeled linear kinetic energy greater than a) 500 kJ, b) 1000 kJ, and c) 1500 kJ. Maps produced using an aggregate of 500 model trajectories. The velocity data from the Nicolum quarry tests was used to demonstrate the validity of the 3D model for determining velocities in hard rock quarries. In this case, since there was no previous 3D calibration performed the 2D model parameters were used as the start point for the calibration. This calibration showed a slight reduction in the roughness relative to the 2D model was needed to produce results that matched the observations well. This makes intuitive sense, as the 3D model allowing for lateral movements will lengthen the true trajectory relative to the path of steepest descent. However, this observation is the opposite of the Vaujany model, where the 196 3D model required higher roughness values than the 2D model, even with the shape factor applied. This suggests that while the 2D calibrations provide a useful starting point for a 3D calibration, simple relationships between the two may not exist for roughness. No benefit was observed from applying the shape factor to the quarry model. This may be a result of the quarry having steep walls with rough edges resulting from uncontrolled blasting, with the rough surface having a much greater effect than the shape of the block itself. For an open talus slope where the boulders are large relative to most of the surface features, the effect of the particle shape may become more pronounced. 197 Chapter 8: Conclusion The recent trend in computer modeling of rock falls has been towards making increasingly complex physically based models. This has been made possible by the great increases in computer processor speeds that have occurred since the earliest rock fall simulation codes were written. The model presented here has returned to a simpler, lumped mass particle representation, with stochastic elements to achieve the variability that is observed in actual rock fall events. The characterization of rock fall sites in engineering practice generally gives a limited amount of information. The level of detail required for the model inputs has been selected with the amount of site information typically available in mind. The speed of modern processors can be more effectively utilized by creating statistically significant populations of simple model results as opposed to making a fewer model runs of a more complex – and more poorly constrained – model. Creating a simple model was one objective; however it was necessary to not oversimplify the problem. This model has been developed with the specific objective of being able to provide estimates of the kinematic properties of rock fall events with a sufficient level of accuracy to produce a safe design, while not being excessively conservative. In order to achieve this, a physical basis was established for all the key empirical parameters: the restitution factors, the rebound process, the shape factor and the stopping criterion. While attempting to find exact solutions for the mechanics of the various impact processes may not be practical, these theoretical constructs are very useful to form the basis of the empirical models. 198 8.1 Applicability of the PIERRE 2 Model This thesis has outlined a relatively simple model that uses a lumped mass, spherical particle representation with a stochastic surface roughness angle, restitution factors defined by hyperbolic functions, and a stochastic particle shape factor. This model has three key hypotheses:  Using a single, stochastic variation of the local slope angle, referred to as the roughness angle, to represent irregularities of both the particle and the slope within an idealized collinear impact model can effectively replicate the effects of non-collinear impacts of irregular particles.  The hyperbolic restitution factors provide a simple way of representing the variability in energy conservation for different impact conditions.  Varying the particle radius at the impact point can represent the complex interactions between translational and rotational motions for non-collinear impacts. The roughness angle hypothesis was tested by calibrating multiple models using runout, velocity and jump height data. The calibrations showed that general model parameters based on a simple description of the material can be used to predict runout distributions and kinematic properties. The roughness angle has been described using a normal distribution truncated so that negative values are excluded for numerical stability. The extreme values that can be produced by this distribution have been shown to improve the model performance relative to a representation using a uniform distribution for the roughness. A linear size correction has been included in the model, which was calibrated using runout versus mass distributions. The addition of this size correction has increased the likelihood of large rocks having longer runout distances, which is consistent with the observed behavior. 199 The term restitution coefficient, which is commonly used in practice, is replaced in this work with the term restitution factor. This is meant to recognize that the conservation of energy or momentum in an impact is not simply a material property. Established principles in impact mechanics, as well as experimental results involving impacts with controlled geometries and impact conditions on known substrate materials, were used to construct a conceptual framework for the restitution factors. From this conceptual framework, normal and tangential restitution factors that are a hyperbolic function of rock diameter multiplied by the normal velocity squared have been defined. The use of hyperbolic functions for the normal and tangential restitution means that the restitution factors for any impact condition can be defined by a single reference point on the hyperbolic curve, keeping with the objective of a simple model. The validity of the restitution factors as defined here has been demonstrated through the calibration process. The assumption that the impacts are collinear, thus conservation of momentum during the impacts can be described using the equations of Goldsmith (1960), simplifies the mathematics of the problem greatly compared to a rigorous rigid body solution. It has been demonstrated here that combining the roughness angle with the hyperbolic restitution factors in this collinear impact framework can effectively represent more complicated interactions. This model is able to recreate observed relationships between the impact incidence angle and the normal restitution, including the phenomenon of normal restitution factors of greater than 1 at shallow incident angles. The model is able to do this as a by-product of the calculation. Stochastic variations can be observed in all terms of the rebound velocity equations, consistent with more complicated stochastic models. To create more realistic estimates of the angular velocity, the shape factor was introduced. The hypothesis that varying the radius at the impact point could improve the angular 200 velocity simulation was tested using the data with the ratio of rotational to linear kinetic energy from the Ehime test site. Theoretically, the cast spheres should not require a shape factor and cubes should have a shape factor ranging from 0.81 to 1.4, if the process can be represented using only the block geometry. The values fit through calibration did not perfectly agree with the theoretical values; however the two ranges were close enough to conclude the geometry based shape factor provides a reasonable first estimate. Calibration using the natural rock dataset showed an improvement in the prediction of angular velocity, while still utilizing using a sphere representation. The model defined in two dimensions can be easily extended into three dimensions. This requires having a slope defined by a surface as opposed to a linear profile, and defining two roughness angles, one in the direction of the incident trajectory and other 90° to the first. Then, by allowing the block to move with 3 degrees of freedom of linear motion and 2 degrees of freedom of rotational motion, the model becomes fully 3D. The assumption of collinear impacts also permits the decoupling of the two components of the translational motion at impact (i.e. in the longitudinal and transverse directions). It has been shown that features developed and more rigorously calibrated in 2D can be successfully applied to 3D with limited additional calibration work. The extension of the model into three dimensions allows for the explicit representation of spatial topographic variations in the model. This allows for the calculation of the spatial distribution of hazards in an area. The ability of the 2D model to produce appropriate values for the design of flexible barriers was demonstrated in Chapter 4. The application of the 3D model to creating hazard maps was demonstrated in Chapter 5. The model was also used for the determination of hazard zones within quarries in Chapter 6. In Chapter 7, the updated 3D model was shown to be able to 201 reproduce results similar to those of the 2D model for design inputs, and provide additional information on the linear extent of protective measures. This demonstrates that the model performs well with limited site information, and satisfies the objective of showing the model is robust from an engineering design standpoint. The best fit parameters found for talus, exposed hard rock and quarry floors are summarized in Table 8.1 for both the 2D and 3D versions of the model. It can be seen the only difference between the 2D and 3D best fit parameters is the splitting of the roughness into 2 independent angles. Table 8.1: Best-fit model parameters for both the 2D and 3D models. Parameter Talus Rock Face Quarry Floor Normal deformation energy, E0.5,n 15 5 5 Tangential deformation energy, E0.95,t 50 50 50 2D Roughness scale value, θscale 0.6 – 0.7 0.65 0.35 3D Longitudinal roughness scale value, θscale 1.0 – 1.5 0.45 - 3D Tangential roughness scale value, ψscale 0.9 – 1.0 0.3 - Shape factor for natural rocks 1.0 – 1.5 - - - Indicates parameters which were not calibrated 8.2 Model Limitations This model does not accurately predict every minute detail of rock fall motion, nor was that the stated goal of the research. The assumed impact conditions of a spherical particle having a collinear impact with a plane surface is a significant simplification of the actual impact conditions of an irregular rock on a natural slope. One of the implications of this can be seen in the simulated distributions of the jump heights for Vaujany and the first impact locations for the Austrian quarry sites, where the model is not able to closely reproduce the observed 202 distributions. However, as was demonstrated, through careful interpretation and use of the data, a good estimate of the flight height can be obtained. Fragmentation can be a significant source of energy loss during impacts, especially during impacts of hard rock on hard rock. The fragmentation process can also produce small, high velocity particles that can do significant damage. This model does not address fragmentation during impact. For impacts involving hard rock slopes, some of the effects of fragmentation may be observed with lower deformation energy values required to match runout distributions. The model has not been rigorously tested against cases with high degrees of fragmentation, as a result its performance for such cases is unknown. The shape factor used provides a good estimate of the average angular velocity; however it still over-predicts the range of values. This is likely another consequence of the assumed collinear impact process. With the knowledge that the model will tend to overestimate the range, a designer can use that information when choosing a design angular kinetic energy. As with any computer model, the quality of the results will be highly dependent on the knowledge, skill and judgement of the modeller. Understanding these limitations is essential to creating results that can be applied successfully in practice. 8.3 Future Research The relatively simple mathematics of this model makes this an ideal method for large scale hazard mapping. Many of the existing GIS integrated rockfall prediction tools intended for large scale hazard mapping, such as STONE or Rockfall Analyst, use over-simplified representations of the impact process, and are often calibrated by matching the distribution of run out paths only. This model provides an opportunity to perform more realistic simulations of rock 203 fall trajectories and kinematics, while maintaining the ability to perform a large number of simulations in a reasonable time frame. There has been a significant amount of work done on the study of fracture mechanics applied to underground excavations; however the fracture of rocks during rock fall impacts remains poorly understood. Due to the complexity of the problem, this could be another area where a simplified physical model combined with stochastic elements could provide useful engineering tools. The model described here is highly sensitive to the roughness angle. Two distributions, uniform and truncated normal, have been tested, but other distributions could be used. A broader search of roughness angle distributions may result in improved model performance. A linear reduction in the roughness scale value has been developed here to model the ability of larger rocks to over-ride surficial irregularities and achieve longer runout distances. More rigorous testing of this assumed linear relationship could be performed to optimize it further. Another important area for future investigation is in regards to the particle shape. Intuitively, more irregular particles are likely to have more erratic trajectories relative to particles with a more regular geometry. Defining two roughness distributions, one for the substrate roughness and the other for the particle could help to represent this. The roughness angle and particle shape factor are currently considered independent random variables. Considering correlation of these two may improve the representation of extreme events, such as a rock impacting on its long axis and achieving a high bounce height. Another phenomenon not considered in this model are gyroscopic effects, where centripetal forces can stabilize a rotating block, allowing it to achieve an abnormally long runout. This would require the shape factor to be correlated with the velocity, with increasing velocity tending to result in a larger shape factor. 204 One of the hindrances in improving the representation of shape and rotational effects in rock fall modeling is the poor understanding of rotational motions associated with rock falls. The relatively small subset of rock fall experiments with detailed kinematic information becomes even smaller when looking at studies of rotational motions as well. Most protection systems are tested and designed for linear impacts only. The Ehime dataset shows that the angular kinetic energy can be up to 40% of the linear kinetic energy, which means that nearly 30% of the falling rock’s total energy is simply neglected for most design methods. Careful study of rotational motions, the factors affecting rotational motions, and integration of this information into design practice remain an under-explored area of rock fall engineering. There are still many questions to be answered regarding the behavior of rock falls, and this will remain an active area of research for some time to come. This model has built on the current state of knowledge in the field in an attempt to create a useful engineering tool. As new information becomes available, the methodologies presented here will still be useful as a means to evaluate new model features or further optimize the model. 205 Bibliography Abbruzzese, J. M. and Labiouse, V. 2014. New Cadanav methodology for quantitative rock fall hazard assessment and zoning at the local scale. Landslides, 11: 551-564. Agliardi, F. and Crosta, G. 2003. High resolution three-dimensional numerical modelling of rockfalls. International Journal of Rock Mechanics and Mining Sciences, 40: 455–471. Andrew, R., Hume, H., Bartingale, R., Rock, A., Zhang, R. 2012. CRSP-3D user’s manual, Colorado rockfall simulation program. Federal Highway Administration, Central Federal Lands Highway Division, Lakewood, CO. FHWA-CFL/TD-12-007. Ashayer, P. 2007. Application of rigid body impact mechanics and discrete element modeling to rockfall simulation. Ph.D. thesis, Graduate Department of Civil Engineering, Toronto, ON. Asteriou, P., Saroglou, H., Tsiambaos, G. 2012. Geotechnical and kinematic parameters affecting the coefficients of restitution for rock fall analysis. International Journal of Rock Mechanics & Mining Sciences, 54: 103-113. Azimi, C. and Desvarreux, P. 1977. Calcul de chutes de blocs et vérification sur modèle réduit, Internal Technical Report, ADRGT. Azzoni, A., la Barbera, G., Zaninetti, A. 1995. Analysis and prediction of rockfalls using a mathematical model. International Journal of Rock Mechanics & Mining Sciences, 32: 709-724. Berger, F., and Dorren, L. 2006. Objective comparison of rockfall models using real size experimental data. In: Disaster mitigation of debris flows, slope failures and landslides. Edited by H. Marui. Universal Academy Press, Inc, Tokyo, Japan, pp. 245–252. 206 Bourrier, F. 2008. Modélisation de l’impact d’un bloc rocheux sur un terrain naturel, application à la trajectographie des chutes de blocs. Ph.D. thesis, Institut Polytechnique de Grenoble, FR. Bourrier, F., Dorren, L., Nicot, F., Berger, F., Darve, F. 2009a. Toward objective rockfall trajectory simulation using a stochastic impact model. Geomorphology, 110: 68-79. Bourrier, F., Eckert, N., Nicot, F., and Darve, F. 2009b. Bayesian stochastic modeling of a rock bouncing on a coarse soil. Natural Hazards and Earth Systems Sciences, 9, 831-846. Bourrier, F., Berger, F., Tardif, P., Dorren, L., Hungr, O. 2012. Rockfall rebound: comparison of detailed field experiments and alternative modelling approaches. Earth Surface Processes and Landforms, 37: 656-665. Bourrier, F. and Hungr, O. 2012. Rockfall dynamics: A critical review of collision and rebound models. In Rockfall Engineering. Editied by S. Lambert and F. Nicot. John Wiley & Sons, Inc., Hoboken, NJ, USA. pp 175-209. Bozzolo, D. and Pamini, R. 1986. Simulation of rockfalls down a valley side. Acta Mechanica, 63: 113–130. Bruckno, B. 2010. Baseline rockfall rates and rockfall protection in Virginia, Ph.D. thesis, University of Nebraska, Lincoln, NB. Buzzi, O., Giacomini, A., Spadari, M., Fityus, S. 2011. Numerical modeling of a rock fall mesh perforation upon impact. Proceedings of the 13th International Conference of the IACMAG 2011. Sydney, Austrailia: 1141 – 1146. CEMAGREF. 2003. Réalisation d’un test d’étalonnage des modèles de trajectographie en utilisant des données provenant d’expérimentations grandeur nature. Département 207 Gestion des Territoires, Unité de recherche Ecosystèmes Paysages Montagnards, Grenoble, FR. Chau, K. T., Wong, R. H. C., Wu, J. J. 2002. Coefficient of restitution and rotational motions of rockfall impacts. International Journal of Rock Mechanics & Mining Sciences, 39: 69-77. Cleveland, W. S., Grosse, E., Shyu, W. M. 1992. Local Regression Models, In Statistical models in S. Edited by J. M. Chambers and T. J. Hastie. Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, California. Descoeudres, F., and Zimmermann, T. 1987. Three dimensional dynamic calculation of rock-falls. In: Proceedings of the Sixth International Congress of Rock Mechanics, Montreal, Canada, 337-342. Dorren, L.K.A., Berger, F., Le Hir, C., Mermin, E., Tardif, P. 2005. Mechanisms, effects and management implications of rockfall in forests. Forest Ecology and Management, 215(1–3): 183–195 Dorren L.K.A., 2012. Rockyfor3D (v5.1) revealed – transparent description of the complete 3D rockfall model, ecorisQ paper. http://www.ecorisq.org. Duncan, J. M. 2000. Factor of safety and reliability in geotechnical engineering. Journal of Geotechnical and Geoenvironmental Engineering, 126: 307-316. Evans, S. G. and Hungr, O. 1993. The assessment or rockfall hazard at the base of talus slopes. Canadian Geotechnical Journal, 30: 620 – 636. Falcetta, J.L. 1985. Un nouveau modèle de calcul de trajectoires de blocs rocheux. Revue Française de Géotechnique, 30: 11–17. 208 Ferrari, F., Giani, G. P., Apuani, T. 2013. Towards the comprehension of rockfall motion, with the aid of in situ tests. Italian Journal of Engineering Geology and Environment, DOI: 10.4408/IJEGE.2013-06.B-13. Forrestal, M. J. and Luk, V. K. 1992. Penetration into soil targets. International Journal of Impact Engineering, 12(3): 427-444. Glover, J.; Schweizer, A.; Christen, M.; Gerber, W.; Leine, R.; Bartelt, P., 2012: Numerical investigation of the influence of rock shape on rockfall trajectory. [Abstract] Geophys. Res. Abstr. 14: EGU2012-11022-1. Giacomini, A., Buzzi, O., Renard, B., Giani, G. P. 2009. Experimental studies on fragmentation of rock falls on impact with rock surfaces. International Journal of Rock Mechanics & Mining Sciences, 46: 708-715. Goldsmith, W. 1960. Impact: The theory and physical behaviour of colliding solids. Edward Arnold (Publishers) Ltd, London, pp. 397. GraphPad Software Inc. Interpreting results: Kolmogorov-Smironov test. Accessed online, 10 April 2014: http://www.graphpad.com/guides/prism/6/statistics/index.htm?interpreting_results_kolmogorov-smirnov_test.htm Guzzetti, F., Crosta, G., Detti, R., and Agliardi, F. 2002. STONE: a computer program for the three dimensional simulation of rock-falls. Computers & Geosciences, 28: 1079–1093. Heidenreich, B. 2004. Small- and half-scale experimental studies of rockfall impacts on sandy slopes, Ph.D. thesis, École Polytechnique Fédérale de Lausanne, CH. Hungr, O., and Evans, S.G. 1989. Engineering aspect of rockfall hazards in Canada, Geological Survey of Canada. 209 Hungr, O., Evans, S. G., Hazzard, J. 1999. Magnitude and frequency of rock falls and rock slides along the main transportation corridors of southwestern British Columbia. Canadian Geotechnical Journal, 36: 244 – 238. Hungr, O., Picarelli, L. and Leroueil, S. 2014. The Varnes classification of landslides-an update. Landslides, 11:167-194. Jaboyedoff, M., Balillifard, F., Philippossian, F., and Rouiller, J. D. 2004. Assessing fracture occurrence using “weighted fracturing density”: a step towards estimating rock instability hazard. Natural Hazards and Earth Systems Sciences, 4: 83–93. Jaboyedoff, M., Dudt, J. P., and Labiouse, V. 2005. An attempt to refine rockfall hazard zoning based on the kinetic energy, frequency and fragmentation degree. Natural Hazards and Earth Systems Sciences, 5: 621–632. Jones, C. L., Higgins, J. D., Andrew, R. D. 2000. Colorado rockfall simulation program, Version 4.0. Colorado Department of Transportation, Denver CO. CDOT-SYMB-CGS-99-1. Kobayashi, Y., Harp, E. L. Kagawa, T. 1990. Simulation of rockfall triggered by earthquakes. Rock Mechanics and Rock Engineering, 23: 1-20. Kolenprat, B. 2012. Tagbauarbeitenverordnung. Berg- und Hüttenmännische Monatshefte, 157(4): 160-164 Labiouse, V., and Heidenreich, B. 2009. Half-scale experimental study of rockfall impacts on sandy slopes. Natural Hazards and Earth System Sciences, 9: 1981-2009. Lan, H., Martin, C. D., and Lim, C.H. 2007. RockFall analyst: A GIS extension for three-dimensional and spatially distributed rockfall hazard modeling. Computers & Geosciences, 33: 262-279. 210 Lee, E. M., and Jones, D. K. 2004. Landslide Risk Assessment. Thomas Tellford Publishing, London, England. pp. 454. Leine, R. I., Schweizer, A., Christen, M., Glover, J., Bartelt, P. Gerber, W. Simulation of rockfall trajectories with consideration of rock shape. Multibody System Dynamics, 32: 241-271. Massey, C. I., Gerstenberger, M., McVerry, G., Litchfield, N. 2012. Canterbury earthquakes 2010/11 Port Hills slope stability: Additional assessment of the life-safety risk from rockfalls (boulder rolls). GNS Science Consultancy, Lower Hutt, NZ, Report 2012/214. MATLAB and Statistics Toolbox Release 2011b. MathWorks, Inc., Natick, Massachusetts, United States Peila, D. and Ronco, C. 2009. Technical note: Design of rockfall net fences and the new ETAG 027 European guideline. Natural Hazards and Earth System Sciences, 9: 1291-1298. Pfeiffer, T., and Bowen, T. 1989. Computer simulation of rockfalls. Bulletin of the Association of Engineering Geologists, 26(1): 135–146. Pichler, B., Hellmich, C., Mang, H. A. 2005. Impact of rocks onto gravel, design and evaluation of experiments. International Journal of Impact Engineering, 31: 559-578. Pierson, L., Davis, S., Van Vickle, R. 1990. The rockfall hazard rating system implementation manual. Federal Highway Administration, Washington, DC. FHWA-OR-EG-90-01. Piteau, D. R., and Clayton, R. 1977. Discussion of paper “Computerized design of rock slopes using interactive graphics” by P.A. Cundall, M.D. Voegelle and C. Fairhurst In: Proceedings of the 16th US Symposium on Rock Mechanics, 62-63, 1977. R Core Team. 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://R-project.org. Rocscience, Inc. 2010. Rocfall program downloads, URL www.rocscience.com. 211 Scott, D. W. 2011. Box-Muller transformation. Computational Statistics, 3: 177-179. Spadari, M., Giacomini, A., Buzzi, O., Fityus, S., Giani, G. P. 2011. In situ rock fall testing in New South Wales, Australia. International Journal of Rock Mechanics and Mining Science 49: 84–93. Spadari, M., Giacomini, A., Buzzi, O., Hambleton, H. 2012. Prediction of the bullet effect for rockfall barriers: a scaling approach. Rock Mechanics and Rock Engineering, 45: 131-144. Spadari, M., Kardani, M. De Carteret, R., Giacomini, A., Buzzi, O., Fityus, S. 2013. Statistical evaluation of rockfall energy ranges for different geological settings of New South Wales, Australia. Engineering Geology, 158: 57-65. Spang, R.M, and Rautenstrauch, R.W. 1988. Empirical and mathematical approaches to rockfall protection and their practical applications. In: proceedings of the Fifth International Symposium on Landslides, 1237-1243. Stevens, W. D. 1998. Rocfall: a tool for probabilistic analysis, design of remedial measures, and prediction of rockfalls. M. A. Sc. Thesis, Graduate Department of Civil Engineering, Toronto, ON. Stronge, W. J. 2000. Impacts mechanics. Cambridge University Press, Cambridge. Turner, A. K., and Duffy, J. D. 2012. Modelling and prediction of rockfall. In Rockfall: characterization and control. Edited by A. K. Turner, and R. L. Schuster. Transportation Research Board, National Academy of Sciences, Washington DC. pp. 334-406. Ushiro, T., Kusumoto, M., Shinohara, S., Kinoshita, K. 2006. An experimental study related to rock fall movement mechanism (In Japanese). Journal of Japanese Society of Civil Engineers, 62(2): 377-386. 212 Venables, W. N. and Ripley, B. D. 2002. Modern applied statistics with S, Fourth Edition. Springer, New York, NY. Vijayakumar S., Yacoub T., Ranjram M. and Curran, J.H. 2012. Effect of rockfall shape on normal coefficient of restitution. In: Proceedings of the 46th US Rock Mechanics / Geomechanics Symposium, Chicago, IL. Volkwein, A., Schellenberg, K., Labiouse, V., Agliardi, F., Berger, F., Bourrier, F., Dorren, L. K. A., Gerber, W., Jaboyedoff, M. 2011. Rockfall characteriszation and structural protection – a review. Natural Hazards and Earth Systems Sciences, 11: 2617-2651. Volkwein, A. and Klette, J., 2014. Semi-automatic determination of rockfall trajectories. Sensors, 14, 18187-18210; doi:10.3390/s141018187 Wyllie, D.C. 2014. Calibration of rock fall modeling parameters. International Journal of Rock Mechanics & Mining Sciences, 67: 170-180. Zhang, Z. X., Kou, S. Q., Jiang, L. G., Lindqvist, P. A. 2000. Effects of loading rate on rock fracture: fracture characteristics and energy partitioning. International Journal of Rock Mechanics and Mining Sciences, 37: 745-762. 213 Appendices 214 Appendix A Roughness Angle Distribution Calibration A.1 Minimum Number of Simulations to Reduce Variability The stochastic roughness angle included in the model means that the results obtained from the same input parameters will vary from simulation to simulation. The model was run with the same input parameters then the distributions of the results were compared. With a sufficient large number of simulations, the distributions should approximately converge on average, allowing more confidence in the comparison of one set of simulation results to the next. This calibration was performed using the Vaujany and Preg Quarry Section 1 datasets. The Vaujany dataset had 100 boulders grouped into 16 mass classes. The Vaujany model was run 5 times with the same input parameters to see the variation from run to run. The results of the 5 runs were then averaged. The model was then run an additional 5 times to determine if the results would have a similar range and result in similar averages. Averaging the runout distributions did not have a significant affect, i.e. the goodness of fit of any individual run was not significantly different than the average of multiple runs. However, the variation at the screen locations was much greater, especially in the predicted maximum values. Results are shown below for jump heights, Table A.1, and velocities, Table A.2. 215 Table A.1: Comparison of jump heights for the Vaujany model, effect of the number of simulations. Run 1 – 5 trials Run 2 – 5 trials Range for Individual Simulations Simulations Averaged Range for Individual Simulations Simulations Averaged Screen 1 Mean 1.15 – 1.37 1.24 1.13 – 1.44 1.23 Max 3.78 – 4.75 4.09 3.29 – 4.65 4.20 Screen 2 Mean 1.27 – 1.39 1.32 1.22 – 1.40 1.32 Max 3.84 – 4.96 4.53 3.92 – 5.25 4.54 Table A.2: Comparison of velocities for the Vaujany model, effect of the number of simulations. Run 1 – 5 trials Run 2 – 5 trials Range for Individual Simulations Simulations Averaged Range for Individual Simulations Simulations Averaged Screen 1 Mean 11.9 – 12.3 12.1 11.9 – 12.7 12.3 Max 21.8 – 25.2 23.2 20.9 – 25.0 22.6 Screen 2 Mean 13.3 – 13.8 13.7 12.8 – 13.6 13.2 Max 23.5 – 30.6 25.9 23.0 – 27.2 24.8 The range in the results for individual simulations is significant, especially when it comes to the prediction of the maximum velocity. The average values for the mean jump heights and velocities show very good agreement when 5 simulations are aggregated. There is still some 216 variation in the predicted maxima; however, the stochastic nature of the model does mean the predicted maxima will be very sensitive to single outliers. An aggregate of 5 individual simulations was selected as an acceptable minimum number of throws for the Vaujany model. This process was repeated with the Preg Quarry Section 1 model. This dataset had 25 boulders with their individual masses recorded. Initially, the results of 5 model simulations were averaged. It was found that there was still a significant amount of scatter in the predicted mean values. The number of throws per boulder was increased to 10, and the means showed much better agreement, as shown in Table A.3. There is a significant amount of scatter in the maximum runout predicted, even with ten simulations aggregated. This underscores the impact of single outliers on the interpretation of stochastic models. Table A.3: Preg Quarry results, effect of the number of simulations. Run 1 – 10 trials Run 2 – 10 trials Run 3 – 10 trials Simulations Averaged Simulations Averaged Simulations Averaged First Impact Mean 2.03 1.94 1.89 Max 4.83 3.46 3.83 Runout Mean 5.77 5.57 5.60 Max 20.3 30 25.0 Despite the variation in the maximum runout values, ten simulations were used as the minimum for a representative calibration for the Preg Quarry simulations. 217 A.2 Input Mean and Standard Deviation for Truncated Normal Distribution The Vaujany dataset was used in the calibration of the roughness angle distribution. The effect of varying the input mean, μ, and standard deviation, σ, was tested. A summary of the variation in all the distributions tested is provided in Table A.4. The means of the four roughness angle distributions shown are very similar. Increasing the value of the standard deviation does have the effect of giving the distribution more scatter. Table A.4: Summary of roughness angle distributions used in calibration. Calculated Distribution standard deviation 0.15 0.25 0.35 0.45 Distribution Mean 0.3 Mean angle - - 19.5 21.9 Max angle - - 56.7 60.6 0.4 Mean angle 19.3 20.1 22.0 22.8 Max angle 39.8 53.1 59.2 61.8 0.5 Mean angle 23.5 23.6 24.7 26.7 Max angle 45.7 52.1 59.8 63.6 The effect that the variation in the roughness distribution had on the model outputs was also examined. The first output examined was the runout distribution, shown in Figure A.1. The differences in this comparison are much more pronounced than those from the distributions of the roughness angles, with increases in both μ and σ having an effect of decreasing the runout. 218 Figure A.1: Comparison of runout distributions with varying roughness distribution inputs. Similar trends are observed for the jump height distribution and velocity distribution at Screen 1. Increasing the mean and the standard deviation lowers the mean jump height, Figure A.2, and velocity, Figure A.3. Using a higher σ value results in more scatter about the mean for both outputs. 219 Figure A.2: Jump height distributions at Screen 1 with varying roughness distribution inputs. Figure A.3: Velocity distributions at Screen 1 with varying roughness distribution inputs. 220 A.3 Size Dependent Roughness Reduction The effect of the roughness reduction equations were tested on the Vaujany model. When the runout versus mass was plotted, shown in Figure A.4, the no roughness case showed a slight tendency for larger rocks to stop sooner, When the roughness reduction parameters were increased to m = 0.3 and b = 1.36 the trend was reversed, with a slight tendency for larger rocks to have a longer runout. Figure A.4: Effect of different roughness reduction on the runout distribution for the Vaujany model with the runout averaged for each rock size. The dataset for the Vaujany experiments did not indicate the mass corresponding to each of the runout distances, thus the relationships shown in Figure A.4 cannot be compared directly to the experimental results. 221 As desired, the effect on the size dependent roughness scaling on total runout distribution, jump height and velocity distributions at Screens 1 and 2 was negligible. This is demonstrated in Figure A.5, Figure A.6, and Figure A.7 which show the effect of different roughness reduction functions while using the same reference deformation energy and roughness scale values. Figure A.5: Effect of the roughness reduction function on the runout distribution. 222 Figure A.6: Effect of the roughness reduction function on the Screen 1 jump height distribution. Figure A.7: Effect of the roughness reduction function on the Screen 1 velocity distribution. 223 Appendix B Stopping Criterion The stopping criterion was calibrated using the total duration from release to stopping for the boulders. A linear regression was done on the observed data points and the 95% prediction interval (PI) was calculated. The stopping conditions were set relative to the total energy head, which is the linear plus rotational kinetic energy divided by the boulder mass. Four values multiplied by the block radius were tested, 0.1R, 0.25R, 0.5R and 0.75R. The simulation data was plotted against the observed data and the regression lines to see how well the two agreed for the different stopping conditions. The percentage of the data exceeding the upper bound of the 95% prediction limit was calculated. If the model is reproducing the observations well there should be approximately 2.5% of the points plotting above this upper bound. 224 Figure B.1: Runout duration versus distance for Preg Section 1 for 4 different stopping criteria compared to the observed data points and the linear regression and prediction interval lines calculated from the observed points. 225 Figure B.2: Runout duration versus distance for Preg Section 3 for 4 different stopping criteria compared to the observed data points and the linear regression and prediction interval lines calculated from the observed points. 226 Figure B.3: Runout duration versus distance for Preg Section 4 for 4 different stopping criteria compared to the observed data points and the linear regression and prediction interval lines calculated from the observed points. The simulation values were compared to the upper bound of the 95% prediction interval to determine what percentage exceeded that value. The simulated durations relative to the observed varied from section to section, with the durations tending to be largest for Section1 and 227 smallest for Section 4. The percentage of the simulated values for each section that exceeded the 95% prediction interval is summarized in Table B.1. Table B.1: Comparison of simulation results to the 95% prediction interval. Case Percentage of Simulations Exceeding 95% Prediction Interval By Stopping Criterion 0.1*R 0.25*R 0.5*R 0.75*R Section 1 34% 6% 0% 0% Section 3 2.4% 0% 0% 0% Section 4 0% 0% 0% 0% Average 13% 2% 0% 0% The effect of the stopping criteria was also checked relative to the predicted runout distributions. As can be seen in the example from Section 1, shown in Figure B.4, stopping criteria did not have a noticeable effect on the runout distributions, or the calculated 95th percentile runout. This indicates that main effect of the different stopping criteria is in the elimination of small bounces near the end of the trajectory. 228 Figure B.4: Effect of the stopping criteria on the simulated runout distributions for Preg Section 1. The P95 value for each is shown by the dashed vertical line. 229 Appendix C Ehime Additional Calibration Plots The following plots show the ratio of simulated linear to angular kinetic energy as well as the envelopes given by Ushiro et al. (2006) for the mean and upper and lower bounds for this ratio. Figure C.1: Ehime sphere dataset: a) without the inclusion of a shape factor, and b) with a shape factor uniformly distributed between 1.0 and 1.2. 230 Figure C.2: Ehime natural rock dataset, a) without the inclusion of a shape factor, and b) with a shape factor uniformly distributed between 1.0 and 1.5. The following plots show: a) the observed and simulated points plotted against elevation for the dataset given, b) the distribution of residuals about the loess curve fit to the observed data. c) distribution of residuals from the simulation relative to the loess curve fit to the observed data. The solid and dashed lines on b) and c) indicate the mean and ± 2 standard deviations on the residual distributions. Figure C.3: Plots for incident velocities vs elevation for the Ehime sphere dataset. 231 Figure C.4: Plots for rebound velocities vs elevation for the Ehime sphere dataset. Figure C.5: Plots for rebound angular velocities vs elevation for the Ehime sphere dataset. 232 Figure C.6: Plots for incident velocities vs elevation for the Ehime cube dataset. Figure C.7: Plots for rebound velocities vs elevation for the Ehime cube dataset. 233 Figure C.8: Plots for rebound angular velocities vs elevation for the Ehime cube dataset. Figure C.9: Plots for incident velocities vs elevation for the Ehime natural rock dataset. 234 Figure C.10: Plots for rebound velocities vs elevation for the Ehime natural rock dataset. 235 Appendix D Statistical Functions for Reliability Calculations In order to calculate the reliability of the model results shown in Section 4.7, a method of assigning probabilities to values not within the observed dataset was required. Visual inspection of the distributions shown in Figure 4.41 and Figure 4.42 indicate the data is skewed. Log-normal and Weibull distributions were fit to the observed distributions for Screens 1 and 2 using the “fitdistr” function in R (R Core Team, 2013, Venables and Ripley, 2002). These two distributions were selected because they have a skewed shape, and are commonly used in practice. A one-sample KS test was used to evaluate how well each of these distributions matches the observed data, summarized in Table D.1. If the KS test indicates a low test statistic value, D, with a high significance level, p-value, the hypothesis that the observed data set tested belongs to the distribution fit to it cannot be rejected. Table D.1: KS test results for log-normal and Weibull distributions fit to the Vaujany screen data. KS Test Results Log-normal Weibull Test statistic D P-value Test statistic D P-value Screen 1 Jump Height 0.298 9.88E-8 0.274 1.28E-6 Velocity 0.108 0.214 0.109 0.206 Energy 0.168 0.00937 0.129 0.0840 Screen 2 Jump Height 0.251 1.38E-5 0.235 6.19E-5 Velocity 0.060 0.887 0.109 0.216 Energy 0.108 0.223 0.109 0.215 236 Comparing two distributions, they show relatively similar D values and p-values. Neither test results in a high significance level for the jump height distribution, indicating a low confidence in the fit. The upper tail of the kinetic energy distribution of the blocks is the primary concern for assessing the impact capacity of a protection structure, so this was examined to determine the statistical distribution used in the reliability calculations. The kinetic energy cumulative distribution was plotted with both the log-normal and Weibull distributions to visually compare the goodness of fit for the upper 10% of the distribution, shown in Figure D.1. Figure D.1: Comparison of log-normal and Weibull probability functions for the upper 10% of observed kinetic energy at Screens 1 and 2. It can be seen that the Weibull function is a better fit to the observed data for the upper 10% of the observed kinetic energies; as a result, the Weibull distribution was selected for the 237 reliability calculations. For consistency, the Weibull distribution was also used for the jump height reliability calculations, despite the low significance level of the KS-test. After finding the SEL, MEL and design interception height at the Screen 1 and 2 locations, the Weibull function fit to the observed kinetic energy and jump height distributions was used to predict the probability of an impact exceeding these design values. A Weibull distribution is described using a shape coefficient, λ, and a scale coefficient, k. The probability of an impact exceeding the design conditions, dc, also known as the reliability, is designated R, is calculated using Equation D.1. Equation D.1: 𝑅 = 1 − 𝑒−(𝑑𝑐λ)𝑘 The values of λ and k obtained from the distribution fitting process used for the reliability calculations are summarized in Table D.2. Table D.2: Summary of shape and scale parameters for Weibull distributions fit to the observed data. λ k Screen 1 Jump Height 1.55 1.37 Kinetic Energy 258.59 1.72 Screen 2 Jump Height 1.81 1.23 Kinetic Energy 246.77 1.42
Thesis/Dissertation
eng
Geological Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivs 2.5 Canada
Graduate
Computer modeling of rock fall : improvement and calibration of the PIERRE 2 stochastic rock fall simulation program
Text