Back analysis of fault-slip in burst prone environment

Journal of Applied Geophysics 134 (2016) 159–171

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Back analysis of fault-slip in burst prone environment Atsushi Sainoki ⁎, Hani S. Mitri Department of Mining and Materials Engineering, McGill University, Montreal H3A 0E8, Canada

a r t i c l e

i n f o

Article history: Received 11 May 2015 Received in revised form 23 June 2016 Accepted 6 September 2016 Available online 10 September 2016 Keywords: Fault-slip Seismic source parameters Peak particle velocity Mining-induced seismicity

a b s t r a c t In deep underground mines, stress re-distribution induced by mining activities could cause fault-slip. Seismic waves arising from fault-slip occasionally induce rock ejection when hitting the boundary of mine openings, and as a result, severe damage could be inflicted. In general, it is difficult to estimate fault-slip-induced ground motion in the vicinity of mine openings because of the complexity of the dynamic response of faults and the presence of geological structures. In this paper, a case study is conducted for a Canadian underground mine, herein called “Mine-A”, which is known for its seismic activities. Using a microseismic database collected from the mine, a back analysis of fault-slip is carried out with mine-wide 3-dimensional numerical modeling. A back analysis is conducted to estimate the physical and mechanical properties of the causative fracture or shear zones. One large seismic event has been selected for the back analysis to detect a fault-slip related seismic event. In the back analysis, the shear zone properties are estimated with respect to moment magnitude of the seismic event and peak particle velocity (PPV) recorded by a strong ground motion sensor. The estimated properties are then validated through comparison with peak ground acceleration recorded by accelerometers. Lastly, ground motion in active mining areas is estimated by conducting dynamic analysis with the estimated values. The present study implies that it would be possible to estimate the magnitude of seismic events that might occur in the near future by applying the estimated properties to the numerical model. Although the case study is conducted for a specific mine, the developed methodology can be equally applied to other mines suffering from fault-slip related seismic events. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Fault-slip that occurs in underground mines occasionally entails intense seismic waves, which could inflict severe damage when hitting the surfaces of nearby mine openings (Ortlepp, 2000). The damage to rockmasses can be predicted on the basis of peak particle velocity (PPV) caused by the seismic waves (Saharan, 2004), seismic moment and seismically radiated energy of the seismic events (Hedley, 1992). PPV is defined as the maximum particle velocity at an arbitrary location in a transmitting medium in the course of the propagation of stress (seismic) waves. It is widely used as a threshold for damage to the rockmass (Brinkmann, 1987) and is associated with the severity of damage caused by mining-induced seismic events (Hedley, 1992). To date, a number of studies have been conducted numerically and analytically to estimate these values. McGarr (2002) proposed a formulation to predict near-field PPV provided that seismic efficiency of mininginduced tremors is less than 6% (McGarr, 1999). In this context, the term “near-field” implies the vicinity of a source location of a seismic event where it is unnecessary to consider seismic wave attenuation or propagation. Note that seismic efficiency is defined as a ratio of seismically radiated energy to the change in potential energy calculated from ⁎ Corresponding author. E-mail address: [email protected] (A. Sainoki).

http://dx.doi.org/10.1016/j.jappgeo.2016.09.009 0926-9851/© 2016 Elsevier B.V. All rights reserved.

(τ1 + τ2)DA/2, where τ1 and τ2 are the stresses before and after an earthquake; A is a fault area loaded to failure by the applied stress τ1; D is a final slip when the shear stress decreases to τ2 (Kanamori, 2001; McGarr, 1994). Back analysis to simulate seismic source parameters has been also carried out by many researchers. Hofmann and Scheepers (2011) simulated the fault-slip area by calibrating the cohesive strength of a causative fault with seismic moment. Potvin et al. (2010) and Sjöberg et al. (2012) conducted back analysis with seismic moment, considering an increase in shear displacements on causative geological structures for each mining stage. Although the seismic moment of fault-slip related seismic events and near-field particle velocity excited by seismic waves might be approximated with the aforementioned methods, it is not still straightforward to estimate the PPV in the vicinity of mine openings for several reasons. First, the formulation proposed by McGarr (2002) is applicable only in the vicinity of source regions; hence ground motion in active mining areas cannot be estimated in the case that the seismic events are located far from the mining areas due to regional unclamping (Castro et al., 2009). Here, regional unclamping denotes stress change, which results in the reduction of normal stress acting on a fault, taking place in an extensive region caused by the extraction of a large amount of ore. Secondly, previous back analyses were carried out under static conditions by means of the classical Mohr-Coulomb criterion, which implies that the results obtained from the back analyses are not validated

160

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

with respect to the dynamic behaviour of rockmass, such as PPV. Furthermore, as reported by Ryder (1988), the cause of fault-slip in underground mines is a sudden stress drop induced by asperity shear on fault surfaces. The classical Mohr-Coulomb criterion cannot take fault surface asperities into account. As the seismic source parameters of fault-slip are strongly dependent on the fault surface properties (Sainoki and Mitri, 2013), it is postulated that the intensity of seismic waves is also dictated by the fault surface properties. Static analyses adopted in previous studies cannot simulate the propagation of the seismic waves nor the ground motion created by the seismic waves. For these reasons, there is still a pressing need to develop a robust methodology to simulate fault-slip that occurs in underground mines, so that the magnitude of particle velocity excited by seismic waves within active mining areas can be estimated. Based on the estimated particle velocity, the selection of appropriate support system to control ground motion due to the seismic waves can be achieved. This is of paramount importance for steady mining production and the safety of mine operators. In the present study, back analysis of a fault-slip related seismic event that took place in “Mine A” is attempted under dynamic conditions, considering asperity shear as the cause of the seismic event. The back analysis is carried out with respect to the physical and mechanical properties of causative fracture/shear zones within which the seismic event took place. Moment magnitude (i.e. seismic moment) estimated from recorded waveforms in the mine and far-field PPV recorded by a strong ground motion sensor are used for the calibration of the fault properties. The estimated values are then validated with ground acceleration recorded by accelerometers. 2. Case study mine The case study mine, herein called “Mine-A” is a base metal operation extracting ore primarily with sublevel stoping method with delayed backfill. There are two mining ore zones: shallow and deep. The upper zone is situated above the 1400 m level and is completely mined out. Current mining activities take place in the deeper zone below the 1400 m level, which consists of two orebodies #1 and #2 striking nearly in the east-west direction. Fig. 1 shows a schematic of the two orebodies on the 1500 m level. The #1 orebody extends from surface to 1700 m level and dips generally at 75° to the south with an average width of 9 m. The #2 orebody extends from 850 m level to 1800 m level and dips sub-parallel to the #1 orebody at 60° to the south. An important geological structure in the area is the shear/fracture zone identified in the vicinity of the orebodies, which is shown in Fig. 1. The width of the shear zone ranges from 15 m to 150 m, the strike is north-west (dip direction is 45°), and the dip is 85° to the east. The rocks in Mine-A are generally classified into five types, namely Norite (NR), Greenstone (GS), Olivine Diabase (OLDI), Massive Sulphide

