UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Influence of a weak layer on depth of rock failure in underground limestone mines in high horizontal… Metcalfe, Joshua Robert 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2016_february_metcalfe_joshua.pdf [ 6.59MB ]
Metadata
JSON: 24-1.0223119.json
JSON-LD: 24-1.0223119-ld.json
RDF/XML (Pretty): 24-1.0223119-rdf.xml
RDF/JSON: 24-1.0223119-rdf.json
Turtle: 24-1.0223119-turtle.txt
N-Triples: 24-1.0223119-rdf-ntriples.txt
Original Record: 24-1.0223119-source.json
Full Text
24-1.0223119-fulltext.txt
Citation
24-1.0223119.ris

Full Text

  INFLUENCE OF A WEAK LAYER ON DEPTH OF ROCK FAILURE IN UNDERGROUND LIMESTONE MINES IN HIGH HORIZONTAL STRESS ENVIRONMENTS  by  Joshua Robert Metcalfe B.Sc in Geological Engineering, Queen’s University, 2007  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)  December 2015   © Joshua Robert Metcalfe, 2015ii  Abstract	High  horizontal  stress  environments  have plagued  the  aggregate  industry  in  shallow  room  and pillar mining operations, particularly  in  the northeastern United States.   With population density  increasing around  Southern  Ontario  and  environmental  regulations  becoming  more  stringent,  it  appears  that Ontario’s aggregate industry could be looking underground in the near future as the source of material to  meet  the  ever‐increasing  demand.  Given  that  the  near  surface  horizontal  stress  conditions  in Southern Ontario are uniquely high (H/V = 4‐6 to depths of 200m; Lo, 1978),  it can be expected that slabbing and buckling failures observed in similar mining operations with lower stress regimes in the U.S. will be exacerbated  in Southern Ontario.    In an effort  to be proactive with  this expected geotechnical design  issue,  a  distinct  element  analysis  using UDEC was  carried  out  to  understand  the mechanisms driving  failure  in  stratified  rock  environments  under  high  horizontal  stress  conditions  as  well  as  to observe the impact of the high horizontal stress on the maximum depth of failure in the roof where the roof is composed of limestone interbedded with a weak shale layer. Accordingly,  116 models  representing  a  variation  of  rock mass  conditions  subjected  to  stress  ratios (H/V) ranging from 1 to 4 were simulated.  A Voronoi tessellation was used to represent the intact rock mass  directly  above  the  excavation  so  that  the  failure  profile  through  intact  rock  could  be  explicitly modelled. Key conclusions from the modelling were as follows: 1) A shale layer within a defined distance from the roof of the excavation could increase the depth of  failure  to  three  times  what  would  normally  be  estimated  for  stratified  rock  masses  of limestone only. 2) Failure  driven  by  the  influence  of  a weak  shale  interbed  occurs  through  diagonal  fracturing, supporting an experimental conclusion published by Stimpson and Ahmed (1992), as opposed to slabbing  or  buckling.  Slabbing  and  buckling  are  the  common  failure mechanisms  in  stratified rock masses without weak interbeds. Therefore, it is critical to understand the interbedded nature of the rock mass comprising the roof over a mine opening so that a proper ground support design can be developed. iii  Preface	The topic of the research presented in this thesis was developed in collaboration with Golder Associates based  on  industrial  challenges  observed  by  the  author  and  work  colleagues  during  his  previous employment with  the  company.   The  research program  for  the  thesis was designed,  carried out and analyzed  solely  by  the  author.    The  case  studies  discussed  in  Chapter  5 were  facilitated  by  Golder Associates as well.   iv  Table	of	Contents	Abstract ......................................................................................................................................................... ii Preface ......................................................................................................................................................... iii Table of Contents ......................................................................................................................................... iv List of Tables .............................................................................................................................................. viii List of Figures ............................................................................................................................................... ix List of Symbols ............................................................................................................................................. xi Acknowledgements ..................................................................................................................................... xii CHAPTER 1  INTRODUCTION ..................................................................................................................... 1 1.1  Thesis Structure ............................................................................................................................ 3 CHAPTER 2  LIMESTONE AND THE HISTORY OF ITS MINING .................................................................... 5 2.1  Limestone, Its Uses and the Importance of Its Impurities ............................................................ 5 2.2  Mining Methodology .................................................................................................................... 6 2.3  Underground Mining as a Viable Option for Limestone Extraction .............................................. 8 CHAPTER 3  LIMESTONE MINING POTENTIAL IN SOUTHERN ONTARIO ................................................. 10 3.1  Regional Geology ........................................................................................................................ 10 3.2  Current State of Limestone Extraction in Southern Ontario ...................................................... 14 3.3  Potential Target Formations for Limestone and Dolostone Mining ........................................... 16 CHAPTER 4  MINING CONDITIONS IN SOUTHERN ONTARIO .................................................................. 19 4.1  Southern Ontario In‐Situ Stress Regime ..................................................................................... 19 4.2  Structural Characteristics of Southern Ontario Sedimentary Sequences ................................... 23 4.3  Geotechnical Experience from the US Aggregate Mining Industry ............................................ 24 CHAPTER 5  FAILURE MECHANISMS ....................................................................................................... 28 5.1  Types of Rock Failure .................................................................................................................. 29 5.1.1  Stress Controlled Failure ..................................................................................................... 30 v  5.1.2  Structurally Controlled Failure ............................................................................................ 31 5.1.3  Compound Failure Mechanisms ......................................................................................... 31 5.2  Modes of Failure and Failure Mechanisms in Stratified Rocks ................................................... 32 5.2.1  The Arching Phenomenon ................................................................................................... 35 5.2.2  Slabbing Failure ................................................................................................................... 36 5.2.3  Buckling Failure ................................................................................................................... 37 5.2.4  Crushing Failure .................................................................................................................. 38 5.2.5  Diagonal Fracturing Failure ................................................................................................. 39 5.2.6  Effects of High Horizontal Stresses ..................................................................................... 40 5.3  Failure Mechanisms from In‐Situ Observations .......................................................................... 41 5.3.1  Mine A ................................................................................................................................. 41 5.3.2  Mine B ................................................................................................................................. 44 5.4  Summary ..................................................................................................................................... 46 CHAPTER 6  NUMERICAL MODELLING .................................................................................................... 47 6.1  Selection of a Suitable Model Method ....................................................................................... 47 6.2  UDEC ........................................................................................................................................... 49 6.2.1  Voronoi Tessellation Generator .......................................................................................... 50 6.3  Two Dimensional Versus Three Dimensional Modelling ............................................................ 51 CHAPTER 7  MICRO‐MECHANICAL PROPERTIES CALIBRATION............................................................... 53 7.1  Intact Rock Mass Properties of Limestone and Shale in Southern Ontario ................................ 54 7.2  Equivalent Mohr‐Coulomb Failure Criterion Parameters ........................................................... 55 7.3  Developing a Properties Range using a UCS Simulation ............................................................. 56 7.3.1  UCS Modelling Conditions ................................................................................................... 57 7.3.2  Strength Overestimation by Using Elastic Blocks ................................................................ 58 7.3.3  Measuring the Numerically‐Simulated UCS ........................................................................ 59 7.3.4  Monitoring Shear and Tensile Fracture Development ........................................................ 60 vi  7.3.5  UCS Testing Results ............................................................................................................. 63 7.3.6  Relationship between Voronoi Edge Length and Failure Mechanism Bias ........................ 65 7.3.7  Dependency of Strength Properties on Voronoi Edge Length ............................................ 66 7.3.8  Suitable Baseline of Voronoi Contact Properties ................................................................ 67 7.4  Calibration Conclusions ............................................................................................................... 70 CHAPTER 8  ENGINEERING MODEL ......................................................................................................... 73 8.1  Model Geometry ......................................................................................................................... 73 8.2  Boundary Conditions ................................................................................................................... 76 8.3  Voronoi Edge Length ................................................................................................................... 76 8.4  Rock Mass and Joint Properties for Numerical Modelling .......................................................... 77 8.5  Variable Parameters.................................................................................................................... 78 8.5.1  Stress Field .......................................................................................................................... 78 8.5.2  Bedding Spacing .................................................................................................................. 80 8.5.3  Distance of Shale Layer from Excavation Roof ................................................................... 82 8.6  Conditional Assumptions ............................................................................................................ 83 8.6.1  Influence of Water .............................................................................................................. 83 8.6.2  Pillar Stability ...................................................................................................................... 83 8.7  Discretization of Discrete Blocks ................................................................................................. 83 8.8  Establishing an Equilibrated Model ............................................................................................ 84 8.9  Validation Model against Current Literature .............................................................................. 86 CHAPTER 9  NUMERICAL MODELLING RESULTS ..................................................................................... 89 9.1  Base Case Models without a Shale Layer .................................................................................... 89 9.2  Shale Layer Models ..................................................................................................................... 91 9.2.1  Shale Layer Modelling Observations ................................................................................... 94 9.3  Shale Influence Boundaries ......................................................................................................... 94 9.4  Mechanisms and Modes of Failure ............................................................................................. 96 vii  9.4.1  Shale Layer Influence on Stress Redistribution ................................................................... 98 9.5  Limitations of Numerical Modelling Results ............................................................................. 100 CHAPTER 10  CONCLUSIONS AND RECOMMENDATIONS ................................................................... 102 10.1  Influence of a Weak  Layer on  the Depth of Failure  in Stratified Rock within High Horizontal Stress Environments ............................................................................................................................. 102 10.2  Key Contributions ...................................................................................................................... 103 10.3  Possible Future Research .......................................................................................................... 103 10.3.1  Snaking Failures ................................................................................................................ 103 10.3.2  Ground Support Design Guidelines ................................................................................... 104 10.3.3  Strength Differential Requirements .................................................................................. 105 Bibliography .............................................................................................................................................. 106 Appendices ................................................................................................................................................ 113 Appendix A      Summary of Roof Failure Numerical Model Results ..................................................... 113 Appendix B      Summary of UCS Test Numerical Modelling Results ..................................................... 117 Appendix C      Computer Code Generated Using the FISH Language .................................................. 122 C.1  Roof Failure in Limestone Mine Simulation Code ................................................................. 122 C.2  UCS Test Simulation Code ..................................................................................................... 130   viii  List	of	Tables	Table 1.  Quarried Formations of Southern Ontario ................................................................................... 15 Table 2.  Summary of Stress Ratios by Depth below Ground Surface ........................................................ 21 Table 3.  Summary of Underground Mining Dimensions from US Limestone Mines ................................. 26 Table 4.  Classifications of Roof Failure during Survey (Esterhuizen et al., 2007) ...................................... 27 Table 5.  Overview of Modelling Method Groups....................................................................................... 48 Table 6.  Range of Intact Rock Mass Properties for Target Ordovician Formations ................................... 54 Table 7.  Corresponding Mohr‐Coulomb Properties from RocLab ............................................................. 55 Table 8.  Model Elastic Properties .............................................................................................................. 58 Table 9.  Constant Voronoi Contact Properties (Referred to as Joint Properties in UDEC) ........................ 58 Table 10.  Corresponding Voronoi Contact Properties Established from UCS Simulations ........................ 70 Table 11.   Summary of Underground Mining Dimensions from US Limestone Mines  (Esterhuizen et al., 2007) ........................................................................................................................................................... 73 Table 12.  Selected Micro‐Mechanical Properties for Voronoi Contacts .................................................... 76 Table 13.  Rock Mass and Joint Properties Used in Numerical Modelling .................................................. 77 Table 14.  Absolute Stress Values for each Stress Ratio Simulated ............................................................ 79 Table 15. Bedding Spacing and Shale Layer Thickness for Numerical Modelling Simulations ................... 82 Table 16.  Summary of Main Failure Mode Classification........................................................................... 96 Table 17.  Summary of Shale Layer Modelling Results ............................................................................. 113 Table 18.  Summary of UCT Test Numerical Modelling Results ................................................................ 117  ix  List	of	Figures	Figure 1.  Schematic of a Typical Room and Pillar Mining Operation (Hamrin, 2001).................................. 7 Figure 2.  Principal Paleozoic and Mesozoic Tectonic Elements of Ontario (Johnson et al., 1992) ............ 11 Figure 3.  Paleozoic Geology of Southern Ontario (Birchard et al., 2004) .................................................. 13 Figure  4.    Regional  In‐Situ  Stress  Ratios  as  a  Function  of  Depth:  (a)  H/V,  (b)  h/V,  and  (c)  H/h (Nuclear Waste Management Organization and AECOM Canada Ltd., 2011) ............................................ 22 Figure 5.  In‐Situ Stress Measurements at Darlington Nuclear Power Plant (Haimson & Lee, 1980) ........ 23 Figure 6.   Structural Data from across Southern Ontario: (a) Kincardine, ON (Cruden, 2011), (b) Lincoln, ON (Gartner Lee Limited, 1996) .................................................................................................................. 24 Figure 7.  North American In‐Situ Stress Map (adapted from Heidbach et al., 2010) ................................ 25 Figure 8.   General Rock Failure Mechanisms  for Underground Excavation  (adapted  from Martin et. al, 1999) ........................................................................................................................................................... 29 Figure 9.  Voussoir Beam Criteria ................................................................................................................ 33 Figure 10.  Specific Failure Modes in Stratified Rock Masses ..................................................................... 33 Figure 11.  Conditions Leading to Different Failure Modes in Stratified Rock Masses ............................... 34 Figure 12.  Compression Arches Composed of (a) Single Voussoir Beam, (b) Stack of Voussoir Beams .... 36 Figure 13.  Slabbing Failure ......................................................................................................................... 37 Figure 14.  Buckling Failure ......................................................................................................................... 38 Figure 15.  Crushing Failure......................................................................................................................... 38 Figure 16.  Diagonal Fracturing Failure ....................................................................................................... 39 Figure 17.  Stress Concentration and Stress Shadows around a Roof Fall (Iannacchione et al., 2001) ...... 41 Figure 18.  Large Roof Fall which has been Supported with Steel Straps to Prevent Further Failure ........ 42 Figure 19.  Weak Interbeds Controlling Failure Profile ............................................................................... 43 Figure 20.  Second Large Roof Fall in Mine A .............................................................................................. 43 Figure 21.  Large Roof Fall in Mine B ........................................................................................................... 45 Figure 22.  Second Large Roof Fall in Mine B .............................................................................................. 46 Figure 23.   Visual of Directional Bias due  to  Iteration Number:  (a)  Iteration Number = 1,  (b)  Iteration Number = 100, and (c) Iteration Number = 10,000. ................................................................................... 51 Figure 24.  Example of UCS Test Model for Voronoi Contact Property Calibration (Voronoi Edge Length = 0.005; 4% of sample length). ...................................................................................................................... 56 Figure 25.  Axial Stress Measurement Locations along Platen‐Sample Interface ....................................... 59 x  Figure 26.  Number of Cycles before Failure .............................................................................................. 60 Figure 27.  Rejected Failures: (a) UCS Test #32, (b) UCS Test #40 .............................................................. 62 Figure 28.  UDEC UCS Failure Simulation Compared to Laboratory Failed Core Samples .......................... 63 Figure 29.  UCS Test #47 Failure Monitoring: (a) Axial Stress ‐ Axial Strain Curve, (b) Contacts in Failure 64 Figure 30.  Dependency of Mechanism Initiating Cracking on Voronoi Edge Length Ratio ....................... 65 Figure 31.  Effect on Strength by Voronoi Edge Length .............................................................................. 67 Figure 32.  Contact Tensile Strength Required to Initiate Suitable Tensile Cracking .................................. 68 Figure 33.  Sensitivity Analyses: (a) Tensile Strength, and (b) Cohesion .................................................... 69 Figure 34.  Sensitivity Analysis on Friction Angle ........................................................................................ 70 Figure 35.  Geometry and Boundary Dimensions for the Numerical Modelling ........................................ 74 Figure 36.  Voronoi Tessellation Regions .................................................................................................... 75 Figure 37.    Summary of Ordovician  Limestone Bedding  and  Interbed  Literature Review:  a)  Limestone Bedding Spacing, b) Shale Interbed Thickness ............................................................................................ 81 Figure 38.  Unbalanced Forces Indicating Model Instability Negligible ...................................................... 85 Figure  39.   Displacement Vectors  Showing Movement  Less  Than  5mm  (Within Acceptable  Range  for Model Equilibrium) ..................................................................................................................................... 85 Figure  40.    Span  Charts  for  Preventing  Buckling  (Snap‐Thru)  and  Crushing  Failure  (Hutchinson  & Diederichs, 1996) ........................................................................................................................................ 87 Figure 41.  Example of the Progression of Buckling Failure: (a) Failure in Motion, (b) Complete Failure. . 88 Figure 42.  Base Case Simulations (without Shale Layer) ........................................................................... 89 Figure 43.  Empirical Ground Support Design Chart (Grimstad and Barton, 1993) .................................... 90 Figure 44.  Modelling Results for k=1 ......................................................................................................... 91 Figure 45.  Modelling Results for k=2 ......................................................................................................... 92 Figure 46.  Modeling Results for k=3 .......................................................................................................... 92 Figure 47.  Modelling Results for k=4 ......................................................................................................... 93 Figure 48.  Shale Layer Influence Chart ...................................................................................................... 95 Figure 49.   Examples of Failures Observed: (a) Failure to Shale Layer by Diagonal Fracturing, (b) Failure with No Shale Layer Present by Buckling, (c) Failure Not Extending to Shale Layer by Buckling ............... 97 Figure 50. Conditions Leading to Different Failure Modes in Stratified Rock Masses ................................ 98 Figure 51.  Stress Distribution around Room: (a) Stable Arch (No Shale Layer), (b) Shale Layer Failure with Stable Arch Attempting to Form Below Shale, (c) Shale Layer Failure ....................................................... 99 Figure 52.  Snaking Failure Observed in US Limestone Mine (Esterhuizen et al., 2008) .......................... 104 xi  List	of	Symbols	c  Cohesion   E  Young’s Modulus G  Shear Modulus K  Bulk Modulus k  In‐Situ Stress Ratio kN  Joint Normal Stiffness kS  Joint Shear Stiffness  Internal Friction Angle  Specific Weight   Poisson’s Ratio   Density 1  Major Principal Stress 2  Intermediate Principal Stress 3  Minor Principal Stress C  Unconfined Compressive Strength H  Major Horizontal Stress h  Minor Horizontal Stress T   Tensile Strength V  Vertical Stress  xii  Acknowledgements	This  research was  funded by  the National  Sciences  and  Engineering Research Council of Canada  and sponsored by Golder Associates Limited. Over the time it took me to carry out my research and finally put my ideas into word form, I have met a myriad of people  that, whether  they know  it or not, have  influenced  the outcome of  this  thesis and  I would like to acknowledge them: Dr. Erik Eberhardt and Dr. Oldrich Hungr for their constructive feedback and patience.  The staff of the Department of Earth, Ocean and Atmospheric Sciences at UBC for helping me through the technicalities of being a graduate  student. My  friends and  colleagues at Golder Associates  for being  flexible  in my work commitments  to provide me with  the  time off  to carry out my  research as well as coordinating learning opportunities for me at operating room and pillar mines in the United States.  I do not think any other  company would have been  so generous during  the  challenging economic  climate  in  the mining industry.  I am incredibly indebted to Laura Mason, Aynsley Neufeld, Audra Vair, Geidy Baldeon, Mika McKinnon, Lindsay Callan, Neda Zangeneh, Lindsay Tetreault and Kathi Unglert.   I am sure that none of you know this because  I am not one  to use my words  to express emotion but all of you at one  time or another were my biggest motivation to push on with my research and remind me that I am capable of finishing this thesis. Of course, I would  like to thank my family and second family – Will and Cathy Pitman for their endless support, even though  I am sure  that at some point each one of them had severe doubts that  I would ever graduate.   Whether  it was nagging, reading me the “riot” act,  lending an ear or simply being the whiteboard in which I could untangle the rat’s nest my brain was in so that I could move forward, I am extremely grateful to have you all as a support network. Finally,  I would  like  to  acknowledge  all of  the  guys who  came  in  and out of my  life  throughout  this adventure.  You provided for one heck of a distraction but also made these years the toughest learning experience of my life which I hope never to endure again.  However, Steve Dmytrasz, I truly appreciate all your love and support as I gallivanted across Canada for many years to complete this thesis.  It is too bad that you did not get to be part of the end of the journey you supported at the start. 1  CHAPTER	1 INTRODUCTION	Limestone is a highly valued resource required for cement and aggregate production. Preference is given to mining  it  from  surface  (quarrying) as  this  is generally  the most cost effective option. Even  though, worldwide, underground mining of  limestone already contributes  to  the supply of aggregate, building stone  and  lime,  Canada  has  yet  to  develop  its  first  underground  limestone  mine  because  of  the abundance of quality limestone and dolostone resources available from surface; particularly in Southern Ontario around the great  lakes.   However, with  increasing population growth, the depletion of readily accessible  resources, and  societal concerns  for preserving environmentally  sensitive areas  resulting  in stricter  regulations  for  the  opening  of  new  quarry  operations,  underground  mining  is  now  being considered as an alternative option for accessing and developing new limestone sources.   Limestone deposits are located near existing or past water bodies (nature of deposition), which is where 90% of  the population of Ontario  lives  (Ontario Geological Survey, 1992).   The Greater Toronto Area (GTA)  limestone  and  dolostone  industry  has  been  the most  impacted  because  it  holds  the  largest aggregate  and  construction  materials  markets  whilst  population  densities  and  environmental regulations prevent quarry expansions and/or new quarry permitting close to these markets.  For a low value, bulk commodity such as  limestone, transportation costs represent a major part of the product’s end  cost  and  therefore  distance  to market  has  a  significant  influence  on  the  economic  feasibility  of aggregate production.   Consequently, the greater transportation costs for supplying these materials to the GTA market  result  in higher material prices, and  therefore, as documented by Thompson  (1982), local underground mines in this area have potential to be more economical than surface quarries which are up to three times the distance from the market. Underground mining  of  stratified  deposits,  such  as  limestone  and  dolostone,  is  typically  carried  out using  the  room  and  pillar mining method  and  presents  different  engineering  challenges  than  those common in quarrying.  Experience from the United States of America (US) suggests that roof failures via slabbing  or  buckling  constitute  the most  prominent  hazard  and  stability  concern  (Esterhuizen  et  al., 2008).  This is related to the geology of limestone and other sedimentary rock deposits, which typically respond to applied stress like beams.  Much research has been completed to date on the development of a stress arch in homogenous, stratified rock masses and Voussoir beams (Diederichs, 1999; Ran et al., 1994; Beer & Meek, 1982; Evans, 1941)  to assist with establishing suitable ground support  for mining operations by understanding the probable depth of roof  failure  in these environments.   However, the 2  inclusion of a weak  layer  (creating a heterogeneous, stratified  rock mass) seems  to govern  rock mass failure mechanisms  resulting  in  a  depth  of  failure  beyond  the  anticipated  depth  of  failure  predicted using  the  stress  arch  theory. As  such,  ground  support  approaches  in  these environments  are  at best moderately successful  in preventing roof  failures as observed  in the US.   This degree of uncertainty  is not acceptable/desirable from the perspective of safety and resource isolation; large failure can prevent access  to  high  quality  stone.    Additionally,  geotechnical  misunderstanding  resulting  in  on‐going remediation  to  deal  with  ground  failures  significantly  increases  the  operating  costs  for  the mining operation and decreases the feasibility of the project. Exacerbating  this  situation  in  Southern  Ontario  is  that  the  area  is  known  to  have  abnormally  high horizontal  in‐situ  stresses.  The  relative  alignment  of  these  high  stresses  and  the weak  layer  in  the stratified rock mass increases the potential for excessive slabbing/overbreak (depth of failure) where the thickness  between  the  roof  and weak  layer  is  insufficient  to  allow  for  proper  redistribution  of  the induced stresses onto adjacent pillars.   The National  Institute  for Occupational  Safety  and Health  (NIOSH) based  in Pennsylvania has  funded numerous observational studies over the past 20 years in an effort to reduce the number of roof failures in  US  limestone  mines  under  high  horizontal  stress  conditions.    The  outcome  of  these  studies (Esterhuizen  et  al.,  2008;  Esterhuizen  et  al.,  2007;  Iannacchione  et  al.,  2001) was  a  series  of  design criteria  for  developing  a  room  and  pillar mine  so  that  the  impact  of  the  high  horizontal  stresses  is minimized by allowing suitable redistribution of stresses: 1) Increase  the  total  proportion  of  headings  driven  parallel  to  the  direction  of  the maximum horizontal stress. 2) Reduce  the number and  size of  cross‐cuts driven perpendicular  to  the direction of maximum horizontal stress. 3) Drive  cross‐cuts  into  existing  headings  instead  of  having  the  heading  and  cross‐cut  faces meeting at an intersection. 4) Offset  cross‐cuts  to  create  only  three‐way  intersections  which  are  more  stable  in  stressed environments, and so that  if  failure does develop  in the cross‐cut then  it  is truncated by a rib pillar. 5) Maintain  a wedge‐shaped mining  front  parallel  to  the  direction  of  the maximum  horizontal stress. 3  This design criteria, although having proven effective at reducing  the number of roof  failures  in  these conditions, is solely based on stress redistribution and does not account for the coupled impact of stress and  structure.   As  such,  roof  failures are  still occurring  characterized by  failure depths beyond  those expected which  indicates that a better understanding of the mechanisms generated from this coupled relationship of stress and structure in these environments is necessary.  This  thesis presents  the  research and numerical modelling completed  to understand  the  impacts of a weak  layer within  a  stratified  rock mass  (such  as  an  interbedded  limestone  and  shale  sedimentary sequence)  on  the  maximum  depth  of  failure  in  high  horizontal  stress  environments.    Based  on observations from the US where a smaller high horizontal stress regime exists than in Southern Ontario, the depth of  failure  truncates  at  the weak  layer  above  the  failure depth predicted using  stress  arch theory.  However, it is expected that at some depth the weak layer will be too far from the excavation proving  inconsequential  to  roof stability and conventional stress arch  theory will control  the depth of failure.  By  understanding  the  governing mechanisms  leading  to  increased  depth  of  failure  and  the distance  in  which  a  weak  layer  proves  inconsequential,  hopefully  more  effective  ground  support programs,  from a  safety and  cost perspective,  can be designed and  implemented when underground limestone and dolostone mining imminently begins in Southern Ontario.   1.1 Thesis	Structure	This thesis has been structured so that the thesis flows from a review of current practices, theories and specifics regarding the ground conditions in Southern Ontario in Chapters 2 through 6, and progresses to the author’s research methodology, results, conclusions and research limitations in Chapters 7 through 10.  Chapter 2 provides a general background of limestone mining discussing its purpose in the industry and how  it  has  changed  over  the  years,  the  mining  methodology,  and  limestone’s  viability  as  an underground resource.  This background is intended to highlight the importance of limestone within the aggregate  and  construction  industries  worldwide  and  hopefully  establish  an  appreciation  for  the probable direction of limestone mining in the future. Chapter 3 focusses on the potential for  limestone mining  in Southern Ontario, outlining the geological setting, current state of the surface mining (quarrying) industry and the potential target formations for underground limestone mining in the province.   4  Chapter  4  details  the mining  conditions  in  Southern Ontario  providing  the  background  geotechnical information  (in‐situ  stress  and  rock mass  structure)  in which  limestone mines  are most  likely  to  be constructed within.    This  information  sets  the  basis  for  the  conditions  represented  in  the  numerical modelling portion of this research.   Finally, international experience in which the aggregate industry of Southern Ontario  can draw upon  in preparation  for exploiting underground  limestone and dolostone resources. Chapter 5 discusses failure mechanisms, transitioning from an elementary classification of failure modes (stress  and  structurally  controlled  failures)  to  specific  modes  of  failure  observed  in  stratified  rock masses.  The effects of high horizontal stress on failure modes and the phenomenon of arching are also visited.    The  chapter  concludes with  two  case  studies  that  summarize  the  failures  observed  by  the author in limestone mines under high horizontal stresses in the northeastern United States. Chapter  6 provides  a background on  numerical modelling  and  discusses  the  selection of  the distinct element approach using the commercial software, UDEC, to carry out the research. Chapter  7  highlights  the  importance  of  model  calibration,  particularly  pertaining  to  rock  mass properties, and steps through the process of establishing suitable micro‐mechanical properties required for the Voronoi tessellation employed in the numerical modelling.  Chapter 8 details the engineering model and methodology for assessing the influence of a weak layer on the depth of failure under variable in‐situ stress ratios and bedding thicknesses. Chapter 9 summarizes the results from the numerical modelling and discusses trends observed at each bedding thickness and in‐situ stress ratio.  Limitations to the research are also detailed at the end of this chapter. Chapter 10 concludes the thesis by discussing the key outcomes and contributions of the research to the field of study.  Potential applications of the research findings to the industry are also listed.  Finally, this chapter closes with recommendations for additional research which could be developed from the work conducted as part of this thesis. 5  CHAPTER	2 LIMESTONE	AND	THE	HISTORY	OF	ITS	MINING	Underground limestone mining became popular in England during the Industrial Revolution because the limestone was  close  to  or  directly  underlying  seams  of  coal,  iron  ore  and  fireclays  that were  being heavily exploited as raw materials to make iron (Brook, 1991).  Although extraction of sedimentary rocks has been  recorded back  to Neolithic  flint mines over 10,000 years ago    (Stocks, 1979),  the  Industrial Revolution marked the first boom of  large scale room and pillar mining for  limestone.   The first North American underground limestone mine began in Kansas City, US in the late 1890’s producing aggregate as  a  construction material.    This  industry  eventually  took  off,  particularly  in  the  northeastern  states (West  Virginia,  Pennsylvania,  Michigan,  Ohio)  and  by  2002  there  were  83  operating  underground limestone mines producing 60Mt of material  for  the aggregate  and  construction  industries  in  the US  (Freas  et  al.,  2006).   Additionally,  production  grew  from  a  few  thousand  tonnes  per  year  to  several million tonnes per year (Stocks, 1979).  Many other countries worldwide have at some time employed the room and pillar mining method for sedimentary rock extraction; most notably are  large operations  in Scandinavian countries and  in Japan (Stocks,  1979).    However,  an  economic  situation  which  favours  underground  mining  over  surface quarrying  for  limestone has  yet  to present  itself  in Canada  and  therefore no underground  limestone mines have ever been developed in the country.  This chapter discusses the role of  limestone  in the construction  industry, the natural characteristics of limestone which control its feasibility as a mineral resource, room and pillar mining as the most common mining method for extracting the ore, and the studies conducted globally  indicating that underground mining of limestone is becoming an economically viable option once again. 2.1 Limestone,	Its	Uses	and	the	Importance	of	Its	Impurities		Carbonate rocks comprise the basis of the construction industry producing aggregate, cement, lime and building stone.   In the US,  it  is estimated that carbonate rocks comprise 75% of all stone quarried, and only  next  to  sand  and  gravel,  they  are  produced  in  the  greatest  quantities  out  of  the  mineral commodities (Carr et al., 1994).  In industry, limestone (CaCO3) and dolostone (CaCO3‐MgCO3) composed mostly of  the mineral  calcite  (CaCO3)  represent  the majority of  the  carbonate  rocks exploited due  to their desirable physical and chemical properties.   6  Limestone and dolostone are sedimentary rocks formed by the deposition of material within bodies of water.    Accordingly,  these  deposits  are  generally  substantial  in  surface  area  and  flat‐lying.    Most material  is  of  a  biologic  origin  (shells,  fecal material,  algae  bi‐products)  and  is  deposited  in  shallow marine  environments, where  the  energy  during  deposition  can  be  quite  variable  (Carr  et  al.,  1994).  Deposits  in  high  energy  environments  are  generally  composed  of  high  purity  carbonate material  of larger grain size and are more economic  than deposits  formed  from  low energy sedimentation, which often  are  diluted  by  fine‐grained,  non‐carbonate  impurities.  These  impurities  not  only  affect  the economic viability of a deposit, but also substantially impact the rock mass quality and the ability of the rock mass to be self‐supporting during mining. The most  common  impurities  are  clay minerals  –  kaolinite,  illite,  chlorite,  smectite  and mixed‐lattice types  –  precipitated  in  laminae  or  thin  partings  (due  to  their  silica  tetrahedral  and  alumina  and/or magnesia octahedral structure) (Carr et al., 1994).  These mineral concentrations are commonly referred to  as  shale  interbeds  amongst  the  carbonate  rocks  and  due  to  their  relatively  weaker  engineering properties form a natural plane of weakness within the rock mass.   Consequently, shale  interbeds are often the origin of rock mass instability on many scales when mining carbonate deposits.  Chert  is  another  common mineral  impurity.    Chert  is  formed  when  a  sedimentary  rock  undergoes mineral replacement as a result of diagenesis.   Therefore,  it tends to concentrate  in  lenses or nodules, instead of layers.  Although they can act as the origin of local stress concentrations and expedite intact rock  failure, chert  inclusions  rarely pose  the same engineering concerns as shale  layers during mining because they are so localized. In  summary,  not  all  limestone  or  dolostone  deposits  are  a  resource.    Their  physical  and  chemical composition have  to be desirable  for aggregate, cement,  lime and building  stone production and any impurities with the rock mass not only degrade the resource but also have potential to increase the risk of geotechnical failure.  2.2 Mining	Methodology	The selection of a mining methodology for extraction of a mineral involves understanding the nature of the  deposit  and  the  production  requirements1  to  ensure  a  profit  can  be made.    From  an  economic                                                             1  Production  requirements  in  this  context  is meant  to  include  all  costs  incurred  by  the mining  company  from construction  of  the  mine  to  closing  of  the  mine  (including  extraction  of  the  mineral,  processing  costs, transportation costs and mine operating and maintenance costs).   7  perspective,  the optimum mining methodology achieves  the highest  revenue  from mineral extraction and minimizes capital and operating expenditures required to extract the resource.  Since limestone and dolostone are low‐value commodities, they need to be mined in large volumes to be profitable.  Due to the high capital and operating costs of underground mining compared to surface mining, surface mining (quarrying) would be the preferred method of extraction for limestone and dolostone, assuming a near surface high purity carbonate formations existed near market. However, due  to a  variety of  restrictions and market  requirements, underground mining may be  the only economic option for some markets.  As mentioned, limestone and dolostone need to be mined in large  volumes  to  be  profitable.    Therefore, most  underground  limestone mines  are  room  and  pillar operations because  the extraction  ratio  is generally  the highest using  this method when mining  large volumes from laterally extensive, flat‐lying deposits; where deposits are not laterally extensive but thick, long hole stoping is a preferred option even though it’s at a lower extraction ratio (Brown et al., 2010).  Headings  and  cross‐drifts  in  room  and pillar mines  are often developed using  several benches  to  an overall excavation height of up to 30m  (Figure 1); the excavation height  is generally controlled by the thickness of the deposit and/or stability limitations.  Figure 1.  Schematic of a Typical Room and Pillar Mining Operation (Hamrin, 2001) 8   It is common for a mine to develop adjacent headings in the same orientation and regularly blast cross‐drifts between the development headings to maintain adequate ventilation of the workings.   With this mining method which  has multiple working  faces,  stress  redistribution  in  the  rock mass  around  the excavations  (pillars  and  roof)  is  continuous  and  if  excavations  are  sequenced  improperly  such  that unsupported large spans or stress concentrations are created, roof or pillar failure can occur risking the safety of workers and potentially comprising the feasibility of the mining operation. 2.3 Underground	Mining	as	a	Viable	Option	for	Limestone	Extraction		Within Canada, limestone is currently only mined from surface to meet the aggregate and construction materials  demand.    However,  many  resource  and  feasibility  assessments  for  subsurface  limestone extraction  have  been  conducted  for  the  region  of  Southern  Ontario  by  request  of  the  Ontario government  since  the  mid‐1970s  (Proctor  &  Redfern,  1974;  Acres  Consulting  Services  Ltd.,  1976; Thompson, 1982; Grass and Dupak, 1986; Planning Initiatives Ltd. and Associates, 1992; Shinobe, 1997).  These studies had  fluctuant conclusions between marginally economical and marginally uneconomical depending on market prices at the time of the study, resource location and the study assumptions; the primary resource target for these studies was under Lake Ontario and in the Niagara Escapement, both located at close proximity to the GTA.  Canada is not the only country monitoring the economic potential of limestone mining since the 1970s. Even  though  the United  Kingdom’s  (UK)  limestone mining  declined  throughout  the  20th  century  and quarrying operations dominated  as  the  iron  industry  subsided  and  the  aggregate demand  increased, over  the past 35 years,  there has been a growing  interest  in extraction of  limestone via underground mining methods once again  (Stocks, 1979).    In 2010, the UK government funded a study to assess the feasibility of the underground mining of aggregates within their country (Brown, et al., 2010).  This study was initiated because it was recognized that the surface hard rock resources suitable for aggregate are in the northern and western parts of England while the primary market demand is from the area around London and the south‐eastern and eastern regions of the country.   After considering  factors  including cost models  for aggregate production, stone processing, haulage of product to market, environmental impact mitigation, health and safety, decommissioning and restoration, it was concluded that although the capital costs were 1.33 – 1.65 times higher than surface operations of comparable size, the cost per tonne of aggregate production by underground mining methods (£13.03 – £13.93 per tonne)  is within the range of those costs reported by surface quarries (£10.95 – £16.48 per tonne) (Brown, et al., 2010).  9  However,  the  study  also  concluded  that  underground  mines  almost  always  have  a  greater  carbon footprint than surface quarries on the basis of CO2, eq/tonne produced, even though they can be much closer  to  the market, because  the underground mines  require  continuous  ventilation.   As developed countries  implement a  form of carbon  tax  to encourage greenhouse gas  reductions, having a greater carbon footprint could tip the scales back in the favour of surface quarrying. Feasibility studies  from many other countries worldwide arriving at similar conclusions as Canada and the UK provide a clear  indication that the global aggregate  industry  is preparing to tread underground once  again,  but  on  a  larger  scale.  Every  study  cited  the  driving  factors  as  stricter  environmental regulations  and  a deficit  in  available high quality  aggregate deposits  close  to market. Greece’s  study (Benardos, et al., 2001) looked beyond the immediate environmental and cost benefit to state that the underground space will be of critical importance in the future. Publications  released  for  the US aggregate  industry also  cite environmental  restrictions and  resource sterilization  due  to  population  growth  as  causing  a  shortage  in  limestone  based  products,  such  as Portland  cement  (US  Geological  Survey,  2012).    This  shortage  is  causing  construction  delays  and requiring an increase in limestone product imports from countries including Canada.  Given that the US remains the world’s top producer of aggregate in terms of number of operating mines and total tonnage of aggregate produced  (Benardos et al., 2001; US Geological Survey, 1999) and  they are struggling  to meet the aggregate industry demand because of mining restrictions, the concept of Canada tapping into underground limestone resources is truly imminent.  10  CHAPTER	3 LIMESTONE	MINING	POTENTIAL	IN	SOUTHERN	ONTARIO	Southern Ontario has been the focus of feasibility studies for  limestone mining within Canada because some of  the  thickest  limestone  formations  exist  there.   Also,  Southern Ontario  is densely populated containing  over  30%  of  the  country’s  population  (Government  of  Canada,  2015)  and  the GTA  is  the fastest  growing  region  of  Ontario  (Ministry  of  Finance,  2011).    Consequently  Southern  Ontario represents the largest domestic aggregate and construction materials market.  These two factors make Southern Ontario the most  likely region  in Canada  in which underground mining of  limestone could be economically feasible. This  chapter  discusses  the  geological  events  that  lead  to  the  sedimentary  sequences  in  Southern Ontario, the characteristics of the sedimentary sequences as a result of the depositional environment, the current state of limestone extraction in Southern Ontario and the specific limestone and dolostone formations that have potential to be mined using underground mining techniques. 3.1 Regional	Geology	Shallow epicontinental seas situated over Southern Ontario in the Paleozoic and Mesozoic eras resulted in  an  abundant  sedimentary  sequence  of  carbonates,  siliciclastics  and  evaporites  over  Precambrian rocks of  the Canadian  Shield  (Johnson  et  al., 1992; Derry  et  al., 1989).    Two major  regional  tectonic events during these eras have given the sedimentary  layers unique characteristics across the province.  Primarily two basins – the Appalachian Basin and Michigan Basin – as well as two arches – the Frontenac Arch  and Algonquin Arch  – were  formed  tilting  the  initially  flat‐lying  sedimentary  deposits  to  dip  at about 4 to 9m/km towards the basins (to the west north of the Algonquin arch and to the south of the Algonquin  arch)  (Johnson  et  al.,  1992;  Carmichael &  Smith,  1978).    As  shown  in  Figure  2,  Southern Ontario  is centered on  the Algonquin Axis and  therefore sedimentary deposits  tend  to be  thicker and better  preserved  towards  the  basins  where  deeper  waters  protected  the  deposits  from  the  many erosional events since the  late Cambrian as a result of eustatic sea  level change.   Most notably at the beginning of  the Middle Ordovician much of  the  sediments of Southern Ontario were eroded  leaving little evidence of  late Cambrian  and  early Ordovician deposits, however,  these  sequences have been preserved in the basins.  For perspective, the sedimentary deposits are preserved up to 1525m thick in some areas of Southern Ontario, and up to 13,000m thick in the Appalachian Basin.  Uniquely, although the sedimentary sequences were deposited around major tectonic events, they show  little evidence of structural disturbance;  this  is  important  for  geotechnical design purposes because  adverse  structural 11  conditions often are the deciding  factor when determining  if a deposit  is economical via underground mining for low‐yielding commodities.  Figure 2.  Principal Paleozoic and Mesozoic Tectonic Elements of Ontario (Johnson et al., 1992) 12  The tilted nature of the deposits and wide‐spread erosion has resulted in surficial geology mimicking the stratigraphic sequence of the Southern Ontario sedimentary system with successively younger rock units outcropping  from  the  northeast  to  the  southwest  of  the  region  (Birchard  et  al.  2004),  as  shown  on Figure  3.    The  surficial  geology  formations  across  Southern  Ontario  are  of  Ordovician,  Silurian  and Devonian age (359 – 488Ma) with the Ordovician exposed around the Kingston area (near the Frontenac Arch) and the Devonian exposed in the Detroit‐Windsor area (near the Michigan Basin).  13    Figure 3.  Paleozoic Geology of Southern Ontario (Birchard et al., 2004)  14  3.2 Current	State	of	Limestone	Extraction	in	Southern	Ontario	There are approximately 50 active quarrying operations  in Southern Ontario which supply much of the regions demand for construction materials and aggregates; no underground mines currently exist.  Table 1 summarizes the formations in Southern Ontario which have been or are currently being quarried, their thicknesses,  and  their  respective  products  (Derry,  1989).    Due  to  more  stringent  government regulations,  several  of  these  operations  have  been  denied  expansion  permits  and  few  prospective quarry sites are being granted surface mining rights.   As a result, many of the  limestone producers are seeking new opportunities to extract stone since the market demand continues to increase.   Because  of  the  tilt  and  lateral  persistence  of  the  sedimentary  formations  across  Southern  Ontario (particularly the Ordovician  limestone units which are present at surface  in Ottawa and exhibit similar characteristics when drilled at depth 500km away  in Niagara Falls),  limestone producers are  looking at the potential of extracting similar rock formations by underground mining methods at suitable locations to  the  southeast  of  their  current  operations.    Most  notably,  limestone  producers  with  quarrying operations  in  between  the  towns  of  Lindsay  and  Bowmanville  who  mine  the  upper  Ordovician formations are  studying  the  feasibility of using underground mining  to extract  the middle Ordovician formations from properties closer to the major construction market (Greater Toronto Area), since they are anticipating only greater difficulties with obtaining quarry expansion permits in the future.  Although this sort of study may seem a bit extreme given that underground mining is generally about three times the cost of surface mining, transportation costs are constantly increasing and therefore the reduction in distance  to  market  (significantly  reducing  transportation  costs)  is  making  the  underground  mining option cost comparable with that for production and transport from traditional surface quarries located further from market. In summary, quarries  in Southern Ontario are already actively searching alternative ore bodies  (either below  their  current  quarries  or  to  the  southwest where  the  same  ore  body  is  at  a  suitable  depth underground or closer to market) with the potential for extracting limestone via underground methods.  As the aggregate  industry  looks underground to meet the  limestone demand, the quarries of Southern Ontario provide key  insight  into  the desirable geological  formations which are  likely  to  represent  the target  formations  for  underground  extraction  in  the  future  since  most  formations  exist  at  depth somewhere else in the province.   15  Table 1.  Quarried Formations of Southern Ontario                                                             2 (Wolfe, 1993) 3 (Derry, 1989) 4 (Wolfe, 1993) 5 (Derry, 1989) 6 (Derry, 1989) 7 (Wolfe, 1993) 8 (Derry, 1989) 9 (Wolfe, 1993) 10 (Wolfe, 1993) 11 (Derry, 1989) 12 (Wolfe, 1993) Formation  Age  Geology  Products  Thickness (m) Dundee  Middle Devonian  Limestone  Crushed stone, aggregate2, cement, rock dust, armour stone3  35 – 45 Lucas  Middle Devonian  Limestone/Dolostone Aggregate, lime, cement, fillers, flux Anderdon Mb. – steel, cement and chemical industries 40 – 75 Amherstburg  Middle Devonian  Limestone/Dolostone Aggregate, lime, cement, fillers, flux Formosa Reef Limestone – lime, cement, flux, building stone4 40 – 60 (< 15) Onondaga  Middle Devonian  Limestone Crushed stone5 30 – 33.5Bois Blanc  Lower Devonian  Limestone/Dolostone Crushed stone6 3 – 50Bertie  Upper Silurian  Dolostone Crushed stone7 < 14Guelph  Middle Silurian  Dolostone Lime, crushed stone, aggregate, flux  4 – 100Lockport  Middle Silurian  Dolostone Crushed stone, lime, aggregate; Eramosa Mb. – building and landscape stone, flux Up to 46 (< 20) Amabel  Middle Silurian  Dolostone Crushed stone, agricultural lime, armour stone8 Wiarton/Colpoy Bay Mb – dimension stone Up to 38 (< 25) Decew  Middle Silurian  Dolostone  Cement9 ‐mined only in combination with the Lockport Fm.  Up to 4 St. Edmund  Middle Silurian  Dolostone Aggregate, crushed stone10  Up to 25 Manitoulin  Lower Silurian  Dolostone  Crushed stone, aggregate, building stone11  Up to 25 Lindsay  Middle Ordovician  Limestone  Crushed stone, aggregate and cement  Up to 67 Verulam  Middle Ordovician  Limestone  Crushed stone, aggregate and cement  32 – 65 Bobcaygeon  Middle Ordovician  Limestone Crushed stone, lime, building stone12, aggregate  7 – 87 Gull River  Middle Ordovician  Limestone Crushed stone, building stone, lime, aggregate  7.5 – 136 16  3.3 Potential	Target	Formations	for	Limestone	and	Dolostone	Mining	In Southern Ontario, the succession of the Gull River, Bobcaygeon, Verulam and Lindsay Formations of the Middle Ordovician represent  the most desirable targets  for underground  limestone mining due to their thickness and chemical composition.  The deposits range from coarse‐grained bioclastic carbonates to  carbonate mudstones with  interbedded  calcareous  and  non‐calcareous  shales  (Ontario Geological Survey, 1992).   These units formed through continuous deposition as water  levels generally decreased over the area (which was a marine reef or shelf at the time) resulting  in gradational contacts between the formations.  Consequently, field identification of individual formations is difficult which could prove a  challenge  for mine design.    The  interbedded  shales not only  impact  the quality of  the deposit but based on industry experience have major implications on the structural stability of the rock mass around underground excavations. The  upper  Lindsay  Formation  gradates  into  highly  bituminous  shales  in  central  Southern  Ontario (graphical region extending from Toronto‐Oshawa, east to Lake Huron, and north to the Bruce Peninsula and Manitoulin  Island)  (Ontario Geological  Survey,  1992).    In  the  literature,  these  shales  have  been specifically  classified  as  the  Collingwood Member  of  the  Lindsay  Formation  or  the  lower  gradational contact of the Blue Mountain Formation with the Lindsay Formation.  These shales represent the upper limit of potentially feasible limestone units of the Ordovician age.  From a geotechnical perspective, the shales of the upper Ordovician, although a poor quality rock mass to design within and may have some environmental  concerns  due  to  their  bitumen  content,  are  an  effective water  barrier  and minimize vertical  transmissivity  of  fluid.    Therefore,  situating  a  mine  at  a  sufficient  distance  below  these formations could reduce the  influence of groundwater  in the mine; water management  is a significant cost in underground mining. The  Silurian  age  represents  a  gradation  from  interbedded  sandstones  and  dolostones with  shales  to interbedded  dolomitic  limestones  with  shales  to  evaporites  and  dolostones  in  succession.    Many Formations  (Manitoulin,  St.  Edmund,  Fossil Hill,  Reynales,  and  Bertie) within  the  Silurian  have  been quarried in the past and/or continue to produce today, primarily for aggregate or crushed stone, but the sequence of  the Amabel,  Lockport and Guelph Formations of  the Middle Silurian exhibit  the greatest potential for further exploitation through underground mining.   The Amabel, Lockport and Guelph Formations were deposited as water levels decreased over the region creating a shallow shelf, and reef‐interreef and barrier depositional environment.  The resulting deposits 17  are  generally  thinly bedded  to massive,  coarse‐grained,  fossiliferous dolostones,  containing  a  varying degree of bitumen.  However, it should be noted that the deposits of the Lockport Formation are more complex than the Amabel Formation, although  laterally equivalent13, due to the deeper water  levels  in the Appalachian Basin.  The more complex deposition conditions of the Lockport Formation resulted in higher content of chert and shale (Goat  Island member) making the Formation a more difficult mining horizon than the Amabel Formation.  Overall, the thicknesses and purity of these dolostone Formations (specifically  the Wiarton/Colpoy Bay Member of  the Amabel Formation and Eramosa Member of  the Lockport  Formation)  make  them  desirable  for  a  variety  of  applications  including  building  stone, dimension stone, aggregate and crushed stone, and consequently potential underground mining targets. The  onset  of  the  Acadian Orogeny  in  the  latter  part  of  the Middle  Silurian  altered  the  depositional environment  over  Southern Ontario  by  restricting  circulation  in  the  inland  seas.    Consequently,  the Guelph Formation was overlain by a sequence of evaporites and dolostones (Salina Formation).  Due to the low vertical permeability of the Salina Formation, the Guelph and Amabel/Lockport Formations act as  an  aquifer.    Currently,  these  formations  are  the  main  water  supply  for  the  cities  of  Kitchener, Waterloo,  Cambridge  and  Guelph.    As  such,  regulations  have  been  implemented  restricting  the exploitation of these Formations for any activity which has potential to disrupt or contaminate the water supply to these cities. In 2011, The Highland Companies tested these regulations by submitting an application to the Ontario government  to  develop  a  limestone  quarry  in  the  Township  of  Melancthon  (approximately  75km northwest of the GTA) targeting the Amabel Formation.  Even though the proposed quarry location is an equivalent  75km  northeast  of  Guelph  (the  closest  city  being  supplied  with  water  by  the  Amabel Formation),  cooperation  with  the  province  and  the  local  community  could  never  be  gained  citing concerns  about  the  potential  for  water  contamination  (VanDyken,  2011).  The  company  ultimately withdrew  its application  in  late 2012.   Both supporters and protestors of  the project have highlighted the key fact that the application was not officially denied by the government, so potential to target this formation in the future may exist (National Farmers Union, 2015; D'aliesio & Howlett, 2012); in the near future though, the Amabel Formation appears off‐limits for mineral extraction within a feasible distance from the GTA market.                                                              13 The Lockport Formation  is primarily  located within the Appalachian Basin.   The Amabel Formation  is primarily located within the Michigan Basin.  The two deposits are geographically separated by the topographical high of the Algonquin Arch.   However,  it  is  unclear  as  to whether  these  units were  laterally  gradational  at  one  time  and greater erosion over the Algonquin Arch eventually separated them. 18  Although  the  limestone  and  dolostone  formations  of  the  Devonian  age  (Bois  Blanc,  Onondaga, Amherstburg, Lucas, Dundee) are currently quarried extensively for aggregate, building stone, lime and cement in the southwest of Ontario, where the deposits have not been eroded, they are shallow (<200 m)  and  therefore  it  would  be  more  economical  to  mine  these  units  via  surface  extraction  where permitting allows.  Additional detailed investigation into the feasibility of subsurface extraction of these units  might  find  it marginally  economical  near  Sarnia  where  the  formation  depths  and  population density are  slightly greater, however,  for  the purpose of  the numerical analysis  in  this  thesis,  specific cases related to the underground extraction of Devonian age formations will not be modelled.   Depending on the location of market demand and existing infrastructure (such as an existing port on the Great Lakes),  the discussed  formations  represent potentially viable  limestone or dolostone  resources.  However,  as  summarized  by  Lee  and  White  (1993),  the  Gull  River  of  the  Ordovician  Formation  is considered  the most  viable mining  horizon  for  underground  extraction  because  of  the  purity  of  the limestone and the structural competency of the rock mass.  The numerical modelling for this thesis will be conceptually based on  the extraction of  the Gull River Formation below  the GTA  (an approximate depth of 150m below ground surface). 19  CHAPTER	4 MINING	CONDITIONS	IN	SOUTHERN	ONTARIO	Currently,  underground mining  in  Southern  Ontario  is  for  the  extraction  of  salt  and  gypsum.  Both minerals  are  relatively  soft  in  comparison  to  limestone  and  therefore  present  a  unique  set  of geotechnical  design  challenges  than would be  anticipated  for  limestone.  To predict  the  geotechnical challenges that may face an underground limestone mine, the mining conditions need to be understood.  From  a  geotechnical  perspective,  the  in‐situ  stress  regime,  structural  nature  of  the  rock mass  and strength of  the  rock mass are  the  key  aspects  for  input  into a mine design  and  for geotechnical  risk assessment; hydrogeology  is another major aspect but  this  field of  study has been omitted  from  this research because it does not contribute to the numerical analysis.  In knowing the range in these mining conditions,  an  engineer  can  be  proactive  in  assessing  probably  of  failure  via  statistical  or  empirical analysis for the many potential failure modes, and consequently develop an appropriate ground control program (including ground support) to mitigate geotechnical risk. This chapter discusses the in‐situ stress regime in Southern Ontario, focusing on the GTA, as well as the regional  structural  characteristics  of  the  sedimentary  sequence  in  the  same  area.    The  rock  mass strength characteristics and specific structural characteristics are detailed  in Chapters 7 and 8.   Finally, published  experience  from  the  US  aggregate  mining  industry  within  similar  mining  conditions  as anticipated for Southern Ontario is reviewed.  Typical room and pillar mining geometries and a synopsis of the geotechnical stability issues plaguing the US limestone mines are presented.  4.1 Southern	Ontario	In‐Situ	Stress	Regime	The  existence  of  high  horizontal  stress  in  the  sedimentary  rocks  of  Southern Ontario  has  been well documented (Lee, 1981; Lo, 1978).  Squeezing observations made during the construction of the wheel pits of the Canadian Niagara and Toronto Power Plants in 1903 first triggered the notion of high in‐situ stresses.  In the years following, heave of quarry floors, identification of “pop up” structures, buckling of canal floors, movement of bridge foundations and cracking of tunnel walls during or after construction provided  further  evidence.    But  it wasn’t  until  the  1970’s  that  the  uniquely  high  horizontal  stresses relative  to  the  vertical  stresses  at  shallow  depths  were  characterized  for  the  southern  part  of  the province through in situ stress measurements performed in various limestone and shale formations (Lo, 1978). 20  The  high  horizontal  in‐situ  stresses  are  a  function  of  the  current  tectonic  plate motions  and  glacial unloading.    Ontario  is  situated  on  the  North  American  plate  which  is  currently moving  in  a WSW direction as  it  is pushed away  from  the Mid‐Atlantic  ridge creating a  regional horizontal  in‐situ  stress regime along an ENE trend  (Nuclear Waste Management Organization and AECOM Canada Ltd., 2011; Esterhuizen et al., 2007).    In  the past,  continental glaciers over  the area  increased  the vertical  stress resulting  in a more  isostatic  stress  field.   As glacial unloading occurred,  the  relative difference  in  the horizontal and vertical stresses increased leaving the uniquely high horizontal stress (relative to vertical stress) conditions measured today. Locally, the magnitude and orientation of the in‐situ stress regime can deviate from the regional average due to rock mass quality and stiffness; stresses accumulate in massive, stiff rigid rock masses (granite or limestone) but not in weak, softer ones (soft shale or salt).  Another possibility resulting in a unique local stress regime relative to the regional one  is that the horizontal stresses could have been relieved over geological  time  by  local  events  such  as  outcropping  or  folding  (Iannacchione  et  al.,  2001);  pop‐up structures in Southern Ontario limestones are a common sign of local stress relief.  As part of the studies conducted to assess the risk of developing a deep geologic repository for low and intermediate  nuclear waste  storage  at  the  Bruce Nuclear  Power  Plant  in  Southern Ontario  (Nuclear Waste Management Organization and AECOM Canada Ltd., 2011), published stress testing data from 26 sites within the Appalachian and Michigan basins which were plotted to empirically determine the stress ratios between the maximum and minimum horizontal stresses and the vertical stress  (Figure 4).   Key relationships from the in‐situ stress data are as follows: 1) At depths of approximately 200m and 800m,  there are distinct decreases  in  the horizontal  to vertical stress ratios. 2) The maximum horizontal to minimum horizontal stress ratio remains nearly the same regardless of depth. 3) Near surface (<50m deep) stress ratios can be very high (k > 6) because of local concentrations of high horizontal stress are being compared to the low vertical stress induced by the weight of the bedrock.  These uniquely high stress ratios were disregarded from this research because the minimum depth for mining was considered to be 100m.  The stress ratios as a function of depth are summarized in Table 2. 21  Table 2.  Summary of Stress Ratios by Depth below Ground Surface Depth  H/V  h/V  H/h 0m – 200m  3.8  2.0  1.9 200m – 800m  2.0  1.1  1.8 >800m    1.4  0.8  1.8  Stress measurements taken at the Darlington Nuclear Power Plant, which is located about 30km east of the Greater Toronto Area and adjacent to one of the larger surface quarries in Ontario, concluded that the principal  in‐situ stress  in the Ordovician  limestone formations  is horizontal and oriented at ENE as shown in Figure 5  (Haimson & Lee, 1980).  Additionally, the horizontal stresses are relatively constant at approximately 13MPa  for  the maximum horizontal  stress within  the Ordovician  limestone  formations and not dependent on depth. 22     Figure 4.  Regional In‐Situ Stress Ratios as a Function of Depth: (a) H/V, (b) h/V, and (c) H/h (Nuclear Waste Management Organization and AECOM Canada Ltd., 2011) 23    Figure 5.  In‐Situ Stress Measurements at Darlington Nuclear Power Plant (Haimson & Lee, 1980) 4.2 Structural	Characteristics	of	Southern	Ontario	Sedimentary	Sequences	Sedimentary sequences  in Ontario often are characterized by a minimum of  three major discontinuity sets: one  flat‐lying set representing  the bedding plane, and  two  relatively orthogonal sub‐vertical sets postulated to be the result of tensile fracturing as  isostatic rebound occurred (Goudie, 2006).   Bedding across Southern Ontario dips at approximately 0.5 degree towards the southwest.  The strike of the sub‐vertical  sets  can  vary  significantly on  a  local  scale but  regionally one  set  aligns with  the direction of maximum horizontal stress and the other aligns with the direction of the minimum horizontal stress.  A unanimous  conclusion  explaining  the  “coincidental”  alignment  of  the  fracture  network  with  the maximum and minimum horizontal stress directions has yet to be reached, even though many scientists believe  the  current  in‐situ  stress  direction  has  regionally  been  constant  since  the  Paleozoic  and therefore  could  have  resulted  in  the  orthogonal  fracture  network mapped  across  Southern  Ontario (Gross  &  Engelder,  1991;  Holst,  1982).    Geotechnical  data  collected  from  limestone  formations  at 24  various  locations  in  Southern  Ontario  is  presented  in  Figure  6  and  shows  this  common  structural relationship between bedding and two orthogonal joint sets in sedimentary rock formations.    Figure  6.    Structural Data  from  across  Southern Ontario:  (a)  Kincardine, ON  (Cruden,  2011),  (b)  Lincoln, ON (Gartner Lee Limited, 1996) 4.3 Geotechnical	Experience	from	the	US	Aggregate	Mining	Industry	Room  and  pillar  aggregate mining  has  flourished  in  the US  over  the  past  century.   Aggregate mines located in the northeastern US exhibit similar stress conditions to those in Southern Ontario – near the Greater  Toronto  Area  (as  shown  in  Figure  7)  –  and  therefore  these mines  provide  a  good  basis  for predicting  the mining and geotechnical  challenges  that are  likely  to be encountered  in  the  inevitable Southern Ontario aggregate mining operations.  25   Figure 7.  North American In‐Situ Stress Map (adapted from Heidbach et al., 2010) After  recognizing  similar  roof  stability  issues  in underground  limestone mines across  the Eastern and Midwestern USA, Esterhuizen et al. (2007) conducted a survey of 34 mines collecting data pertaining to roof spans, rock mass properties, support practices and roof  instabilities as part of a National  Institute for Occupational Safety and Health research project.  From the survey, there are three conclusions that provide  useful  insight  into  limestone  mining  under  high  horizontal  stress  conditions  that  can  be translated back to the predictive research being conducted as part of this thesis for mining in Southern Ontario. Firstly, Table 3 summarizes the mining dimensions measured during the survey which provide a basis for the excavation sizes likely to be employed in Southern Ontario.    26  Table 3.  Summary of Underground Mining Dimensions from US Limestone Mines Parameter  Literature Average  Minimum  Maximum Mining Height (m)  11.6  4.8  38.0 Pillar Width (m)  13.8  4.6  28.6 Room Width (m)  13.5  9.1  16.8 Intersection Diagonal (m)  21.7  16.1  29.6 Secondly, small failures were observed in 30 of the 34 mines surveyed.  Large failures were observed in 19 of  the 30 mines  that had small  failures.   These statistics  indicate  that  the majority of mines under high horizontal stress  in conditions similar to those expected  in Southern Ontario are still not properly designed  to  prevent  large  failures  and  consequently  the  associated  risk  to  human  life  would  be considered unacceptable in today’s mining industry.   Small and large failures were defined as rock falls and roof falls respectively based on their failure characteristics, as shown in Table 4. Roof  falls  are  the  target  of  this  thesis’  research.    Esterhuizen  et  al.  (2007)  concluded  that  the main factors contributing to large roof falls are horizontal stress, large joints and insufficient thickness of the immediate  roof beam.   Of  the  roof  falls,  stress was  the main contributing  factor  in 36% of  them and these  failures were observed at a variety of depths.   Beam  failures  constituted 28% of  roof  falls and always truncated at the boundary between  limestone and some overlying weak band or parting plane.  Block failures governed by large discontinuities extending across entire rooms marked 21% of roof falls.  Caving failures which represented only 15% of roof falls were all related to the collapse of weak shale exposed  in  erosion  channels  or  progressive  failure  of weak  roof  rocks.   According  to  these  roof  fall observations, weak  layers or shales contributed  to a minimum of 43% of  failures  (not  including  those classified  as  stress  failures which may  have  had  a weak  layer  or  shale  influence  that  could  not  be observed)  in  limestone mining environments and  therefore  these weak  layer  influences are critical  to understand  for  safe  mine  design  in  the  future.    In  the  author’s  opinion,  Esterhuizen  et  al.  (2007) adequately portray  the detrimental  risks  to a mining  company  associated with  roof  falls  through  the following statement,  “although large roof falls only make up a small percentage of total roof exposure, their potential impact on safety and mine operations can be very significant. Most cases of large roof falls require barricading‐off or abandonment of the affected entry.”  27  Based on  the observation noted  in  this statement  that  roof  falls  typically  result  in a  truncated mining header, one can easily infer that loss of equipment or life would occur if these entities were involved in the failure.  Table 4.  Classifications of Roof Failure during Survey (Esterhuizen et al., 2007) Size of Failure  Type of Failure  Description Small Rock Falls  Isolated rock fragments less than about 1m in length Slabs  Thin slabs caused by weathering or stress spalling,  less than about 30cm in length and about 25mm thick Blocks  Blocky  rock  fragments  caused  by  the  intersection  of  joint  planes, blasting fractures, bedding or stress fractures Beams  Stepped roof or brow formed by fall propagation to bedding plane Large Roof Falls  Falls  larger  than  about  1m  in  length,  typically  consisting  of multiple rock fragments Blocks  Large discontinuities and joints associated with fall Beams  Bedded layers in the roof fail under gravity loading Stress  Horizontal stress‐related shearing and buckling of roof beds Caving  Fall  caused  by  progressive  spalling,  blocky  roof  or weathering  of weak strata The failure mechanisms leading to the rock and roof fall categories listed by Esterhuizen et al. (2007) are thoroughly discussed in Chapter 5. Lastly, 25 of the 34 mines employed a mine design that accounted for a minimum thickness of limestone beam  in  the  immediate  roof.   Mines  that  preserved  a  competent  roof  beam  thickness  of  2.25m  or greater were able  to mine without  regular  support, compared  to  those with average competent  roof beam thickness of 1.3m required regular ground support in the roof. This experience from aggregate mining in the Eastern and Midwestern US provides valuable insight into the potential geotechnical issues that engineers will need to account for when completing mine designs for similar mining conditions  in Southern Ontario.   However,  the maximum  in‐situ horizontal stress  in the US  aggregate mines  typically  ranges  between  2  to  3  times  the  in‐situ  vertical  stress  (Agapito & Gilbride, 2002).   Therefore, the geotechnical  issues observed may be exacerbated under the Southern Ontario  stress  regime where  the maximum  in‐situ  horizontal  stress  can  be  3  to  4  times  the  in‐situ vertical stress.  28  CHAPTER	5 FAILURE	MECHANISMS	Rocks, and  the  conditions under which  they  fail, have been  the  subject of  study  since as  far back as Ancient  Egyptian  times.   More  coordinated  efforts  targeted  at  understanding  the  driving  forces  and complex mechanisms  leading to rock failure have dominated engineering publications since the 1960s, and reiterative analyses efforts continue as new technologies allow for more realistic representations of failure environments and  loading conditions  in experimental studies.   Understanding of these complex relationships  remains  very  difficult  because many  of  the  rock mass  properties  or  loading  conditions which appear to drive failure tend to have a coupled functionality, and subsequently are dependent on at  least  one  other  variable  in  the  system  when  governing  failure.    Consequently,  outlining  of  the boundaries that define the ranges in which each failure mechanism initiates is environment specific, and many failure mechanisms can initiate simultaneously.  Over the years of experimentation, it has become clear  that  stress  acting  on  a  rock mass  and  structure within  the  rock mass  in  their  varying  degrees govern  the  mechanisms  defining  the  different  modes  of  failure  possible  around  an  underground excavation  (Figure  8).    Numerous modes  of  failure  exist  because many  factors  (physical  properties, homogeneity  and  anisotropy  of  the  rock mass,  discontinuity  properties,  free  face  geometry,  in‐situ stress, etc.) affect the degree of  impact  in which the  in‐situ stress has on the rock mass.   Therefore,  in order to avoid the many forms of rock failure when performing engineering design, it is critical to have a thorough  understanding  of  the  rock  failure mechanisms  pertaining  to  the  environment  in which  the design is for. 29   Figure 8.  General Rock Failure Mechanisms for Underground Excavation (adapted from Martin et. al, 1999) This chapter begins by outlining the general modes of stress and structurally controlled failure as well as failure  modes  governed  by  both  stress  and  structural  mechanisms  interdependently.    Following,  a comprehensive review of the common mechanisms and modes of failure observed  in stratified rock  is given, highlighting the factors that have the greatest impact on the mechanisms of failure observed. 5.1 Types	of	Rock	Failure	Rock failures are often described as either being stress controlled or structurally controlled based on the mechanism which appeared to drive failure.  However, rock failure in reality is rarely that simple and a unique  combination  of  applied  stress  and  rock  mass  structure  must  exist  to  result  in  the  failure 30  mechanisms  observed.    This  section  discusses  the  spectrum  of  failure mechanisms  from  dominantly stress  controlled  to  dominantly  structurally  controlled,  as well  as  those with  interdependent  failure relationship. 5.1.1 Stress	Controlled	Failure	In  homogeneous  rocks,  failure  initiates  when  stress  being  applied  to  them  exceeds  either  their compressive, tensile or shear strength.  Tang and Hudson (2010) even suggest that at a micro‐scale level, intact  rock  failure  initiation  is only a  function of  the  tensile and  shear  components;  since  the  tensile strength  of  rocks  is  about  10%  of  its  compressive  strength,  C,  (Jumikis,  1979)  the  indirect  tensile stresses  induced  in  a  compressive  stress  environment  cause  the  rock  to  fail  in  tension  before  the compressive strength  is exceeded.    In‐situ rock failure observations  indicate that a combination of the three  forms  of  stress  application  characterize  the mode  of  failure  observed.    For  example,  during  a uniaxial  compression  test  on  a  homogeneous  rock  sample,  failure  usually  initiates  in  tension  in  the middle of  the  sample.   Once  the  tensile  fractures  coalesce under  increased  loading,  failure becomes shear dominated until the sample yields at its maximum unconfined compressive strength.   Rock masses  typically  fail  in  two distinct manners, brittle or non‐brittle, depending on  the  intact  rock strength characteristics, and  the stress conditions acting on  them.   Brittle  failure  through  intact  rocks commonly occurs around excavations  in strong, massive rock masses (RMR14 > 75) under moderate to high in‐situ stress conditions (σ115/σC > 0.15); i.e. Modes 4 and 7 shown in Figure 8.  During brittle failure, instability often occurs  in  the  form of axial  splitting, and  spalling or  “bursting”  (Li & Nordlund, 1993; Stacey & Page, 1986).   As Martin (1997) has characterized, these failures are typically governed by the tensile strength of the intact rock and intact rock fracturing (i.e. spalling) initiates as a result of induced tensile stresses at grain‐scale when  the major principal stress  (1) equates to about 40% of the rock’s compressive strength.   Non‐brittle failure around excavations  in massive rock masses under  low stress conditions  (σ1/σC  <  0.15);  i.e. Mode  1  in  Figure  8,  are  not  common  as  the  rock  typically  responds elastically  to  the stress and  remains stable.    In certain cases, where  the  rock mass  is weak and under                                                             14  RMR  (Rock Mass  Rating)  is  a  rock mass  classification  system  develop  by  Bieniawski  (1976) which  has  been refined and modified over the year; notable amendments are by Bieniawski (1989) and Laubscher (1990).  RMR is a summation  process  that  assigns  a  quantitative  value  to  represent  the  quality  of  a  rock  mass.    The  RMR classification  takes  into  account  intact  rock  strength,  rock  quality  designation  (RQD),  joint  spacing,  joint characteristics (roughness, alteration and aperture) and in‐situ water conditions. 15  In‐situ  stresses  acting on  a  rock mass  are  analyzed using  the magnitudes  and orientations of  their principal stress: 1 – major principal  stress, 2 –  intermediate principal  stress, 3 – minor principal  stress.   The principal stresses are mutually orthogonal and characterize the three‐dimensional stress field acting on a rock mass.  31  high stress conditions (σ1/σC < 0.4), the rock mass will continuously deform over time in the excavation in  the  form  of  squeezing  or  swelling,  similar  to  Mode  9  in  Figure  8.  Therefore,  a  comprehensive understanding of the intact rock strength properties (, E, c) and in‐situ stress regime (1, 2, and 3) is essential for assessing potential for instability in stress controlled environments. 5.1.2 Structurally	Controlled	Failure	Discontinuities from joints to large scale faults separate the rock mass into wedges or blocks of varying sizes, and interrupt the continuous nature of the mass by introducing boundaries with unique properties between the blocks changing the kinematics of failure.  Structurally controlled failures are governed by the  strength  properties  of  the  discontinuities.    These  failures  occur  in  the  form  of  gravity  or  sliding failures;  i.e. Modes 2 and 3  in Figure 8. The orientations of the discontinuities will determine whether gravity  or  sliding  failure  is  possible.    The  joint  tensile  strength  governs  failure  initiation  for  gravity failures whereas  the  shear  strength properties  (effective  angle of  friction  and  joint  cohesion)  govern stability for sliding failures (Stacey & Page, 1986).  When  the nature of  rock mass  is heavily  fractured, disturbed or non‐cohesive,  failure  typically occurs along particle boundaries when the cohesive and frictional properties between the particles of the rock mass  are  exceeded.    Accordingly,  failure  is  governed  by  the  shear  strength  properties  of  the discontinuities.  However, contrary to moderately fractured rock masses, heavily fractured rock masses often respond  to stresses similar  to weak, homogenous rock masses.   For example, under high  in‐situ stress conditions (σ1/σC > 0.40), the small blocks composing the rock mass will continuously move and bridge over time in the excavation in the form of what is observed as squeezing; i.e. Mode 9 in Figure 8.  For  low and moderate  in‐situ stress conditions (σ1/σC < 0.40), failure occurs  in the form of unravelling; i.e. Modes 3 and 6 in Figure 8.   5.1.3 Compound	Failure	Mechanisms	Much of the complexity associated with understanding the failure mechanisms acting in a rock mass has to do with  the nature of  the  rock.   A  rock mass  is  rarely an  ideal homogeneous,  isotropic and elastic solid.    The  degree  of  heterogeneity,  anisotropy,  confinement,  block  geometry,  boundaries  and  pre‐existing  defects  affects  the  mechanisms  and  modes  of  failure  which  are  able  to  develop.    The heterogeneities  in  a  rock mass  can be  a  function of differential  solidification processes,  stratification differences  or  differential metamorphic  processes  as  the  rock  formed.    Because  of  the  randomness usually associated with heterogeneities in a rock mass and spacing of discontinuity sets, it is common for 32  both structurally controlled and stress controlled modes of  failure  to be operating simultaneously.    In these cases, the mechanisms  leading to  failure are a combination of  localized  fracturing of  intact rock and an exceedance of the strength properties of a discontinuity.  This is illustrated in Modes 5 and 8 in Figure 8 where in a moderately fractured environment (50 < RMR < 75) under moderate to high in‐situ stress conditions (σ1/σC > 0.15), failure will be a combination of  localized brittle  intact rock failure and sliding or tensile failure along the discontinuity.   Due  to  the  nature  of  deposition  of  interbedded  limestone  and  shale  rock  masses,  semi‐persistent occurrence of  sub‐vertical  jointing and  relatively weak discontinuity  strength properties  compared  to the  intact  rock  strength properties of  the  limestone,  the modes of  failure documented  for  room and pillar limestone mines are the result of compound failure mechanisms. 5.2 Modes	of	Failure	and	Failure	Mechanisms	in	Stratified	Rocks	Stratified  rock masses,  such  as  interbedded  limestone  and  shale,  are  geometrically  characterized  by planar and persistent bedding which has low or zero tensile strength in the direction orthogonal to the bedding plane and low shear strength on the surfaces compared to the intact rock strength properties of the  limestone.   Consequently, upon excavating a room for the extraction of  limestone, the  intact rock between  the  bedding  planes  overlying  the  excavation  becomes  a  self‐supporting  “beam”,  and  two different mechanisms of failure may initiate (Brady & Brown, 1985): 1) Surface  spalling  and  internal  fracturing  of  the  intact  rock  comprising  the  beam  or  beams overlying the excavation (particularly at the ends and centre of the beam), and 2) Separation of  the bedding planes  due  to  the weight of  the  beam or  indirect  tensile  stresses induced by compressive  loads  (such as  in a horizontally dominated  in‐situ  stress environment (Jumikis, 1979)) and shearing along the bedding planes as the beam or beams deflect  into the excavation. Depending on the strength of the intact rock relative to the in‐situ stress conditions and the geometry of the  beam  (thickness  to  span  or  slenderness  ratio),  both mechanisms  of  failure  are  likely  to  initiate resulting in a compound rock failure.  In cases such as Southern Ontario where the cross‐cutting  joints are orthogonal to the bedding planes within a tolerance of one third to one half of the effective friction angle of the joints (Figure 9), then the 33  rock mass  is considered  to be a composition of a  series of Voussoir beams16 of which Voussoir beam theory can be applied to predict failure in homogenous rock masses (Ran et al., 1994).    Figure 9.  Voussoir Beam Criteria The  modes  of  failure  associated  with  Voussoir  beams  are  slabbing,  buckling,  crushing  or  diagonal fracturing (Figure 10).  These failure modes fall under Modes 2, 5 and 8 in Figure 8 depending on the in‐situ stress conditions.    Figure 10.  Specific Failure Modes in Stratified Rock Masses                                                              16  The  Voussoir  beam  label  for  a  beam  of  rock  spanning  an  excavation was  adopted  from  the  Voussoir  arch considered  in masonry structures because of similar relationships between vertical deflection,  lateral  thrust and beam stability observed during Evans (1941) investigations of roof deformation mechanics.  34  Voussoir beam theory indicates that the maximum deflection that a beam can withstand prior to failure is a function of the specific weight (), Young’s modulus (E), uniaxial compressive strength (UCS), span and  thickness  of  the  beam  (Diederichs &  Kaiser,  1999;  Stacey &  Page,  1986;  Brady &  Brown,  1985; Sterling,  1980).  Following  a  series  of  experiments  using  a  base  friction  machine,  Goodman  (1989) extrapolated Voussoir beam theory for a single beam to calculate the stresses and deflections acting on the lower beam in a multiple beam failing system.  Goodman’s equations were primarily functions of the beams’ geometries and their individual  and E.  Figure 11 highlights the relative conditions under which each mode of  failure  is dominant.   Division  lines between  failure modes  in Figure 11 are hypothetical and based on published qualitative observations; they have not been validated with quantitative data.  In general  though, where  the slenderness ratio  is high  (beam  thickness  is  less  than 1/10 of  the span), shearing governs failure.  In contrast, where the slenderness ratio is low (beam thickness is greater than 1/2 of the span), tensile fracturing – direct from beam deflection and indirect from compressive forces – governs failure.    Figure 11.  Conditions Leading to Different Failure Modes in Stratified Rock Masses  35  It should be noted  that  in cases where additional moderately dipping  (30 – 45 degree dip)  rock mass structure exists, Voussoir beam theory does not apply because the stability of the beam is compromised by the moderately dipping structure; slip along these structures is likely to occur resulting in premature shear failure (Ran et al., 1994).  Additionally, the Voussoir analogue is not applicable for poor rock mass conditions where there are more than 3 joint sets or the RQD is less than 50  (Diederichs & Kaiser, 1999).  For the purpose of this research, only a rock mass composed of classic Voussoir beams was considered. The following sections discuss the conditions which result  in each of the Voussoir beam failure modes observed  in  stratified  rock masses  and  how  those  conditions  contribute  to  the  development  of  the failure mechanisms that define each mode of failure. 5.2.1 The	Arching	Phenomenon	At a certain beam thickness to span ratio of a Voussoir beam, a phenomenon called “arching” occurs.  Arching is the process of stress redistribution in the beam to its abutments when the beam is suspended under  its own weight.    If  the beam’s  thickness  is enough  to allow  for  complete  redistribution of  the stresses across  the beam’s span  to  its abutments  then stability can be achieved naturally.   Brady and Brown  (1985)  mathematically  characterized  this  phenomenon  by  calculating  the  stress  distribution curve within a single deflecting beam. Their work has since been further developed by Diederichs and Kaiser (1999) so that an assessment of the stability of a Voussoir beam can be made quantitatively.  The stress distribution curve in a deflecting beam is a function of the beam’s deflection as well as the lateral thrust generated by the force of the horizontal stress. Although arching within a single beam has been mathematically explained  in recent years, the arching phenomenon was discovered back in 1885 by Henri Fayol, a French mining engineer, while investigating the behavior of stacks of beams spanning a simple support system by observing  the deflection of  the lowest beam as successive beams were added on top (Brady & Brown, 1985).  Fayol concluded that at a certain thickness of stack, the load applied by adding successive beams loaded onto the stack does not influence the deflection of the  lowest beam.   Having no mathematical proof  for this phenomenon, he could only postulate that the additional  load was transferred  laterally to the beam supports  instead of onto the lowest beam.  The mathematics explaining arching in a stack of Voussoir beams – an accurate representation of  a  rock  arch  formed  above  a mining  excavation  in  a  stratified  rock mass  –  remains estimated at best due to the number of complex  interactions between the various beams forming the 36  arch.   Figure 12 present a visual  interpretation of how arching occurs within a single Voussoir beam as well as a stack of Voussoir beams.  Figure 12.  Compression Arches Composed of (a) Single Voussoir Beam, (b) Stack of Voussoir Beams  The  research  being  conducted  as  part  of  this  thesis  does  not  require  the  understanding  of  the mathematical  relationships pertaining  to Voussoir beams. However,  an  appreciation  for  the Voussoir beam  properties  directly  contributing  to  the  development  of  a  multiple  beam  rock  arch  –  span, thickness, Young’s modulus, specific weight and joint friction, cohesion and tensile strength – should be had  since  the arch controls  the depth of  failure  in  the  stratified  roof  rock mass over an underground excavation and  the purpose of  this research  is  to assess how a shale  layer  influences  the ability of an interbedded limestone‐shale roof to create an arch and become self‐supporting. 5.2.2 Slabbing	Failure	Slabbing  failure  occurs  when  the  beam  thickness  to  span  ratio  is  low  to moderate  and  the  in‐situ horizontal stress is low.  Under these conditions, the weight of the beam exceeds the tensile strength of the bedding plane with  the overlying  strata allowing  separation  from  it and causing  shear  stresses  to develop  vertically  at  the  abutments of  the  self‐supporting beam.    If  the  load of  the beam or beams exceeds the rock mass shear strength, the beam will fail at the abutments, as shown in Figure 13.  With very  thin  beams  (thickness  <  1/100  of  the  span),  slabbing  failure  still  occurs  but  shear  stress 37  concentrations  tend  to  form  at  local  heterogeneities  within  the  slab  instead  of  at  the  abutments.  Accordingly, very thin beams often fail in smaller slabs but have more frequent occurrences of failures.  If  the rock mass  is much weaker  than  the weight of  the slabs, slabbing  failure can progress over  time with successively higher beams failing resulting in a “chimney”‐style failure.  In cases where the mining excavations are near surface, the chimney failure can create surface subsidence or a sinkhole.  This style of failure and corresponding sinkholes have plagued the UK  in the recent past due to standard mining practices from the late 1800’s and early 1900’s for extracting high quality limestone beneath weak coal seams which left near surface mine workings unsupported and open after mine closure  (Brook, 1991).  Figure 13.  Slabbing Failure  5.2.3 Buckling	Failure	Buckling failure commonly occurs when the bedding thickness is small in comparison to the beam length and the in‐situ horizontal stress is moderate to high, or if the beam is jointed.  Under these conditions, stress  redistribution within  the beam as  shearing occurs along  the bedding planes  results  in  localized stress concentrations generating a new plane of weakness across the beam.  For intact beams, fractures initiate in tension at the localized stress concentration points, which is often a result of a heterogeneity within  the  rock mass,  and  shear  forces  lead  to  coalescence  of  the micro‐fractures  into  a  diagonal fracture  across  the  beam  (Tang & Hudson,  2010).    For  jointed  rock masses,  stresses  concentrate  on existing discontinuities that extend across the beam.  As shearing initiates along the diagonal fracture or existing discontinuity, induced lateral thrust results in tensile failure in the beam abutments.  Ultimately, shearing  continues  along  the  diagonal  fracture  or  existing  discontinuity  and  corresponding  tensile fracturing propagation continues in the abutments as the moment increases leading to inward failure of the beam  into the excavation  (Figure 14).    In  thin, strong  (high modulus) beams, buckling can happen rapidly under its own weight (Diederichs, 1999). 38   Figure 14.  Buckling Failure 5.2.4 Crushing	Failure	Crushing failure commonly occurs when the bedding thickness is moderate to high in comparison to the beam  length and the  in‐situ horizontal stress  is high.   Under these conditions, failure  is governed by a combination of tensile and shear failure due to induced lateral thrusts from the beam’s deflection (Tang & Hudson, 2010).   As  the beam undergoes  flexural bending,  tensile stresses develop mid‐beam on  its underside while  conjugate  compressive  stresses  concentrate on  the  top  side,  as  shown  in  Figure 15.  Simultaneously,  tensile  stresses  develop  in  the  abutments  at  the  upper  portion  of  the  beam while compressive stresses concentrate on the lower side.  Stress in the mid‐span is approximately one half of that in the abutments and therefore the abutments fail first.  Since rock is weak in tension – estimated at only about 10% of  its  compressive  strength  (Jumikis, 1979) –  fracturing  initiates at  the  top of  the abutments and propagates downwards;  this  initial  fracturing  is  invisible  to an observer underground.  Once  the  abutments  begin  to  fail,  the  beam  becomes  essentially  self‐supporting  which  leads  to subsequent tensile fracturing at the lower mid‐span and crushing at the upper mid‐span and lower side of the abutments as stresses are redistributed within the beam.  Figure 15.  Crushing Failure 39  Diederichs  and  Kaiser  (1999)  developed  a mathematical  relationship  to  characterize  crushing  failure.  They  determined  that  fracturing  initiates when  beam  deflection  is  approximately  10%  of  the  beam thickness and ultimate  failure  is  imminent when beam deflection reaches 25% of  the beam  thickness.  These  conclusions  have  strong  practical  application  as  alarm  limits when monitoring  deformation  in underground excavations; particularly against crushing failure where  initial fracturing  is  invisible to the eye. 5.2.5 Diagonal	Fracturing	Failure	Diagonal  fracturing  failure commonly occurs when  the bedding  thickness  is high  in comparison  to  the beam  length and  the  in‐situ horizontal stress  is high.   Under  these conditions, when  the underground excavation  is  created  and  a  stable  compression  arch  can  nearly  develop  to  resist  beam  deflection, discrete tensile fracturing occurs from the lower abutments to the upper mid‐span following the inclined stress trajectories in the beam – parallel to the compression arch (Stimpson & Ahmed, 1992).  Once the diagonal fractures reduce the ability of the beam to be self‐supporting, an arched portion of the beam will fail under gravity into the excavation, as shown in Figure 16.  Figure 16.  Diagonal Fracturing Failure Stimpson and Ahmed (1992) characterized diagonal fracturing failure by conducting physical model tests in  a  laboratory  and  reproducing  the  fracturing  propagation  using  finite  element  numerical  analysis.  During their study, they also concluded that this mechanism may be  important where weak or broken material exists above the beam.  The results of this thesis will be able to support or refute Stimpson and Ahmed’s conclusion since the shale interbed will be modelled much weaker than the limestone beam. 40  5.2.6 Effects	of	High	Horizontal	Stresses	High horizontal stresses alone do not necessarily  indicate  increased risk of failure.   The geometric and boundary conditions of the structure of the rock mass, the strength of the intact rock and the stiffness of structures control the effects that the high horizontal stresses have on the rock mass.   In massive rock masses, the horizontal stresses provide confinement which has proven through experimental evidence to significantly strengthen the rock mass (Tang & Hudson, 2010).   Both the peak and residual frictional strengths  increase with  increasing confining pressure.    In  rock masses characterized by  joint sets with dips greater  than 60 degrees and a dip direction normal  to  the direction of  the maximum horizontal stress, the horizontal stresses strengthen the rock mass by increasing the normal forces acting to resist shearing  along  the  joint  planes  creating  a  “clamping”  effect  and  ultimately  resisting  slabbing  failure (Diederichs, 1999).    In  stratified  rock masses, characterized by  relatively  flat‐lying bedding  (dip angles less  than  30  degrees),  stress  concentrations  along  the  bedding  planes  are  exacerbated,  encouraging greater  extents  of  shear  failure  laterally  and  tensile  failure  vertically.    Consequently,  slabbing  and buckling failures can develop  in rock masses which under normal stress conditions (H/V≈1) would be considered  stable.    Furthermore,  in  underground  excavations  such  as  room  and  pillar mines,  stress redistribution is constantly occurring as mine development persists generating stress concentrations and stress shadows within the roof beams and pillars (Figure 17).  Large roof falls observed in room and pillar mines  often  have  an  elliptical  shape  with  the  long  axis  oriented  perpendicular  to  the  direction  of maximum horizontal stress.   Based on  this observation,  Iannacchione et al.  (2001) recommended that mines under high horizontal stress conditions adjust their mine layouts and orient their headings parallel to  the  direction  of maximum  horizontal  stress.    Therefore,  high  horizontal  stresses  can  significantly increase the risk of failure in stratified rock where bedding planes are inherently weaker and horizontal. 41   Figure 17.  Stress Concentration and Stress Shadows around a Roof Fall (Iannacchione et al., 2001) 5.3 Failure	Mechanisms	from	In‐Situ	Observations	Since the commencement of this research, the author has had the opportunity to visit two underground limestone mines  in the US which experience high horizontal stress conditions and have had  large roof failures; one  in West Virginia  (Mine A) and one  in Pennsylvania  (Mine B).   For confidentiality reasons, specific  details  of  the mines  had  to  be withheld  from  this  thesis,  but  the  observations  can  provide valuable  insight  into  the  failure mechanisms  at work  and depths of  failure under  these  conditions  – which are expected to be similar to those present in Southern Ontario. 5.3.1 Mine	A	Mine  A  was  an  underground  limestone  mine  for  high  quality  stone  and  aggregate  which  was interbedded with bands of competent shale,  friable shale and marker beds of soft, weak green shale.  Horizontal stresses at the mine were postulated at approximately 1.5 – 2 times the vertical.   The mine had recurring  issues with  large roof falls, particularly at  intersections where spans were greater.   They did  not  have  any  system  for monitoring  deflections  other  than  visual  inspection.    If  deflection was visually observed by chance, the mine would use steel straps to attempt to contain the roof fall.  Often 42  the roof fall presented no warning signs visible to the human eye since the mine workings were 12m in height and lighting in the mine was provided by head lamps or lights on mobile equipment; exceptions were at a working  face or maintenance area where  lighting was brighter.   Following a  roof  fall,  steel straps were installed to prevent further failure, as shown in Figure 18.  Figure 18.  Large Roof Fall which has been Supported with Steel Straps to Prevent Further Failure A supported roof failure was further observed to look for signs of its failure mechanisms.  It was evident that the failure was comprised of a series of Voussoir beam failures, even though the large mass failed at once, which resulted in a stepped failure profile.  The fractures across the beams were steep resembling those that are formed by the process of diagonal fracturing.   However this  is unusual according to the current understanding of  failure  in  stratified  rock masses because  in  this case  the  individual Voussoir beams were thin to moderate in thickness (thickness to span ratio equals 1/40 to 1/20) and the in‐situ stress considered moderate which  typically would have resulted  in buckling or slabbing  failure.   Upon closer  inspection  of  the  steps,  it  could  be  seen  that  they  corresponded  to  bedding  planes  or  shale interbeds and the failure ultimately terminated at a weak green shale marker bed (Figure 19).  Depth of failure was estimated at approximately 2.5m.  43   Figure 19.  Weak Interbeds Controlling Failure Profile A  second  failure  in  the mine was  observed  as  shown  in  Figure  20.    This  failure  also  occurred  at  an intersection where the span was greater, and  it presented a similar stepped  failure profile as the  first failure.  However, where the first failure was more circular in nature, this failure was distinctly elliptical.  Failure in the shortened axis orientation was observed to be prominently stepped on one side whereas failure  in  the elongated axis was observed  to be  steep with no obvious  steps.    Since  the  failure was unsupported,  a  closer  observation  of  the  steps  within  it  was  not  possible.    Depth  of  failure  was estimated at approximately 3m for this roof fall.    Figure 20.  Second Large Roof Fall in Mine A 44  The steep fractures are characteristic of diagonal fracturing which is not surprising if the failure was one 3m  thick competent beam  (beam  thickness  is approximately 1/5  failure span).   However,  the stepped failure  is  characteristic  of  a  series  of  slabbing  or  buckling  failures  where  individual  thinner  (beam thickness is approximately 1/20 failure span) competent beams fail.  Assuming the failure truncated on a weak, shale interbed similar to the first failure observed, it would appear that the weak shale layer fails before  the more competent  interbeds closer  to  the excavation and  therefore  the  failure  initiates as a thick beam at the final depth of failure; as opposed to a series of thin beams successively failing upwards to the final depth of failure.  5.3.2 Mine	B	Mine  B was  an  underground  limestone mine  primarily  for  high  quality  lime  production.    Horizontal stresses at the mine were measured multiple times with varying results, but are considered to average between 2 – 3 times the vertical.  The limestone formation mined was massive except for the upper 2m which was a transitional contact with the overlying shale formation.   The shale formation was weaker than the limestone formation and therefore the transitional contact was comprised of weak and strong interbeds.  The strata undulated across the mine with dips ranging from 5 to 15 degrees.  The mine plan was to follow the limestone – shale contact and use this interface as the excavation roof.  Unfortunately the undulation of the strata was impractical to predict and it was not uncommon for the interface to be lost while driving a header or drift.   When  the  interface was  lost and an  intersection was  created by blasting a cross‐drift into an adjacent development header, near its mining face, a large roof fall would occur, as shown in Figure 21.  The mine regularly installed “angels” – a monitoring instrument installed in the roof which releases a red visual (often a piece of flagging) when separation of the bedding in the roof is detected – as a warning system for unstable roof conditions.  However, since the large roof falls occurred  during development,  the  angels were only useful  for  indicating  if  a  failure was progressing laterally.  If an angel was triggered, then the mine would isolate that heading, drift or intersection. 45   Figure 21.  Large Roof Fall in Mine B Figure 21 and Figure 22 were two large roof falls observed by the author during a site visit.  Although the failures  could  only  be  observed  from  the  floor,  the  stepped  failure  profile  and  elliptical  shape were clearly evident  indicating that these falls were a series of Voussoir beam failures similar to the failures observed  in Mine A.   Since stress  testing had been completed at  the mine and  the orientation of  the principal stresses was generally known, it was confirmed that the elongated axis of the roof fall aligned with the minor horizontal in‐situ stress orientation.  Also, like Mine A, the failure mechanisms appear to be  a  combination  of  diagonal  fracturing  and  slabbing  or  buckling.    The  roof  fall  shown  in  Figure  21 truncated at the overlying shale contact; the overlying shale  is dark grey to black  in colour making this observation visible from the floor.  The truncating layer for the roof fall shown in Figure 22 could not as easily be discerned, however,  the  thick  light  coloured  step  (a  thick  limestone bed)  suggests  that  the truncation was  likely within  the  transition  zone.   Most  notably,  both  roof  falls  had  been  previously supported by 8ft (2.4m)  long rock bolts and the roof falls exceeded the supported depth.   Accordingly, depth of failure was estimated to be approximately 4.5m and 4m respectively for the two roof falls. 46   Figure 22.  Second Large Roof Fall in Mine B The observations made at Mine B highlight the need to understand the mechanisms and depth of failure in these environments so that adequate ground support can be installed to prevent major roof falls and safeguard personnel.  5.4 Summary	Rock  failure  and  the  complex  relationships  of  the mechanisms  driving  it  continue  to  challenge  the scientific community when it comes to precision engineering; unfortunately rock failure is not as simple as  a  stress  controlled  or  a  structurally  controlled  failure. When  excavating  in  stratified  rock masses, failure  modes  and  mechanisms,  although  generically  characterized  by  Voussoir  beams  theory,  are uniquely defined by  the beam geometries,  rock mass  strengths  (or difference of  strength) and  in‐situ stresses acting on the beam.  These conditions can alter the amount of deflection before a beam will fail together with  the depth of  failure.   Therefore,  it  is  important  to understand  the many variables  that contribute  both  through  stress  and  structurally,  as well  as  their  interdependencies,  before  one  can assess the potential for failure in stratified rock masses. 47  CHAPTER	6 NUMERICAL	MODELLING	Numerical modelling methods have been employed in the fields of civil and mining engineering since the 1960s (Mukherjee, 2003) as a tool for better understanding of the mechanisms of failure in rock and soil on  both  small  and  large  scales.    The  natural  environment  is  often more  complex  than what  can  be practically  represented  by  the mesh  and  equations  governing  the  simulation  of  failure  in  numerical modelling.   Given  this  limitation,  a myriad  of mathematical methods  have  been  developed  over  the years which adopt certain simplifications best suited for specific environments.   Each method  is based on unique closed‐form solutions or partial differential equations (PDEs) to most accurately represent the targeted  environment  of  application;  for  example,  representation  of  the  ground  conditions  as  a continuum  or  discontinuum.    Simply,  depending  on  the  soil  or  rock  strength,  and  presence  of discontinuities, heterogeneities, anisotropy and plasticity  in  the rock mass, certain modelling methods are more appropriate to use than others.  This chapter provides a brief overview of numerical modelling methods, introduces UDEC and discusses why it was selected for conducting the numerical modelling carried out in this research. 6.1 Selection	of	a	Suitable	Model	Method	Numerical modelling methods generally differ in how the way that blocks or elements within blocks are allowed to interact with one another while the modelled material deforms leading to failure or explicit movement during failure.  The most common types of numerical modelling methods can be grouped as continuum,  discontinuum  or  a  hybrid  of  the  two.    The  fundamental  concepts,  advantages  and disadvantages of each method are summarized in Table 5. 48  Table 5.  Overview of Modelling Method Groups   Continuum Methods Discontinuum Methods Hybrid MethodsFundamental Concept Represents the system as a continuous material through interconnected elements. Discontinuities are implicitly represented by degraded material properties through a homogenization process.1 Best suited for scenarios where joints have little influence on the failure mechanisms. Represents the system as a set of interacting blocks with unique contact strength properties governing block movement.  Discontinuities are explicitly represented, but may also be implicit within larger blocks.  Best suited for scenarios where joints influence the failure mechanism. Represents the system as a combination of interacting blocks that are further divisible.  Intact blocks are assigned a brittle fracture criterion. Best suited for scenarios involving brittle fracturing and fragmentation. Advantages  Geometrically simple andgenerally fast computing because adjacent elements remain connected during deformation and yield making computational processes streamlined. Separation, displacement and rotation of the blocks is permitted.  Discontinuities can be explicitly represented using laboratory joint properties to show directional weakness. Combines the continuum and discontinuum methods to reduce limitations of each method independently. Disadvantages  Separation, displacement and rotation of regions are not permitted; solution becomes unstable with large deformations.  Continuity of adjacent boundaries must be maintained.  Requires more input parameters to be calibrated and constrained (i.e., both intact blocks and joint properties).  More computationally demanding as block contacts and interactions must be tracked.  Individual blocks are continuous and cannot be fractured during failure. Potential for continuity and compatibility issues at interfaces between regions of different methods; particularly where adjacent elements employ different assumptions. Common Methods  Boundary Element Method  Finite Element Method  Finite Difference Method  Distinct Element Method  Discrete Element Method  Particle Flow Method  Boundary Element / Finite Element Methods  Distinct Element / Finite Element Methods 1Discontinuities  can  be  explicitly modelled  as  small  strain  contacts  or  as  narrow  zones  of  elements  assigned properties  to  simulate  joint  strength  if desired; however,  still no  separation, displacement or  rotation of  those boundaries can occur. 49  In an effort  to verify his  theories pertaining  to bedded  roof  rock behavior over mined  spans, Sterling (1980) conducted a series of laboratory experiments.  A key conclusion from these experiments was that roof beds cannot be simulated by continuous, elastic beams or plates since their behaviour is dominated by  the  rock mass blocks created by natural cross  joints or  induced  transverse  fractures.   The author’s observational  experience  in  US  limestone mines  as  discussed  in  Section  5.3  concurs with  Sterling’s conclusion because the author observed that the failure mechanisms in room and pillar limestone mines subject  to  high  horizontal  stresses  have  a  distinct  and  important  structural  component  controlling failure.   The role of structure  in  the  failures observed was evident because  the stepped  failure profile indicated  that a  release of blocks on bedding planes had occurred extending  to an ultimate depth of failure controlled by a notably weaker shale interbed.  Accordingly, a discontinuum method was chosen to  numerically model  the  anticipated  limestone mining  conditions  for  Southern Ontario  so  that  the structural  component  of  the  failure  mechanism  could  be  adequately  simulated.    In  particular,  the distinct  element  method  using  the  commercial  software  UDEC  was  selected  because  this  method permits modelling  of  the  block  separation,  translation  and  rotational movements  that  occur  during beam failure.   6.2 UDEC	The Universal Distinct  Element Code  (UDEC)  is  a  two dimensional discontinuum numerical modelling code which models  the  quasi‐static  or  dynamic  response  to  loads  applied  to  a  simulated  rock mass containing multiple,  intersecting  discontinuity  planes  (Itasca  Consulting Group,  Inc.,  2015).    The  rock mass  is  represented  as  an  assembly  of  discrete  blocks  (either  rigid  or  deformable)  and  the discontinuities are treated as boundary conditions between blocks.  Material strength and stiffness are applied to both the intact blocks as well as the contacts between the blocks themselves.  UDEC assesses failure within blocks against defined failure criteria and calculates displacements by converting stresses acting on each block  relative  to adjacent blocks  through  the bulk and  shear moduli  to  the associated strains.  Deformations occur through both squeezing and bending of blocks, as well as through shearing and  separation  of  the  blocks  along  their  contacts.    Additionally  in  UDEC,  the  blocks  are  free  to delaminate and separate from adjacent blocks during failure.   UDEC has a built‐in scripting language called FISH which allows the user to customize or automate many of  the programs operations.   An  automatic Voronoi  tessellation  generator  is  also  available  to  create multidirectional contacts that can be bonded together to simulate  intact rock fracturing.   Of particular value, the Voronoi network of potential intact rock fracture pathways can be superimposed with a joint 50  network representing the geological structures present in the rock mass.  UDEC was specifically chosen for conducting this research because of these features. Since  investigating  the  depth  of  failure  was  going  to  involve  iterative  modelling  under  various geometrical  and  in‐situ  stress  conditions,  the  automation  ability  of  built‐in  scripting  increased  the efficiency of the modelling by reducing the number of manual adjustments required to carry out a series of  simulations.    Secondly,  since  the observed  failures  in  room and pillar  limestone mines under high horizontal  stress  conditions  show  stepped  failures  indicating  a  combination  of  structural  and  stress controlled failure, the Voronoi tessellation would allow for the simulation of intact rock failure between the bedding planes as observed in the US limestone mines and expected under Southern Ontario mining conditions. The  FISH  script written  to  semi‐automate  the modelling  is  presented  and  discussed  in  Appendix  C.  Details of the Voronoi tessellation generator and a discussion of its use in the modelling are presented in Section 6.2.1 and in Chapters 7 and 8 respectively.   6.2.1 Voronoi	Tessellation	Generator	The Voronoi tessellation generator creates a series of randomly oriented breaks  in an existing block to form polygons ranging from perfectly shaped and evenly sized hexagons to irregularly shaped and sized blocks.    The  polygons  are  formed  by  using  the  Voronoi  algorithm  to  randomly  distribute  points throughout  a  defined  region  (which  act  as  the  centroids  to  the  polygons),  then  constructing  lines bisecting the distance between each point and its adjacent points (Itasca Consulting Group, Inc., 2015).  When  the  contacts  are  given  suitable  micro‐mechanical  intact  rock  properties,  degrading  to  joint properties after  fracture,  the Voronoi  tessellation  is useful  for  simulating  the brittle  fracture of  intact rock.   The only  criteria maintained  is  the  fracture  length  (edge  length).   The user defines an average edge  length  for  the Voronoi  joints  in which  the program ensures  that  this average  length  is achieved over the region in which the Voronoi tessellation is used.   The user can also control the degree of irregularity of the polygons by setting the iteration number.  The iteration number controls an iteration procedure which shifts the polygon centroids resulting in variable lengths and orientations of the bisecting lines (Itasca Consulting Group, Inc., 2015).  To obtain a uniform distribution of hexagons over the region being model using the Voronoi tessellation, a higher  iteration number is used (greater than 1000).  Although, utilizing regularly shaped polygons creates a directional bias in simulating the fractures because the majority of edges are aligned with the upward direction of 51  the model.    By  using  a  low  iteration  number  (such  as  1),  the  polygons  are  comprised  of  edges  in  a variation  of  directions  which  reduces  fracture  generation  bias  created  by  the  use  of  the  Voronoi tessellation generator.   The concept of directional bias due to  iteration number  is shown  in Figure 23.  Accordingly, the modelling carried out as part of this research, utilized a  low  iteration number so that fracturing could develop more “naturally”.        Figure 23.  Visual of Directional Bias due to Iteration Number: (a) Iteration Number = 1, (b) Iteration Number = 100, and (c) Iteration Number = 10,000.  Finally, where the Voronoi tessellation generated contacts are used to simulate fracturing within intact rock, the Voronoi contacts cannot simply be assigned rock properties because these are size dependent and  can underrepresent  the  strength  and  stiffness of  intact  rock.    Therefore,  they must be  assigned micro‐mechanical properties which are calibrated to reproduce the strength and stiffness of the  intact rock being modelled (Wang et al., 2013).  The micro‐mechanical joint calibration process is discussed in detail in Chapter 7. 6.3 Two	Dimensional	Versus	Three	Dimensional	Modelling	A key question in numerical modelling is whether a two‐dimensional plane strain assumption is valid or whether  the  problem  being  analyzed  is  truly  three  dimensional.    For  the  context  of  observing  the (a)  (b) (c)52  increased maximum depth of failure as a result of high horizontal stresses  in a stratified environment, the greatest  impacts of the high stress can be captured  in a  linear alignment with the direction of the high stress.  As a result, it is not necessary to have the added complexity and computing demands of a three dimensional modelling program (such as 3DEC ‐ UDEC’s three dimensional modelling companion) to  understand  the  operating mechanisms  leading  to  the  change  in  ultimate  depth  of  failure  in  high horizontal  stress  environments.    By  opting  to  employ  UDEC  for  this modelling,  a  relatively  detailed Voronoi tessellation  (described  in Section 6.2.1) could be used without simplification. This would have not been possible  in 3DEC because  the computing demands required  to handle a Voronoi  tessellation with  the  equivalent  amount  of  detail  in  the  third  dimension  would  have  exceeded  computing capabilities at the time of modelling. 53  CHAPTER	7 MICRO‐MECHANICAL	PROPERTIES	CALIBRATION	In order to monitor the generation of the brittle fracture driven mechanisms developed in the haunches of  the beams  that  lead  to  the beam  failure observed  in‐situ,  the expected  failure area  (plus  suitable additional area to eliminate boundary effects) was modelled using a Voronoi tessellation.  By using the Voronoi  tessellation  in  this  region,  intact  rock  failure  can  be  explicitly  modelled  by  allowing displacements  and  separations  along  the  Voronoi  polygons  representing  the  limestone  beam  and interbedded shale layer.  However, invoking a synthetic rock mass using the Voronoi tessellation comes with  its challenges;  the  contact properties given  to  the edges of  the Voronoi polygons must  simulate intact  rock  and  respective  strength  properties  to  achieve  an  unbiased/valid  representation.   Many papers  (Wang  et  al.,  2013;  Scholtes  et  al.,  2011;  Koyama &  Lanru,  2007)  have  reported  that Mohr‐Coulomb strength properties (cohesion, friction and tensile strength) derived from laboratory testing of intact rocks when used directly as the contact cohesion, contact friction and contact tensile strength in a numerical model analysis do not accurately  simulate an equivalent  intact  rock  strength;  the  synthetic rock mass would perform weaker  than desired.   Therefore  the  laboratory derived strength properties need  to be  calibrated  to  find  the  equivalent micro‐mechanical  strength properties  so  that  the  intact synthetic  rock  strength  can  be  properly  modelled  with  the  Voronoi  tessellation.    Determining  the equivalent micro‐mechanical strength properties  to be used  in  the Voronoi  tessellation  is not a  trivial task since the amount of adjustment of these properties  is directly related to the Voronoi edge  length used  in  the model.    The  edge  length  is  user  defined  and  based  on  the  scale  of  the model  and  the mechanisms  that are  the  target of  the modelling exercise.   Needless  to  say,  the user must know  the modelling capacity of the computer/software and have a basic understanding of the scale and type of movement  required  to  be  observed  to  achieve  the  goals  of  the  modelling  in  order  to  select  an appropriate  Voronoi  edge  length  and  begin  determining  the  equivalent  micro‐mechanical  strength properties.   This process  can be  iterative and  time  consuming depending on  the  level of preciseness required for the numerical model. This  chapter  outlines  the  unconfined  compressive  strength  (UCS)  test  simulations  conducted  to determine  equivalent  micro‐mechanical  strength  properties  to  be  used  as  the  Voronoi  contact properties in the numerical modelling analysis.  In addition, it will discuss the importance of obtaining a suitable  Voronoi  polygon  size  (via  edge  length  input)  and  the  associated  adjustment  of  contact properties required to develop a valid synthetic rock mass model in which the strength (and associated failure mechanisms) of an intact rock beam is represented adequately by the Voronoi tessellation. 54  7.1 Intact	Rock	Mass	Properties	of	Limestone	and	Shale	in	Southern	Ontario	A  review  of  published  rock  mass  properties  for  the  Ordovician  limestone  formations  in  Southern Ontario, the target units for underground limestone mining in the region, was conducted to establish a range  of  strength  and  deformation  properties  (UCS,  T,  c,  E)  to  be  represented  in  the  numerical modelling  analysis.    The  Lindsay,  Bobcaygeon  and  Gull  River  Formations  exhibit  similar  strength properties (average UCS of 110MPa) whereas the Verulam Formation, as might be expected given that it is characterized by thin limestone strata interbedded by equally thick shale, exhibits noticeably weaker strength properties (average UCS of 23MPa) (Lo, 1989).  Due to the high percentage of shale content in the Verulam Formation  resulting  in a more uniform  strength distribution between  the  limestone and shale  layers,  and  therefore  outside  the  focus  of  this  research  (i.e.  the  influence  of  a weak  layer  is absent), the Verulam strength properties were not be considered.  In addition, the Cobourg Formation is the Michigan Basin equivalent nomenclature of  the Appalachian Basin  Lindsay  Formation  (Damjanac, 2008; Ontario Geological Survey, 1992).   Since a  low  to  intermediate  level nuclear waste  repository  is planned  for  construction  in  the  Cobourg  Formation  in  Southern Ontario  about  220km  northwest  of Toronto,  considerable  amount  of  public  domain  strength  data  exists  for  this  limestone  unit  and therefore its values were also incorporated.   In  summary,  for  this  study,  a  range  of  the  reported  Lindsay,  Bobcaygeon,  Gull  River  and  Cobourg Formations  values  were  analyzed  to  establish  a  representative  range  of  intact  rock  mass  strength properties for the most economical mineable limestone units.  These intact properties will be the target range for identifying correlating joint properties to be used in the Voronoi tessellation of the numerical analysis.  The established average, minimum and maximum intact rock mass strength properties for the limestone are presented in Table 6.       Table 6.  Range of Intact Rock Mass Properties for Target Ordovician Formations  Strength  properties  for  the  weak  shale  layer  are  a  range  of  public  domain  values  published  for laboratory  testing  conducted  on  the  Ordovician  shale  units  (Queenston,  Georgian  Bay  and  Blue Mountain Formations).  Although an erosional contact exists between the middle Ordovician limestones and upper Ordovician shales, it is expected that a range of these shale strength properties will suitably   UCS (MPa) Tensile Strength, T (MPa) Young’s Modulus, E (GPa) Poisson’s Ratio,  Limestone  110 (72 ‐ 143)  6 (3 ‐10)  50 (32 ‐ 63)  0.19 (0.12 – 0.26) Shale  30 (20 ‐ 44)  10 (5 ‐ 14)  12 (2 ‐ 27)  0.10 (0.02 – 0.17) 55  represent those of the interbedded shale layers for the purpose of this study. The established average, minimum and maximum intact rock mass strength properties for the shale are also presented in Table 6.   7.2 Equivalent	Mohr‐Coulomb	Failure	Criterion	Parameters	Numerical modelling programs assess rock failure by comparing the stress applied to an element in the model to a user‐defined failure criterion for the material.  Hoek‐Brown and Mohr‐Coulomb are the most common  failure  criteria  for  intact  rock  available  in most  numerical modeling  software.    A  standard approach  in rock engineering  is to use measured  intact rock mass strength properties from  laboratory testing  and  observable  rock mass  classification  parameters  to  define  the  rock mass  strength.    This process involves carrying out a series of mathematical calculations derived from empirical relationships which  are  unique  to  the  failure  criterion  targeted.    RocLab,  a  software  distributed  by  RocScience, streamlines this process for engineers because  it accepts the user defined strength  inputs and outputs the  equivalent  Hoek‐Brown  and  Mohr‐Coulomb  intact  rock  failure  criterion  parameters.    The relationships used are based on those published by Hoek et al. (2002). Accordingly, intact rock mass properties compiled from literature (presented in Table 6) were input into RocLab to establish baseline cohesion and friction angle values; the basis of the Mohr‐Coulomb failure criterion.   The methodology RocLab  follows requires  the  input parameters  to be  in  the  form of Hoek‐Brown developed rock mass classification properties  (UCS, GSI17, mi18, D19).   The  inputted Hoek‐Brown properties and calculated equivalent Mohr‐Coulomb failure criterion parameters are shown in Table 7. Table 7.  Corresponding Mohr‐Coulomb Properties from RocLab   Hoek‐Brown Properties  Failure Envelope Range Mohr‐Coulomb Failure Parameters UCS (MPa)  GSI  mi  D 3 max20 (MPa)  Friction Angle,  Cohesion, c (MPa) Limestone  110  76  10  0.7  2.06  52˚  3.06 Shale  30  50  6  0  2.06  37˚  0.58                                                             17  Geological  Strength  Index  (GSI)  is  a  rock mass  characterization  system  developed  to  assist with  converting reliable rock mass property data into failure criterion input data for numerical analysis (Marinos et al., 2007).  For good quality rock masses (RMR > 40), GSI is considered equivalent to Bieniawski’s RMR.  18 mi  is an  intact material constant analogous  to  frictional strength which can be derived  from  triaxial  testing of intact rock (Eberhardt, 2012). 19 D is a disturbance factor which accounts for blast damage and stress relaxation (Eberhardt, 2012). 20 3 max. was calculated using a tunnel depth of 150m and a unit weight of 26 kN/m3. 56  The  Mohr‐Coulomb  failure  parameters  were  used  to  define  the  intact  rock  mass  strength  of  the limestone  throughout  the  numerical modeling  calibration  exercise  to  establish  corresponding micro‐mechanical contact properties for the Voronoi tessellation contacts. 7.3 Developing	a	Properties	Range	using	a	UCS	Simulation	The UCS of  the many  limestone  formations  in Southern Ontario  is quite well defined  from  laboratory tests conducted for the design of large quarrying operations, building foundations, nuclear facilities and waste  repositories.   Much  of  this  data  is  in  the  public  domain,  particularly  from  the  nuclear waste repository project.    Subsequently,  a model  representing  a  rock  sample  (51mm dia.) with  a  length  to width ratio of 2.5 being compressed at a rate of 0.02 m/s – a standard UCS laboratory test as shown in Figure 24 – was used to  identify a range of Voronoi edge  lengths and corresponding adjusted Voronoi contact  properties  to  represent  the  synthetic  rock mass  comprising  the  rock  beam  over  the mine excavation.   The  steel  platens were modelled  at  2.5  times  the width  of  the  rock  sample  to  prevent  axial  stress concentrations at the edges of the rock sample and ensure even loading on the rock‐platen interface.  Figure 24.  Example of UCS Test Model for Voronoi Contact Property Calibration (Voronoi Edge Length = 0.005; 4% of sample length).  57   Given that UDEC is a two‐dimensional program, this simulation is actually representing an infinitely long beam with cross‐sectional dimensions shown in Figure 24 and not a cylindrical core sample that is used in  laboratory  UCS  experiments.    Therefore  identifying  exact  Voronoi  contact  properties  to mimic  a specific UCS value is extraneous and the best that can be achieved from this study is to obtain a range in model  properties which  correlate  to  the  range  in  typical  limestone  strength  properties  presented  in Table 6 and Table 7.  To  determine  suitable  adjusted  Voronoi  contact  properties  for  numerical  analysis,  a  systematic approach was carried out  increasing the contact friction, contact cohesion and contact tensile strength values in a variety of combinations until the simulated UCS fell within the range of laboratory measured UCS.    The maximum  axial  strength  detected  at  the  rock‐platen  interface  in  the UCS  simulation was compared  to  the  reported  UCS  data  available  to  determine  if  the  tested  Voronoi  properties  of  the synthetic  rock mass  simulated  intact  rock  strength properties  from  laboratory  testing.   The  following sections  present  the model  set‐up,  systematic  approach  to modelling,  a  thorough  discussion  of  the results  and  concludes  with  the  equivalent  range  in  Voronoi  contact  properties  to  be  used  in  the limestone mining numerical study.   7.3.1 UCS	Modelling	Conditions	Figure 24 displays the dimensions of the UCS simulation where a rock specimen is situated between two steel platens; the bottom platen is fixed in both the x‐ and y‐ directions and the upper platen is fixed in the x‐direction but has a constant y‐velocity of ‐0.02m/s; this rate of compression aligns with the loading rate required to meet ASTM International D7012‐14 standard (2014) for conducting UCS tests on intact rock specimens (Park, 2004).  The contacts between the steel platens and the rock specimen are fixed in the x‐direction  to eliminate sliding along  that  interface which could be encouraged simply by Voronoi geometry  bias.    If  such  slippage  occurred  in  a  laboratory  UCS  test,  the  sample  failure  would  be considered an  invalid  test and  that  strength measurement discarded.   Both  the  steel platens and  the rock  mass  were  considered  deformable,  elastic  blocks  and  modelled  with  the  strength  properties presented in Table 8. 58  Table 8.  Model Elastic Properties   Steel Platens  Intact Rock Density ()  7750 kg/m3  2650 kg/m3 Bulk Modulus (K)21  160 GPa  30.6 GPa Shear Modulus (G)18  79.0 GPa  22.9 GPa The edges of the Voronoi polygons representing potential fracture pathways through  intact rock were given the strength properties presented in Table 9.  These values remained constant while the Voronoi contact properties  (referred  to as  joint properties  in UDEC) were varied  to determine  suitable micro‐mechanical  properties.    These  constants were  considered  appropriate  for  intact  rock  and  the  fresh surfaces created during the failure of the intact rock in the simulated UCS test. Table 9.  Constant Voronoi Contact Properties (Referred to as Joint Properties in UDEC) Voronoi Property  Constant Value Joint Normal Stiffness (kN)  1600 GPa/m Joint Shear Stiffness (kS)  160 GPa/m Joint Residual Friction  31˚ Joint Residual Cohesion  0 MPa Joint Residual Tensile Strength  0 MPa Joint Dilation  0˚ 7.3.2 Strength	Overestimation	by	Using	Elastic	Blocks	Elastic properties were chosen for the Voronoi polygons because the purpose of  introducing a Voronoi tessellation  into the modelling was to allow for explicit failure of the  intact rock mass through contact separation of adjacent blocks.  It should be emphasized that failure through brittle fracture is restricted to the Voronoi contacts only and not through the Voronoi polygons themselves.  Unlike plastic blocks in UDEC, elastic blocks do not yield internally nor deform appreciably.   As a result, at greater edge length ratios,  there  is  a higher possibility of  “pseudo‐bridging”22 which will  cause  the maximum  axial  stress required to ultimately fail the simulated intact rock to be overestimated.                                                                 21  The  bulk  and  shear  moduli  are  calculated  from  the  Young’s  Modulus  (E)  and  Poisson’s  ratio  ()  using  the  formulas:                 and   22 Pseudo‐bridging is the terminology used here because the bridging is a result of acute‐angled corners of blocks within  the rock mass preventing movement of much  larger blocks.    In‐situ  these acute‐angled corners would be crushed.  However, given that fractures cannot be generated through blocks in UDEC, this ‘pseudo‐bridging’ is able to develop.   59  7.3.3 Measuring	the	Numerically‐Simulated	UCS	Axial  stress versus axial  strain plots were generated  to monitor  the  sample’s  strength as  it deformed with the maximum axial stress representing the UCS of the sample.   The axial stress was calculated  in the numerical simulation by averaging  the axial stress measured at 23 elements across  the  top of  the sample as shown in Figure 25.  Figure 25.  Axial Stress Measurement Locations along Platen‐Sample Interface Every  500  cycles  the  axial  strain was  calculated  by measuring  the  total  vertical  displacement  of  the upper  platen  and  dividing  by  the  length  of  the  sample  as well  as  the  corresponding  axial  stress  to produce stress‐strain curves for each failure simulation.  As shown in Figure 26, most failures occur after 500,000 cycles regardless of Voronoi edge  length and therefore calculating axial strain and axial stress every 500 cycles provides more  than sufficient detail  to accurately  identify  the numerically‐computed UCS to one decimal place. 60              Figure 26.  Number of Cycles before Failure 7.3.4 Monitoring	Shear	and	Tensile	Fracture	Development	Failure  through  intact  rock  initiates  from  local  concentrations  of  tensile  stresses  exceeding  the  local tensile strength creating micro‐fractures which through the continued application of the load and stress redistribution  result  in  coalescence of  the micro‐fractures  into macro‐fractures  along which  shearing occurs (Tang & Hudson, 2010).   A correct selection of Voronoi properties should allow for this fracture generation  process  to materialize.    Ergo  the  numbers  of  tensile  fractures  and  shear  fractures  as  a percentage of total contacts created by the Voronoi tessellation were monitored to ensure that tensile fractures dominated  in  the early stages of  the sample being compressively  loaded, and as  the sample approached its ultimate strength, the shearing mechanism dominated.   The  tensile and shear  fractures were determined by having UDEC scan every Voronoi edge contact at 500‐cycle  intervals  and  count  the  number  of  contacts  that  have  normal  displacements  greater  than 0.1mm for tensile fractures and shear displacements greater than 0.1mm for shear fractures.  Tang and Hudson (2010) suggest that cracks initiate when the stresses on the intact rock sample reach 30‐40% of the sample strength in discontinuity‐free, homogenous rock.  Therefore, monitoring of the development of shear and tensile fractures as the sample approaches failure was carried out by plotting the tensile and shear fracture counts against the axial stress as a percentage of the numerically‐devised UCS.  61  The  percentage  of  the  UCS  at  which  the  sample  begins  to  develop  fractures  and  the  relationship between the numbers of stress‐induced shear versus tensile fractures were used as criteria for accepting or  rejecting synthetic  rock mass property  results.   For example,  if a set of Voronoi contact properties resulted in a UCS value within the range of actual laboratory UCS values but tensile fractures initiated at 20% of  the UCS,  as  shown  in  Figure 27a,  then  those properties were  rejected  as  suitable equivalent model  input properties.   Similarly, as  shown  in Figure 27b,  if  fracture  initiation was  tensile dominant around  at  about  30%  of  the  UCS  value  but  the  UCS  value  did  not  fall  within  the  range  of  actual laboratory UCS values, then those properties were also rejected.  62    Figure 27.  Rejected Failures: (a) UCS Test #32, (b) UCS Test #40 63  7.3.5 UCS	Testing	Results	A total of 77 numerical simulations were carried out with varying combinations of Voronoi edge lengths, contact  friction  angles,  contact  cohesions,  and  contact  tensile  strengths  to  establish  suitable micro‐mechanical  properties.    Each  simulation was monitored  to  ensure  that  failure  occurred  through  the modelled core sample, similar to a valid laboratory UCS test failure (Figure 28).  Axial stress – axial strain curves were produced to ensure a relatively  linear slope occurred prior to failure (Figure 29a).   Finally, the  percentage  of  peak  strength  in which  fracturing  initiated  in  the  sample  as well  as  the  type  of fracturing  (shear or  tensile) were monitored  to ensure  failure within  the  sample developed occurred accordingly to Tang and Hudson  (2010) extensive research on fracture propagation within  intact rocks (Figure 29b).  Figure 28.  UDEC UCS Failure Simulation Compared to Laboratory Failed Core Samples 64      Figure 29.  UCS Test #47 Failure Monitoring: (a) Axial Stress ‐ Axial Strain Curve, (b) Contacts in Failure 65  7.3.6 Relationship	between	Voronoi	Edge	Length	and	Failure	Mechanism	Bias	Intuitively,  by  using  a  greater  edge  length  (relative  to  the model  size), UDEC  perceives  the  Voronoi polygons as wedges moving amongst one another  in which  shear  strength properties  influence  crack initiation and failure.   At smaller edge lengths (relative to the model size), UDEC perceives the Voronoi polygons as  individual grains  in which  the Voronoi contacts  initially  fail  in  tension until a macro‐scale discontinuity develops  and  shearing has  a dominant  influence on  failure.   The  latter  is  the preferred mechanism  of  failure  for  the  purpose  of  establishing  Voronoi  contact  properties  that  correspond  to intact rock properties.  However, the smaller the Voronoi polygons are, the greater the demands are on the modelling computer both in computing power and memory.  Accordingly, numerous simulations were conducted using a variety of properties over a range of edge lengths to determine the maximum Voronoi edge length that could be used while tensile fracturing still initiated cracking within the sample between 30 – 40% of the peak strength.   The range of properties assessed  covered  expected  and maximum  adjusted  values  for  each  of  friction  angle,  cohesion  and tensile strength.   Figure 30.  Dependency of Mechanism Initiating Cracking on Voronoi Edge Length Ratio  66  Figure 30 shows the results of the simulations indicating that at ratios of edge length to beam length23 of less  than 8%24,  tensile  stresses  initiate  cracking, and at greater  than 17.5%,  shearing  is  the dominant mechanism for crack initiation.  This boundary appears to exist for all cases independent of the Voronoi joint properties.    Edge  length  to beam  length  ratios  from 8%  to 17.5%  represent  a  transition  region where the contacts dynamics (tensile separation and shearing) are mutually dominant.   In summary,  it can be postulated that a ratio of Voronoi edge length to a critical length measurement, in this case beam length,  in a numerical model  should always be  less  than 8% when  the desired mechanism  is  tension dominated.  The orange highlighted region in Figure 30 represents the Voronoi joint properties and edge length  combinations  that  meet  the  criteria  of  crack  initiation  occurring  between  30‐40%  of  the maximum axial strength. 7.3.7 Dependency	of	Strength	Properties	on	Voronoi	Edge	Length	Given that Voronoi edge length has a distinct impact on the mechanism of failure, numerous simulations were conducted  to assess  the  sensitivity of  the maximum axial  strength  to  the edge  length using  the same range of Voronoi contact properties as  inputs.   As shown  in Figure 31, a bi‐modal strength curve appears to exist as the edge  length ratio  increases.   The changes  in the slope of the curves that define strength  increase  closely  resemble  the  boundaries  interpreted  in  Section  7.3.5  existing  between governing failure mechanism within the model (shaded areas in Figure 30). For  the  range of properties  tested, a Voronoi edge  length  that  is 7% of  the beam  thickness would be equivalent  to  an  intact UCS  between  80  –  160MPa  (more  appropriately  closer  to  100MPa  since  the higher  strength  curves  are  those  representing maximum  friction  angle  and  cohesion  values) which  is within a typical range of UCS values for Ordovician limestone units in Southern Ontario.                                                             23  In the discussion about Voronoi edge  length ratio  (edge  length to beam  length ratio), beam  length  is the UCS sample  length  in  the  micro‐mechanical  joint  property  analysis.  Beam  length  was  chosen  as  the  normalizing dimension because ultimately  in  the depth of  failure  study,  the high horizontal  stress compresses  the  synthetic rock mass on the end of the Voussoir beam.   24 Since length of UCS sample is 128mm, Voronoi edge length is 0.01m for an edge length to beam length ratio of 8%.  67   Figure 31.  Effect on Strength by Voronoi Edge Length However, based on  the  curves, a Voronoi edge  length  to beam  length  ratio of 4%  should be used  to represent purely tensile failure.  This is evident since the strength profiles are similar at ratios less than that indicating that the tensile strength is controlling failure – the desired mechanism of failure.  Beyond 4%, the Voronoi contact property combinations which used the maximum values for friction angle and cohesion  showed  that  these  properties  significantly  contributed  to  the  strength  of  the  rock  mass indicating that the Voronoi tessellation is sensitive to shear strength properties above this length ratio.  However, at a 4% ratio, the UCS simulated is between 60 – 100MPa which is a lower range in strength than  exhibited  by  typical  Ordovician  limestone  units  in  Southern  Ontario  (as  shown  in  Table  6).  Therefore, additional simulations at a 4% ratio needed to be conducted (primarily increasing the tensile strength  input)  to  identify  suitable Voronoi  contact property  values which  correspond  to UCS  values between 72 ‐ 143MPa.  7.3.8 Suitable	Baseline	of	Voronoi	Contact	Properties	In order to establish a suitable range in Voronoi contact properties that are equivalent to the target UCS value  range  for  the Ordovician  limestone units of Southern Ontario, an  iterative process of  increasing the  tensile  strength,  friction  angle,  and  cohesive  strength  inputs was  completed.    Since  the Voronoi blocks represent an intact rock mass to be fractured in tension, it is imperative that the tensile strength is the limiting rock mass strength property.  Figure 32 (adapted from Figure 30 to show the target crack initiation  range, 30 – 40% of 72 – 143MPa)  shows  that within  the  tension‐dominated  range  (Voronoi 68  edge length ratios < 8%), the tensile strength is the limiting input as desired.  This is evident because out of  the  combinations  simulated  only  three  traces  existed  representing  the  three  different  tensile strengths tested, regardless of the cohesion and friction angle.  However, it is clearly shown even though all  combinations of properties are  capable of producing  the  strengths equivalent  to  those within  the lower range of the target zone (crack initiation pressure of 22 – 33MPa, light orange region), only those with tensile strengths of 25MPa meet the target range for average strength (crack initiation pressure of 33 – 44MPa, dark orange region).   Therefore, a contact tensile strength of 25MPa was considered the minimum tensile strength for further numerical analysis.    Figure 32.  Contact Tensile Strength Required to Initiate Suitable Tensile Cracking A series of additional simulations at a Voronoi edge  length ratio of 4% (0.005m) were carried out as a form  of  sensitivity  analysis  to  refine  a  range  in  contact  cohesion,  contact  friction  angle  and  contact tensile strength for that edge length ratio.  Since the combination of T = 25MPa, c = 50MPa and  = 50° met the lower strength criteria for acceptable equivalent Voronoi contact properties, they were used as the basis for refine a suitable range.   Initially, the tensile strength was varied from 9MPa to 55MPa while maintaining a constant cohesion of 50MPa and friction angle of 50° to determine how much the tensile strength could be  increased while still having the failure initiate in tension. As shown in Figure 33a, tensile strength can only be increased 69  to approximately 30MPa with  the constants before  failure  initiates  in  shear  instead of  tension, which indicates  that  tensile  strength  is as high as  it can be unless  the  shear properties  (c, ) are  increased.  However, peak strength values are still at the very low range of acceptable UCS values.   Secondly, the cohesion was varied from 50MPa to 110MPa while maintaining a friction angle of 50° and tensile strength of 25 MPa.  As shown in Figure 33b, additional rock mass strength was achieved through the increase in the cohesion and all failure initiated in tension. However, at cohesion values greater than 80MPa,  failure  initiation  began  below  30%  of  the  peak  strength  so  for  a  tensile  strength  of  25MPa, cohesion cannot be increased greater than 80MPa and additional strength must be gained by increasing the friction angle.  Figure 33.  Sensitivity Analyses: (a) Tensile Strength, and (b) Cohesion 70  Lastly, the friction angle was varied from 30° to 70° while maintaining constant cohesion of 80MPa and tensile  strength of 25MPa.   As  shown  in  Figure 34,  a marginally  strength  increase occurred with  the increase in friction angle.  However, friction angles above approximately 54° resulted in failure initiation below 30% of the peak strength.   Figure 34.  Sensitivity Analysis on Friction Angle This  iterative process was  continued  for a Voronoi edge  length  ratio of 4%  (0.005m) and  the  refined range of baseline Voronoi  contact properties  representing  intact  rock mass values  for  the Ordovician limestone formations of Southern Ontario are presented in Table 10.                Table 10.  Corresponding Voronoi Contact Properties Established from UCS Simulations   Friction Angle, ° Cohesion, c (MPa) Tensile Strength, T (MPa) Limestone  45 ‐ 55  80 ‐ 110  25 – 30 Results from all the UCS tests numerically modelled using the various combinations of contact friction angle,  contact  cohesion  and  contact  tensile  strength  to  establish  suitable micro‐mechanical  contact properties for the synthetic rock mass are summarized in Appendix B. 7.4 Calibration	Conclusions	From the 77 successful UCS test simulations completed to calibrate the micro‐mechanical properties for use with the Voronoi tessellation generator in UDEC, some key conclusions became evident: 71  1) Micro‐mechanical properties utilized  in  synthetic  rock mass  simulations must be calibrated  to obtain a realistic model; using  laboratory strength values directly for the Voronoi contacts will underrepresent the modelled intact rock strength.   2) Numerous  combinations  of  micro‐mechanical  properties  can  achieve  the  desire  intact  rock strength.   Without  understanding  the mechanisms  of  intact  rock  failure  and monitoring  the mechanism  of  failure  during  the  calibration,  it  is  possible  to  select  a  combination  of micro‐mechanical  properties which  do  not  accurately  represent  intact  rock  failure.    Consequently, failure  mechanisms  in  future  modelling  exercises  employing  these  “calibrated”  micro‐mechanical properties will also be inaccurate. 3) The edge length of the Voronoi contact has a significant impact on the intact rock strength and the mechanisms of  fracture  initiation  that  can  reasonably develop.    Edge  length  ratios  (edge length of  the Voronoi  contact normalized by  a  length under  axial  compression)  less  than 8% perform  like grains within a rock mass and  tensile strength  initiates  failure.    In addition, small edge  length  ratios  appear  to  require  greater  inflation  of  laboratory  strength  properties  to achieve  an  equivalent  synthetic  intact  rock  strength.   On  the  other  hand,  edge  length  ratios greater  than 17.5% perform  like a  collection of wedges  (or a blocky  rock mass) and  shearing mechanisms dominate failure.  Accordingly, large edge lengths should not be used to represent intact  rock  failure.    If  large  edge  lengths  are  used  to  represent  blocky  ground,  the  Voronoi contact properties are anticipated to be similar to laboratory strength properties based on peak strength values observed during the calibration exercise.   Since  large edge  lengths were not of interest  for  this  study and not  investigated  further,  conclusion about  them are  conceptual at best. 4) For failure through  intact rock where fracture  initiation occurs  in tension, an edge  length ratio less than 8% should be used.   Contrarily, for intact rock failure where fracture mechanisms are predominantly shear, edge length ratios greater than 17.5% should be employed. 5) Edge  length ratios between 8% and 17.5% do not appear  to have any usefulness  in numerical modelling.  This range in edge lengths seems to be a transition zone where intact rock strength values are extremely sensitive  to any change  in  input properties, which  is not practical  for an engineer performing analysis since input properties contain an element of uncertainty. 72  The conclusions from the micro‐mechanical property calibration were considered when developing the engineer  modelling  for  assessing  the  influence  of  a  weak  layer  on  roof  failure  in  room  and  pillar limestone mines. 73  CHAPTER	8 ENGINEERING	MODEL	The mechanisms contributing  to  the depth of  failure  in  the  roof of a  room and pillar  limestone mine were understood by developing  a  geometrically  accurate  scale model  in UDEC of  a  typical  room.    In doing so, the scale effects of the numerical model can be minimized resulting  in an  increased  level of confidence that the model is an accurate representation of in‐situ conditions. This chapter discusses the numerical model setup  in UDEC  including geometries, boundary conditions, rock  mass  and  joint  properties,  variable  and  assumed  input  parameters  as  well  as  the  modelling methodology  employed  to  establish  a  suitable understanding of  the  influence of  a  competent, weak layer on  the depth of  failure  in  the  roof within variable stress  regimes and where  the  limestone  is of different bedding thicknesses. 8.1 Model	Geometry	A numerical model representing the geometries of the anticipated room conditions in a limestone mine was  developed  using  ranges  of  in‐situ  stress  conditions  and  bedding  thicknesses  likely  to  be encountered within Southern Ontario limestone mines.  Since UDEC is a two‐dimensional software, and observed  failures  reported  in  the  literature  and  seen  by  the  author  were most  commonly  running diagonally at  four‐way  intersections, a  cross‐section of  the diagonal of a  room  in a  typical  room and pillar limestone mine was modelled.  The dimensions considered typical for these mines were based on the  research  published  by  Esterhuizen  et  al.  (2007) which  summarized  the mining  geometries  of  34 limestone mines  located within  the Eastern and Midwestern United States  (Table 11).   Accordingly, a 22m wide by 12m high excavation was selected for this research, as shown in Figure 35. Table 11.  Summary of Underground Mining Dimensions from US Limestone Mines (Esterhuizen et al., 2007) Parameter  Literature Average  Minimum  Maximum Mining Height (m)  11.6  4.8  38.0 Intersection Diagonal (m)  21.7  16.1  29.6 The external boundary  for  the model was extended  five  times  the  room  size  away  from  the  area of interest; i.e. 110m wide (5 times the room width, 22m) on each side of the room and 60m high (5 times the  room  height,  12m)  below  and  above  the  room.    At  this  distance  from  the  excavation,  external boundary  effects  impacting  the  stress  distributions  and  resulting  strain  around  the  room  will  be negligible.  Since limestone mines will be at depths greater than 60m and the in‐situ horizontal stresses 74  are dominant over the gravity driven vertical stresses within the mining environment, it is not necessary to  represent  the ground surface distance  from  the excavation  for  the upper external boundary  in  the model. Within the external boundary, the model has been separated  into 3 distinct areas of mesh generation, as shown  in Figure 35, so that computational demands can be managed within the model and a more detailed analysis  can be  conducted  in areas of greater  interest.   For  the purpose of  this  study, mesh discretization includes the number of discontinuities represented within a comparatively sized area and the element size within the block.  Figure 35.  Geometry and Boundary Dimensions for the Numerical Modelling The coarsest mesh corresponds  to  the area of  least  interest and provides a buffer  zone between  the external boundary and the area where stress redistribution is expected to occur (shown as “Outer Area” in Figure 35).  Within the Outer Area, only bedding at a 2m spacing is modelled.  The purpose of this was so that any regional stress redistribution along bedding planes as a result of the high horizontal stresses could occur.  Since this is not an area of primary interest, the bedding spacing was kept constant at 2m for all simulations; 2m spacing is twice that of the largest bedding spacing modelled during the analysis. 75  The “Inner Area” as shown in Figure 35 is the area in which stress redistribution as a result of the in‐situ stresses will influence the mechanisms around the excavation leading to failure.  Within the Inner Area, the bedding spacing is varied according to the simulation being modelled.  Secondly, a randomly spaced, semi‐persistent sub‐vertical joint set is modelled as would be expected in most sedimentary formations. It  is  also  important  to  note  that  the  Inner  Area  provides  a  transitional  change  from  the  complex geometry modelled  in  the  immediate  roof  rock  above  the  excavation  and  the  Outer  Area  so  that boundary effects remain negligible between the areas of high and low geometric complexity.   Figure 36.  Voronoi Tessellation Regions The densest mesh  corresponds  to  the area of greatest  interest  to  the  research and  incorporates  the immediate  roof and haunches  (shown as “Roof Voronoi”  in Figure 35).   To explicitly  show  intact  rock failure within the strata “beams” and the depth of failure in the roof, this area has been modelled using both the Voronoi tessellation generator and a discontinuity set representing the bedding thickness.  The Roof Voronoi area extends 11m above the roof of the room, 3m along the sidewalls and 5m wide of the sidewalls of the room.  To ensure that artificial boundary effects would not be generated by having small Voronoi blocks adjacent  to  large  tabular blocks,  the Roof Voronoi area was partitioned  into 7  regions 76  (shown  in  Figure  36)  in which  the Voronoi  edge  lengths were  gradually  increased  in  size.    Since  the research  is  focused  on  high  horizontal  in‐situ  stress  environments,  a  gradual  transition  was  not considered necessary in the vertical direction.  Lastly, a shale  layer of constant  thickness was  included at some distance  into  the roof  for simulations involving a shale layer.  The shale layer extended across all areas of the model (as shown in Figure 36) to ensure that no unnecessary boundary effects were introduced.  8.2 Boundary	Conditions	A  few controls were  implemented on  the external boundary of  the model  in order  to simulate  in‐situ conditions as much as possible.   Firstly, the base of the model was restricted such that  it cannot move vertically.   This was  required so  that when gravity was  included  in  the model,  the model  remained  in place.   All other boundaries were  free  to move both horizontally and vertically.   Secondly,  the  in‐situ horizontal and vertical  stresses were applied  to  the external boundary  to match  those applied  to  the elements within the model.  Finally, although the vertical stress increase due to gravity was considered negligible compared to the high in‐situ stresses, it was applied so that detached blocks in the roof would fall towards the excavation void, either bridging or collapsing into the opening.  The controls implemented on the model are shown in Figure 35: F – Fixed velocity boundary, O – Free to move boundary, S – applied stress boundary. 8.3 Voronoi	Edge	Length	From  the micro‐mechanical property  calibration  completed  (discussed  in Chapter  7),  a Voronoi  edge length equivalent to 4% of the beam length in the roof was selected to be used in the model.  Since the micro‐mechanical property calibration was carried out via a UCS test simulation which compresses the synthetic rock mass along  its  long axis, and the  in‐situ stress regime being modelled will compress the strata parallel to the bedding, the  length of the beam above the excavation was considered a suitable normalizing dimension.  Accordingly, for a span of 22m, the Voronoi edge length utilized was 0.88m. Table 12.  Selected Micro‐Mechanical Properties for Voronoi Contacts   Friction Angle, ° Cohesion, c (MPa) Tensile Strength, T (MPa) Calibrated Range of Properties  45 ‐ 55  80 ‐ 110  25 ‐ 30 Selected Properties  45  80  25 77  Table 12 shows the range in Voronoi contact properties suitable for representing intact rock mass values for  the  Ordovician  limestone  units  of  Southern  Ontario  as  established  from  the  micro‐mechanical calibration.  8.4 Rock	Mass	and	Joint	Properties	for	Numerical	Modelling	The  limestone blocks  in the Outer Area and  Inner Area and shale  layer were modelled as plastic solids with  the  ability  to  internally  fail  and  have  their movements  governed  by  residual  properties.    The limestone  blocks  in  the  Roof Voronoi  area,  however, were modelled  elastically  because  the Voronoi joints provide an explicit fracture plane for failure to occur in this rock mass.  The elastic and plastic rock mass  properties  given  to  the  limestone  and  shale  blocks  were  selected  as  the  average  rock mass properties  summarized  from  literature  (discussed  in  Section  7.1)  and  are  presented  in  Table  13.  Residual cohesion and residual tensile strength values were considered 0 because cohesion and tensile strength usually provide no significant resistance to failure once they have been exceeded.    Table 13.  Rock Mass and Joint Properties Used in Numerical Modelling  UCS (MPa) Density,  (kg/m3 ) Young’s Modulus. E (GPa) Poisson’s Ratio,  Bulk Modulus, K (GPa) Shear Modulus, G (GPa) Friction Angle,  Residual Friction Angle Cohesion, c (MPa) Residual Cohesion (MPa) Tensile Strength, T (MPa)Residual Tensile Strength (MPa) Normal Stiffness (GPa/m) Shear Stiffness (GPa/m) Block Properties Limestone  110  2680  36.0  0.19 19.4 15.1 51° 31°  2.16 0  ‐0.74  0  ‐  ‐ Shale  32  2600  17.7  0.09 7.2  8.1  37° 31°  0.65 0  ‐0.12  0  ‐  ‐ Joint Properties Bedding  ‐  ‐  ‐  ‐  ‐  ‐  31° 31°  0.57 0  0.95  0  160  16 Sub‐vertical Set  ‐  ‐  ‐  ‐  ‐  ‐  45° 31°  0.57 0  0.95  0  160  16 Micro‐Mechanical Voronoi Contact Properties Limestone   ‐  ‐  ‐  ‐  ‐  ‐  45° 31°  80  0  25  0  160  16 Itasca Consulting Group (2008) summarized a series of direct shear tests conducted on bedding planes of Ordovician  limestones  for  the Ontario  Power Generations’ Deep Geologic Repository  for  Low  and Intermediate Level Waste project.   Direct shear tests provide values of the friction angle and cohesion for the discontinuities tested.  The bedding plane properties for the modelling were chosen to mimic the 78  Itasca  reported  friction  angle  and  cohesion  values.    Based  on  these  values,  the  associated  tensile strength for the bedding planes was calculated using the Coulomb slip law,   c + NtanWhere   ‐ shear stress   c = cohesion   N = normal stress (or negative tensile strength)    = friction angle  Very  little  laboratory  joint  strength data was  located  for  sub‐vertical  joint  sets within  the Ordovician limestone  formations  and  therefore  these  discontinuities  were  given  similar  cohesion  and  tensile strength as the bedding planes.   However based on the author’s experience, sub‐vertical  joints sets  in these  environments  often  have more  irregular  joint  surfaces  and  therefore  for modelling  the  joint friction angle was  increased  for  this set  to make  these discontinuities more  resistant  to  shear  failure.  The bedding plane and sub‐vertical joint properties used for the simulations are presented in Table 13. The lower bounds of the range of micro‐mechanical properties considered reasonable to represent the intact  properties  of  limestone were  selected  for  the modelling  so  that  an  element  of  conservatism would exist in the results since there is no reasonable method to measure the factor of safety when the mode  of  failure  is  complex  such  as  roof  failures  in  stratified  deposits  under  high  horizontal  stress conditions.    Stiffness  and  residual  strength  properties  were  chosen  the  same  as  the  other  joint properties because once intact rock mass failure has occurred, failure resistance is primarily governed by the residual friction angle.  The micro‐mechanical Voronoi joint properties used for the simulations are presented in Table 13. 8.5 Variable	Parameters	The in‐situ stress ratio, bedding spacing within the limestone and the distance of the shale layer in the roof above the excavation were varied to represent the most probable range of mining conditions for Southern Ontario.  The range in each parameter modelled as part of understanding the influence of the shale layer on the depth of failure in the roof are discussed in the following sections. 8.5.1 Stress	Field	Within  a  horizontally  dominated  stress  regime,  there  can  be  significant  local  variation  in  the  in‐situ stress ratio, k (horizontal stress, H: vertical stress, V), as a result of local stress relief features or simply 79  changes  in  depth  causing  an  increase  in  the  vertical  stress  (Esterhuizen  et  al.,  2007; Dolinar,  2003).  Therefore,  the  research conducted simulated various stress  regimes centred about  the  regional stress field discussed in Section 4.1.  In recapitulation, within 200m of surface, the maximum horizontal stress averages 3.8 times the vertical stress.  Below 200m depth, the average horizontal to vertical stress ratio halves  to  2.0.    Based  on  the  geology  of  Southern  Ontario,  the  Ordovician  limestones  with mining potential are situated between 100m and 300m below surface.   Therefore, horizontal to vertical stress ratios, k, of 1H:1V, 2H:1V, 3H:1V and 4H:1V were simulated during the numerical modelling.  A ratio of k=1H:1V was simulated for two reasons: 1) to compare results with those published in literature under isotropic stress conditions, and 2) to present a continuous understanding of the  influence of the shale layer  from  isotropic  (no  dominant  stress  orientation)  to  anisotropic  (horizontally  dominant  stress orientation) stress conditions. The absolute value  for  the vertical stress used  in  the modelling was calculated  for a mine situated at 150m depth and using a rock mass density of 2650kg/m3 for the overlying bedrock.  The vertical stress was calculated using the formula, V = density ( x gravity (g) x depth below surface (h) A depth of 150m was  chosen as  the  constant  for  stress  calculations because  literature  suggests  that stress ratios where the horizontal stress is greater than 3 times the vertical stress exist primarily above 200m  depth  and  the  Bobcaygeon  Formation  (Ordovician  limestone  formation  in  the  eastern  part  of Southern Ontario considered an economic source of aggregate using underground mining techniques for extraction)  is situated approximately 150m below ground surface.   Accordingly,  this depth equated  to calculated  stress  values  considered  reasonable  for many  anticipated mining  conditions  in  Southern Ontario.   The absolute horizontal stress values were calculated based on  the desired stress  ratio  for  the model and  the vertical stress value.   The absolute stress values simulated  for each  stress  ratio are shown  in Table 14. Table 14.  Absolute Stress Values for each Stress Ratio Simulated Stress Ratio, k  Vertical Stress , V, (MPa)  Horizontal Stress, H, (MPa) 1H:1V  4.0 MPa  4.0 MPa 2H:1V  4.0 MPa  8.0 MPa 3H:1V  4.0 MPa  12.0 MPa 4H:1V  4.0 MPa  16.0 MPa 80  As  a  form  of  calibration,  the  calculated  absolute  stress  values were  compared with  absolute  values obtained by Haimson and Lee (1980) from in‐situ stress measurements conducted in a deep borehole in Bowmanville, Ontario using hydraulic fracturing methods, where the stress ratio  is between k=2.5 and k=3.5.   Results of  the  in‐situ measurements  indicated  that  the maximum horizontal  stress within  the Ordovician limestone formations varied between 10.6 and 15.4MPa, which is in good agreement for the respective stress ratios simulated. 8.5.2 Bedding	Spacing	The bedding spacing  in  the strata above  the excavation equals  the  thickness of beam  in beam  theory calculations when assessing beam deflection limits.  As discussed in Section 5.2, beam theory states that a thicker beam has better flexural resistance than a narrow beam.   Subsequently, a greater  limestone bedding spacing is hypothesized to naturally reduce the influence that the shale layer has on the depth of failure in the roof.   Shale, on the other hand, has a bedding spacing which  is often on the scale of millimeters; although a shale  interbed may be thicker.   However, as an  interbed, shale does not form a competent beam that would provide any  resistance  to beam  failure.   Therefore,  for  the purpose of  this  research,  the  shale layer was modelled as a homogenous interbed without any internal bedding.  Also, since the shale layer is relatively much weaker than the limestone and both literature and the author’s observations of these types of  roof  failures  suggest  that  the entire  shale  interbed  fails  regardless of  its  thickness,  the  shale layer was modelled as having a constant thickness instead of a range of thicknesses.   In  order  to  identify  a  suitable  range  of  limestone  bedding  spacing  and  average  thickness  of  shale interbed  anticipated  within  the  Ordovician  limestone  formations  of  Southern  Ontario,  a  literature review was conducted.  During the review, bedding spacing and interbed thickness were compiled into groups – very thin (0.1 – 5cm), thin (5 – 20cm), medium (20 – 50cm), thick (50 – 100cm), massive (100 – 300cm), very massive  (>300 cm) – representing their respective dimensions. The results of the review are compiled in Figure 37.     81    Figure 37.   Summary of Ordovician Limestone Bedding and  Interbed Literature Review: a) Limestone Bedding Spacing, b) Shale Interbed Thickness 82  From  the  107  limestone  and  31  shale  data,  a  range  of  bedding  spacing  and  the  average  interbed thickness were selected for the numerical modelling simulations, as shown in Table 15.  The range was selected  by  assuming  the  percentage  of  cases  per  formation  having  the  same  spacing/thickness correlates with  the  probability  that  the mining  conditions  in  Southern Ontario will  exhibit  the  same bedding spacings and interbed thickness. Table 15. Bedding Spacing and Shale Layer Thickness for Numerical Modelling Simulations  Limestone Bedding Spacings  0.1m, 0.2m, 0.5m, 1.0m Shale Average Interbed Thickness  0.2m  8.5.3 Distance	of	Shale	Layer	from	Excavation	Roof	Initially, baseline simulations were conducted without a shale layer to establish the minimum depths of failure for each stress regime and bedding spacing combination.  Subsequent modelling then included a shale  layer  situated  at  varying distances  above  the  roof of  the  excavation,  starting  at  2m  above  the excavation  and  increasing  in  1m  intervals  to  8m  above  the  excavation.    A  one  meter  interval  for adjusting the shale layer between models was selected to match the largest bedding spacing modelled.  Since  the primary goal of  the  research  is  to understand how  failure mechanisms change as  the  shale layer impacts the system and ultimately influences the depth of failure, as well as being cognizant of the limitations of the modelling, knowing the precise distance from the roof at which the shale layer ceases to influence the depth of failure is not necessary; or reasonably valid for that matter.   This research did not include numerical analysis of the influence the shale would have within 2m of the excavation  because  from  the  author’s  experience,  in  high  horizontal  stress  environments,  the combination of blast damage with a weak layer in close proximity to the roof would result  in excessive overbreak during excavation.   Secondly, most mining operations utilize a mechanic scaler which would easily scale down loose slabs impacted by a weak layer within 2m of the roof.  Finally, standard ground support for intersections involves a minimum of systematic bolting using 2.4m long rock bolts (and mesh for  permanent  travelways)  which  in  most  cases  would  be  sufficient  to  retain  any  slabs  remaining following excavation and scaling procedures.     83  8.6 Conditional	Assumptions	Other environmental conditions were assumed negligible  for  the purpose of conducting  this  research because of two main reasons: 1) they are not mentioned in the literature as factors associated with high horizontal stress mining failure, or 2) they have  such  a  strong  influence on  failure  that  they  effectively  change  the mechanisms of failure preventing the influence of the shale layer from being detectable.   These conditions and their assumptions are discussed in the following sections. 8.6.1 Influence	of	Water	It  is well‐known  that water  has  detrimental  effects  on  the  rock mass  quality  and  joint  conditions  in sedimentary environments; particularly its ability to cause shales to swell and wash out.  However, it is assumed  that once mining has commenced  that water would be controlled by drainage and pumping systems and therefore any flow transmitting through the roof strata would be minimal, if any, and have negligible impact on failure compared to the impacts of the high stress acting on the bedding planes.   Effects from changing air moisture content (humidity) are not negligible and they are discussed in detail in Section 9.5. 8.6.2 Pillar	Stability	Although  pillar  stability  is  a  concern  for  fully  benched  room  and  pillar  limestone mining  operations where  horizontal  stresses  are  less  dominant  (or  anisotropically  dominant),  this  aspect  was  not addressed  in  this  research.   The height of  the excavation was chosen  to align with  the average  room height  from  the  survey  conducted by Esterhuizen et al.  (2007) of 34 aggregate mines  in Eastern and Midwestern US.  It is assumed that pillars of this height would remain competent over time and that the horizontally dominant stress  regimes modelled mitigate pillar stress by concentrating  the stress along the relatively horizontal bedding planes and creating stress shadows in the centre of the pillars. 8.7 Discretization	of	Discrete	Blocks	The discrete blocks within the model were made deformable by discretizing them into zones using both the quadrilateral and  triangular meshing  function  in UDEC.   Movement of a  zone  relative  to adjacent zones within a block will be calculated using finite‐difference theory.  The quadrilateral meshing function automatically  generates  diagonally  opposed  triangular  zones  in  blocks with  greater  than  3  corners.  Where possible, blocks were discretized using  the quadrilateral meshing  function because plastic  flow 84  calculations are  improved with this style of discretization (Itasca Consulting Group,  Inc., 2015).   Blocks which could not be discretized using the quadrilateral meshing function because they were defined by only 3 corners were discretized using the triangular meshing function.  The triangular meshing function will work for all block shapes. Regardless of the meshing function, an edge length characterizing the size of the discretized zone must be  defined  by  the  user.    A  fine mesh,  created  by  using  a  small  edge  length, would  allow  for more detailed plastic yielding and deformation  representation within  the model.   However, a  finer mesh  is more computationally demanding than a coarse mesh, created by using a larger edge length.  Since the numerical model developed  for  this  research employs  the Voronoi  tessellation  generator  resulting  in computationally  intense  calculations,  a  coarse  mesh  that  would  properly  discretize  all  blocks  was chosen.   To optimize  the discretization within  the model, coarser meshes were used at distal  regions from  the  room  excavation  and  finer meshes were  used  in  the  immediate  vicinity  of  the  excavation where movement is anticipated. 8.8 Establishing	an	Equilibrated	Model	Prior to commencing any numerical simulations, the model must be brought to an initial equilibrium or static state under the boundary and  in‐situ stress conditions.   In an effort to expedite this process, the model was first run through 30,000 cycles with all discrete blocks being modelled as elastic so that no permanent  plastic  deformations  would  be  present  before mining  of  the  excavation  was  simulated.  Secondly,  the  desired  plastic  properties  were  subsequently  assigned  to  the  discrete  blocks  after initialization  (as discussed  in Section 8.4) and  the model was allowed  to equilibrate  for an additional 40,000 cycles. After  the  combined  70,000  cycles,  the  number  of  unbalanced  forces  had  exponentially  decreased towards  an  asymptotic  value  considered  negligible  for  modelling  (as  shown  in  Figure  38)  and  the displacement vectors depicting total displacement indicated that movements due to model initialization were  less than 5mm  (as shown  in Figure 39).   Accordingly,  it was concluded that the  initialized model had  successfully  reached an equilibrated  state and  that  the model was  suitable  for continued use  for conducting the required numerical simulations. 85   Figure 38.  Unbalanced Forces Indicating Model Instability Negligible  Figure  39.    Displacement  Vectors  Showing Movement  Less  Than  5mm  (Within  Acceptable  Range  for Model Equilibrium) 86  8.9 Validation	Model	against	Current	Literature	A  final validation of the base model was conducted by simulating a homogenous  limestone rock mass under  isotropic  in‐situ  stress conditions  for a variation of bedding  thicknesses  (0.1m, 0.2m, 0.5m and 1.0m)  and  comparing  the  results  to  the minimum  thickness of bedding predicted by Hutchinson  and Diederichs  (1996)  in which  stability would  be  achieved.  Hutchinson  and  Diederichs  (1996)  prepared charts  based  on  snap‐thru  (buckling)  and  crushing  modes  of  failure  which  correlate  thickness  of lamination, rock mass modulus and UCS to a maximum stable span, as shown  in Figure 40.   Since the UCS of the  limestone represented  is 110MPa, under  isotropic  in‐situ stress conditions, crushing failure should not occur according to Hutchinson and Diederichs (1996).   However, for a span of 22m and the rock mass modulus  of  36GPa  – which  remained  constant  throughout  the  research  ‐  the  anticipated minimum bedding thickness for stability was back‐calculated to be 0.67m; buckling failure  is predicted when modelling thinner bedding thicknesses than this. 87   Figure 40.  Span Charts for Preventing Buckling (Snap‐Thru) and Crushing Failure (Hutchinson & Diederichs, 1996) 88  The  equilibrated model  under  isotropic  stress  conditions  showed  failure  in  all models.    For models where bedding spacing was less than 0.5m, multi‐beam failure occurred to a maximum depth of failure of 0.4m.   Where the bedding spacing was 0.5m or greater, failure did not occur.    In all models, failure occurred by buckling  (snap‐thru) as  shown  in Figure 41. These  results are  considered  consistent with Hutchinson and Diederichs (1996) failure chart.   Figure 41.  Example of the Progression of Buckling Failure: (a) Failure in Motion, (b) Complete Failure.  89  CHAPTER	9 NUMERICAL	MODELLING	RESULTS	A total of 132 models representing various bedding thicknesses,  in‐situ stress regimes and distance of shale  layer  from  the  roof of  the excavation were conducted  to understand  the  influence  that a weak layer has on the maximum depth of failure over rooms in room and pillar limestone mines. This chapter presents the results of the numerical modelling and highlights observed trends in the data unique  to each  stress  regime and bedding  thickness. Also  the mechanisms  that govern  failure as  the shale layer influences the depth of failure are discussed.   The collected data for each of the 132 simulations is presented in Appendix A. 9.1 Base	Case	Models	without	a	Shale	Layer	Twenty  simulations were  completed  for  the  range  in  bedding  thicknesses  and  in‐situ  stress  regimes modelled to act as the basis from which the shale  layer models can be compared to.   Figure 42 shows that depth of  failure  is under 1m except under unusually high horizontal stress conditions  (k ≥ 3).   All failures are less than 2m in depth.  Figure 42.  Base Case Simulations (without Shale Layer) The  failure  depths  presented  in  Figure  42  align with  the majority  of  failures  observed  in  limestone mines.    This  is  likely  why  the  industry  is  having  such  difficulty  understanding  large  roof  falls; 90  conservatism  in  ground  support  design  charts  (including  bolt  length  formulas)  recommend  ground support that will adequately support failures  in all  in‐situ stress environments.   For example,  if ground support  is designed using  industry  standard design  curves published by Grimstad  and Barton  (1993), then  for  a  conservative  range  in  rock mass quality  for  the Ordovician  limestone  formations  (poor  to good rock), reinforcement between systematic bolting and 50 – 90mm of fibre reinforced shotcrete with bolting  is  recommended  (Figure  43);  these  recommendations  assume  that  the  mine  intersections modelled are permanent openings.    Figure 43.  Empirical Ground Support Design Chart (Grimstad and Barton, 1993) Bolt length design in industry is guided by Barton et al. (1980) formula,  For mine  intersections of 22m  span,  the  recommended bolt  length  is 4.1m, which  is 2.1 beyond  the greatest depth of  failure.   For practicality reasons, mines rarely use 4.1m  long bolts due to challenges with  installation and  limited usefulness  throughout  the mining operation;  they often use 2.4m or 3m long  rock bolts  and decrease  the pattern  spacing  to make up  for  the  shorter  length. Therefore even systematic  bolting with  commonly  used  bolt  lengths,  provides  adequate  ground  support  in  the  vast majority of mining conditions.   As such, large roof falls appear to be anomalies which are unpredictable 91  to  the  industry.    However, when  both  stress  and  structure  (interbedded weak  layer)  combined  are analyzed, large roof falls are not as anomalous under certain mining conditions.  9.2 Shale	Layer	Models	A  total  of  112  simulations were  carried  out  to  observe  the  influence  that  a weak  layer,  such  as  an interbedded  shale,  has  on  the  depth  of  failure  as well  as  the  failure mechanisms  leading  to  failure.  Figures 44 through 47 present the results of the modelling by stress regime and comparing the distance between the roof and the shale layer to the depth of failure.  These graphs were further interpreted to identify the zone of shale influence and zone of no shale influence.   In general, the simulated data shows that beyond some critical distance that the shale layer is from the roof  of  the  excavation,  it  does  not  influence  the  depth  of  failure.    Also,  as  the  bedding  thickness decreases and stress ratio increased, the shale had a greater influence on the depth of failure.                Figure 44.  Modelling Results for k=1 At a stress ratio of 1 (Figure 44), for the range  in mining conditions modelled, the  influence of shale  is minimal and typical ground support design would have a high probability of preventing any  large roof falls. At a stress ratio of 2 (Figure 45), the influence of the shale becomes quite prominent for thinly bedded strata extending the depth of failure beyond the length of typical ground support installed in limestone 92  mines.  Under these stress conditions, a shale layer within 5m of the roof has potential to influence roof stability.  Although, bedding thicknesses of 0.5m and 1m remain minimally influenced.                Figure 45.  Modelling Results for k=2                 Figure 46.  Modeling Results for k=3 93  At a stress ratio of 3 (Figure 46), the entire range of bedding thicknesses modelled are influenced by the weak  layer when the shale  is  less than 4m  from the roof of the excavation.   Thinly bedded strata can have  failure depths exceeding 6m  from the roof of the excavation  if the conditions are right.   Ground support design  strategies need  to  change  in  these environments  to mitigate  large  roof  falls;  thinly  to moderately bedded strata will need to be supported with closely spaced cable bolts or a combination of cable bolts and shotcrete as opposed to systematic rock bolting.  Finally, at a  stress  ratio of 4  (Figure 47),  the entire  range of bedding  is  influenced by  the weak  layer when the shale is <5m from the roof.  Thinly bedded strata can have failure depths exceeding 7m from the  roof of  the excavation.   Ground support design strategies  for all bedding  thicknesses will need  to change  to mitigate  large  roof  falls;  a  combination of  cable bolts  and  shotcrete will  likely need  to be required as opposed to systematic rock bolting.  Figure 47.  Modelling Results for k=4 Zone of No Shale Influence Zone of Shale Influence 94  9.2.1 Shale	Layer	Modelling	Observations	Throughout the numerical modelling, a few general observations were made that were not intuitive and as such may have some significance in interpreting the data.  Primarily, the depth of failure in the zone of no shale influence during the shale layer modelling only occasionally aligned with the depth of failure under the same mining conditions during the base case modelling; the depths of failure were similar but not identical.  This observation provides an indication that even though the shale layer does not appear to have an influence on the depth of failure, there may indeed be some minor influence.  From another perspective, the small variance (less than 0.20cm) between the depths of failure may simply be a result of model boundary conditions with the introduction of elements of differing properties through the area of interest. Secondly,  at  the  lower  stress  regimes,  it  was  common  for  the  depth  of  failure  between  bedding thicknesses of 0.1 and 0.2 to be similar, contrary to their significant differences at higher stress regimes.  This  could  be  an  indication  that  the  thickness  of  beam  required  to  establish  a  compression  arch  is directly  proportional  to  the  in‐situ  stress  regime;  high  stress  regimes  require  greater  thicknesses  of bedding to achieve stability.  If this observation has merit, then in‐situ stress conditions will need to be incorporated  into  the  mathematics  defining  compression  arches  (Diederichs  and  Kaiser,  1999)  and defining stable beam thicknesses (Hutchinson and Diederichs, 1996). 9.3 Shale	Influence	Boundaries	A maximum depth of failure was interpreted for each stress regime and bedding thickness combination from  the  simulated  data  to  develop  a  series  of  curves  defining  the  boundary  in which  a  shale  layer significantly influences the depth of failure in the roof.  These curves are shown in red in the Shale Layer Influence Chart  (Figure 48)  in addition  to  the  curves  representing  the depth of  failure  curves when a shale layer does not influence the depth of failure (shown in purple). 95   Figure 48.  Shale Layer Influence Chart This  chart was developed  as  a  risk management  tool  to  assist  geotechnical  engineers with design of limestone mines under high horizontal  stress  conditions.    For  a  given  set of  in‐situ  stress  and  strata structural conditions (bedding thickness), the geotechnical engineer can use this chart to determine the critical depth of competent limestone in which they must ensure exists over the roof of an intersection in  order  to mitigate  the  risk  of  large  roof  falls  and minimize  ground  support  costs.  For  example  (as shown by the blue dashed lines in Figure 48), if the limestone bedding thickness is 0.5m and the in‐situ stress  ratio  (k)  is 3,  then a shale  layer within a distance of 5m  from  the  roof will cause  failure  to  the depth of the shale layer (i.e. if the shale layer is 4m from the roof than depth of failure would be 4m).  If there is no shale layer within 5m, then the depth of failure expected will be approximately 1m.  In order to  minimize  ground  support  requirements  in  this  case,  a  geotechnical  engineer  would  conduct investigation to determine if a weak shale interbed exists within 5m of the planned mining horizon, and if  so,  would make  recommendations  for  adjustment  to  the mine  plan  or  design  a  suitable  ground support system for the conditions. 96  9.4 Mechanisms	and	Modes	of	Failure	The main failure in each mining condition simulated was classified based on the 4 primary failure modes in  stratified  rock  (slabbing,  buckling  crushing  and  diagonal  fracturing).    This  data  was  gathered  as justification to support or reject Stimpson and Ahmed’s (1992) conclusion that weaker  layers  influence the  failure  mechanisms  within  stratified  rock  masses.    Furthermore,  it  will  confirm  the  author’s observation  from  large roof  falls  in US  limestone mines characterized by  thinly  to moderately bedded strata;  in  these  environments,  the  author  postulated  that  the  failure  profile  resembled  diagonal fracturing as opposed to a series of buckling or slabbing failures as current literature suggests should be the mechanisms of failure. A summary of the main failure classification is provided in Table 16.  These results strongly indicate that buckling failure  is the dominant mode of failure when a shale  layer does not have an  influence on the failure.    In  all  cases when  a  failure occurred when  a  shale  layer was present but  the  failure did not extend to the shale layer, buckling was the main failure mode.  Slabbing was only observed at low stress ratios (k = 0.5, 1) where a shale  layer was not present.   The most notable statistic though  is that  in all cases when failure extended to the shale layer, failure was by diagonal fracturing.  This indicates that a weak  layer  inclusion  in  a  rock  mass  not  only  influences  the  depth  of  failure  in  the  roof  but  also influences  the mechanisms  that  lead  to  failure, which  fully  supports  Stimpson  and  Ahmed’s  (1992) conclusions.  Table 16.  Summary of Main Failure Mode Classification Condition Failure Mode ‐ # of observed Total Simulations Slabbing Buckling Crushing Diagonal Fracturing No Failure No shale layer  5  11  0  0  4  20 Failure extended to shale layer  0  0  0  42  0  42 Failure did not extend to shale layer  0  47  0  0  23  70 Examples of the dominant failure mode observed for the each condition are shown in Figure 49. 97     Figure 49.  Examples of Failures Observed: (a) Failure to Shale Layer by Diagonal Fracturing, (b) Failure with No Shale Layer Present by Buckling, (c) Failure Not Extending to Shale Layer by Buckling 98  It should be noted that there could be a bias by using the Voronoi tessellation which results in buckling occurring  as  opposed  to  slabbing  or  crushing.    Since  UDEC  does  not  allow  for  fracture  generation through  blocks  and  the Voronoi  tessellation  provides  angular  failure  planes  through  intact  rock,  any crushing failures would present as buckling failures and some slabbing failures can present as diagonal fracturing failures.   However, since slabbing and diagonal fracturing have distinctly opposite conditions under which  their mechanisms  develop,  as  shown  in  Figure  50,  confusion  between  these modes  of failure would  only  exist  near  their  boundary  (beam  thickness  to  span  ratio  is moderate  and  stress conditions relative to the rock strength is also moderate).  Figure 50. Conditions Leading to Different Failure Modes in Stratified Rock Masses 9.4.1 Shale	Layer	Influence	on	Stress	Redistribution	The  principal  stress  distribution  around  the  room  was  observed  to  understand  the  role  of  stress redistribution  on  the  mode  of  failure  and  the  development  of  a  stability  arch.    Intuitively,  stress concentrations primarily occurred  along bedding planes  and  the  relatively weaker  shale  layer.    In  all cases, stability was achieved by developing an arch through multiple beams. However, between bedding planes, stress redistribution occurred differently depending on the mode of failure. Figure 51 shows examples of different stress redistribution patterns observed in achieving a stable arch. 99     Figure 51.  Stress Distribution around Room: (a) Stable Arch (No Shale Layer), (b) Shale Layer Failure with Stable Arch Attempting to Form Below Shale, (c) Shale Layer Failure 100  Figure 51a shows that when a stable arch developed over the excavation, stress vectors are shallowly angled.  In contrast, Figure 51c shows that when failure extended to the shale layer, stress vectors were steeply angled  from  the abutments  to  the  shale  layer. Also,  the  stress vectors aligned with  the steep failure profile observed during diagonal  fracturing  failures.    It appears  that  the  influence of  the  shale layer on stress redistribution is the reason that all failures extending to the shale layer failed by diagonal fracturing.  Finally, Figure 51b shows conditions where a stable arch was attempting to form below the shale layer, but ultimately was unable to and failure occurred to the depth of the shale layer via diagonal fracturing.  This stress distribution pattern was interpreted as failure conditions in which the shale layer was close to the critical distance from the roof in which the shale layer would not influence the depth of failure. 9.5 Limitations	of	Numerical	Modelling	Results	The major  limitation of the research  is that the modelling needed to be simplified to two dimensional analysis so  that  the Voronoi  tessellation could be used  to observe  the  failure mechanics  that develop through  intact  rock.   As  a  result  any  influences  to  stability  in  the  third dimension,  such  as  clamping effects that would be induced as a result of the secondary principal stress also being horizontal, would have  been  neglected  during  the  study.   As Agapito  and Gilbride  (2002)  concluded  in  their  research, anisotropic  stress  conditions,  where  the  maximum  horizontal  stress  was  greater  than  3  times  the minimum horizontal stress, can result in roof failures depending on the orientations of the workings.  In addition, they concluded that the shear strengths of the rock mass are often  lower  in anisotropic rock masses than isotropic ones.  Although the measured stress regime in Southern Ontario indicates that the maximum  to minimum horizontal  stress  ratio  is approximately 2,  there  could be  some merit  in using three dimensional software to understand highly anisotropic horizontal stress conditions (maximum to minimum horizontal stress greater than 3) to confirm the conclusion of the two dimensional modelling.  Secondly, the  influence of adjacent mine workings has not been captured  in this research.   Room and pillar mines have adjacent mine workings  in every direction so these workings will control  local stress concentrations  and  shadows.    Both  stress  concentrations  and  stress  shadows  have  potential  to exacerbate the influence of the shale layer and reduce it.  Thirdly, water (including air moisture) and the element of time were not included within the modelling.  Each simulation was run until unstable forces within the model became negligible.  This is reasonable for the  study  performed.    However,  carbonate minerals  are  known  to  degrade  through  reactions with 101  water.   Mine B visited by  the author had  reported  increased  roof  falls during  times of high humidity within the mine (summer months); so much that at one point they injected a water sealant into the rock mass over a critical  intersection which was exhibiting signs of progressive roof failure to  inhibit further failure.   At  the  time  of  the  site  visit,  they  had  concluded  that  it  indeed  slowed  the  development  of failure.  But the injection process had cost them over a million dollars to implement.   In  addition,  by  not  having  a  time  reference  for  the modelled  failures,  the  length  of  time  between underground excavation and roof failure is unknown and therefore it is difficult to understand the risk of these roof falls.  If a roof fall occurs during blasting, then the risk to people and equipment is minimal.  There  is  just a  financial  implication with overbreak;  this assumes  that  the roof stabilizes  following  the initial failure.  But if the roof fall occurs 3 – 12 months later, then the risk that people and equipment are working under the unstable ground conditions is much higher.  It  would  be  naïve  to  consciously  omit  the  potential  impacts  of  water  and  time,  if  computationally possible to include them in the modelling. Finally, the research only focused on a narrow range of structural orientation and rock mass and  joint properties.  These input parameters have substantial influence on the performance of the model.  As a result, the research conclusions are only valid for a specific rock mass environment.  For example, if the bedding  was  dipping  at  10  degrees,  stress  redistribution  would  occur  differently  in  the  roof  and mechanisms of failure would in turn develop uniquely.  Buckling failure is further encouraged in dipping strata so it can be anticipated that an even greater percentage of failures would be classified as buckling failure than were concluded in this research. 102  CHAPTER	10 CONCLUSIONS	AND	RECOMMENDATIONS	This chapter provides a final summary of the research conducted as well as outlines recommendations for future study. 10.1 Influence	 of	 a	 Weak	 Layer	 on	 the	 Depth	 of	 Failure	 in	 Stratified	 Rock	within	High	Horizontal	Stress	Environments	Intuitively, high horizontal stresses are associated with clamping effects which result in greater stability with greater stress.   However,  this  research shows  that  in stratified environments  interbedded with a weak  layer, there  is greater  instability with greater horizontal stress.   This research concluded that the maximum depth of failure when a weak layer is present could increase as much as four times compared to  that  for a homogenous  rock mass.   The comparison  in depth of  failure  from  rock masses with and without the influence of a weak layer is summarized in Figure 48. Another significant conclusion is that all failures observed where the shale layer has an influence on the depth of failure were via diagonal fracturing. This observation agrees with Stimpson and Ahmed’s (1992) conclusions  from  their  laboratory  and  numerical  study, which  used  discrete  crack  propagation  finite element  software.    It appears  that  the  shale  layer having much weaker  strength properties  than  the bedding planes between the limestones causes the roof between the excavation and the shale layer to act  as  one  thick  beam  instead  of  multiple  thinner  beams.  This  often  results  in  separation  of  the shale/limestone  interface  before  separation  of  the  limestone  bedding  planes  and  subsequent  stress redistribution causing diagonal fracturing from the haunches to the shale layer causing failure. Finally, in most cases, failure that extended to the shale layer truncated at the shale layer.  It is possible that failure would progress beyond the shale when the shale layer is less than 2m from the roof because a compression arch providing stabilization in the new roof beam is less likely to be able to form over the remaining span.  Although  not  heavily  discussed  in  this  thesis, many  of  the  research  conclusions  are  dependent  on  a significant  strength  differential  between  the  limestone  and  the  shale  (weak  layer).    The  strength differential  is the primary contributor to the  large  increase  in depth of failure as well as the change  in preferential failure mechanism to diagonal fracturing.  The strength differential allows for unusually high stress concentrations to develop  in the roof rock as well as  it does not allow for stress transfer across the weak boundary.   103  10.2 Key	Contributions	The key contribution  from  this  research  is  the Shale Layer  Influence Chart  (Figure 48) which provides guidance  for  ground  support  and mine  design.    Also,  it  is  the  author’s  hope  that  this  dissertation presents  justifiable  reason  to change  the philosophy  for ground  support design  in  these mines under high horizontal conditions in the following ways: 1) There  needs  to  be  a  greater  importance  put  towards  identifying  the  frequency  of  shale interbeds where the shale is significantly weaker relative to the limestone.   2) The range in bedding thickness within the limestone needs to be understood, as opposed to just using the rock mass classification used in empirical ground support charts.   3) Due  to  the  depth  of  failure,  cable  bolts  should  be  considered  for  intersections  where  the diagonal spans are greater as opposed to rock bolts because standard 2.4m or 3.0m rock bolts do not provide enough anchorage depth  to  support  the weight of  the  slabs generated under these conditions.   Ultimately,  the  best  approach  to  ground  control would  be  to  isolate  the  required  depth  of massive limestone  such  that  a  suitable  compression  arch  can  form  and  ground  support  requirements  can be minimized.    The  Shale  Layer  Influence  Chart  (Figure  48)  developed  from  the  research  can  provide guidance for the thickness of competent beam that should be isolated to achieve this goal. 10.3 Possible	Future	Research	The conclusions of  this  research could be  further explored  to develop a greater understanding of  the influence  of  a  weak  layer  within  high  horizontal  stress  environments  from  a  three‐dimensional perspective  as  computing  power  increases,  and  to  establish  guidelines  for  ground  support  for underground  aggregate  mines  with  high  in‐situ  stresses  and  interbedded  weak  layers.    Finally  the strength  differential  conditions  between  the  shale  interbed  and  limestone  in  order  for  weak  layer influences  to  occur  can  be  investigated.    The  following  sections  discuss  these  potential  research directions.   10.3.1 Snaking	Failures	In  some  aggregate mining  operations  in  the US,  snaking  failures  connecting multiple  roof  falls  have formed suggesting that the failures simulated as part of this research can extend three dimensionally to distances of hundreds of meters.  Figure 53 shows an example of a snaking failure which occurred as a series of  failure over  a 70 day period of  time.   The  interesting part  about  these  failures  is  that  they progress opposite of the mining direction until they connect with a prior roof fall. 104   Figure 52.  Snaking Failure Observed in US Limestone Mine (Esterhuizen et al., 2008) At this point, snaking failures have been attempted to be explained by simple three‐dimensional stress redistribution models (Esterhuizen et al., 2008).  However, since observations published in the literature and supported by  the author concluded  that many  roof  falls  truncate at a weak  layer such as a shale interbed, it may be possible that the lateral continuity of the weak layer allows for the migration of the roof failures over time as stresses do redistribute after a roof fall or a mining heading  is extended.   As such, there could be some merit  in using three‐dimensional discontinuous modelling software to back‐analyze a snaking failure. 10.3.2 Ground	Support	Design	Guidelines	From this research, the understanding of the increase in depth of failure as a result of the influence of a shale  layer can be used  to establish guidelines  for ground support  in mines subject  to high horizontal stress conditions.  Even though the results of this research are mostly conceptual, the trends relating to how the depth of failure could be used to provide guidance on the length and type of bolts that should be used  in a ground support package  for the mine at  intersections, or the thickness of  limestone that should  be  isolated  between  the  roof  and  a  known weak  layer  so  that  a  less  robust  ground  support 105  package is required.  It should be noted that it is not just the length of bolt being optimized by isolating a competent beam of rock in the roof but also the robustness and frequency of ground support because with greater possible depths of failure comes the risk of larger volume slabs or wedges that need to be supported.    10.3.3 Strength	Differential	Requirements	The definition of a weak  layer  is not trivial.   The relative difference  in rock mass strength between the host rock and weak  layer controls whether there  is an  increase  in the depth of failure.    If the strength differential  is not beyond  a  critical  threshold,  then  separation on  the weak  layer prior  to  intact  rock failure is unlikely to occur.  Consequently if separation of the bedding does not occur first, then failure will be  in  the  form of a  series of  single beam buckling  failures as opposed  to a multi‐beam diagonal failure.    Therefore,  a  future  research  topic  could  be  to  investigate  the  critical  strength  differential required between  the shale and  limestone such  that  the  influence of  the weak  layer on  the depth of failure  is  important.   From an engineering perspective,  the merit  in  investigating  the  critical  strength differential  requirement  is  to understand under which  rock mass conditions  the  interbeds need  to be considered independently in the mine and ground support design. An associated  research  topic may be  to  investigate  the minimum  thickness of  shale  required  for  it  to perform  as  an  independent weak  layer.    Intuitively,  the  strength differential  requirement  is  likely  an integral part of determining the minimum shale  thickness requirements; a thin, very weak shale  likely has the same influence as a thick moderately weak shale. 106  Bibliography	Acres  Consulting  Services  Ltd.  (1976).  Underground  Aggregate  Mining  Preliminary  Feasibility  Study. Niagara Falls. Agapito,  J., & Gilbride, L.  (2002). Horizontal Stresses as  Indicators of Roof Stability. 2002 SME Annual Meeting. Phoenix. ASTM  International.  (2014).  ASTM  D7012‐14,  Standard  Test  Methods  for  Compressive  Strength  and Elastic Moduli of Intact Rock Core Specimens under Varying States of Stress and Temperatures. West Conshohocken: ASTM International. Barton, N., Loset, F., Lien, R., & Lunde, J. (1980). Application of the Q‐System in Design Decisions. In M. Bergman (Ed.), Subsurface Space (Vol. 2, pp. 553‐561). New York: Pergamon. Beer, G., & Meek, J. (1982). Design Curves for Roofs and Hangingwalls  in Bedded Rock. Transactions of the Institution of Mining and Metallurgy based on Voussoir Beam and Plate Solutions, A18‐A22. Benardos,  A.,  Kaliampakos,  D.,  Prousiotis,  J., Mavrikos,  A.,  &  Skoparantzos,  K.  (2001).  Underground Aggregate Mining  in Athens: A Promising  Investment Plan. Tunnelling and Underground Space Technology(16), 323‐329. Bieniawski, Z. (1976). Rock Mass Classification in Rock Engineering. In Z. Bieniawski, & A. Balkema (Ed.), Symposium on Exploration for Rock Engineering, (pp. 97‐106). Rotterdam. Bieniawski, Z. (1989). Engineering Rock Mass Classifications. New York: John Wiley & Sons. Birchard, M., Rutka, M., & Brunton, F.  (2004). Lithofacies and Geochemistry of  the Lucas Formation  in the  Subsurface  of  Southwestern  Ontario:  A  High‐Purity  Limestone  and  Potential  High‐Purity Dolostone Resource. Ontario Geological Survey, Open File Report 6137, 180p. Brady,  B., &  Brown,  E.  (1985).  Rock Mechanics  for Underground Mining  (3rd  ed.).  The Netherlands: Springer. Brook, D. (1991). Abandoned Limestone Mines in the West Midlands of England ‐ A Strategy for Action. Fourth International Symposium on Land Subsidence. No. 200, pp. 215 ‐ 223. IAHS. 107  Brown, T. J., Coggan, J. S., Evans, D. J., Foster, P. J., Hewitt, J., Kruyswijk, J. B., . . . Steadman, E. J. (2010). Underground Mining  of  Aggregates Main  Report  ‐  Assess  the  Feasibility  of  the Underground Mining of Aggregates. ASRP Project No. 7. Carmichael, T., & Smith, G. (1978). Nuclear Power Plants, Underground Space and Engineering Geology. In A. Currie, & W. Mackasey (Ed.), Toronto  '78: Field Trips (pp. 106‐113). Toronto: The Carswell Printing Company. Carr, D. D., Rooney,  L.  F., &  Freas, R.  C.  (1994).  Limestone  and Dolomite.  In D.  Carr  (Ed.),  Industrial Minerals and Rocks (6th ed., pp. 605‐629). Society for Mining, Metallurgy, and Exploration. Cruden,  A.  (2011).  OPG's  Deep  Geologic  Repository  for  Low  &  Intermediate  Level  Waste:  Outcrop Fracture Mapping. NWMO DGR‐TR‐2011‐43. D'aliesio, R., & Howlett, K. (2012, November 21). Critics Celebrate Surprise End of Mega Quarry north of Toronto.  The  Globe  and  Mail.  Retrieved  from http://www.theglobeandmail.com/news/national/critics‐celebrate‐surprise‐end‐of‐mega‐quarry‐north‐of‐toronto/article5520026/ Derry, Michener, Booth, Wahl, & Ontario Geological  Survey.  (1989).  Limestone  Industries  of Ontario, Volume  1  ‐ Geology,  Properties  and  Economics. Ontario Ministry  of Natural  Resources,  Land Management Branch, 158p. Diederichs, M.  (1999).  Instability  of  Hard  Rockmasses:  The  Role  of  Tensile  Damage  and  Relaxation. University of Waterloo, Waterloo. Diederichs, M., & Kaiser, P.  (1999).  Stability of  large  excavations  in  laminated hard  rock masses:  the voussoir analogue revisited. Internation Journal of Rock Mechanics and Mining Sciences, 36, 97‐117. Dolinar, D. R.  (2003). Variation of Horizontal Stresses and Strains  in Mines  in Bedded Deposits  in  the Eastern  and Midwestern United  States.  Proceedings  of  the  22nd  International  Conference  on Ground Control in Mining. Morgantown, VW. Eberhardt, E.  (2012). The Hoek‐Brown Failure Criterion. Rock Mechanics and Rock Engineering, 45(6), 981‐988. 108  Economopoulos, J., Koronakis, N., & Sofianos, A. (1993). Roof Failure Mechanisms in Greek Underground Bauxite  Mines.  In  L.  Ribeiro,  E.  Sousa,  &  N.  Grossmann  (Ed.),  Eurock  '93:  Safety  and Environmental Issues in Rock Engineering. 1, pp. 795‐802. Lisbon: A.A. Balkema. Esterhuizen,  G.  S.,  Dolinar,  D.  R.,  &  Ellenberger,  J.  L.  (2007).  Observations  and  evaluation  of  Floor Benching Effects on Pillar Stability  in U.S. Limestone Mines. Proceedings of  the 1st Canada‐US Rock Mechanics Symposium. Vancouver, BC. Esterhuizen, G. S., Dolinar, D. R., & Iannacchione, A. T. (2008). Field Observations and Numerical Studies of  Horizontal  Stress  Effects  on  Roof  Stability  in  U.S.  Limestone Mines.  Journal  of  the  South African Institute of Mining and Metallurgy, 108(6), 345 ‐ 352. Esterhuizen, G.  S., Dolinar, D.  R.,  Ellenberger,  J.  L.,  Prosser,  L.  J., &  Iannacchione, A.  T.  (2007).  Roof Stability  Issues  in Underground Limestone Mines  in the United States. Proceedings of the 26th International Conference on Ground Control in Mining, (pp. 320 ‐ 327). Morgantown, WV. Evans, W.  (1941).  The  Strength  of Undermined  Strata.  Transactions  of  the  Institution  of Mining  and Metallurgy, 50, 475‐500. Freas, R., Hayden,  J., & Pryor  Jr., C.  (2006). Limestone and Dolomite.  In  Industrial Minerals and Rocks: Commodities, Markets and Uses  (7 ed., pp. 581‐597). Littleton: Society for Mining, Metallurgy, and Exploration. Gartner Lee Limited. (1996). Regional Geologic Model: Smithville Phase IV. Smithville Phase IV Bedrock Remediation Program. Goodman, R. (1989). Introduction to Rock Mechanics (2nd ed.). John Wiley & Sons. Goudie, A.  (Ed.).  (2006).  Encyclopedia of Geomorphology  (Vol. 1 & 2). New York: Taylor &  Francis  e‐Library. Government  of  Canada.  (2015,  02  11).  Population  of  Census  Metropolitan  Areas.  Retrieved  from Statistics  Canada:  http://www.statcan.gc.ca/tables‐tableaux/sum‐som/l01/cst01/demo05a‐eng.htm Grass, J., & Dupak, D. (1986). Underground Aggregate Mining and Potential Uses of Underground Space ‐ Report 86058. Toronto: Ontario Hydro, Geotechnical and Hydraulic Engineering Dept. 109  Grimstad, E., & Barton, N. (1993). Updating the Q‐System for NMT. International Symposium on Sprayed Concrete:  Modern  Use  of  Wet  Mix  Sprayed  Concrete  for  Underground  Support  (pp.  46‐66). Fagernes: Norwegian Concrete Association. Gross, M., & Engelder, T. (1991). A Case for Neotectonic joints along the Niagara Escarpment. Tectonics, 10, 631‐641. Haimson,  B.  C.,  &  Lee,  C.  F.  (1980).  Hydrofracturing  Stress  Determinations  at  Darlington,  Ontario. Underground  Rock  Engineering,  13th  Canadian  Symposium  on  Rock  Mechanics,  Canadian Institute of Mining and Metallurgy Special Vol. 22, pp. 42‐50. Montreal. Hamrin, H. (2001). Chapter 1: Underground Mining Methods and Applications.  In W. A. Hustrulid, & R. Bullock (Eds.), Underground Mining Methods: Engineering Fundamentals and International Case Studies (pp. 3 ‐ 14). Littleton, Colorado: Society for Mining, Metallurgy and Exploration. Heidbach, O., Tingay, M., Barth, A., Reinecker,  J., Kurfeβ, D., & Müller, B.  (2010). Global Crustal Stress Pattern  based  on  the World  Stress Map  Database  release  2008.  Tectonophysics,  482,  3‐15. doi:10.1016/j.tecto.2009.07.023 Hoek,  E.,  Carranza‐Torres,  C.,  &  Corkum,  B.  (2002).  Hoek‐Brown  failure  criterion  ‐  2002  Edition. Proceedings from NARMS‐TAC Conference, 1, pp. 267‐273. Toronto. Holst, T. (1982). Regional Jointing in the Michigan Basin. Geology, 10, 971‐974. Hutchinson,  D.,  &  Diederichs,  M.  (1996).  Cablebolting  in  Underground  Mines.  Richmond:  BiTech Publishers Ltd. Iannacchione, A. T., Marshall, T. E., Burke, L., Melville, R., & Litsenberger, J. (2001). Safer Mine Layouts for Underground Stone Mines subjected  to Excessive Levels of Horizontal Stress  . SME Annual Meeting. Denver, CO. Itasca Consulting Group.  (2008). OPG's Deep Geologic Repository  for Low &  Intermediate Level Waste: Phase I Long‐Term Cavern Stability.  Itasca Consulting Group, Inc. (2015). UDEC Version 6.0: Distinct‐Element Modeling of Jointed and Blocky Material in 2D. 110  Johnson, M., Armstrong, D., Sanford, B., Telford, P., & Rutka, M. (1992). Paleozoic and Mesozoic Geology of Ontario. In P. Thurston, H. Williams, R. Sutcliffe, & G. Stott (Eds.), Geology of Ontario, Ontario Geological Survey, Special Volume 4, Part 2 (pp. 907‐1008). Jumikis, A. R. (1979). Rock Mechanics. Rockport, MA: Trans Tech Publications. Kendorski, F. S. (2002). Understanding and Solving Roof Control Problems  in Stone Mines. SME Annual Meeting. Phoenix, AZ. Koyama, T., & Lanru, J. (2007). Effects of Model Scale and Particle Size on Micro‐Mechanical Properties and  Failure  Processes  of  Rocks  ‐  A  Particle Mechanics  Approach.  Engineering  Analysis  with Boundary Elements, 31(5), 458 ‐ 472. Lam,  T.,  Martin,  D.,  &  McCreath,  D.  (2007).  Characterising  the  Geomechanics  Properties  of  the Sedimentary Rocks for the DGR Excavations. Canadian Geotechnical Conference, (pp. 636 ‐ 644). Ottawa. Laubscher, D. (1990). A Geomechanics Classification System for the Rating of Rock Mass in Mine Design. Journal of the South African Institute of Mining and Metallurgy, 90(86), 257‐273. Lee, C.  (1981).  In  Situ  Stress Measurements  in  Southern Ontario. Proc. 22nd US  Symposium on Rock Mechanics (pp. 465 ‐ 472). Massachusetts: Rock Mechanics from Research to Application. Lee,  C.,  &  White,  O.  (1993).  Engineering  Geology  in  Underground  Aggregate  Mining  and  Space Utilisation. Engineering Geology, 35(3‐4), 247‐257. Li, C., & Nordlund, E. (1993). A model for simulation of the failure of brittle rocks. In L. Ribeiro, E. Sousa, & N. Grossman  (Ed.), Eurock  '93: Safety and Environmental  Issues  in Rock Engineering. 1, pp. 143‐150. Lisbon: A.A. Balkema. Lo, K. (1978). Regional Distribution of In Situ Horizontal Stresses in Rocks of Southern Ontario. Canadian Geotechnical Journal(15), 371 ‐ 381. Lo, K.  (1989). Recent advances  in design and evaluation of performance of underground structures  in rocks. Tunnelling and Underground Space Technology, 4(2), 171‐183. 111  Marinos, P., Marinos, V., & Hoek, E. (2007). The Geological Strength Index (GSI): A Characterization Tool for Assessing Engineering Properties  for Rock Masses.  In M. Romana, A. Perucho, & C. Olalla (Ed.), ISRM Workshop W1 (pp. 13‐21). Madrid: Taylor and Francis Group. Martin,  C.  (1997).  The  Effect  of  Cohesion  Loss  and  Stress  Path  on  Brittle  Rock  Strength.  Canadian Geotechnical Journal, 34(5), 698‐725. Martin, C., Kaiser, P., & McCreath, D. (1999). Hoek‐Brown parameters for predicting the depth of brittle failure around tunnels. Canadian Geotechnical Journal, 36(1), 136‐151. Ministry  of  Finance.  (2011). Ontario  Population  Projections  Update:  2010  ‐  2036 Ontario  and  Its  49 Census Divisions. Government of Ontario. Mottahed, P., & Szeki, A. (1982). The Collapse of Room and Pillar Workings in a Shaley Gypsum Mine due to Dynamic Loading.  In  I. Farmer (Ed.), Symposium on Strata Mechanics. 32, pp. 260‐263. New York: Elsevier Scientific Publishing Company. Mukherjee,  S.  (2003).  Boundary  Element  Methods  in  Solid  Mechanics  ‐  a  Tribute  to  Frank  Rizzo. Electronic Journal of Boundary Elements, 1(1), 47‐55. National  Farmers  Union.  (2015,  November  10).  Impact  of  Water  Taking  Operations  in  Wellington County.  North  Dufferin  Agricultural  and  Community  Taskforce.  Retrieved  from http://www.ndact.com/ Nuclear  Waste  Management  Organization  and  AECOM  Canada  Ltd.  (2011).  OPG's  Deep  Geologic Repository  for  Low &  Intermediate  Level Waste:  Regional Geomechanics  ‐  Southern Ontario. NWMO DGR‐TR‐2011‐13. Park, E. M. (2004). Numerical Simulation of the Mechanical Behavior of Discontinuous Rock Masses. In R. H. Y. Shimizu  (Ed.), Numerical Modeling  in Micromechanics via Particle Methods  (pp. 85  ‐ 91). London: Taylor and Francis Group. Planning Initiatives Ltd. and Associates. (1992). Aggregate Resources of Southern Ontario. A State of the Resource Study. Toronto: Ontario Ministry of Natural Resources. Proctor & Redfern  Ltd.  (1974). Mineral Aggregates  Study  ‐  Central Ontario  Planning  Region. Ontario Ministry of Natural Resources. 112  Ran, J., Passaris, E., & Mottahed, P. (1994). Shear Sliding Failure of the Jointed Roof  in Laminated Rock Mass. Rock Mechanics and Rock Engineering, 27(4), 235‐251. Scholtes, L., Donze, F.‐V., & Khanal, M.  (2011). Scale Effects on Strength of Geomaterials, Case Study: Coal. Journal of the Mechanics and Physics of Solids, 59(5), 1131 ‐ 1146. Shinobe, A.  (1997). Economics of Underground Conversion  in an Operating Limestone Mine. Montreal: University of McGill. Stacey, T., & Page, C. (1986). Practical Handbook for Underground Rock Mechanics. Federal Republic of Germany: Trans Tech Publications. Sterling, R. (1980). The Ultimate Load Behaviour of Laterally Constrained Rock Beams. Proceedings of the 21st US Symposium on Rock Mechanics, 533‐542. Stimpson, B., & Ahmed, M. (1992). Failure of a  linear Voussoir arch: a  laboratory and numerical study. Canadian Geotechnical Journal, 29(2), 188‐194. Stocks,  J.  (1979). The Production of Aggregates by Underground Mining Methods  ‐ Development of a Simulation Technique for Prefeasibility Analysis. Minerals and the Environment, 1(4), 148 ‐ 160. Tang, C., & Hudson, J. A. (2010). Rock Failure Mechanisms: Explained and Illustrated. London: CRC Press. Thompson, L. G. (1982). Underground Aggregate Mining. Ministry of Natural Resources. US Geological Survey. (1999). Natural Agreggates ‐ Foundation of America's Future.  US Geological Survey. (2012). Limestone ‐ A Crucial and Versatile Industrial Mineral Commodity.  VanDyken,  R.  (2011, December). Melancthon Mega‐Quarry: Water  Concerns  and  Environmental  Risk Assessments  impede Plans  for a Mega‐Quarry near  the GTA. Canadian Geographic. Retrieved from http://www.canadiangeographic.ca/magazine/dec11/ontario_mega_quarry.asp Wang,  H.,  Xu,  C.,  Li,  T.,  Liu,  S.,  &  Peng,  N.  (2013).  The  Relationships  between Macro‐  and Micro‐Mechanical  Parameters  for  Modeling  Diamond  Rock  Cutting  using  DEM.  3rd  International FLAC/DEM Symposium. Hangzhou, China. Wolfe, R.  (1993). An  inventory of  inactive quarries  in  the Paleozoic  limestone and dolostone  strata of Ontario (Vol. Open File Report 5863). (O. G. Survey, Ed.) Queen's Printer of Ontario. 113  Appendices	These appendices provide  summary  tables  for each of  the numerical models  for both  the  calibration exercise and the shale layer influence modelling.  The generic forms of the UDEC codes used to carry out the modelling are provided as well.   Appendix	A						Summary	of	Roof	Failure	Numerical	Model	Results	This appendix provides a summary of the 132 numerical simulations of the varying combinations of  in‐situ stress, bedding thickness and distance the shale  layer  is from the roof of the excavation to assess the  influence  that  a  weak  layer  has  on  the  depth  of  failure  under  horizontally  dominant  stress conditions. The conditions modelled, depth of failure and main mode of failure are provided in Table 17. Table 17.  Summary of Shale Layer Modelling Results Test ID  Distance of Shale from Roof Bedding Thickness (m) Depth of Failure(m above roof)  Primary Mode of Failure Stress Ratio = 0.5H:1V Model_116  No shale  0.1  0.2  Slabbing Model_115  No shale  0.2  0.2  Slabbing Model_114  No shale  0.5  0  No failure Model_13  No shale  1  0  No failure Stress Ratio = 1H:1V Model_104  No shale  0.1  0.4  Slabbing Model_103  No shale  0.2  0.2  Slabbing Model_102  No shale  0.5  0  No failure Model_12  No shale  1  0  No failure Model_120  2  0.1  2  Diagonal fracturing Model_119  2  0.2  2  Diagonal fracturing Model_118  2  0.5  0  No failure Model_117  2  1  0  No failure Model_42  3  0.1  0.4  Buckling Model_30  3  0.2  0.2  Buckling Model_18  3  0.5  0  No failure Model_10  3  1  0  No failure Model_86  4  0.1  0.4  Buckling Model_85  4  0.2  0.2  Buckling Model_84  4  0.5  0  No failure Model_83  4  1  0  No failure Model_46  5  0.1  0.4  Buckling Model_34  5  0.2  0.2  Buckling Model_22  5  0.5  0  No failure 114  Test ID  Distance of Shale from Roof Bedding Thickness (m) Depth of Failure(m above roof)  Primary Mode of Failure Model_11  5  1  0  No failure Model_82  6  0.1  0.4  Buckling Model_81  6  0.2  0.2  Buckling Model_80  6  0.5  0  No failure Model_79  6  1  0  No failure Model_50  7  0.1  0.4  Buckling Model_38  7  0.2  0.2  Buckling Model_26  7  0.5  0  No failure Model_14  7  1  0  No failure Model_101  8  0.1  0.4  Buckling Model_100  8  0.2  0.2  Buckling Model_99  8  0.5  0  No failure Model_98  8  1  0  No failure Stress Ratio = 2H:1V Model_107  No shale  0.1  0.6  Buckling Model_106  No shale  0.2  0.6  Buckling Model_105  No shale  0.5  0.5  Buckling Model_3  No shale  1  0  No failure Model_124  2  0.1  2  Diagonal fracturing Model_123  2  0.2  2  Diagonal fracturing Model_122  2  0.5  2  Diagonal fracturing Model_121  2  1  2  Diagonal fracturing Model_43  3  0.1  3  Diagonal fracturing Model_31  3  0.2  3  Diagonal fracturing Model_19  3  0.5  0.5  Buckling Model_8  3  1  0  No failure Model_62  4  0.1  4  Diagonal fracturing Model_61  4  0.2  0.6  Buckling Model_78  4  0.5  0.5  Buckling Model_77  4  1  0  No failure Model_47  5  0.1  0.7  Buckling Model_35  5  0.2  0.6  Buckling Model_23  5  0.5  0  No failure Model_9  5  1  0  No failure Model_76  6  0.1  0.7  Buckling Model_75  6  0.2  0.6  Buckling Model_74  6  0.5  0  No failure Model_73  6  1  0  No failure Model_51  7  0.1  0.7  Buckling 115  Test ID  Distance of Shale from Roof Bedding Thickness (m) Depth of Failure(m above roof)  Primary Mode of Failure Model_39  7  0.2  0.6  Buckling Model_27  7  0.5  0  No failure Model_15  7  1  0  No failure Model_97  8  0.1  0.7  Buckling Model_96  8  0.2  0.6  Buckling Model_95  8  0.5  0  No failure Model_94  8  1  0  No failure Stress Ratio = 3H:1V Model_110  No shale  0.1  1.4  Buckling Model_109  No shale  0.2  1.3  Buckling Model_108  No shale  0.5  1  Buckling Model_2  No shale  1  1  Buckling Model_128  2  0.1  2  Diagonal fracturing Model_127  2  0.2  2  Diagonal fracturing Model_126  2  0.5  2  Diagonal fracturing Model_125  2  1  2  Diagonal fracturing Model_44  3  0.1  3  Diagonal fracturing Model_32  3  0.2  3  Diagonal fracturing Model_20  3  0.5  3  Diagonal fracturing Model_4  3  1  3  Diagonal fracturing Model_72  4  0.1  4  Diagonal fracturing Model_71  4  0.2  4  Diagonal fracturing Model_70  4  0.5  4  Diagonal fracturing Model_60  4  1  1  Buckling Model_48  5  0.1  5  Diagonal fracturing Model_36  5  0.2  5  Diagonal fracturing Model_24  5  0.5  2.5  Buckling Model_5  5  1  1  Buckling Model_59  6  0.1  6  Diagonal fracturing Model_58  6  0.2  2.8  Buckling Model_69  6  0.5  2.5  Buckling Model_68  6  1  1  Buckling Model_52  7  0.1  3  Buckling Model_40  7  0.2  2.8  Buckling Model_28  7  0.5  2.5  Buckling Model_16  7  1  1  Buckling Model_93  8  0.1  3  Buckling Model_92  8  0.2  2.8  Buckling Model_91  8  0.5  2.5  Buckling 116  Test ID  Distance of Shale from Roof Bedding Thickness (m) Depth of Failure(m above roof)  Primary Mode of Failure Model_90  8  1  1  Buckling Stress Ratio = 4H:1V Model_113  No shale  0.1  1.9  Buckling Model_112  No shale  0.2  1.7  Buckling Model_111  No shale  0.5  1.5  Buckling Model_1  No shale  1  1.4  Buckling Model_132  2  0.1  2  Diagonal fracturing Model_131  2  0.2  2  Diagonal fracturing Model_130  2  0.5  2  Diagonal fracturing Model_129  2  1  2  Diagonal fracturing Model_45  3  0.1  3  Diagonal fracturing Model_33  3  0.2  3  Diagonal fracturing Model_21  3  0.5  3  Diagonal fracturing Model_7  3  1  3  Diagonal fracturing Model_67  4  0.1  4  Diagonal fracturing Model_66  4  0.2  4  Diagonal fracturing Model_65  4  0.5  4  Diagonal fracturing Model_54  4  1  4  Diagonal fracturing Model_49  5  0.1  5  Diagonal fracturing Model_37  5  0.2  5  Diagonal fracturing Model_25  5  0.5  5  Diagonal fracturing Model_6  5  1  1.4  Buckling Model_64  6  0.1  6  Diagonal fracturing Model_56  6  0.2  6  Diagonal fracturing Model_55  6  0.5  2.5  Buckling Model_63  6  1  1.4  Buckling Model_53  7  0.1  7  Diagonal fracturing Model_41  7  0.2  2.8  Buckling Model_29  7  0.5  2.5  Buckling Model_17  7  1  1.4  Buckling Model_57  8  0.1  3  Buckling Model_89  8  0.2  2.8  Buckling Model_88  8  0.5  2.5  Buckling Model_87  8  1  1.4  Buckling 117  Appendix	B						Summary	of	UCS	Test	Numerical	Modelling	Results	This appendix provides a summary of the 77 UCS test models conducted to establish suitable micro‐mechanical properties for the roof failure numerical study.  The input parameters, peak strength, number of cycles to failure and tensile and shear failure initiation values are provided in Table 18. Table 18.  Summary of UCT Test Numerical Modelling Results Test ID Input Parameters Max. Axial Stress/Peak Strength (MPa) # of Cycles to Failure # of Joint Contacts Tensile Failure Initiation (% of strength) Shear Failure Initiation (% of strength) Point when Shear Contacts exceed Tensile Contacts (% of strength) Voronoi Edge Length Joint Cohesion, c (MPa) Joint Tensile Strength, T (MPa) Joint Friction,  (°) Residual Joint Friction,  R (°) Bulk Modulus, K (GPa) Shear Modulus, G (GPa) 26  0.001  50  9  50  31  30.56  22.92  out of memory  ‐  ‐  ‐  ‐  ‐ 17  0.0025  45  15  50  31  30.56  22.92  47.5  2300000  5390  38  39  48 25  0.0025  50  9  50  31  30.56  22.92  42.7  2040000  5346  39  41  50 28  0.0025  80  9  60  31  30.56  22.92  45.1  2220000  5400  30  32  40 52  0.0025  110  9  50  31  30.56  22.92  46.2  2265000  5605  35  38  45 53  0.0025  50  25  50  31  30.56  22.92  49.4  2270000  5640  44  40  0 54  0.0025  50  9  70  31  30.56  22.92  41.2  2080000  5585  40  42  52 61  0.0025  110  25  50  31  30.56  22.92  58.5  3150000  5638  38  34  0 16  0.005  45  15  50  31  30.56  22.92  67.8  1184000  1619  34  34  45 24  0.005  50  9  50  31  30.56  22.92  64.8  1330000  1616  29  32  49 44  0.005  110  9  50  31  30.56  22.92  82  1760000  1777  23  25  40 47  0.005  50  25  50  31  30.56  22.92  72.1  1250000  1752  37  40  46 51  0.005  50  9  70  31  30.56  22.92  74.4  1385000  1736  25  28  46 118  Test ID Input Parameters Max. Axial Stress/Peak Strength (MPa) # of Cycles to Failure # of Joint Contacts Tensile Failure Initiation (% of strength) Shear Failure Initiation (% of strength) Point when Shear Contacts exceed Tensile Contacts (% of strength) Voronoi Edge Length Joint Cohesion, c (MPa) Joint Tensile Strength, T (MPa) Joint Friction,  (°) Residual Joint Friction,  R (°) Bulk Modulus, K (GPa) Shear Modulus, G (GPa) 60  0.005  110  25  50  31  30.56  22.92  100  1910000  1779  27  28  33 62  0.005  50  35  50  31  30.56  22.92  78.7  1360000  1754  40  39  0 63  0.005  50  45  50  31  30.56  22.92  70.4  1292000  1758  50  46  0 64  0.005  50  55  50  31  30.56  22.92  70.4  1292000  1758  50  46  0 65  0.005  50  25  30  31  30.56  22.92  76.6  1350000  1736  35  37  43 66  0.005  50  25  40  31  30.56  22.92  76.3  1400000  1757  35  37  43 67  0.005  50  25  60  31  30.56  22.92  81.1  1520000  1763  33  35  40 68  0.005  50  25  70  31  30.56  22.92  75.2  1540000  1756  31  35  42 69  0.005  80  25  30  31  30.56  22.92  86.7  1670000  1764  31  33  38 70  0.005  80  25  40  31  30.56  22.92  86.6  1482000  1756  31  33  38 71  0.005  80  25  50  31  30.56  22.92  88.4  1620000  1765  30  32  37 72  0.005  80  25  60  31  30.56  22.92  92.1  1750000  1767  29  31  36 73  0.005  80  25  70  31  30.56  22.92  94.6  1835000  1783  28  30  35 74  0.005  80  35  30  31  30.56  22.92  94.3  1797000  1776  34  33  0 75  0.005  80  35  40  31  30.56  22.92  94.3  1891000  1751  34  33  0 76  0.005  80  35  50  31  30.56  22.92  102.4  1997000  1786  31  30  0 77  0.005  80  35  60  31  30.56  22.92  96.5  1880000  1784  33  32  0 78  0.005  80  35  70  31  30.56  22.92  95  1925000  1777  30  31  34 79  0.005  110  35  30  31  30.56  22.92  102.1  1940000  1785  31  30.5  0 80  0.005  110  35  40  31  30.56  22.92  102.7  1917000  1744  31  30.5  0 119  Test ID Input Parameters Max. Axial Stress/Peak Strength (MPa) # of Cycles to Failure # of Joint Contacts Tensile Failure Initiation (% of strength) Shear Failure Initiation (% of strength) Point when Shear Contacts exceed Tensile Contacts (% of strength) Voronoi Edge Length Joint Cohesion, c (MPa) Joint Tensile Strength, T (MPa) Joint Friction,  (°) Residual Joint Friction,  R (°) Bulk Modulus, K (GPa) Shear Modulus, G (GPa) 81  0.005  110  35  50  31  30.56  22.92  101.9  1915000  1782  31  30.5  0 82  0.005  110  35  60  31  30.56  22.92  102.8  1934000  1782  31  30  0 83  0.005  110  35  70  31  30.56  22.92  106.8  2210000  1788  30  29  0 84  0.005  140  35  60  31  30.56  22.92  118.9  2366000  1795  26.5  26  0 85  0.005  170  35  60  31  30.56  22.92  121.5  2428000  1793  28.5  25.5  0 86  0.005  80  30  60  31  30.56  22.92  91.7  1834000  1773  33.5  32  0 87  0.005  110  30  60  31  30.56  22.92  100.7  2040000  1783  30.5  29  0 88  0.005  140  30  60  31  30.56  22.92  112.1  2250000  1788  28  27  0 89  0.005  170  30  60  31  30.56  22.92  113.5  2266000  1792  27  26  0 90  0.005  260  30  60  31  30.56  22.92  113.5  2266000  1793  27  26  0 91  0.005  80  30  40  31  30.56  22.92  94.2  1860000  1740  33  31  0 14  0.01  45  15  50  31  30.56  22.92  115.3  697000  502  26  26  32 18  0.01  45  9  50  31  30.56  22.92  100.2  575000  496  26  28  36 19  0.01  50  9  50  31  30.56  22.92  112.8  737000  491  23  25  32 27  0.01  80  9  60  31  30.56  22.92  150.6  1117000  509  17  19  24 29  0.01  50  25  50  31  30.56  22.92  98.8  573500  488  36  33  0 30  0.01  80  25  50  31  30.56  22.92  138  850000  503  25  24  0 31  0.01  110  25  50  31  30.56  22.92  169.5  925000  506  21  19  0 32  0.01  50  3  50  31  30.56  22.92  107.5  637000  489  23  25  31 33  0.01  50  25  60  31  30.56  22.92  129.7  820000  487  26  24  0 120  Test ID Input Parameters Max. Axial Stress/Peak Strength (MPa) # of Cycles to Failure # of Joint Contacts Tensile Failure Initiation (% of strength) Shear Failure Initiation (% of strength) Point when Shear Contacts exceed Tensile Contacts (% of strength) Voronoi Edge Length Joint Cohesion, c (MPa) Joint Tensile Strength, T (MPa) Joint Friction,  (°) Residual Joint Friction,  R (°) Bulk Modulus, K (GPa) Shear Modulus, G (GPa) 34  0.01  50  25  70  31  30.56  22.92  181.3  1244000  509  17  17  22 35  0.01  50  15  50  31  30.56  22.92  124.5  756000  499  24  24  30 36  0.01  80  9  50  31  30.56  22.92  130  713000  503  20  22  28 37  0.01  50  25  40  31  30.56  22.92  97.2  533000  488  36  33  0 38  0.01  50  25  30  31  30.56  22.92  95.4  483000  493  37  34  0 39  0.01  110  9  50  31  30.56  22.92  163  919000  506  16  17  22 40  0.01  20  9  50  31  30.56  22.92  67.7  619000  477  38  41  53 41  0.01  20  25  50  31  30.56  22.92  83.5  928000  482  42  42  0 50  0.01  50  9  70  31  30.56  22.92  185.6  1303000  508  14  15  20 55  0.01  50  9  30  31  30.56  22.92  94  408000  492  28  30  39 56  0.01  50  9  40  31  30.56  22.92  94.1  409000  491  28  30  39 57  0.01  50  9  60  31  30.56  22.92  116.8  723000  487  22  24  31 20  0.02  45  15  50  31  30.56  22.92  305.8  1012000  140  8.5  9.2  10 22  0.02  50  9  50  31  30.56  22.92  311.7  925000  137  6.6  6.6  0 43  0.02  110  9  50  31  30.56  22.92  343.2  975000  138  6  6  0 46  0.02  50  25  50  31  30.56  22.92  331.2  904000  139  9  9  9 49  0.02  50  9  70  31  30.56  22.92  373.8  973000  138  6  6  0 59  0.02  110  25  50  31  30.56  22.92  238.2  603000  141  12  12  0 21  0.04  45  15  50  31  30.56  22.92  445.4  440000  49  13  12  0 23  0.04  50  9  50  31  30.56  22.92  444.8  950000  47  11  8  0 121  Test ID Input Parameters Max. Axial Stress/Peak Strength (MPa) # of Cycles to Failure # of Joint Contacts Tensile Failure Initiation (% of strength) Shear Failure Initiation (% of strength) Point when Shear Contacts exceed Tensile Contacts (% of strength) Voronoi Edge Length Joint Cohesion, c (MPa) Joint Tensile Strength, T (MPa) Joint Friction,  (°) Residual Joint Friction,  R (°) Bulk Modulus, K (GPa) Shear Modulus, G (GPa) 42  0.04  110  9  50  31  30.56  22.92  444.6  850000  49  11  8  0 45  0.04  50  25  50  31  30.56  22.92  445.4  900000  49  14  12  0 48  0.04  50  9  70  31  30.56  22.92  461  800000  50  10  8  0 58  0.04  110  25  50  31  30.56  22.92  444.5  800000  49  14  12  0   122  Appendix	C						Computer	Code	Generated	Using	the	FISH	Language	Coding using the FISH language within UDEC was required to adequately observe the displacements and stresses within the model and understand the mechanisms  leading to failure.   As part of this research, UCS  test  simulations were  completed  to  ensure  that micro‐fracturing was  simulated  properly when using  the  Voronoi  tessellation,  and modelling  of  beam  failure  in  the  roof  of  classic  room  and  pillar limestone mines  in high horizontal stress environments was conducted   to understand the  impact of a weak layer on the depth of failure.  This appendix discusses the codes and FISH language used during the research to simplify and automate the collection of data so that the numerical modelling was effective from both a time and computational perspective. C.1	Roof	Failure	in	Limestone	Mine	Simulation	Code	In  an  effort  to  semi‐automate  the  generation  and  re‐generation  of  the  model  geometry,  material properties  (including  the  automatic  calculation  of  equivalent Mohr‐Coulomb  properties  from  Hoek‐Brown properties), voronoi regions and  joint properties and  locations of history points  for quantifying deflections  as  variations  in  the  bedding  thickness,  shale  distance  from  the  excavation  and  stress conditions were  studied, a code was written using UDEC’s FISH  language.   The  following datafile was called into UDEC once the input parameters were defined by the user under the section of code titled, “INPUT PARAMETERS” to carry out the basic numerical modelling for this research. ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ new ro 0.001 set ovtol 0.1 ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;INPUT PARAMETERS  define inputs h=12 w=22 inner_joint_spacing = 1 outer_joint_spacing = 2 beam_bedding_spacing = 1 beam_thickness = 3 shale_thickness = 0.2 impact_zone_above_shale = 8  ;ROCK MASS PROPERTIES  123  ;LIMESTONE L_intact_UCS = 110e6 L_GSI = 70 L_mi = 10 L_Disturbance_Factor = 0.7 L_Modulus = 36.04e9 L_res_fric = 31 L_res_coh = 0 L_res_t = 0 L_dil = 0 L_poisson = 0.19  L_dens = 2.68E3 vor_jkn = 1.60e11 vor_jks = 1.60e10  ;SHALE S_intact_UCS = 30e6 S_GSI = 50 S_mi = 6 S_Disturbance_Factor = 0 S_Modulus = 17.65e9 S_res_fric = 31 S_res_coh = 0 S_res_t = 0 S_dil = 0 S_poisson = 0.09 S_dens = 2.6e6  ;JOINT PROPERTIES sh_stiff = 1.60e10 nor_stiff = 1.60e11 bed_fric = 31 bed_coh = 0.57e6 bed_tens = 0.95e6 bed_rfric = 31 vor_jfric = 45 vor_jcoh = 80e6 vor_jten = 25e6  ;STRESSES hstress = ‐12e6 vstress = ‐4e6  end inputs ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;CALCULATIONS  124  define fish g1 = 7*h g2 = 7*w g3 = 3*w g4 = 3*h g5 = 4*w g6 = 4*h g7 = 3*w‐3*h g8 = 4*w+3*h g9 = 1*h g10 = 6*h j1 = 5*inner_joint_spacing j2 = 1*inner_joint_spacing vor1 = g3 ‐ 1.5*inner_joint_spacing  vor2 = g5 + 1.5*inner_joint_spacing vor3 = g6 ‐ 1*inner_joint_spacing vor4 = g3 ‐ 3*inner_joint_spacing vor5 = g5 + 3*inner_joint_spacing vor6 = g6 ‐ 3*inner_joint_spacing vor7 = g3 ‐ 5*inner_joint_spacing vor8 = g5 + 5*inner_joint_spacing end fish  ;SHALE LAYER define shale lower_shale_contact = g6+beam_thickness*inner_joint_spacing upper_shale_contact = lower_shale_contact+shale_thickness g19 = lower_shale_contact+impact_zone_above_shale*inner_joint_spacing end shale ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;GEOMETRY bl 0.0,0.0 0,g1 g2,g1 g2,0.0  ;ROOM DEFINITION cr g3,g4 g5,g4 cr g3,g4 g3,g6 join cr g3,g6 g5,g6 cr g5,g4 g5,g6 join  ;INNER BOUNDARY DEFINITION cr g7,g9 g7,g10 cr g8,g10 g8,g9  ;JOINT REGION ID jreg id 1 g7,g9 g7,g10 g8,g10 g8,g9 jreg id 2 vor7,g6 vor7,lower_shale_contact vor8,lower_shale_contact vor8,g6 125  jreg id 3 vor1,vor3 vor1,g19 vor2,g19 vor2,vor3 jreg id 4 vor7,lower_shale_contact vor7,g19 vor8,g19 vor8,lower_shale_contact jreg id 5 vor4,vor3 vor4,g19 vor1,g19 vor1,vor3 jreg id 6 vor4,vor6 vor4,vor3 g3,vor3 g3,vor6 jreg id 7 vor2,vor3 vor2,g19 vor5,g19 vor5,vor3 jreg id 8 g5,vor6 g5,vor3 vor5,vor3 vor5,vor6 jreg id 9 vor7,vor6 vor7,g19 vor4,g19 vor4,vor6 jreg id 10 vor5,vor6 vor5,g19 vor8,g19 vor8,vor6 jreg id 11 g7,lower_shale_contact g7,upper_shale_contact g8,upper_shale_contact g8,lower_shale_contact  ;BEDDING DEFINITION jset a 0 g 0 s outer_joint_spacing o 0,g6 jset a 0 g 0 s inner_joint_spacing o 0,g6 range jreg 1 jdelete  ;RANDOM SUBVERTICAL JOINTS jset a 90 5 g j1 j2 s 3 2 t 2.9 o g7,g9 range jreg 1  jdelete  cr vor1,vor3 vor1,g19  cr vor2,vor3 vor2,g19 cr vor4,vor6 vor4,g19  cr vor5,vor6 vor5,g19  cr vor7,vor6 vor7,g19  cr vor8,vor6 vor8,g19   ;VORONOI BEAM AND ABUTMENTS vo e 0.88 i 1 ro 0.001 seed 850925 range jreg 3 vo e 1.2 i 1 ro 0.001 range jreg 5 vo e 1.2 i 1 ro 0.001 range jreg 6 vo e 1.2 i 1 ro 0.001 range jreg 7 vo e 1.2 i 1 ro 0.001 range jreg 8 vo e 2 i 1 ro 0.001 range jreg 9 vo e 2 i 1 ro 0.001 range jreg 10 jdelete  ;SHALE LAYER cr g7,upper_shale_contact g8,upper_shale_contact  ;MONITORING JOINTS (ALL HAVE PROPERTIES OF INTACT ROCK) jset a 0 g 0 s beam_bedding_spacing o vor7,g6 range jreg 2 jset a 0 g 0 s beam_bedding_spacing o vor7,lower_shale_contact range jreg 4  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;MESHING gen quad 0.5 range jreg 3 gen edge 0.2 range jreg 3 126  gen quad 1 range jreg 5 gen edge .5 range jreg 5 gen quad 1 range jreg 6 gen edge .5 range jreg 6 gen quad 1 range jreg 7 gen edge .5 range jreg 7 gen quad 1 range jreg 8 gen edge .5 range jreg 8 gen quad 1.5 range jreg 9 gen edge .8 range jreg 9 gen quad 1.5 range jreg 10 gen edge .8 range jreg 10 gen quad 1.5 range jreg 1 gen edge .8 range jreg 1 gen quad 2 gen edge 1  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;MATERIAL PROPERTIES  ;ROCLAB PROPERTIES  define Roclab ;LIMESTONE L_mb = L_mi*exp((float(L_GSI)‐100)/(28‐14*L_Disturbance_Factor)) L_s = exp((float(L_GSI)‐100)/(9‐3*L_Disturbance_Factor)) L_a = 0.5+(exp(‐float(L_GSI)/15)‐exp(‐20/3))/6 L_sig_c = L_intact_UCS*L_s^L_a L_RL_t = L_s*L_intact_UCS/L_mb Temp1 = (L_mb+4*L_s‐L_a*(L_mb‐8*L_s))*(L_mb/4+L_s)^(L_a‐1) L_sig_cm = L_intact_UCS*(Temp1)/(2*(1+L_a)*(2+L_a)) L_sig3max = L_sig_cm*0.47*(L_sig_cm/(‐hstress))^(‐0.94) L_sig3n = L_sig3max/L_intact_UCS Temp2 = 6*L_a*L_mb*(L_s+L_mb*L_sig3n)^(L_a‐1) Temp3 = (Temp2)/(2*(1+L_a)*(2+L_a)+6*L_a*L_mb*(L_s+L_mb*L_sig3n)^(L_a‐1)) L_RL_phi = asin(Temp3)*180/pi Temp4 = (1+(6*L_a*L_mb*(L_s+L_mb*L_sig3n)^(L_a‐1)))/((1+L_a)*(2+L_a)) GAB = (1+2*L_a)*L_s+(1‐L_a)*L_mb*L_sig3n Temp5 = L_intact_UCS*(GAB)*(L_s+L_mb*L_sig3n)^(L_a‐1) L_RL_coh = (Temp5)/((1+L_a)*(2+L_a)*sqrt(Temp4)) L_Bulk_Mod = L_Modulus/(3*(1‐2*L_poisson)) L_Shear_Mod = L_Modulus/(2*(1+L_poisson))  ;SHALE S_mb = S_mi*exp((float(S_GSI)‐100)/(28‐14*S_Disturbance_Factor)) S_s = exp((float(S_GSI)‐100)/(9‐3*S_Disturbance_Factor)) S_a = 0.5+(exp(‐float(S_GSI)/15)‐exp(‐20/3))/6 S_sig_c = S_intact_UCS*S_s^S_a 127  S_RS_t = S_s*S_intact_UCS/S_mb Temp6 = (S_mb+4*S_s‐S_a*(S_mb‐8*S_s))*(S_mb/4+S_s)^(S_a‐1) S_sig_cm = S_intact_UCS*(Temp6)/(2*(1+S_a)*(2+S_a)) S_sig3max = S_sig_cm*0.47*(S_sig_cm/(‐hstress))^(‐0.94) S_sig3n = S_sig3max/S_intact_UCS Temp7 = 6*S_a*S_mb*(S_s+S_mb*S_sig3n)^(S_a‐1) Temp8 = (Temp7)/(2*(1+S_a)*(2+S_a)+6*S_a*S_mb*(S_s+S_mb*S_sig3n)^(S_a‐1)) S_RS_phi = asin(Temp8)*180/pi Temp9 = (1+(6*S_a*S_mb*(S_s+S_mb*S_sig3n)^(S_a‐1)))/((1+S_a)*(2+S_a)) SED = (1+2*S_a)*S_s+(1‐S_a)*S_mb*S_sig3n Temp10 = S_intact_UCS*(SED)*(S_s+S_mb*S_sig3n)^(S_a‐1) S_RS_coh = (Temp10)/((1+S_a)*(2+S_a)*sqrt(Temp9)) S_Bulk_Mod = S_Modulus/(3*(1‐2*S_poisson)) S_Shear_Mod = S_Modulus/(2*(1+S_poisson)) end Roclab  ;ASSIGNING ELASTIC MATERIAL PROPERTIES TO ROCK MASS group zone 'block' zone model elastic density L_dens bulk L_Bulk_Mod sh L_Shear_Mod range group 'block' group zone 'shale' range g7 g8 lower_shale_contact upper_shale_contact zone model elastic density S_dens bulk S_Bulk_Mod sh S_Shear_Mod range group 'shale'  ;ASSIGNING JOINT PROPERTIES group joint 'bedding' joint model  residual  jks  sh_stiff  jkn  nor_stiff  jfric  bed_fric  jcoh  bed_coh  jrf  bed_rfric  jresc  0  jdil  0  jt bed_tens jrt 0 range group 'bedding' set jcondf residual jks sh_stiff jkn nor_stiff jfric L_res_fric jcoh L_res_coh jdil 0 jt L_res_t jrf L_res_fric  ;ASSIGNING JOINT PROPERTIES TO VORONOI JOINTS EQUAL TO ROCK MASS PROPERTIES group joint 'voronoi above excavation mass' range jreg 3 angle 5,175  joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi above excavation mass' group joint 'voronoi abument mass 1' range jreg 5 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 1' group joint 'voronoi abument mass 2' range jreg 6 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 2' group joint 'voronoi abument mass 3' range jreg 7 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 3' group joint 'voronoi abument mass 4' range jreg 8 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 4' group joint 'voronoi abument mass 5' range jreg 9 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 5' 128  group joint 'voronoi abument mass 6' range jreg 10 angle 5,175 joint model  residual  jks  vor_jks  jkn  vor_jkn  jfric  vor_jfric  jcoh  vor_jcoh  jt  vor_jten  jrf  L_res_fric  jresc L_res_coh jrt L_res_t range group 'voronoi abument mass 6' group joint 'voronoi shale' range jreg 11 joint model  residual  jks vor_jks  jkn vor_jkn  jfric S_RS_phi  jcoh S_RS_coh  jten S_RS_t  jrf 0  jresc 0  jrt 0 range group 'voronoi shale'  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;BOUNDARY CONDITIONS  set gravity=0.0 ‐9.81 define xx g13 = g1‐0.1 g14 = g1+0.1 g15 = g2‐0.1 g16 = g2+0.1 end xx boundary stress hstress,0,0 range ‐0.1,0.1 0,g1 boundary stress hstress,0,0 range g15,g16 0,g1 boundary stress 0,0,vstress range 0,g2 g13,g14 boundary yvelocity 0 range 0,g2 ‐0.1,0.1 insitu stress hstress,0,vstress szz hstress  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;HISTORY POINTS  ;DEFINITION OF HISTORY POINT IN BEAM CENTRE define yy g17 = 3.5*w‐0.05 g_hist = g6+0.05 number_of_bedding = (lower_shale_contact‐g6)/beam_bedding_spacing loop n (1,number_of_bedding) g18 = g_hist + (n‐1)*beam_bedding_spacing command hist xdisp g17,g18 hist ydisp g17,g18 endcommand end_loop end yy  ;DEFINITION OF HISTORY POINT IN SHALE CENTRE define history_shale hist2 = lower_shale_contact+.01 end history_shale hist ydisp g17,hist2 129  hist xdisp g17,hist2  ;ARRAY OF HISTORY POINTS IN BEAM AND ABUTMENTS define hist_array horizontal_divisions = 10 vertical_divisions = 4 uppercontact = lower_shale_contact‐0.05 lowercontact = g6+0.05 loop m (‐1,horizontal_divisions) hist_pointx = g3+(m+1)*(w/horizontal_divisions) loop n (2,vertical_divisions) hist_pointy = g6+(n‐1)*(float(beam_thickness)/vertical_divisions) command hist ydisp hist_pointx,hist_pointy endcommand end_loop command hist ydisp hist_pointx,lowercontact hist ydisp hist_pointx,uppercontact endcommand end_loop end hist_array  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;MODEL EQUILIBRIUM  cycle 30000  save 2015_Model_ID_elastic.sav  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;PLASTIC PROPERTIES  zone model mohr density L_dens bulk L_Bulk_Mod sh L_Shear_Mod dil L_dil fric L_RL_phi coh L_RL_coh ten L_RL_t range group 'block' zone model elastic density L_dens bulk L_Bulk_Mod sh L_Shear_Mod range jreg 3 zone  model  mohr  density  S_dens  bulk  S_Bulk_Mod  sh  S_Shear_Mod  dil  S_dil  fric  S_RS_phi  coh S_RS_coh ten S_RS_t range group 'shale'  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;MODEL EQUILIBRIUM 2  cycle 40000  save 2015_Model_ID_plastic.sav  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ 130   ;CHANGE PROPERTIES HERE IF NEEDBE (SEE PROPERTY CHANGE DAT.FILE)  ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;EXCAVATION OF ROOM  reset disp jdisp delete range g3,g5 g4,g6 ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;STEPPING  cycle 100000 save 2015_Model_ID_100000.sav; This stepping code would be defined to repeat as necessary by the User ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ C.2	 UCS	Test	Simulation	Code	The purpose of conducting the UCS test simulations was to determine suitable Voronoi joint properties for  a  defined  Voronoi  edge  length  such  that  the  intact  rock mass  had  an  unconfined  compressive strength comparable to those results measured during actual laboratory testing while ensuring that the failure  through  the  simulated  sample  was  initially  tensile  dominated  and  failure  occurred  at approximately 30 – 40% of the unconfined compressive strength.  In order to validate the sample failure and  properly  understand  the  impacts  of  the  failure  mechanisms  on  the  sample  strength,  it  was imperative to produce an axial stress versus axial strain curve as the sample was loaded and ultimately failed.   This was completed by using FISH  language  in UDEC to calculate the average axial stress acting across  the  compressing  platen  –  sample  interface  and  the  corresponding  axial  strain  at  500  cycle intervals.    Since UDEC  does  not  have  an  explicit  time  function  in which  the  loading  velocity  can  be calculated  into  axial  strain,  the  axial  strain  corresponding  to  each  500  cycles  was  calculated  by measuring  the  vertical  displacement  of  the  compressing  platen  and  dividing  it  by  the  length  of  the sample.   So  that  the measurement process could be automated during  the  simulation, FISH  language was used to loop the loading procedure based on a 500 cycle loop and record the axial stress and axial strain values per loop in a table, which could be plotted following the completion of the simulation.  The FISH  language  used  to  code  the  automated  measurement  of  the  axial  stress/axial  strain  curve  is presented  in  the  datafile  code  below  under  “DEVELOPMENT  OF  THE  AXIAL  STRESS/AXIAL  STRAIN CURVE”. 131  A  similar  loop  sequence was  carried  out  to  observe  the  development  of  tensile  and  shear  fractures during the failure of the simulated sample.  In this loop, FISH language was used to measure the amount of  shear and normal displacement on all Voronoi edge  lengths  in  the model at 500 cycle  increments.  Displacements were considered  tensile  fractures  in  two ways: 1) normal displacement on  the contact was greater  than 0.1mm  regardless of  shear displacement along  the  same  joint  (as  long as  the  shear displacement  did  not  exceed  the  normal  displacement,  2)  normal  displacement  on  the  contact was greater  than 0.1mm and shear displacement was  less  than 0.1mm.    In doing so,  the  fractures  in pure tensile  failure  at  the  cycle  of  measurement  could  be  observed;  this  is  of  interest  because  proper simulation micro‐fracture  generation  in  tension  using  the  Voronoi  tessellation  was  the  goal  of  this calibration  exercise.    For  overall  tensile  fracture  count,  it was  assumed  that  as  long  as  the  normal displacement exceeded  the  shear displacement  then  the  fracture has a  tensile origin.   Contacts with shear displacements  greater  than  the normal displacements were  considered  to have  a  shear origin.  These  fracture  counts were  recorded  in  tables with  the  corresponding axial  stress  for each 500  cycle measurement  so  that  they could be plotted  following  the  simulation  to determine  the percentage of peak unconfined compressive strength in which fracturing initiated and the dominant fracture initiation mechanism.   The FISH  language used to code the automated counting of shear and tensile fractures at 500  cycle  increments  during  the UCS  test  simulation  is  presented  in  the  datafile  code  below  under   “COUNTING SHEAR AND TENSILE FRACTURES”.   ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ new ro 0.0001 set ovtol 0.01 ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;GEOMETRY, VORONOI AND MESHING  bl ‐0.06505 ‐0.084 ‐0.06505 0.084 0.06505 0.084 0.06505 ‐0.084 cr ‐0.06505 0.064 0.06505 0.064 cr ‐0.06505 ‐0.064 0.06505 ‐0.064 cr ‐0.02505 ‐0.064 ‐0.02505 0.064 cr 0.02505 ‐0.064 0.02505 0.064 132  jdelete delete bl range ‐0.05,‐0.03 ‐0.01,0.01  delete bl range 0.03,0.05 ‐0.01,0.01 vo e 0.005 i 1 ro 0.0001 seed 050582 range ‐0.02505,0.02505 ‐0.064,0.064 gen quad 0.02 range ‐0.02505,0.02505 ‐0.064,0.064  gen edge 0.01 range ‐0.02505,0.02505 ‐0.064,0.064  table 1 ‐0.02505,‐0.064 ‐0.02505,0.064 0.02505,0.064 0.02505,‐0.064 ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;MATERIAL AND JOINT PROPERTIES  ;ASSIGNING ELASTIC MATERIAL PROPERTIES TO ROCK MASS group zone 'block' range ‐0.02505,0.02505 ‐0.064,0.064 zone model elastic density 2.65e3 bulk 3.056e10 sh 2.292e10 range group 'block'  ;ASSIGNING ELASTIC MATERIAL PROPERTIES TO PLATENS group 'platens' range outside table 1 gen quad 1 range group 'platens' zone model elastic density 7.750e3 bulk 1.60e11 sh 7.9e10 range group 'platens'  ;ASSIGNING JOINT PROPERTIES TO VORONOI JOINTS EQUAL TO ROCK MASS PROPERTIES group joint 'rockmass' joint model  residual  jks 1.60e11  jkn 1.60e12  jfric 50  jcoh 80e6  jt 25e6  jrf 31  jresc 0  jrt 0  range group 'rockmass' set jcondf residual jks 1.60e11 jkn 1.60e12 jfric 31 jcoh 0 jdil 0 jt 0 jrf 31 ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;BOUNDARY CONDITIONS (FIX BOTTOM PLACE AND COMPRESS UPPER PLATEN) bound yvel=‐0.02 range bl 342 bound yvel=0.0 range bl 202 bound xvel=0.0 range group 'platens'   133  def temp ntab = 1 end temp  ;AUTOSAVE FUNCTION def autosave name1 = 'UCS1_'+cycID+'.sav' end ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;DEVELOPMENT OF THE AXIAL STRESS/AXIAL STRAIN CURVE  ;DEFINE STRESS‐STRAIN CURVE def slip_load ntab = ntab + 1 tot_str = 0.0 n_z = 0.0 x_z = 0.0 loop n (1,23) x_z = float(n) Temp11 = ‐0.02505+x_z/455 iz = z_near(Temp11,0.064) tot_str = tot_str+z_syy(iz) n_z = n_z +1 endloop  ;AXIAL STRAIN MEASUREMENT agp_ydisp = gp_near(0,0.06401) astrain_ydisp = gp_ydis(agp_ydisp) ax_strain = ‐astrain_ydisp/0.128   134  ;AXIAL STRESS MEASUREMENT ave_astress = ‐tot_str/n_z  ;RECORDING MEASUREMENTS xtable (2,ntab) = ax_strain ytable (2,ntab) = ave_astress end table 2 (0,0) ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;COUNTING SHEAR AND TENSILE FRACTURES  define failurecount ten1_count = 0 ten2_count = 0 sh_count = 0 tot_count = 0 peak_strength = 138e6 Temp12 = contact_head loop while Temp12 # 0 tot_count = tot_count + 1 if c_ndis(Temp12) > 1e‐4 ten2_count = ten2_count + 1 if c_sdis(Temp12) < 1e‐4 ten1_count = ten1_count + 1 endif endif if c_ndis(Temp12) < 1e‐4 if c_sdis(Temp12) > 1e‐4 sh_count = sh_count + 1 else if c_sdis(Temp12) < ‐1e‐4 sh_count = sh_count + 1 135  endif  endif endif norm_ten1 = float(ten1_count)/tot_count norm_ten2 = float(ten2_count)/tot_count norm_sh = float(sh_count)/tot_count norm_astress = ave_astress/peak_strength xtable (3,ntab) = norm_astress ytable (3,ntab) = norm_ten1 xtable (4,ntab) = norm_astress ytable (4,ntab) = norm_ten2 xtable (5,ntab) = norm_astress ytable (5,ntab) = norm_sh Temp12 = c_next(Temp12) tbs = table_size(3) endloop end table 3 (0,0) table 4 (0,0) table 5 (0,0) ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ;APPLY LOAD def loop_load loop h (1,1700) cycID = h*500 command cycle 500 endcommand ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ loop sk1 (1,20) if h = sk1*100 autosave 136  command save name1 endcommand endif end_loop  slip_load failurecount end_loop end  loop_load ;‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0223119/manifest

Comment

Related Items