(MASU), and Metasediment (MTSD). Norite and Metasediment are generally located in the footwall and hanging wall, respectively. Greenstone appears to be a subdivision of the more massive rock within the metasediments of the hanging wall, and the orebodies consist of MASU. The mechanical properties of the rock masses are shown in a later section. 3. Microseismic analysis Microseismic activity in Mine-A has been recorded for more than 5 years with a microseismic monitoring system that consists of 29 uniaxial accelerometers, 5 triaxial accelerometers and one strong ground motion sensor (4.5 Hz geophone). The accelerometers installed in the mine have a sensitivity of 30 V/g and a clipping Voltage of 3.9 Vare, whereby acceleration up to 1.3 m/s2 can be measured as a maximum. The strong ground motion (SGM) sensor is, on the other hand, installed near the ground surface and capable of capturing large seismic events with low frequencies, of which moment magnitude ranges from 0 to 3. By utilizing these types of instruments, both of small and large seismic events can be detected, and seismic source parameters for each seismic event are calculated on the basis of waveforms recorded by the microseismic monitoring system. The calculations are conducted with ESG software (ESG, 2011), which is a specialized seismic data analysis software widely employed for various kinds of field applications related to induced seismicity. The overview of the monitoring system being used at the mine and an actual case study is summarized by Urbancic and Trifu (2000). In addition, a number of case studies have been undertaken with monitoring systems similar to that installed in this mine, not only for understanding seismicity in underground mines but also for investigating induced seismicity related to long-term production, hydraulic fracturing, and CO2 sequestration all over the world (Li et al., 2014; Maxwell and Urbancic, 2001; Trifu and Shumila, 2010). There are a number of sources for seismic events in underground mines, such as production blasts, generation and nucleation of extension fractures, and propagation of shear rupture. The seismic data analysis software being used at the mine is capable of distinguishing blastinduced seismic events from the others with recorded waveforms. However, it is still necessary to ascertain the seismic events caused by a shear movement along a discontinuity. This is referred to as faultslip related seismic event in this paper. The back analysis conducted in the present study aims at simulating a fault-slip related seismic event and estimating the properties of shear/fracture zones where the seismic event took place. To do so, the present study adopts the magnitude time history (MTH) analysis and time between events (TBE)-rate chart proposed by Hudyma and Beneteau (2012). The fundamental idea behind the methods is that seismic events related to rockmass shearing and a slip along a fault take place without having a strong correlation with blasting activities because the direct triggering mechanism of fault-

Fig. 1. Schematic illustration of the orebodies and intersecting shear zone at 1500 m below surface.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

slip related seismic events is an instantaneous stress drop caused by the shearing of asperities on a fault surface. Such stress drop does not necessarily occur immediately after stress re-distribution due to mining activities. Although it is true that blast-induced stress waves could induce fault-slip (Sainoki and Mitri, 2014d), it occurs only if a fault is located in the vicinity of the blasts and is under a critical condition prior to the blast, i.e., the shear stress acting on the fault almost reaches its maximum shear strength. Detailed descriptions of MTH analysis and TBErate method are provided in the literature (Hudyma and Beneteau, 2012). With the methods, the sources of seismic events can be identified for each cluster of seismic events that occur in the same area. Hence, as the first step to apply the methods, seismically active zones in Mine-A are divided into 6 regions as shown in Fig. 2. Table 1 shows the number of seismic events used for the analyses of microseismic database. TBE-rate chart is then created for the 6 regions, which is shown in Fig. 3. According to Hudyma and Beneteau (2012), seismic events with a low TBE-rate indicates that those are induced by stress re-adjustments that occur immediately after mine developments and mining activities. On the other hand, seismic activities with a high TBE-rate can be considered to be associated with rockmass shearing and fracturing. It is found from the figure that clusters of seismic events that occurred in regions 5 and 6 show comparatively high TBE-rates for the whole range of event magnitude. Note that a TBE-rate corresponds to the slope of the lines depicted in Fig. 3, that is, for instance, the near horizontal line in a range of magnitude from −2.0 to −1.0 for region 2 represents a low TBE-rate in Fig. 3. Thus, the regions are further divided into small zones to carry out more detailed analysis as shown in Fig. 4. MTH analysis and TBE-rate chart are applied to the additionally divided zones. Fig. 5 shows the results of magnitude time history analysis performed for regions 5-1 and 5-4. The increase in the cumulative number of seismic events at a nearly constant rate shown in Fig. 5(a) indicates that the seismic events are related to fracturing of rockmass and slip along geological discontinuities, i.e. fault-slip (Hudyma and Beneteau, 2012). On the other hand, Fig. 5(b) shows the stepwise increase, indicating that the seismic events that took place in region 5-4 are strongly related to blasting activity. The same criteria are applied to the results obtained from MTH analyses performed for the other subregions. This paper presents only part of the results. Based on the results of the MTH analyses in conjunction with TBE-rate charts, it is deduced that the seismic events that occurred in regions 5-1, 5-2, and 6-1 are strongly associated with fault-slip (rock fracturing and slip along geological discontinuities). Considering the fact that the three include the shear zone, the results of the analyses appear to be quite reasonable. Based on that, a large seismic event that occurred in region 5-1 in the past has been chosen for the back analysis. Table 2 shows the location and source parameters of the seismic event, which have been estimated on the basis of recorded waveforms, using the specialized seismic waveform analysis software described above.

Fig. 2. Divided regions for MTH analysis and TBE-rate chart.

161

Table 1 Number of seismic events used for MTH and TBE-rate analyses for each region shown in Fig. 2. 1 Number of seismic events

2

3

4

5

6

Total

43,811 66,223 32,922 17,819 91,019 12,612 356,469

4. Numerical model of Mine-A A 3D mine-wide model encompassing the major geological structures is constructed by means of FLAC3D code, which employs an explicit finite difference method based on continuum mechanics. Using the numerical model, static and dynamic analyses are performed in the present study. The static analysis aims to simulate the pre-mining stress state in the mine as well as the stress state at that time when the targeted seismic event took place. To be more precise, mining activities are simulated from the beginning of the mining sequence in the orebodies, and the static analysis is continued until simulating the situation of mining activities immediately after the occurrence of the targeted seismic event. The static simulation of the pre-mining stress state is carried out while applying tractions to the model outer boundaries (Sainoki and Mitri, 2014c; Shnorhokian et al., 2014). As a result, realization of stress variations resulting from the difference in the stiffness of the rockmasses is achieved. As suggested by Bewick et al. (2009), the simulation method becomes more crucial when the simulation of insitu stress state of a mine with complex geological conditions is attempted. The dynamic analysis is conducted to simulate the seismic event shown in Table 2 and to estimate the physical and mechanical properties of the causative shear zone with back analysis. In this section, first, the geometry of the mine-wide model is shown. Then, the simulation techniques employed as well as the geotechnical properties of the rock masses in the mine are shown. 4.1. Numerical model descriptions Fig. 6 shows the geometry of the mine and its geology. As shown, a very large model has been created in order to generate pre-mining stress state and to prevent model boundaries from affecting analysis results. The width, height, and length of the model are 1220 m, 610 m, and 1220 m, respectively. The shear zone shown in Fig. 1 corresponds to that indicated in Fig. 6 and is modeled while assigning a constitutive model and mechanical properties different from those of the surrounding

Fig. 3. Time-between-event (TBE)-rate chart. In this figure, regions 1 to 6 correspond to those indicated in Fig. 2.

162

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

approximately 90 times difference in zone edge length between coarse and dense discretization. Note that the source region is composed of a number of densely discretized zones representing a failure plane. The total number of zones and grid points in the model are 1146310 and 1186080, respectively. 4.2. Numerical simulation of the pre-mining stress state in the mine

Fig. 4. Subregions further divided from regions 5 and 6 in Fig. 2.

rockmasses to the discretized zones situated within the shear zone in Fig. 1. In this regard, the geometry of the shear zone in Fig. 1 is simplified as it is not comparatively as irregular as those of the other geological structures on the scale considered in the present study. Specifically, the width is approximated to be 60 m, and dip angle and dip direction are assumed to be 85° and 45°, respectively, based on the geological information of the mine. The source region of the seismic event is modeled inside of the shear zone as shown in Fig. 7, presuming that the selected seismic event took place by failure along pre-existing weak planes within the shear/fracture zone. Although, in reality, the orientation of weak planes in the shear zone does not necessarily coincide with that of the shear zone, the dip and dip angle of the source region (failure plane) are assumed to be the same as those for the shear zone. As verified in a later section, the assumption is deemed acceptable as a premise to carry out back analysis. The modeling procedure of the source region is the same as that for the shear zone, i.e., a constitutive model and mechanical properties different from those of the shear zone are assigned to the discretized zones situated in the source region. The dimensions of the source region are determined from the source radius shown in Table 2. As can be seen in Fig. 6, orebodies are located in the center of the model. Stopes are modeled within the orebodies and sequentially extracted in accordance with mining sequences of Mine-A. Careful consideration has been given in generating meshes around the active mining areas, where zones have been densely discretized, while in the vicinity of the model boundaries the number of grid points is reduced in order to decrease calculation time. Specifically, the edge length of a zone located near and in the active mining areas including the source region ranges approximately from 0.6 m to 4 m, while zones near the model outer boundaries have an edge length of 56 m at a maximum. There is

The pre-mining stress state in the mine-wide model is simulated with the boundary traction method (Shnorhokian et al., 2012, 2013) and with the simulation technique allowing for the heterogeneity of the shear stiffness of rockmass within shear zone (Sainoki and Mitri, 2014c). Note that stresses applied to the model outer boundaries are calibrated with the boundary traction method so that simulated stress state coincides with that at a location where stress measurement was conducted at the mine. Accordingly, the stresses applied to the model boundaries are not necessarily consistent with those of far-field regional stresses. The stress measurement was taken after stope extraction had started, but the effect of mining-induced stress is small enough to be considered negligible for the calibration of pre-mining stress state, considering that the measurement was conducted far away from active mining areas. The reason why the two methods are employed to simulate the pre-mining stress state in the mine is as follows. Initially, the simulation of the pre-mining stress state was attempted only with the boundary traction method. The boundary traction method permits simulations of highly stressed regions caused by the difference in stiffness amongst rock masses. As can be seen in the previous studies (Bewick et al., 2009), a stiff rockmass, such as a dyke, is strongly associated with mining-induced seismicity, indicating that it stores high elastic strain energy because of its high stiffness, compared to that of the surrounding rockmass. With an aim to take into account the effect of the complex geological structures in the mine (Fig. 6) in the pre-mining stress state, the preliminary analysis employed the boundary traction method, whereby calibrated stresses are applied to the model outer boundaries. The result obtained from the analysis revealed that although the high stress within the dyke was successfully simulated with the simulation technique, fault-slip potential defined by the difference between the maximum shear strength and the shear stress at the source location of the targeted seismic event takes a negative value with respect to the strike and dip angle of the shear zone. Fault-slip potential is at most − 26 MPa, even if different orientations of the weak plane are considered. It indicates that any fault-slip related seismic events would never happen in the region. Thus, it was found that further consideration with respect to the heterogeneity of rock mass in the shear zone needs to be given. During 2012, 33 large seismic events were detected with the strong ground motion sensor installed near the ground surface. It was then found that the source locations of

Fig. 5. Magnitude time history analysis: (a) region 5-1 and (b) region 5-4.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

163

Table 2 Seismic source parameters of the seismic event selected for back analysis.

Source information

North

East

Depth (m)

Seismic moment (N·m)

Moment magnitude

Source radius (m)

−305

793

1520

1.6 × 1011

1.47

46.7

most large seismic events are situated within the shear zone where fault-slip potential is estimated to be negative according to the numerical analysis not considering heterogeneity. Such consistent finding supports two ideas: (i) the targeted seismic event did take place within the shear zone even after accounting for the location error associated with the events, and (ii) it is necessary to incorporate stiffness heterogeneity of the shear zone in order to produce an in-situ stress state inducing fault-slip in the shear zone. Sainoki and Mitri (2014c) demonstrate that fault-slip potential within shear zones can significantly increase at the location where shear stiffness heterogeneity exists within the shear zone. According to Bandis et al. (1983) and Brady and Brown (1993), shear stiffness of rockmasses including joints are strongly dependent upon joint spacing and joint surface roughness. Thus, it is natural that shear stiffness of shear/fracture zones varies with the properties of discontinuities within the shear zone, consequently producing regions with high shear stiffness. The previous study (Sainoki and Mitri, 2014c) numerically confirms that high fault-slip potential is retained inside the region with high shear stiffness. Based on the same concept as the previous work, the present study simulates the heterogeneity of shear stiffness within the shear zone. In order to simulate the difference in shear stiffness, a transversely isotropic model, which was initially developed by Bruggeman (1935) and used to simulate fault zones (Ivins and Lyzenga, 1986), is applied to the shear zone, excluding the source region of the seismic event. For the source region, it is necessary to apply a constitutive model that is capable of simulating fault behaviour during fault-slip. According to Ryder (1988), fault-slip in underground mines is driven by a sudden stress drop due to breakage of asperities on fault surfaces, although stress re-distribution resulting from mining activities is the trigger for the fault-slip. To date, a number of shear strength models that allow for the physical properties of fault surfaces, such as roughness and asperity height, have been developed (Barton et al., 1985; Barton and Choubey, 1977; Ladanyi and Archambault, 1969). Sainoki and Mitri

(2014a) incorporated Barton's shear strength model into a ubiquitous joint model using FLAC3D, a three-dimensional explicit finitedifference program, with C++ programming language and simulated fault-slip induced by stress changed due to stope extraction, considering asperity shear. In the present study, the same constitutive model is applied to zones within the source region in order to simulate fault-slip due to asperity shear. Barton's shear strength model used for the ubiquitous joint model is expressed as:

  JCS þ φb τmax ¼ σ n JRC log σn

ð1Þ

where τmax and σn represent the maximum shear strength and normal stress acting on a joint, respectively; and JRC, JCS and ϕb represent joint roughness coefficient, joint wall compressive strength and basic friction angle (static friction angle), respectively. As can be seen from Eq. 1, the shear strength model can take joint roughness coefficient into consideration. It is also found from the equation that the maximum shear strength derived from Barton's model is always greater than that obtained from Mohr-Coulomb failure criterion when cohesive strength = 0. During dynamic analysis to simulate fault-slip, when shear stress on failure planes exceeds the maximum shear strength determined by Eq. 1, JRC in the equation is changed to 0 since asperities on the failure planes are assumed to be sheared off due to slip. The residual stress arising from the change in JRC is the fault-slip driving force during the dynamic analysis. Regarding the ubiquitous joint model, its weak plane along which slip occurs is assumed to have the same dip and dip angle as those of the shear zone. The seismic data analysis software performs a moment tensor inversion analysis, which gives 77° and 16° as the dip angle and dip direction of a fault plane of the targeted seismic event, respectively. The dip angle is almost the same as that of the shear zone, while the dip direction differs by 29° from that of the assumed failure plane. Although the orientation of a fault plane derived from the

Fig. 6. Numerical model geometry and geology in Mine-A.

164

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

Fig. 7. Plan view of 1500 m level showing source region.

moment tensor inversion might be accepted as a premise for back analysis, the present study does not adopt the estimated values. The reason stems from the fact that a large seismic event caused by fault-slip takes place along a “pre-existing” weak plane (Ortlepp and Stacey, 1994). The slip plane orientation derived from the moment tensor inversion analysis differs from that of the shear zone. Although it is true that orientations of weak planes in the shear zone vary from place to place, it is difficult to verify the correctness of the inversion analysis. In addition, taking into account such deviation makes it more complicated to simulate fault-slip that could occur in the future because the orientation of a weak plane through which the slip takes place needs to be known in advance. It is unrealistic to map all weak planes in the shear zone. It is obviously more practical if fault-slip with seismic source parameters comparable to an actual one can be simulated assuming that the orientation of the failure plane is the same as that of the shear zone. Hence, the present study makes the assumption and discusses the difference between accelerations obtained from fault-slip simulated with backanalyzed properties and those recorded with the seismic monitoring system. 4.3. Rockmass properties As described in the previous section, the rock masses in Mine-A can be classified into five types, namely Norite (NR), Greenstone (GS), Olivine Diabase (OLDI), Massive Sulphide (MASU), and Metasediment (MTSD). The mechanical properties of the rock masses derived from laboratory experiments are listed in Table 3. As mentioned above, the transversely isotropic model is applied to the finite difference zones within the shear zone, excluding the source region of the seismic event, in order to simulate high slip potential condition resulting from the difference in shear stiffness of rockmasses. By means of the constitutive model, shear moduli of rockmasses with respect to a plane with the same dip and dip angle as those for the shear zone are changed to small values compared to those of the surrounding rockmasses. Preliminary numerical analysis was carried out to estimate the shear moduli of rockmasses within the shear zone. The procedure of the numerical analysis is shown in Fig. 8. As indicated in Fig. 8, in the Table 3 Mechanical properties of rockmasses in Mine-A.

Modulus of elasticity (GPa) Poisson's ratio

NR

OLDI

GS

MTSD

MASU

56 0.25

86 0.26

65 0.23

45 0.24

44 0.30

static analysis, pre-mining stress states were first simulated with the boundary traction method. Stopes in the orebodies were then extracted in accordance with mining sequences provided by the mine until the same period in which the seismic event occurred. The purpose of this preliminary analysis is to examine slip potential; hence slip behaviour was not simulated within the source region. A number of models with different shear moduli were generated, and the analysis procedure was carried out for each model. It was found out from the analyses that slip potential of rockmasses within the source region became positive when shear moduli of rockmasses within the shear zone except the source region was one-tenth of that of the rockmasses where it occurs. As reported by Morrow et al. (1982) and Budiansky and O'Connell (1976), cracks and fractures could decrease shear moduli of rockmasses up to one-fifteenth. In addition, Bandis et al. (1983) proposed an equation that relates joint surface roughness with the shear stiffness of the joint, and shear moduli of rockmass are directly affected by the shear stiffness of joints included in the rockmass. With the equation (Bandis et al., 1983), it is readily found that the shear moduli of rockmasses could drastically change, depending on joint properties. Hence, the result obtained from the preliminary analysis is not an unreasonable value. Table 4 shows the shear moduli of rockmasses within the shear zone excluding the source region. For the source region, the ubiquitous joint model with Barton's model is applied. Hence, it is necessary to find static (basic) friction angle of rockmasses with the region. In the present study, basic friction angles shown in Table 5 are estimated on the basis of publications (Alejano et al., 2012; Barton, 1973). With back analysis carried out in the present study, determining the physical and mechanical properties of the causative shear/fracture zones, namely fault surface roughness (JRC) and dynamic friction angle (ϕd), is attempted. First, the shear moduli of rockmasses within the shear zone is determined with the aforementioned methodology (Fig. 8) that considers the shear stiffness heterogeneity between the source region and the rest of the shear zone. Subsequently, additional preliminary analysis is conducted. The aim of the analysis is to estimate a friction reduction factor (FRF) applied to the basic friction angles (ϕb) so that back-analyzed ϕd falls within a reasonable range. The magnitude of fault-slip is determined by two factors: (i) the magnitude of stress drop resulting from the transition from ϕb to ϕd due to asperity shearing and (ii) stress state in the source region when the fault-slip takes place. The two influencing factors are affected by the static and dynamic frictional resistance, that is, there is a combination of JRC and an unrealistically low ϕd, whereby fault-slip with magnitude comparable to that shown in Table 2 is simulated. In fact, it was found out that when the basic friction angles in Table 5 are applied, an unrealistically small ϕd is obtained as a result of back analysis. In order to address the issue, FRF is employed, by which ϕb in Table 5 are multiplied. The procedure of estimating FRF is shown in Fig. 9. As shown in the figure, this analysis is based on the numerical model with the estimated shear stiffness. As the first step, the developed ubiquitous joint model is applied to the source region, setting JRC to zero. This is because when JRC = 0, the in-situ stress state in the source region is most prone to fault-slip; hence back-analyzed ϕd will take the maximum value, compared to other models with higher JRC. If the back-analyzed ϕd is less than a plausible range of ϕd, it is necessary to decrease FRF applied to ϕb in order to bring the back-analyzed ϕd to a reasonable value. As shown in Fig. 9, after performing a static analysis to simulate in-situ stress state, a dynamic analysis is carried out, in which fault-slip is simulated while changing the friction angle from static to dynamic only for zones where the shear stress acting on the assumed fault plane reaches the maximum shear strength as determined by Eq. (1). The dynamic analysis is repeatedly conducted until back-analyzing ϕd whereby fault-slip with magnitude comparable to that shown in Table 2 is simulated. As explained above, if the back-analyzed ϕd is smaller than a reasonable value, FRF is reduced. Then, the static analysis is performed again with reduced ϕb. In this way, FRF is estimated.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

165

Fig. 8. Methodology to generate fault-slip potential at the source location of a seismic event within a weak shear zone.

The plausible range of ϕd used to estimate FRF in Fig. 9 is estimated from previous numerical, analytical and experimental studies. Accordingly, back-analyzed ϕd in a later section with peak particle velocity (PPV) and moment magnitude falls into the estimated range of ϕd. ϕd is strongly dependent on slip rates during fault-slip. According to a frictional resistance test with a modified torsional Kolsky bar performed by Yuana and Prakash (2008), the coefficient of kinetic friction of rockanalog materials ranges from 0.2 to 0.3 when slip rate is between 2 and 5 m/s. The experimentally obtained coefficient of kinetic friction (Yuana and Prakash, 2008) is approximately equivalent to friction angles between 12° and 17°. Note that the friction angles shown in Table 5 are basic friction angles measured under static conditions. According to Dieterich and Kilgore (1996), when slip rates are less than 1 mm/s, the coefficient of friction of granite drops to 0.65, which is approximately equivalent to a friction angle of 33°. McGarr (1991) estimates slip rates of actual fault-slip events that took place in a deep underground mine. The estimated slip rates range from 1 m/s to 5 m/s, which are estimated from far-field S wave pulses by McGarr (1991). Sainoki and Mitri (2014b) simulate fault-slip in dynamic conditions, considering slip-weakening distance defined as the distance over which the shear stress acting on a fault decreases to a residual value during slip. The study shows that slip rates of fault-slip taking place in a deep underground mine vary from 0.6 m/s to 3.3 m/s, depending on the slip-weakening distance. The information on the dynamic frictional resistance and slip rate during fault-slip obtained from previous studies is utilized to constrain the range of back-analyzed ϕd. Consequently, it would be reasonable to assume that the dynamic friction angle of a fault during fault-slip in underground mines is in the vicinity of 20°. The range of ϕd is used to estimate FRF as discussed above and shown in Fig. 9. As a result of iterative analyses, FRF is estimated to 0.8, with which the back-analyzed ϕd falls within the reasonable range when JRC = 0 is assumed, i.e., when FRF is greater than 0.8, back-analyzed

ϕd is less than the reasonable range even if JRC is set to zero, meaning that ϕd is back-analyzed to far less than the plausible range when JRC is set to greater than zero because static shear resistance in the source region increases with JRC. 4.4. Analysis conditions of static and dynamic analysis During the static analysis to simulate stepwise stope extraction as per mining sequence at the mine, the model outer boundaries are constrained. The constrained boundary condition is widely employed for stress analysis with a mine-wide model where the area of interest is sufficiently distant from model outer boundaries (Bewick et al., 2009), so that stress re-distribution caused by the sequential stope extraction is unaffected by the model boundary. As shown in Fig. 6, the orebodies are located far away from the model outer boundaries. For dynamic analysis to simulate fault-slip in the source region, the boundary conditions are changed to viscous in order to prevent the model boundaries from reflecting seismic waves arising from the fault-slip. The viscous boundary condition is based on the use of independent dashpots (mechanical viscous dampers) in the normal and shear directions at the model boundaries (Lysmer and Kuhlemeyer, 1969) in order to effectively absorb the energy of seismic waves, particularly when angles of incidence are greater than 30° (Itasca, 2009). The timestep used for the dynamic analysis is automatically calculated based on the volume of each zone of the model, the P-wave velocity obtained from rockmass mechanical properties and the face area of each zone (Itasca, 2009). Regarding the attenuation of seismic waves, local damping embodied in FLAC3D is applied to the numerical model. The damping system is defined as follows (Itasca, 2009):
F i þ F ðiÞ ¼ M

dvi dt

 ð2Þ

Table 4 Mechanical properties of rockmasses within the shear zone outside source region of the seismic event.

Shear modulus of any planes normal to the weak plane with dip angle and direction of 45° and 85° (GPa) Shear modulus in the plane of isotropy (GPa)

NR

OLDI

GS

MTSD

MASU

2.3 23

3.4 34

2.6 26

1.8 18

1.7 17

166

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

Table 5 Basic friction angles (ϕb) of rockmasses within source regions of the seismic event.

Basic (static) friction angle (°)

OLDI

GS

36

32

where Fi, M, and v represent force acting on a node, nodal mass, and velocity; and F(i) denotes damping force. As shown in Eq. (2), the damping force is directly added to the equation of motion, and then the equation is solved. The damping force, F(i), is given as follows:

F ðiÞ ¼ −πDj F i jsignðvi Þ

D¼

ΔW=W 4π

8 < þ1 if y N 0 signðyÞ ¼ −1 if y b 0 : 0 if y ¼ 0

ð3Þ

ð4Þ

where ΔW and W are the amount of energy removed per cycle during oscillation and the mean kinetic energy at the instant of removal, respectively; and D is a local damping coefficient related the ratio of ΔW to W, respectively. As indicated with Eq. (4), the damping coefficient, D, essentially represents the intensity of damping, i.e., as D increases, more energy is removed per a cycle of oscillation, approaching a critically damped condition where no oscillation occurs. Hence, D is directly associated with a damping ratio defined as an actual damping coefficient

divided by a critical damping coefficient (Itasca, 2009). Damping ratio = 1 corresponds to a critically damped condition. The local damping system embedded in FLAC3D operates by adding or subtracting mass from a grid point at certain times during a cycle of oscillation. Although local damping cannot capture the energy loss properly when waveforms are complicated, it enables frequencyindependent damping. According to ABAQUS (2003), a damping ratio for rock falls into 0.02 to 0.05. Thus, 2% of critical damping is tentatively adopted as an initial value, i.e., 0.02 is substituted to D, which is used to calculate damping force expressed by Eq. (3). An appropriate damping coefficient is estimated later through the back analysis by comparing peak particle velocity (PPV) computed from the dynamic analysis with that recorded by the strong ground motion sensor when the targeted event took place. 5. Back analysis Fig. 10 shows the procedure of the back analysis. The numerical model with the FRF and shear modulus of the weak shear zone estimated in previous sections is used. As discussed, 0.02 (2%) is assumed as a tentative value for the damping coefficient, D, which is estimated later. The first step of the back analysis is to ascertain combinations of JRC and ϕd whereby fault-slip with moment magnitude comparable to that shown in Table 2 is simulated. As ϕd is not related to the static analysis, JRC is assumed first, and a static analysis is performed, in which stepwise stope extraction is carried out after simulating the pre-

Fig. 9. Procedure to estimate the friction reduction factor (FRF) applied to the static friction angles in Table 5.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

167

Fig. 10. Back analysis procedure for JRC, ϕd and D.

mining stress state. The value of assumed JRC is arbitrary within the range from 0 to 20 (Barton, 1973). For each assumed JRC, ϕd is estimated. The mining sequence proceeds until the time, at which the seismic event took place. Using the stress state obtained from the static analysis, dynamic analysis is performed, assuming ϕd. During the dynamic analysis, fault-slip is simulated in the source region by setting JRC to zero and ϕb to ϕd for the zones where shear stress exceeds the maximum shear strength. Changing the two mechanical properties is intended to trigger

fault-slip and has a physical meaning that corresponds to shearing of fault surface asperities and reduction in frictional resistance on the fault surface due to slip. On the basis of relative shear displacement increments obtained from the dynamic analysis, moment magnitude is computed. If the computed moment magnitude does not agree with that of the targeted seismic event shown in Table 2, the dynamic analysis is repeated after modifying ϕd. The present study does not consider changes of ϕd in a fraction of a degree, i.e., ϕd is decreased or increased

168

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

at least by one degree for the calibration. This is because the time required for the calibration significantly increases if small increments less than one degree are taken into account in the calibration process. In this way, combinations of JRC and ϕd with which fault-slip with magnitude comparable to that of the targeted seismic event is simulated are obtained. As shown in Fig. 10, the next step is to conduct dynamic analyses with the obtained combinations of JRC and ϕd while changing D. During the dynamic analyses, velocities of several grid points located outside of the source region are recorded, and the peak particle velocity (PPV) is computed for each location. The computed PPV is used to select the most appropriate combination of JRC, ϕd, and D through a comparison with PPV recorded by the strong ground motion sensor installed at the mine. The procedure to select the most appropriate combination is described in more detail later. After determining the most appropriate combination, the procedure moves to a final phase, where validation of a numerical model with the estimated properties is conducted, using ground peak acceleration recorded by accelerometers at the mine. Note that the peak acceleration used to validate the numerical model is obtained from sensors different from the one providing the data of PPV, i.e., the data of acceleration is completely independent of PPV used for the calibration. 5.1. Calibration with moment magnitude Table 6 shows the back-analyzed properties and moment magnitude computed from numerical models using the estimated friction reduction factor (FRF) and shear stiffness. As shown in the table, ϕd has been determined for JRC of 0, 5, and 10 so that computed moment magnitude becomes comparable to that estimated from recorded waveforms. According to Barton and Choubey (1977), JRC of 20 represents the roughest joint surface generated by tension failure. Hence, it would be advisable to assume JRC within a range from 0 to 15 for a shear zone for the back analysis. In the present study, JRC of 15 was found inappropriate because ϕd is back-analyzed to a value smaller than the plausible range when JRC = 15 is applied. As can be found from Table 6, when JRC = 10, the estimated ϕd is 19°, which is approximately a lower limit of the plausible range. The computed moment magnitude is very sensitive to ϕd. An increase or decrease in ϕd by one degree results in change in moment magnitude between 0.5 and 1.0. Thus, it is indeed possible to back-analyze ϕd more precisely if needed. However, as the primary objective of the present study is to propose a methodology of fault-slip taking into account the dynamic behaviour of rockmass, further calibration of ϕd is not conducted. The computed moment magnitude is sensitive to the dynamic friction angle. It should be noted that the solution is non-unique, that is, it is impossible to determine only one combination of JRC and the dynamic friction angle on the basis of moment magnitude because moment magnitude obtained from numerical analysis can be adjusted by changing the combinations as shown in the table. Hence, further calibration is carried out using PPV obtained from strong ground motion (SGM) sensor.

the numerical model, PPVs are investigated at several observation points located on a line between the locations of the SGM sensor and the hypocenter of the seismic event. Then, PPV at the location of SGM sensor is extrapolated from PPVs at the observation points in the model. It is to be noted, as described in the overall procedure shown in Fig. 10, the damping coefficient, D, is also estimated at this stage. For the purpose, dynamic analyses are carried out to simulate faultslip for the obtained combinations while changing D from 0% to 5%, which is the plausible range of damping for rockmass as discussed in the previous section. Fig. 11 shows PPVs on the line when JRC = 0, ϕd = 23°, and D = 0.02 (2%). Note that Fig. 11 is shown as an example and the same procedure is taken for the other combinations in Table 6 and damping coefficients. As shown in the figure, PPV is investigated at six observation points on the line. In order to extrapolate PPV at the location of SGM sensor, many fitted curves are tentatively generated employing various kinds of functions, such as polynomial and exponential ones. As a result, power functions are found to be the most suitable to simulating seismic wave attenuation within the numerical model. Squared coefficient of correlation of the fitted curve shown in the figure is approximately 0.999. Based on the fitted curve, PPV at the location of SGM sensor can be extrapolated. Likewise, PPV at the location of SGM sensor is extrapolated for JRC = 5 and 10, which follow the same trend, although constants of the fitted exponential function are completely different amongst them. Note that, for each combination of JRC and ϕd, dynamic analyses are performed while changing the damping coefficient. A damping coefficient, D, is directly related to seismic wave attenuation, meaning that the intensity of seismic waves observed in active mining areas away from the source region of the seismic event is strongly dependent on D. Thus, it is crucial to estimate the coefficient. Fig. 12 summarizes the results. As can be seen, PPV is less than an observed value for the whole range of damping coefficient when JRC = 0. The result indicates that the combination of JRC = 0, where ϕd = 23° is not appropriate, although moment magnitude obtained from the model is in a good agreement with that estimated from recorded waveforms. In the same way, although a model with JRC = 10 and 0% damping (D = 0) provides exactly the same PPV as the recorded value, 0% damping is unrealistic. Thus, the combination also appears not to be appropriate. Lastly, for the model with JRC = 5 and 0% damping (D = 0), PPV extrapolated from PPVs in the same way as shown in Fig. 11 is greater than the observed value. When damping is changed to 2% (D = 0.02), extrapolated PPV is less than the observed value, indicating that PPV estimated from numerical analysis should coincide with the observed value when an appropriate damping is selected between 0% and 2% for the combination (JRC = 5 and ϕd = 22°).

5.2. Calibration with peak particle velocity In order to determine the most appropriate combination of JRC and ϕd, PPV is examined with respect to the obtained combinations listed in Table 6. The SGM sensor used to record ground velocity at the mine is installed close to the ground surface which is outside of the generated numerical model. In order to calculate PPV at the appropriate distance in Table 6 Back-analyzed JRC and ϕd. Model No.

JRC

ϕd (°)

Moment magnitude

1 2 3

0 5 10

23 22 19

1.51 1.42 1.50

Fig. 11. PPV on a line between SGM sensor and hypocenter of seismic event.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

169

Fig. 12. Extrapolated PPV at the location of SGM sensor. Fig. 13. Comparison of peak ground acceleration.

5.3. Numerical model validation with back-analyzed properties JRC and ϕd have been estimated at 5 and 22°, respectively, by means of the aforementioned procedure. The next step is validating the estimated properties. As mentioned in the previous section, ground acceleration is recorded in Mine-A with 34 accelerometers installed in the environs of active mining areas, which are included inside of the mine-wide model. The peak ground acceleration recorded by the accelerometers is utilized for the validation by comparing that with results obtained from dynamic analysis using the estimated properties. For the validation, only accelerometers that recorded reasonable peak accelerations are picked from amongst all the accelerometers installed at the mine, based on an assumption that peak accelerations should decrease with the distance from the hypocenter of the targeted event. In fact, there are abnormally low or high peak accelerations, which are not used for the validation. The causes for such anomalies are assumed to be local rockmass conditions, such as a shear zone, densely spaced joints, and fractured rockmass. The maximum allowable acceleration for the installed accelerometers is 1.3 m/s2, which would be a reason why low peak accelerations were observed in the vicinity of the source region when the seismic event took place. Some of the installed uniaxial and triaxial accelerometers recorded only waveforms assumed to be noise at that time, which are inappropriate for the validation. Furthermore, waveforms recorded by several accelerometers installed at locations outside of the mine wide model cannot be used. Based on those facts, peak accelerations in the z-direction (depth direction) obtained from 7 uniaxial accelerometers located between 300 m and 740 m from the source region were used for the validation. Fig. 13 shows peak acceleration obtained from the dynamic analysis using the estimated properties. The peak acceleration is determined on the basis of results for the period of 300 ms after the onset of the dynamic analysis. The time period is determined as follows. S-waves propagate at approximately 3 km/s in the host rock (NR), and the distance between the hypocenter of the targeted seismic event and the location of the farthest accelerometer located in NR is 740 m. Thus, it takes 240 ms for S-waves arising from the fault-slip to reach the observation point. Hence, 300 ms is considered sufficient to investigate peak acceleration at the point. Broken lines in Fig. 13 show peak acceleration in z-direction, which is the same direction as that for the uniaxial accelerometers. As shown in the figures, the dynamic analysis is performed for 0% and 2% damping (D = 0 and 0.02), using the indicated result that appropriate damping appears to fall between 0% and 2% as shown in Fig. 12. Thus the upper and lower boundaries of the peak acceleration can be estimated by conducting dynamic analyses for the two damping

conditions. It is found from Fig. 13 that the observed peak acceleration falls between the upper and lower peak acceleration obtained from the dynamic analysis except the two locations 500 m and 740 m away from the source region. Peak acceleration in the z-direction obtained from the dynamic analysis is smaller than the observed values for all locations. As a reason why the peak acceleration in the z-direction obtained from the dynamic analysis is smaller than the observed values, it is postulated that initial motion of the simulated fault-slip differs from that of the actual seismic event, that is to say, fault-slip took place in directions different from the assumed failure plane during the seismic event in reality. Although the simulated peak ground acceleration in the z-direction is smaller compared to the observed values, these results indicate that reasonable estimation of peak ground acceleration in the vicinity of active mining areas could be performed on the basis of the estimated values since the simulated peak acceleration still falls into the reasonable values. This validation then suggests that the estimated properties and assumptions made for conducting the back analysis are reasonably acceptable. 6. Simulating peak particle velocity in active mining areas Lastly, propagation of seismic waves arising from the simulated seismic event is shown with respect to ground velocity. Although McGarr (2002) proposed the equation to estimate near-field peak particle velocity, the equation is valid only in the vicinity of source regions of seismic events provided that rupture velocity is known. In the present study, fault properties are calibrated in terms of moment magnitude and peak particle velocity, and the calibrated values are then validated with peak ground acceleration. Thus, more reliable ground motion than that estimated from an empirical chart (Wang and Cai, 2015) could be estimated by performing dynamic analysis with the calibrated values. Fig. 14 depicts ground motion induced by propagation of seismic waves arising from the seismic event simulated with the calibrated values. The illustrated area is from 750E to 1050E and from − 150 N to −450 N, which corresponds with the scale shown in Fig. 1. An active mining area, where stopes are extracted and mine openings are located, is also shown in the Fig. 14. The simulated ground motion is considered the upper case when 0% damping is taken for the analysis. 0% damping is unrealistic, but it provides the upper limit of PPV experienced in the active mining area. Instead of applying 0% damping, it is possible to make more precise estimation of the damping coefficient from Fig. 12. Fig. 14 shows seismic waves propagating outward through the rockmasses from the source region and reaching the active mining

170

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

Fig. 14. Ground motion induced by seismic event simulated with calibrated values.

area as time elapses. PPV is widely employed as an indicator to estimate the severity of ground motion and damage to rockmass caused by seismic waves (Diederichs, 2007; Hedley, 1992) because PPV is directly related to kinetic energy induced by seismic waves arising from fault-slip. Dynamic support systems installed in underground mines are intended to absorb the kinetic energy transmitted by seismic waves for the purpose of preventing rock ejection and falls. Hence, estimating PPV is crucial for the design of dynamic support systems. In this regard, Fig. 14 provides valuable information, i.e., based on particle velocity shown in this figure, the intensity of ground motion (PPV) that could occur in the vicinity of the active mining area can be estimated. The estimated PPV can be used in conjunction with empirical charts showing the relation of PPV with damage to rockmass and seismic hazard (Brinkmann, 1987; Diederichs, 2007; Hedley, 1992). In this modeled case, particle velocities of rockmasses within the active mining area reach approximately 0.6 m/s. Calculations of this value with simulations would support the design of appropriate ground support systems. In addition to the simulation of seismic wave propagation, this study sheds light on the estimation of the other seismic source parameters, such as near-field energy and rupture velocity, which were not the focus of previous studies on back analysis of mining-induced fault-slip (Alber and Fritschen, 2011; Hofmann and Scheepers, 2011). As the present study simulates fault-slip in dynamic conditions, not only static stress drop but also near-field and far-field seismically released energy can be estimated from the simulated fault-slip. Using those parameters, it would be possible to estimate seismic efficiency and rupture velocity (McGarr and Fletcher, 2001). The study is recommended as future work.

356,469 seismic events, with magnitude-time-history analyses and time-between-event rate charts. In a 3D numerical model, the shear stiffness of the rockmass within the seismic source region is kept the same as the host rock, and a ubiquitous joint model with Barton's shear strength model is applied to simulate fault-slip. The rest of the shear zone outside the source region of the seismic event is treated as transversely isotropic hence the shear stiffness is decreased with respect to the orientations of assumed planes of weakness. Back analysis is conducted by comparing moment magnitude of the seismic event and peak particle velocity recorded by strong ground motion sensor with results obtained from the dynamic analysis. As a result, surface roughness and a dynamic friction angle of failure plane are estimated with the back analysis. In addition, the estimated properties are validated with peak ground acceleration recorded by seven uniaxial accelerometers. These results indicate that the simulation technique developed and employed for the back analysis could be useful in order to determine the physical and mechanical properties of causative shear/fracture zones in underground mines, which are generally difficult to obtain. By conducting the back analysis for other seismic events within the shear/fracture zones, the physical and mechanical properties of rockmasses within the zones could be statistically estimated. Thus, it would be possible to predict seismic events that could occur in response to planned mining activities. Although the present study is conducted for Mine-A, the developed methodology can be applicable to other mines suffering from seismic events that occur especially in shear/fracture zones. Acknowledgement

7. Conclusion In the present study, back analysis of a fault-slip seismic event that took place within a shear or fracture zone in Mine-A is carried out. The source mechanism of the targeted seismic event is deduced from

This work is financially supported by a grant from the Natural Science and Engineering Research Council of Canada (NSERC) (Grant no. RGPIN 41603-10) - Discovery Grant Program, and McGill University (MEDA Fellowship Program); the authors are grateful for their support.

A. Sainoki, H.S. Mitri / Journal of Applied Geophysics 134 (2016) 159–171

Reference ABAQUS, 2003. ABAQUS online documentation: Ver 6.4-1. Dassault Systemes, France. Alber, M., Fritschen, R., 2011. Rock mechanical analysis of a M1 = 4.0 seismic event induced by mining in the Saar District, Germany. Geophys. J. Int. 186, 359–372. Alejano, L.R., Gonzalez, J., Muralha, J., 2012. Comparison of different techniques of tilt testing and basic friction angle variability assessment. Rock Mech. Rock. Eng. 45 (6), 1023–1035. Bandis, S., Lumsden, A.C., Barton, N.R., 1983. Fundamentals of rock joint deformation. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 20 (6), 249–268. Barton, N., 1973. Review of a New Shear-strength Criterion for Rock Joints Engineering Geology 7. pp. 287–332. Barton, N., Choubey, V., 1977. The shear strength of rock joints in theory and practice. Rock Mech. 10, 1–54. Barton, N., Bandis, S., Bakhtar, K., 1985. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 22 (3), 121–140. Bewick, R.P., Valley, B., Runnalls, S., Whitney, J., Krynicki, Y., 2009. Global approach to managing deep mining hazard. The 3rd CANUS Rock Mechanics Symposium, Toronto. Brady, B.H.G., Brown, E.T., 1993. Rock Mechanics: for Underground Mining. Springer Science & Business Media. Brinkmann, J.R., 1987. Separating shock wave and gas expansion breakage mechanism. 2nd International Symposium on Rock Fragmentation by Blasting, pp. 6–15. Bruggeman, D.A.G., 1935. Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen. Ann. Phys. 416 (7), 636–664. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12 (2), 81–97. Castro, L.A.M., Carter, T.G., Lightfoot, N., 2009. Investigating factors influencing fault-slip in seismically active structures. 3rd CANUS Rock Mechanics Symposium, Toronto. Diederichs, M.S., 2007. The 2003 Canadian Geotechnical Colloquium: mechanistic interpretation and practical application of damage and spalling prediction criteria for deep tunnelling. Can. Geotech. J. 44 (9), 1082–1116. Dieterich, J.H., Kilgore, B., 1996. Implications of fault constitutive properties for earthquake prediction. National Academy of Science, U.S.A., pp. 3787–3794 ESG, 2011. WaveVis. ESG solutions, Kingston, ON, Canada. Hedley, D.G.F., 1992. Rockburst Handbook for Ontario Hard Rock Mines. Ontario Mining Association. Hofmann, G.F., Scheepers, L.J., 2011. Simulating fault slip areas of mining induced seismic tremors using static boundary element numerical modelling. Min. Technol. 120 (1), 53–64. Hudyma, M., Beneteau, D.L., 2012. Time - a key seismic source parameter. In: Hawkes, C. (Ed.), 21st Canadian Rock Mechanics Symposium ROCKENG12 - ROCK ENGINEERING FOR NATURAL RESOURCES. The Canadian Rock Mechanics Association, pp. 311–318 Edited by. Itasca, 2009. FLAC3D - Fast Lagrangian Analysis of Continua. Itasca Consulting Group Inc., U.S.A. Ivins, E.R., Lyzenga, G.A., 1986. Stress patterns in an interplane shear zone: an effective anisotropic model and implications for the transverse ranges, California. Phil. Trans. R. Soc. A 318, 285–347. Kanamori, H., 2001. Energy budget of earthquakes and seismic efficiency. Int. Geogr. 76, 293–305. Ladanyi, B., Archambault, G., 1969. Simulation of shear behaviour of a jointed rock mass. In: Somerton, W.H. (Ed.), The 11th US Symposium on Rock Mechanics. NY, Society of Mining Engineers, 1970, The University of California, Berkeley, California, pp. 105–125 Edited by. Li, B., Dai, F., Xu, N.W., Sha, C., Zhu, Y.G., Shi, Y.L., He, G., 2014. Microseismic monitoring system establishment and its engineering applications of deep-buried underground powerhouse. Chin. J. Rock Mech. Eng. 33, 3375–3384. Lysmer, J., Kuhlemeyer, R.L., 1969. Finite dynamic model for infinite media. J. Eng. Mech. 95 (EM4), 859–877. Maxwell, S.C., Urbancic, T.I., 2001. The role of passive microseismic monitoring in the instrumented oil field. Lead. Edge 20 (6), 636–639.

171

McGarr, A., 1991. Observations constraining near-source ground motion estimated from locally recorded seismograms. J. Geophys. Res. 96 (B10), 16495–16508. McGarr, A., 1994. Some comparisons between mining-induced and laboratory earthquakes. PAGEOPH 142, 467–489. McGarr, A., 1999. On relating apparent stress to the stress causing earthquake fault slip. J. Geophys. Res. 104 (B2), 3003–3011. McGarr, A., 2002. Control of strong ground motion of minig-induced earthquakes by the strength of the seismogenic rock mass. J. South. Afr. Inst. Min. Metall. 225–229. McGarr, A., Fletcher, J.B., 2001. A method for mapping apparent stress and energy radiation applied to the 1994 Northridge earthquake fault zone - revisited. Geophys. Res. Lett. 28 (18), 3529–3532. Morrow, C.A., Shi, L.Q., Byerlee, J.D., 1982. Strain hardening and strength of clay-rich fault gouges. J. Geophys. Res. 87, 6771–6780. Ortlepp, W.D., 2000. Observation of mining-induced faults in an intact rock mass at depth. Int. J. Rock Mech. Min. Sci. 37, 423–426. Ortlepp, W.D., Stacey, T.R., 1994. Rockburst mechanisms in tunnels and shafts. Tunn. Undergr. Space Technol. 9, 59–65. Potvin, Y., Jarufe, J., Wesseloo, J., 2010. Interpretation of seismic data and numerical modelling of fault reactivation at El Teniente, Reservas Norte sector. Min. Technol. 119 (3), 175–181. Ryder, J.A., 1988. Excess shear stress in the assessment of geologically hazardous situations. J. South. Afr. Inst. Min. Metall. 88 (1), 27–39. Saharan, M.R., 2004. Dynamic modelling of rock fracturing by destress blasting. Mining & Materials Engineering. McGill University, Montreal, QC, Canada. Sainoki, A., Mitri, H.S., 2014a. Dynamic modelling of fault slip with Barton's shear strength model. Int. J. Rock Mech. Min. Sci. 67, 155–163. http://dx.doi.org/10.1016/j.ijrmms. 2013.12.023. Sainoki, A., Mitri, H.S., 2014b. Effect of slip-weakening distance on selected seismic source parameters of mining-induced fault-slip. Int. J. Rock Mech. Min. Sci. 73, 115–122. http://dx.doi.org/10.1016/j.ijrmms.2014.09.019. Sainoki, A., Mitri, H.S., 2014c. Methodology for the interpretaiton of fault-slip seismicity in a weak shear zone. J. Appl. Geophys. 110, 126–134. http://dx.doi.org/10.1016/j. jappgeo.2014.09.007. Sainoki, A., Mitri, H.S., 2014d. Numerical simulation of rock mass vibrations induced by nearby production blast. Can. Geotech. J. 51 (11), 1253–1262. http://dx.doi.org/10. 1139/cgj-2013-0480. Sainoki, A., Mitri, H.S., Feng, X.-T., Hudson, J.A., Tan, F., 2013. Comparative study of shear strength models for fault-slip analysis. Sinorock. Taylor & Francis Group, London, Shanghai, pp. 551–556 Edited by. Shnorhokian, S., Mitri, H.S., Thibodeau, D., 2014. Numerical simulation of pre-mining stress field in a heterogeneous rockmass. Int. J. Rock Mech. Min. Sci. 66, 13–18. Shnorhokian, S., Mitri, H.S., Thibodeau, D., Moreau-Verlaan, L., 2012. Analysis of the influence of the mining sequence on a remnant pillar with FLAC3D. 21st Canadian Rock Mechanics Symposium. RockEng 12, Edmonton, AB, Canada, pp. 29–41. Shnorhokian, S., Mitri, H.S., Thibodeau, D., Moreau-Verlaan, L., 2013. Optimizing stope sequence in a diminishing ore pillar: a case study. World Mining Congress, Montreal. Sjöberg, J., Perman, F., Quinteiro, C., Malmgren, L., Dahner-Lindkvist, C., Boskovic, M., 2012. Numerical analysis of alternative mining sequences to minimize potential for fault slip rockbursting. Min. Technol. 121 (4), 226–235. Trifu, C.-I., Shumila, V., 2010. Microseismic monitoring of a controlled collapse in Field II at Ocnele Mari, Romania. Pure Appl. Geophys. 167, 27–42. Urbancic, T.I., Trifu, C., 2000. Recent advances in seismic monitoring technology at Canadian mines. J. Appl. Geophys. 45 (4), 225–237. Wang, X., Cai, M., 2015. Influence of wavelength-to-excavation span ratio on ground motion around deep underground excavation. Tunn. Undergr. Space Technol. 49, 438–453. Yuana, F., Prakash, V., 2008. Use of a modified torsional Kolsky bar to study frictional slip resistance in rock-analog materials at coseismic slip rates. Int. J. Solids Struct. 45, 4247–4263.