Owned by: China Association for Science and Technology
Sponsored by: China Coal Society
Published by: Springer Nature
About issue

The International Journal of Coal Science & Technology is a peer-reviewed open access journal. It focuses on key topics of coal scientific research and mining development, serving as a forum for scientists to present research findings and discuss challenging issues.


Coverage includes original research articles, new developments, case studies and critical reviews in all aspects of scientific and engineering research on coal, coal utilizations and coal mining. Among the broad topics receiving attention are coal geology, geochemistry, geophysics, mineralogy, and petrology; coal mining theory, technology and engineering; coal processing, utilization and conversion; coal mining environment and reclamation and related aspects.


The International Journal of Coal Science & Technology is published with China Coal Society, who also cover the publication costs so authors do not need to pay an article-processing charge.


The journal operates a single-blind peer-review system, where the reviewers are aware of the names and affiliations of the authors, but the reviewer reports provided to authors are anonymous.


  • A forum for new research findings, case studies and discussion of important challenges in coal science and mining development

  • Offers an international perspective on coal geology, coal mining, technology and engineering, coal processing, utilization and conversion, coal mining environment and reclamation and more

  • Published with the China Coal Society

Show More
Editors-in-Chief
Suping Peng, Shimin Liu
Managing Editor
Wanjie Wang
Associate Editors
Bo Hyun Kim, Dongjie Xue, Pedram Roghanchi, Wu Xiao, Zhiqiang Wu
Publishing model
Open Access. Learn about publishing OA with us
Home > Volumes and issues > Volume 10, issue 1

Numerical modelling of spatially and temporally distributed on-fault induced seismicity: implication for seismic hazards

Research Article

Open Access

Published: 27 January 2023

0 Accesses

International Journal of Coal Science & Technology Volume 10, article number 4, (2023)

Abstract

Induced seismicity is strongly related to various engineering projects that cause anthropogenic in-situ stress change at a great depth. Hence, there is a need to estimate and mitigate the associated risks. In the past, various simulation methods have been developed and applied to induced seismicity analysis, but there is still a fundamental difference between simulation results and field observations in terms of the spatial distribution of seismic events and its frequency. The present study aims to develop a method to simulate spatially distributed on-fault seismicity whilst reproducing a complex stress state in the fault zone. Hence, an equivalent continuum model is constructed, based on a discrete fracture network within a fault damage zone, by employing the crack tensor theory. A fault core is simulated at the center of the model as a discontinuous plane. Using the model, a heterogeneous stress state with stress anomalies in the fault zone is first simulated by applying tractions on the model outer boundaries. Subsequently, the effective normal stress on the fault plane is decreased in a stepwise manner to induce slip. The simulation result is validated in terms of the b-value and other seismic source parameters, hence demonstrating that the model can reproduce spatially and temporally distributed on-fault seismicity. Further analysis on the parameters shows the variation of frequency-magnitude distribution before the occurrence of large seismic events. This variation is found to be consistent with field observations, thus suggesting the potential use of this simulation method in evaluating the risk for seismic hazards in various engineering projects.

1.Introduction

Induced seismicity is the instantaneous failure of rock mass and/or pre-existing rock discontinuity that releases stress waves because of anthropogenic in-situ stress disturbance in deep underground. For underground mines, induced seismicity can be a cause for rockbursts and could inflict devastating damage to underground facilities (Ortlepp and Stacey 1994; Ortlepp 2000; Kaiser and Cai 2012; White and Whyatt 1999). Owing to the increasing demand for exploiting mineral resource and energy at great depths, the importance of estimating and mitigating the risks for induced seismicity has been increasing around the world, especially for countries where very deep underground mines are being developed. Hence, it is indispensable to gain a better understanding of the mechanism of induced seismicity and to develop a method to predict and control its occurrence.

To elucidate the mechanism of induced seismicity and examine its relation with mining activity, seismic monitoring has been performed in a number of deep underground mines (Ma et al. 2016, 2018; Urban et al. 2016; Kozłowska et al. 2016; Beer et al. 2017; Liu et al. 2019; Wang et al. 2019). Such studies have shown that the hypocenter of induced seismicity is not necessarily concentrated in the proximity of regions undergoing in-situ stress disturbance due to anthropogenic activities. It has been then recognized that geological structures have a large influence on the occurrence of induced seismicity. For instance, Mgumbwa et al. (2017) shows, based on the result of seismic monitoring in a deep underground mine, that seismic activity becomes intense while forming spatially distributed clusters. Smith (2018) also indicates that there is no clear relation between the spatial distribution of mining-induced seismic events and the location of active mining areas; it is then demonstrated that seismically active regions are situated near and/or within pre-existing geological structures.

Furthermore, field observations of induced seismicity indicate the variation in the spatial distribution and source parameters of induced seismicity prior to the occurrence of intense seismic events (Mutke et al. 2016; Zhang et al. 2014; He et al. 2019; Cai et al. 2018; Jun et al. 2019). Specifically, this includes a decrease in the b-value in the Gutenberg formula, an increase in high-energy events, and an abnormal increase in apparent stress. Hence, many studies highlight the importance of microseismic monitoring and suggest the use of seismic source parameters as precursory information for a large seismic event.

As discussed, it has been recognized that the presence of geological structures is the dominant cause of spatially distributed seismicity taking place away from active mining areas and that characterizing the source parameters helps in evaluating the risk for seismic hazards. It should be noted, however, that mining-induced stress change is insignificant in such regions, implying that there are pre-existing stress anomalies that can be a potential source for induced seismicity. In this regard, Sainoki et al. (2016) investigate the effect of meter-scale stiffness heterogeneity due to a fracture network in a rock mass on the generation of stress concentration. It was then demonstrated that locally stressed regions are generated in the fracture network, leading to the conclusion that such meter-scale stress concentration could be a cause for severe seismic activity often observed in geological structures.

Importantly, although the effect of fracture network on the occurrence of induced and triggered seismicity has been clarified, such studies were performed by incorporating a fracture network into a cube-shaped model with an edge length of 5 m constructed in the framework of the discrete element method. Thus, the method cannot be directly applied to a mine-wide scale numerical model. To overcome this limitation, Sainoki et al. (2021) develop a methodology to reproduce the complex, heterogeneous stress state on a mine-wide scale by constructing an equivalent continuum model considering the elasticity of macro fractures, based on the crack tensor theory developed by Oda (1986).

As described above, a significant advancement has been made in terms of numerical modelling of the complex, heterogeneous in-situ stress state within geological structures. It is then postulated that such a heterogeneous stress state with stress anomalies produces spatially distributed seismic events in underground mines and affects their characteristics. However, most of previous studies simulate only a single or a few seismic events on a fault whilst assuming a stress state simply varying with depth in their numerical models (Kong et al. 2019; Sainoki and Mitri 2014a; Li et al. 2021; Wei et al. 2020). Characteristics of the simulated seismic activity are inconsistent with those measured in the field, especially for the number of seismic events, thus making it difficult to compute and evaluate imperative parameters, such as the b-value.

Considering the discussion above, the present study aims to simulate spatially and temporally distributed seismic events on a fault plane whilst considering a heterogeneous stress state in the fault zone. The simulation result is validated and characterized in terms of the source parameters, such as the b-value, seismically released energy, and seismic moment.

2.Discrete fracture network in a fault damage zone

2.1 Fracture characteristics

A fault zone is generally classified into fault core and fault damage zone (Berg and Skar 2005; Lin et al. 2007; Faulkner et al. 2010) as shown in Fig. 1, although the fault damage zone can be further classified into inner and outer damage zones (Joussineau and Aydin 2007; Choi et al. 2016; Gurrik 2016; Ostermeijer et al. 2020). The fault core is a narrow zone with a width of less than tens of centimeters, where shear displacement is accumulated with the evolution of the fault, whilst the fault damage zone is composed of rock masses with dense fractures and extends in the direction perpendicular to the fault core. Field observations indicate that various factors, such as the shear displacement in a fault core, mechanical properties of the protolith and ambient stress conditions, affect the generation and distribution characteristics of fractures in a fault damage zone (Molli et al. 2010; Mitchell et al. 2011; Faulkner et al. 2011).

Fig. 1
figure 1

Characteristics of fractures in a fault zone

Importantly, although there are various contributing factors, field observations consistently show that the density of fractures in a fault damage zone decays with the distance from the fault core and its characteristic can be approximated with a power law function (Johri et al. 2014; Ostermeijer et al. 2020; Wu et al. 2019). Regarding the coefficient of the power law function related to the maximum fracture density near a fault core, the value significantly varies, depending on faults. Although some field observations report a near-core fracture density more than 100 per meter (Ostermeijer et al. 2020), most of field observations show lesser fracture densities (Ju et al. 2014).

The present study defines meter-scale fractures as macro fracture as shown in Fig. 1a. Importantly, the density of such macro fractures is assumed to follow the fracture distribution characteristics in a fault damage zone described above, whilst for meso and micro fractures their densities are presumed to be constant in the damage zone. This is attributed to the fact that there are not sufficient field observations to characterize populations of such fractures. Then, for the macro fractures, a discrete fracture network is statistically generated to construct an equivalent continuum model, whilst the effect of meso and micro fractures is implicitly considered in calculating equivalent, anisotropic compliance matrices by decreasing the elastic modulus of the host rock as described in Sect. 3. It is to be noted that the implication of meso and micro fractures on the frequency-magnitude distribution of seismic events simulated in the numerical model is discussed in later section.

2.2 Generation of discrete fracture network

To generate a fracture network in a fault damage zone, the following parameters need to be determined: fracture density, fracture length, and fracture orientation. For the fracture density, according to previous studies (Agosta et al. 2010), the following power law function is employed.

$$P_{10} = \alpha d^{\beta }$$
(1)

where, α represents a coefficient related to the maximum fracture density near a fault core; β is a scaling factor; and d is the distance from the fault core. Regarding fracture length, the following cumulative density function is employed (Gutierrez and Youn 2015).

$$F\left( l \right) = \frac{{l_{\min }^{1 - a} - l^{1 - a} }}{{l_{\min }^{1 - a} - l_{\max }^{1 - a} }}$$
(2)

where, a is a scaling exponent dictating the ratio between smaller and larger fractures; lmax and lmin denote the maximum and minimum fracture length simulated, respectively. Lastly, as for fracture orientation, the Fisher function expressed with the following equation is used (Fisher 1996).

$$f\left( \theta \right) = \frac{{\kappa \sin \theta \exp \left( {\kappa \cos \theta } \right)}}{{\exp \left( \kappa \right) - \exp \left( { - \kappa } \right)}}$$
(3)

where, θ denotes the angular deviation from the mean vector; κ is a dispersion factor that determines the degree of the angular deviation. Using these formulae, disk-shaped fractures in a fault damage zone are statistically generated.

The present study uses the following parameters for Eqs. (1) to (3): α = 30, β = 0.8, lmax = 10 m, lmin = 2 m, a = 2.5, and, κ = 40. The verification of these parameters and their influence on the stress state in a fault damage zone are described in the study (Sainoki et al. 2021). Figure 2 illustrates the discrete fracture network (DFN) generated, which is composed of approximately twenty million disk-shaped fractures that statistically follow Eqs. (2) and (3). The fault plane (fault core) is located at the center of the model and runs in the y-direction with a dip angle of 90°.

Fig. 2
figure 2

Discrete fracture network in a fault damage zone

3.Numerical model construction

3.1 Equivalent continuum model

Figure 3 illustrates the procedure to calculate an equivalent continuum model representing the behaviour of the fault zone. As shown in the figure, the first step is to convert of the Young’s modulus of intact rock to the elastic modulus of the rock mass including micro and meso fractures. For the conversion, geological strength index (GSI) proposed by Hoek (2007) is employed as the first approximation to consider the effect of meso and micro fractures. In the present study, micro and meso fractures are defined as those with lengths of two meters or less in the fault damage zone since lmin is set at 2 m in Eq. (2). In this regard, an assumption is made that such micro and meso fractures are randomly oriented and distributed in the fault damage zone, thus producing isotropic elastic behaviour. M in Fig. 3 denotes the equivalent elastic compliance tensor of the rock mass with meso and micro fractures.

Fig. 3
figure 3

Procedure to construct an equivalent continuum model

It is to be noted that although GSI is typically employed to consider macro fractures including those with lengths greater than several meters, GSI can be applied to such meso fractures. In fact, Cai et al. (2004) calculated a rock mass block volume to determine a GSI value whilst assuming an average joint length of 2 m.

Furthermore, this method to produce an equivalent continuum model representing a fault damage zone was first proposed by Sainoki et al. (2021) and then verified based on comparison with discrete element method explicitly modelling the fracture network. According to the study, the relative ratio of the elastic modulus of the rock mass including micro and meso fractures to the stiffness of macro fractures is one of the most dominant factors determining the degree of stress heterogeneity and the magnitude of stress anomalies in the fault damage zone. Hence, this influences the characteristics of on-fault seismic activity. In a later section, the validity of the model input parameters is discussed based on the frequency-magnitude distribution and source parameters of seismic events on the fault plane.

After calculating M, the crack tensor theory developed by Oda (1986) is employed to compute the elastic compliance tensor C equivalent to the elasticity of macro fractures with the following equations.

$${F}_{ij} = \frac{1}{V}\frac{\pi }{4}D^{3} {n}_{i} {n}_{j}$$
(4)
$${F}_{ijkl} = \frac{1}{V}\frac{\pi }{4}D^{3} {n}_{i} {n}_{j} {n}_{k} {n}_{l}$$
(5)
$$\begin{aligned} {C}_{{ijkl}} = & \sum\nolimits_{1}^{{\text{NCR}}} {\left[ {\frac{1}{{K_{\text{n}} D}} - \frac{1}{{K_{\text{s}} D}}} \right]} {F}_{{ijkl}} \\ & + \frac{1}{{4K_{\text{s}} D}}\left( {\delta _{{ik}} {F}_{{jl}} + \delta _{{jk}} {F}_{{il}} + \delta _{{il}} {F}_{{jk}} + \delta _{{jl}} {F}_{{ik}} } \right) \\ \end{aligned}$$
(6)
$${T}_{ijkl} = {C}_{ijkl} + {M}_{ijkl}$$
(7)

where, F and n are the crack tensor and the unit normal vector of the fracture intersecting with the zone, respectively; C, M, and T represent the anisotropic compliance tensor equivalent to the elasticity of fractures, the isotropic compliance tensor of the rock mass, and the anisotropic compliance tensor of fault damage zone, respectively; NCR is the total number of fractures intersecting with the zone; Kn and Ks are shear and normal stiffnesses of the fracture, respectively; δ is the Kronecker’s delta; and D denotes the diameter of the fracture.

Eventually, the equivalent compliance tensor T considering all the effect of macro, meso, and micro fractures on the elastic matrix of the rock mass is input to a zone in the continuum model composed of parallelepiped zones illustrated in Fig. 3. It is to be noted that characteristics of macro fractures, i.e., density, orientation, and length, spatially vary as depicted in Fig. 2. Accordingly, the equivalent elastic compliance tensor T spatially varies in the model as well, i.e., for each zone, T is computed based on the crack tensor theory and the fracture characteristics in the local region, thus generating a heterogeneous continuum model and rendering its mechanical behaviour anisotropic.

3.2 Numerical model constructed

Figure 4 illustrates the numerical model constructed following the procedure shown in Fig. 3. In the mode, the red-colored zones correspond to regions where the macro fractures modelled in Fig. 2 are located, so that the anisotropic compliance tensor is calculated and input to the zone. On the other hand, the grey-colored zones are in regions without macro fractures in Fig. 2. As such, only the compliance tensor M is input to the zones. The fault core running parallel to the y-direction with a dip angle of 90° is located at the center of the model, corresponding to the DFN in Fig. 2. The fault plane is explicitly modelled with interface elements as shown in Fig. 4b.

Fig. 4
figure 4

Numerical model showing zones with anisotropic or isotropic compliance tensor

3.3 Model input parameters for the rock mass

It is well known that fractures are intensely developed in a fault zone composed of low-porosity rocks. Accordingly, as the protolith of the fault zone, granite is assumed. The elastic modulus of the rock mass composed of granite is then calculated with the assumptions that the Young’s of the intact granite (Eintact) ranges from 60 to 70 GPa and that its GSI ranges from 60 to 70, considering the effect of meso and micro fractures. Poisson’s ratio ranges from 0.17 to 0.30 (Li et al. 1999; Malek et al. 2009), according to which, 0.2 is assumed. As the density of granite, 2600 kg/m3 is applied to the model.

As for the fracture stiffnesses Kn and Ks in Eq. (7), the fracture normal and shear stiffnesses are set at 30 GPa and 10 GPa, respectively, as joint normal-to-shear stiffness ratio greater than 1.0 is frequently measured and/or assumed (Matsuki et al. 2009). For more details, the shear stiffnesses is determined based on the study (Bandis et al. 1983), assuming a discontinuity with a unit length, as the fracture stiffness is inversely proportional to a fracture diameter in the crack tensor theory. The normal-to-shear stiffness ratio was assumed to be 3 according to the study (Bandis et al. 1983). It is noteworthy that the assumed fracture stiffnesses fall within a reasonable range as per previous studies (Eftekhari et al. 2014; Konstantinovskaya et al. 2014) including granite and carbonate rocks. The stress dependency of the joint stiffness is not considered in the present study.

3.4 Model input parameters for the interface element representing the fault core

Following the same normal-to-shear stiffness ratio, Kn and Ks for the fault core are assumed to be 1.0 GPa/m and 0.3 GPa/m, respectively. The discrepancy between the stiffness of fractures and that of the fault core is attributed to their scale difference, and the stiffness of the fault assumed falls within a reasonable range used in simulating the mechanical behaviour of faults in underground mines (Sainoki and Mitri 2014b).

The friction angle of the fault is assumed to be 30°, which is considered to be a reasonable value according to previous studies (Barton 1976), and the cohesion is assumed to be zero. Then, the Coulomb’s shear strength model with a slip-weakening law is employed. That is, when the shear stress acting on the fault exceeds the shear strength determined by the Coulomb’s shear strength model, the friction angle is instantaneously decreased from 30° to 20°. It is to be noted that the slip-weakening distance is set to zero. This is because the slip-weakening distance is affected by a number of factors and is difficult to be determined (Marone and Kilgore 1993; Cocco et al. 2009). The validity of the slip-weakening law and the input parameters is discussed based on seismic source parameters later. When the proposed method is applied to an actual case study, the calibration of the slip-weakening distance would be required to simulate seismic events comparable to field observations.

4.Simulation of seismic events on the fault plane

As the first step, an initial stress state is simulated by applying regional stresses to the model external boundaries. It is to be noted that the direction of gravity is rotated at 20° in anticlockwise direction to simulate a fault dipping at 70°. The reason why such a vertical fault is simulated in Fig. 4 is that a complicated procedure needs to be taken to generate a DFN for dipping faults. As shown in Fig. 5a, the magnitude of stress applied to the left boundary differs from that of the right boundary because of the rotated regional stress. The stresses applied to the model boundary correspond to those at a depth of 2000 m and a horizontal-to-vertical stress ratio of 0.7. As for the bottom boundary, the displacement in all the directions is fixed at the four vertices as shown in Fig. 5b to prevent the rigid body movement due to the unbalanced force resulting from the difference in tractions between the left and right boundaries. Except for the four vertices, the displacement on the bottom boundary is constrained only in the z-direction.

Fig. 5
figure 5

Boundary conditions to simulate an initial stress state: a Stresses applied to the model outer boundaries b Displacement constraint condition on the bottom boundary

After simulating the initial stress state, the effective normal stress on the fault plane is gradually decreased in a stepwise manner in the circular region on the fault as shown in Fig. 6a. Figure 6b shows initial pore pressure distribution on the fault that varies with the depth. The bottom and top boundaries of the model are situated 2000 m and 1900 m below the ground surface, respectively. As shown in Fig. 6c, in each analysis stage, the effective normal stress is decreased by 0.5 MPa. Cumulatively, 17.5 MPa of the effective normal stress is decreased. This method assumes that on-fault induced seismicity is typically caused by unclamping of faults, i.e., the effective normal stress reduction due to rock mass excavation, fluid injection, etc.

Fig. 6
figure 6

Effective normal stress reduction to induced slip on the fault plane: a Circular region on the fault, where the effective normal stress is decreased b Initial pore pressure distribution on the fault c Effective normal stress reduction in each analysis stage

5.Results and discussion

5.1 Stress state within the fault damage zone surrounding the fault core

To acquire an in-depth understanding of the mechanical behaviour of the fault plane, it is important to investigate the stress state in the surrounding rock mass. Figure 7 shows the maximum and minimum stress distribution in the surrounding rock mass on vertical and horizontal sections. As can be seen, the heterogeneous stress state with stress anomalies has been successfully simulated, because of the heterogeneity of the anisotropic compliance matrix. Local regions with stress concentration correspond to zones with a relatively large stiffness due to a low fracture density.

Fig. 7
figure 7

Heterogeneous stress distribution within the fault damage zone: a Contour of the maximum compressive stress b Contour of the minimum compressive stress

More specifically, the maximum compressive stress in Fig. 7a locally exceeds 100 MPa, which is approximately twice as large as the stress level expected at the depth, i.e., 2000 m below the ground surface. It is found from Fig. 7a that the magnitude of the maximum compressive stress in the vicinity of the fault core tends to be lower than the surrounding region away from the core on the horizontal section. In contrast, the frequency of highly stressed zones increases in the strip-shaped region parallel to the y-direction on the horizontal section. This is attributed to the large fracture density near the core, thus reducing the rock mass stiffness of the near-fault region.

It is noted that this method to simulate the complex stress state in a fault damage zone was developed by Sainoki et al. (2021) to study seismic activity taking place in geological structures away from active mining areas. According to the study, the degree of stress concentration and heterogeneity is influenced by the ratio of the fracture stiffness substituted to Eq. (6) to the rock mass stiffness excluding the influence of macro fractures, i.e., M in Fig. 3. If the fracture stiffness were infinitely large, the anisotropic compliance tensor T would be equal to M, thus resulting in a homogeneous stress state in the damage zone. In other words, the lower the fracture stiffness, the more heterogeneous stress state is generated in the model. Thus, this result is dependent on the ratio of the host rock stiffness to the fracture stiffness. More detailed discussions on the frequency distribution of the maximum and minimum stresses in the model as well as the influence of all the input parameters are found in the study (Sainoki et al. 2021).

It is crucial to note that this heterogeneous stress state directly affects the stress state on the fault plane shown in Fig. 4b, hence exerting an influence on the occurrence of seismicity on the fault plane. The validity of the model and input parameters is discussed in a later section based on characteristics of seismic activity simulated on the fault plane.

5.2 Shear displacement increment on the fault core

Figure 8 shows shear displacement increments on the fault plane during each analysis stage shown in Fig. 6c. It can be seen from Fig. 8a that when the effective normal stress is cumulatively decreased by 5 MPa, i.e., 10th stage, multiple slip patches dispersedly appear on the fault plane. The total stress at the center of the fault is approximately 35 MPa, corresponding to a depth of 1950 m. Considering this, 5 MPa reduction of the effective normal stress is insignificant. That is, this result implies that a relatively small stress change can activate the fault, which is consistent with the characteristic of field observations in deep underground mines. This slip events at 10th stage would be regarded as so-called triggered seismic events because of the relatively small stress change.

Fig. 8
figure 8

Shear displacement increment on the fault plane during each stage (not cumulative): a 10th stage b 15th stage c 20th stage d 25th stage

Next, Fig. 8b shows the shear displacement increment during the stage where the effective normal stress is cumulatively decreased by 7.5 MPa (15th stage) on the fault plane. As can be seen, the area undergoing shear failure further expands with the increase in the magnitude of the shear displacement, and there seem to be new slip patches produced, where they are present in Fig. 8a. With the further decrease in the effective normal stress, the slip patches coalescence to form a single, large event as shown in Fig. 8c. The figure shows that slip patches connect with each other, causing shear failure in an extensive region. Eventually, fault-slip takes place on almost the entire circular region as indicated in Fig. 8d.

5.3 Determination of seismic events based on stress drop

To further explore the characteristic of the seismic activity during the analysis stages, a seismic event is detected based on the magnitude of shear stress drop at each analysis stage. An Interface element undergoing a shear stress drop larger than a threshold value is considered a seismicity-related patch. The reason why stress drop is adopted to detect a seismic event rather than the magnitude of shear displacement is that an interface element that already underwent shear failure in a previous stage continuously undergoes shear failure in subsequent stages because the stress state is located on the failure envelope. In this case, the shear behaviour may not be related to seismicity, i.e., it does not release seismic activity and can be considered as an aseismic activity. It is well known that instantaneous stress change, i.e., brittle failure, is generally required for seismic energy to be released (Kanamori and Rivera 2006, Kanamori 2001). For this reason, the occurrence of seismicity is determined based on stress drop, and the threshold is tentatively set at 0.40 MPa. In a later section, the threshold is varied to investigate its influence on the characteristics of seismic events detected. It is to be noted that the threshold needs to be greater than 0.19 MPa. This is since the effective stress reduction of 0.50 MPa corresponds to a shear stress drop of 0.18 MPa for interfaces with the residual friction angle (20°). This means that a shear stress drop larger than 0.19 MPa is most likely associated with the transition of friction angle from the initial to the residual value, representing brittle shear failure that could produce seismic energy.

Figure 9 shows seismic events identified based on shear stress drop with a threshold of 0.4 MPa. Specifically, 22, 98, 130, and 45 seismic events are identified for 10th, 15th, 20th, and 25th stages, respectively. It is remarkable that a significant number of seismic events are identified on the fault plane. This is a clear difference between the current simulation method and previous studies that predominantly simulate a single or several seismic events at most. The complex, heterogeneous stress state of the fault damage zone results in heterogeneous shear and normal stress distributions on the fault plane, thereby yielding randomly distributed seismic events. When the stress state is homogenous, a single seismic event is simulated as shown in the study (Sainoki and Mitri 2014a).

Fig. 9
figure 9

Seismic events detected on the fault plane based on stress drop: a 10th stage b 15th stage c 20th stage d 25th stage

It is found that the number of seismic events taking place in each stage increases first and then decreases. Figure 10 summarizes the number of seismic events taking place in each analysis stage. As shown in the figure, the number of seismic events starts to increase in 10th stage, takes the maximum value in 20th and 21st stages, and then drastically decreases. Interestingly, Fig. 9 indicates that although the number of seismic events decreases in stage 25th, each event has a large slip area compared to the early stages, implying that larger events tend to take place in later stages.

Fig. 10
figure 10

Number of seismic events detected at each analysis stage while reducing the effective normal stress on the fault core

5.4 Characterization of seismic source parameters of induced and triggered seismic events

Further analysis is conducted, based on seismic source parameters, for seismic events detected in each analysis stage. Figure 11 shows seismic moment increments calculated by summing seismic moments of the events for each analysis stage. As shown, the seismic moment starts to increase gradually in 12th stage, followed by a drastic increase in 20th stage. Then, it takes the maximum value in 24th stage. Thereafter, the seismic moment starts to decrease.

Fig. 11
figure 11

Seismic moment increment at each analysis stage

It is interesting to note that the stage where the number of seismic events takes the maximum value does not coincide with that of the maximum seismic moment. As shown V number of seismic events take place in 20th stage, but the cumulative seismic moment at the stage is insignificant, compared to the maximum value. To further investigate this discrepancy, the frequency distribution of moment magnitude is shown in Fig. 11. It is found from the figure that in 15th and 20th stages, seismic events with moment magnitude ranging from −1.4 to −1.2 occur the most frequently. Compared to the stages, the frequency of seismic events with magnitudes less than −1.0 considerably decreases in 25th stage, whilst large events with magnitudes greater than 0 take place.

Figure 13 depicts the average of seismic moment for each analysis stage, highlighting the different characteristic of seismic events among the analysis stages. As can be seen, the average of seismic moment is approximately 5.0 × 107 N m and 3.5 × 108 N m in 20th and 25th stages, respectively. This discrepancy corroborates the result shown in Fig. 12.

Fig. 12
figure 12

Frequency distribution of moment magnitude at different analysis stages

In summary, in early stages, seismic events with relatively small magnitudes are caused by a small change in the effective normal stress. Then, the further decrease in the effective stress causes larger events in later stages, but the number of events decreases. Thereafter, both the frequency and magnitude decrease (Fig. 13).

Fig. 13
figure 13

Average of seismic moment for each analysis stage

Figure 14 shows the increment of seismically radiated energy for each analysis stage, calculated by summing the energy of seismic events taking place in each analysis stage. It is to be noted that seismically radiated energy is usually evaluated based on the slip rate of fault patches undergoing shear failure (McGarr and Fletcher 2001) or particle velocities of the surrounding rock mass (Rivera and Kanamori 2005). As the present study performed a quasi-static analysis, neither slip rate nor particle velocity can be accurately evaluated during the analysis. Such that, as the first approximately, seismically released energy is evaluated based on shear stress change during shear failure, using the equation below.

$$E_{\text{R}} = \frac{1}{2}\left( {\tau_{1} + \tau_{2} } \right)DA - A\int\limits_{0}^{D} {\tau \text{d}D}$$
(8)

where, ER denotes seismically radiated energy; τ, τ1, and τ2 are the shear stress of a fault patch, the shear stresses before and after the slip motion, respectively, D is a shear displacement increment during the slip event; A is an area of fault patches undergoing slip. In this simulation, the slip-weakening distance is set at zero, which is equivalent to zero fracture energy and corresponds the ideal model described in the study (Kanamori and Rivera 2006). Although this method cannot evaluate seismically radiated energy accurately, it can be nevertheless used as an index to characterize seismic events and is deemed sufficient considering the purpose of this study.

Fig. 14
figure 14

Energy released during each analysis stage

It is found from Fig. 14 that the variation of seismically radiated energy is almost the same as that of seismic moment in Fig. 11. That is, the energy starts to increase in 12th stage, shows the noticeable increase in 20th stage, takes the maximum value in 23rd to 25th stages, and then decreases. Hence, this result corroborates the discrepancy between the frequency of seismic events and its intensity.

5.5 b-value of the seismic activity on the fault core

As for natural earthquakes, its magnitude-frequency relationship is commonly described with the well-known Gutenberg–Richter (G–R) law (Gutenberg and Richter 1944), which is expressed as below.

$$N = 10^{a - bM}$$
(9)

where, N denotes the number of seismic events with magnitudes equal or greater than M; a and b are constants representing the intensity of seismicity and the relative number of strong earthquakes with respect to weak ones, respectively. As such, a large b-value indicates high frequency of earthquakes with relatively small magnitudes. Previous studies indicate that for natural earthquakes, the b-value generally falls within the range from 0.6 to 1.6 (Scholz 2015; Schorlemmer and Wiemer 2004).

When calculating b-value, it is required to determine the minimum magnitude of complete recording Mc. For the determination, the present study employ the maximum curvature method (Wiemer and Wyss 2000). Although this method is known to underestimate Mc in the case of gradually curved frequency distributions of magnitude, this method is fast and reliable (Woessner and Wiemer 2005), providing the first approximation. It should be noted that Mc is interpreted as the minimum earthquake that can be adequately modelled in the model. That is, Mc is strongly associated with the minimum fracture length of the DFN shown in Fig. 2 and the interface element size representing the fault core in the equivalent continuum model. This is because the fracture and element sizes determine the lower limit of seismic events. In other words, the minimum fault-slip area corresponds to the area of a single interface element.

Figure 15 shows frequency-magnitude distributions for different thresholds to determine seismicity-related fault patches as well as b-values calculated based on Mc estimated with the method. It is found from the figures that despite the difference in the threshold, the estimated b-values fall within a range from 0.9 to 1.0. The range is quite reasonable and comparable to that estimated based on seismic monitoring records in the field (El-Isa and Eaton 2014). Figure 15 further indicates that b-value increases with the threshold. This is reasonable since with the increase in the threshold, smaller events with a low stress drop are excluded in detecting seismicity-related fault-slip, hence increasing proportion of large seismic events.

Fig. 15
figure 15

b-values for different thresholds of stress drop to define a seismic event: a 0.2 MPa b 0.3 MPa c 0.4 MPa

5.6 Relation between shear displacement increment and stress drop and its implication

Figure 16 shows the maximum shear displacement and shear stress drop during each stage. It is worth mentioning that there is no significant discrepancy in the magnitude of stress drop between the early and late stages, compared to the maximum shear displacement. As can be seen from the figure, the maximum stress drop during each stage ranges from 1.6 to 1.8 MPa in the early stages until the effective normal stress reduction reaches 8 MPa. In the later stages, the magnitude increases to more than 3.0 MPa. Although the value of the stress drop increases with the decrease in the effective normal stress, the difference is not significant. This implies that even if the magnitude of the seismic event is relatively small, the locally intense stress drop could produce severe ground vibrations in the vicinity of the fault patch undergoing shear failure.

Fig. 16
figure 16

Maximum shear displacement increment and stress drop at each analysis stage

As demonstrated in Figs. 8 and 9, in early stages, seismic events take place in a locally stressed region due to the presence of pre-existing stress anomalies resulting from the heterogeneous stiffness distribution. Hence, the magnitude of the stress drop associated with slip-weakening behaviour may not be correlated with the magnitude of seismic events, especially for triggered seismic events taking place in early stages. This is because the magnitude of shear displacement increment during a seismic event is affected by the area undergoing slip, i.e., the larger the area becomes, the larger the magnitude of the slip event, as shown in Figs. 9 and 11. For these reasons, triggered seismic events with a relatively small magnitude could produce a large stress drop, which may inflict severe damage to a local region in the vicinity of the hypocenter.

Figure 17 depicts shear stress drop increments on the fault plane, which clearly shows that many slip patches undergo relatively large stress drops in 10th stage, which completely differs from the characteristic of the shear displacement increment shown in Fig. 8. It is also noteworthy that patches with a large stress drop are dispersedly distributed on the fault, irrespective of the progress of the analysis stage, and that locations of slip patches with a large shear displacement increment shown in Fig. 8 do not coincide with those of slip patches with a relatively large stress drop. For instance, as shown in Fig. 8d, during 25th stage, the most intense slip takes place at the center of the circular region, and the magnitude of the slip decreases towards the periphery of the region, whilst patches with a relatively large stress drop are spatially distributed in almost the entire circular region, as shown in Fig. 17d. This implies that the location of severe ground motions does not necessarily coincide with the area near the hypocenter, i.e., locally intense damage could be inflicted in the region located near slip patches undergoing a large stress drop.

Fig. 17
figure 17

Shear stress drop on the fault plane during each stage (not cumulative): a 10th stage b 15th stage c 20th stage d 25th stage

Further study would be required to quantify the intense ground motion caused by a large stress drop in a local area. To achieve this, dynamic analysis needs to be performed instead of the quasi-static analysis in the present study.

6.Validity of the simulation method and the seismic activity simulated on the fault plane

The validity of the numerical model and analysis result is examined from various points of view since this study is not based on an actual case study. First, as expected, the analysis result is dependent on the discretization of the model domain, that is, the zone size. According to the study (Sainoki et al. 2021), when the zone edge length is 0.2 m, the stress distribution in the equivalent continuum model derived from the crack tensor theory is quantitatively consistent with that in a discontinuum model constructed in the framework of discrete element method, where fractures are explicitly modelled with discontinuity planes. The comparison was made in a cubed region with dimensions of 20 m × 20 m × 20 m, thus discretizing the domain into 1,000,000 zones. In the study, further analysis was performed by extending the region to 300 m × 500 m × 300 m. In that case, the effect of the zone size was investigated whilst varying the zone volume size from 3 to 350 m3. It was then concluded that the analysis result sufficiently converges, irrespective of the shape of the zone, when the zone volume is reduced to 3 m3. In the present study, the model domain is a cube with an edge length of 100 m, which is discretized into 8,000,000 zones with an edge length of 0.5 m corresponding to a volume of 0.125 m3. Considering the result of the previous study, the zone edge length and volume should be small enough to obtain an accurate result. Admittedly, as implied in Fig. 15, seismic events with magnitudes less than −2.0 cannot be adequately simulated with the zone size. It should be noted, however, that such extremely small seismic events are not generally recorded in the field (Santis et al. 2019; Bischoff et al. 2010; Leptokaropoulos et al. 2017) due to the limitation of seismic monitoring systems. Considering this, the numerical model used in the present study meets the requirement to simulate induced seismicity on a fault plane subjected to a complex stress state due to the fracture network of fault damage zone.

The b-value representing the magnitude-frequency distribution of seismic events is known to be close to 1.0 when calculated based on the seismicity of the whole Earth (El-Isa and Eaton 2014), although it shows spatial and temporal fluctuations, affected by various factors, such as depth, stress state, and the strength of the crust (Mutke et al. 2016, Urban et al. 2016, Scholz 2015, Popandopoulos and Lukk 2014). Calculations for seismic events on the fault plane corroborate the validity of the numerical model and input parameters. As shown in Fig. 15, irrespective of the threshold, b-values computed for all the seismic events that took place during the whole analysis are close to unity, falling within a reasonable range reported in the aforementioned studies.

Another approach to examine the validity of the numerical model is to compute scaled energy defined as the ratio of seismically radiated energy to seismic moment. Over the past decades, significant efforts have been made in determining the scale energy of natural and induced seismicity (Gibowicz et al. 1991; Jost et al. 1998; Yamada et al. 2007; Venkatarangan et al. 2006; Venkatarangan and Kanamori 2004; Mayeda et al. 2005; Ide and Beroza 2001). All the seismological investigations yield similar ranges from 1.0 × 10–7 to 1.0 × 10–4. It is noteworthy that Ide and Beroza (2001) argue that previous studies missed a portion of seismic energy in their calculations due to the limitation of observable bandwidths. They concluded that if the missed energy is taken into consideration, the scaled energy takes a constant value of 3.0 × 10–5.

Figure 18 shows scaled energy computed for seismic events in each analysis stage. It is remarkable that, irrespective of the analysis stage, the computed scaled energy is close to the value (3.0 × 10–5) suggested by Ide and Beroza (2001), without any calibrations with respect to the model input parameters and boundary conditions. Hence, this result strongly suggests that the numerical model and input parameters can produce seismic events with reasonable seismic source parameters as well as appropriate frequency-magnitude distribution.

Fig. 18
figure 18

Scaled energy computed for seismic events in each analysis stage

These results and discussions validate the model discretization scheme and induced seismicity simulated in the model. Since the numerical model is not based on any case studies, the validation methods described above should be sufficient at this stage. When applying the method in this study to a case study, calibrations of input parameters would be required, e.g., the residual strength of the fault and stiffnesses of the host rock and fractures.

7.Interpretation of induced seismicity on the fault plane

The discussion in the previous section indicates that the seismic activity on the fault plane possesses characteristics similar to those of seismicity in the field. In this section, further characterization and interpretation of the seismic activity are performed. First, as discussed in Sect. 5.4 and shown in Fig. 13, the frequency-magnitude distribution varies with the decrease in the effective normal stress. The figure clearly indicates that as the effective stress decreases, the relative frequency of seismic events with small magnitudes less than 1.0 decreases, accompanied by the increase in the number of large seismic events. This variation in the relative frequency corresponds to the decrease in the b-value, i.e., although the b-value is close to 1.0 when all the seismic events are considered, it varies with the decrease in the effective normal stress.

This tendency is consistent with the variation in b-value related to the occurrence of relatively large seismic events reported from case studies for deep underground mines as well as geothermal projects (Zhang et al. 2014; Barth et al. 2011; Mutke et al. 2016), that is, as suggested by such studies, the b-value may be used as a precursor of major failure causing a large seismic event that could produce strong ground motions. In this regard, it appears that the seismic activity simulated on the fault place reproduces this characteristic. As demonstrated, the number of seismic events takes the maximum value in 20th stage, but the frequency of seismic events with magnitudes ranging from -1.4 to -1.2 is the largest, which would yield a relatively large b-value. This indicates a low risk for the occurrence of large seismic events. Thereafter, the average seismic moment drastically increases as shown in Fig. 12, whereas the number of seismic events start to decrease, thus the b-value decreasing. This may be considered as a precursor of the severe seismic activity entailing large events in 23rd to 25th stages as shown in Fig. 11. These results demonstrate a great potential regarding the applicability of the method developed in this study to the risk evaluation for intense seismicity caused by anthropogenic seismic activities.

The primary objective of the present study is to simulate seismic activity on a fault plane subjected to a complex, heterogeneous stress state within a fault damage zone and to characterize the induced seismicity. Hence, seismic events were induced by simply reducing the effective normal stress in the circular region on the fault. In the future, this method will be applied to actual cases, such as deep underground mines with a major fault and fluid injection related to geothermal energy development, where seismicity is induced by unclamping of faults due to orebody extraction or fluid injection-induced pore pressure increase. Then, the characteristic of seismic events simulated will be compared to seismic monitoring results for further validation and the calibration of model input parameters.

8.Conclusions

The present study applies the crack tensor theory to construct an equivalent continuum model discretized into millions of zones in the framework of the finite difference method, in which each zone has a different anisotropic compliance matrix corresponding to the characteristic of a local fracture network in a fault damage zone. A fault plane is simulated with interface elements at the center of the model. Using the model, a complex, heterogeneous stress state with stress anomalies is first simulated by applying stresses on the model outer boundaries. Thereafter, seismicity is induced by decreasing the effective normal stress on the fault plane by 0.5 MPa in each analysis stage.

It was then demonstrated that the method can simulate induced seismicity composed of hundreds of spatially distributed seismic events on the fault plane. The validity of the simulated seismic events was investigated in terms of the b-value and scaled energy. For both the parameters, the induced seismic activity yielded values comparable to those obtained from seismic monitoring records of natural and induced seismicity in the field. Further analysis in terms of the seismic moment and b-value of the seismic events simulated on the fault plane indicated the variation of frequency-magnitude distribution before the occurrence of large seismic events, i.e., the relative frequency of relatively small seismic events decreases before a large event takes place, hence corresponding to the decrease in the b-value. This variation is consistent with the change in the b-value observed in the field, suggesting that the potential use of this simulation method in evaluating the risk for the occurrence of intense seismicity in various engineering projects, such as deep underground mining and geothermal energy development.

References

[1] Agosta F, Alessandroni M, Antonellini M, Tondi E, Giorgioni M (2010) From fractures to flow: a field-based quantitative analysis of an outcropping carbonate reservoir. Tectonophysics 490:197–213
[2] Bandis S, Lumsden AC, Barton NR (1983) Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech 20:249–268
[3] Barth A, Wenzel F, Langenbruch C (2011) Probability of earthquake occurrence and magnitude estimation in the post shut-in phase of geothermal projects. J Seismolog 17:5–11
[4] Barton N (1976) Rock mechanics review the shear strength of rock and rock joints. Int J Rock Mech Min Sci Geomech Abstr 13:255–279
[5] Beer W, Jalbout A, Riyanto E, Ginting A, Sullivan M, Collins DS (2017) The design, optimisation, and use of the seismic system at the deep and high-stress block cave Deep Mill Level Zone mine. In: Hudyma M, Potvin Y (eds) Undergound mining technology
[6] Berg SS, Skar T (2005) Controls on damage zone asymmetry of a normal fault zone: outcrop analyses of a segment of the Moab fault, SE Utah. J Struct Geol 27:1803–1822
[7] Bischoff M, Cete A, Fritschen R, Meier T (2010) Coal mining induced seismicity in the Ruhr area, Germany. Pure Appl Geophys 167:63–75
[8] Cai M, Kaiser PK, Uno H, Tasaka Y, Minimi M (2004) Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the GSI system. Int J Rock Mech Min Sci 41:3–19
[9] Cai W, Dou L, Zhang M, Cao W, Shi JQ, Feng L (2018) A fuzzy comprehensive evaluation methodology for rock burst forecasting using microseismic monitoring. Tunn Undergr Space Technol 80:232–245
[10] Choi JH, Edwards P, Ko K, Kim YS (2016) Definition and classification of fault damage zones: a review and a new methodological approach. Earth Sci Rev 152:70–87
[11] Cocco M, Tinti E, Marone C, Piatanesi A (2009) Scaling of slip weakening distance with final slip during dynamic earthquake rupture. In: Fukuyama E (ed) Fault-zone properties and earthquake rupture dynamics. Elsevier Academic Press, Cambridge
[12] Eftekhari M, Baghbanan A, Bagherpour R (2014) The effect of fracture patterns on penetration rate of TBM in fractured rock mass using probabilistic numerical approach. Arabian J Geosci 7:5321–5331
[13] El-Isa ZH, Eaton DW (2014) Spatiotemporal variations in the b-value of earthquake magnitude–frequency distributions: classification and causes. Tectonophysics 615–616:1–11
[14] Faulkner DR, Jackson CAL, Lunn RJ, Schlische RW, Shipton ZK, Wibberley CAJ, Withjack MO (2010) A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones. J Struct Geol 32:1557–1575
[15] Faulkner DR, Mitchell TM, Jense E, Cembrano J (2011) Scaling of fault damage zones with displacement and the implications for fault growth processes. J Geophys Res 116:B05403
[16] Fisher NI (1996) Statistical analysis of circular data. Cambridge University Press, Cambridge
[17] Gibowicz SJ, Young RP, Talebi S, Rawlence DJ (1991) Source parameters of seismic events at the Underground Research Laboratory in Manitoba, Canada: scaling relations for events with moment magnitude smaller than -2. Bull Seismol Soc Am 81:1157–1182
[18] Gurrik RO (2016) Damage zone characteristics of an exhumed reservoir-caprock succession. Master, University of OSLO
[19] Gutenberg B, Richter CF (1944) Frequency of earthquakes in California. Bull Seismol Soc Am 34:185–188
[20] Gutierrez M, Youn DJ (2015) Effects of fracture distribution and length scale on the equivalent continuum elastic compliance of fractured rock masses. J Rock Mech Geotech Eng 7:626–637
[21] He S, Song DZ, Li Z, He X, Chen J, Li D, Tian X (2019) Precursor of spatio-temporal evolution law of MS and AE activities for rock burst warning in steeply inclined and extremely thick coal seams under caving mining conditions. Rock Mech Rock Eng 52:2415–2435
[22] Hoek E (2007) Practical rock engineering
[23] Ide S, Beroza G (2001) Does apparent stress vary with earthquake size? Geophys Res Lett 28:3349–3352
[24] Johri M, Dunham EM, Zoback MD, Fang Z (2014) Predicting fault damage zone by modeling dynamic rupture propagation and comparison with field observations. Solid Earth 119:1251–1272
[25] Jost M, Bübelberg T, Jost Ö, Harjes H-P (1998) Source parameters of injection-induced microearthquakes at 9 km depth at the KTB deep drilling site. Bull Seismol Soc Am 88:815–832
[26] Joussineau G, Aydin A (2007) The evolution of the damage zone with fault growth in sandstone and its multiscale characteristics. J Geophys Res 112:B12401
[27] Ju W, Hou D, Zhang B (2014) Insights into the damage zones in fault-bend folds from geomechanical models and field data. Tectonophysics 610:182–194
[28] Jun A, Jianguo N, Pengqi Q, Shang Y, Hefu S (2019) Microseismic monitoring and its precursory parameter of hard roof collapse in longwall faces: a case study. Geomech Eng 17:375–383
[29] Kaiser PK, Cai M (2012) Design of rock support system under rockburst conditions. J Rock Mech Geotech Eng 4:215–227
[30] Kanamori H (2001) Energy budget of earthquakes and seismic efficiency. Int Geophys 76:293–305
[31] Kanamori H, Rivera L (2006) Energy partitioning during an earthquake. In: Abercrombie R, McGarr A, Di Toro G, Kanamori H (eds) Earthquakes: radiated energy and the physics of faulting. Americal Geopysical Union, Washington, D.C
[32] Kong P, Jiang L, Shu J, Wang L (2019) Mining stress distribution and fault-slip behavior: a case study of fault-influenced longwall coal mining. Energies 12:2494
[33] Konstantinovskaya E, Rutqvist J, Malo M (2014) CO2 storage and potential fault instability in the St. Lawrence Lowlands sedimentary basin (Quebec, Canada): insights from coupled reservoir-geomechanical modeling. Int J Greenhouse Gas Control 22:88–110
[34] Kozłowska M, Orlecka-Sikora B, Rudziński Ł, Cielesta S, Mutke G (2016) Atypical evolution of seismicity patterns resulting from the coupled natural, human-induced and coseismic stresses in a longwall coal mining environment. Int J Rock Mech Min Sci 86:5–15
[35] Leptokaropoulos K, Staszek M, Cielesta S, Urban P, Olszewska D, Lizurek G (2017) Time-dependent seismic hazard in Bobrek coal mine, Poland, assuming different magnitude distribution estimations. Acta Geophys 65:493–505
[36] Li HB, Zhao J, Li TJ (1999) Triaxial compression tests on a granite at different strain rates and confining pressures. Int J Rock Mech Min Sci 36:1057–1063
[37] Li Z, Wang C, Shan R, Yuan H, Zhao Y, Wei Y (2021) Study on the influence of the fault dip angle on the stress evolution and slip risk of normal faults in mining. Bull Eng Geol Env 80:3537–3551
[38] Lin A, Maruyama T, Kobayashi K (2007) Tectonic implications of damage zone-related fault-fracture networks revealed in drill core through the Nojima fault, Japan. Tectonophysics 443:161–173
[39] Liu JP, Xu SD, Li YH, Lei G (2019) Analysis of rock mass stability based on mining-induced seismicity: a case study at the Hongtoushan copper mine in China. Rock Mech Rock Eng 52:265–276
[40] Ma X, Westman EC, Fahrman BP, Thibodeau D (2016) Imaging of temporal stress redistribution due to triggered seismicity at a deep nickel mine. Geomech Energy Environ 5:55–64
[41] Ma J, Dineva S, Cesca S, Heimann S (2018) Moment tensor inversion with three-dimensional sensor configuration of mining-induced seismicity (Kiruna mine, Sweden). Geophys J Int 213:2147–2160
[42] Malek F, Suorineni FT, Vasak P (2009) Geomechanics strategies for rockburst management at Vale Inco Creighton Mine. In: Diederichs M, Grasselli G (eds) ROCKENG09, Toronto
[43] Marone C, Kilgore B (1993) Scaling of the critical slip distance for seismic faulting with shear strain in fault zones. Nature 362:618–621
[44] Matsuki K, Nakama S, Sato T (2009) Estimation of regional stress by FEM for a heterogeneous rock mass with a large fault. Int J Rock Mech Min Sci 46:31–50
[45] Mayeda K, Gök R, Walter W, Hofstetter A (2005) Evidence for non-constant energy/moment scaling from coda-derived source spectra. Geophys Res Lett. https://doi.org/10.1029/2005GL022405
[46] McGarr A, Fletcher JB (2001) A method for mapping apparent stress and energy radiation applied to the 1994 Northridge earthquake fault zone—revisited. Geophys Res Lett 28:3529–3532
[47] Mgumbwa J, Page A, Human L, Dunn MJ (2017) Managing a change in rock mass response to mining at the Frog’s Leg underground mine. In: Wesseloo J (ed) Deep mining
[48] Mitchell TM, Ben-Zion Y, Shimamoto T (2011) Pulverized fault rocks and damage asymmetry along the Arima-Takatsuki Tectonic Line, Japan. Earth Planet Sci Lett 308:284–297
[49] Molli G, Cortecci G, Vaselli L, Ottria G, Cortopassi A, Dinelli E, Barbieri M (2010) Fault zone structure and fluid–rock interaction of a high angle normal fault in Carrara marble (NW Tuscany, Italy). J Struct Geol 32:1334–1348
[50] Mutke G, Pierzyna A, Baranski A (2016) b-value as a criterion for the evaluation of rockburst hazard in coal mines. In: The 3rd international symposium on mine safety science and engineering, Montreal, Canada
[51] Oda M (1986) An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour Res 22:1845–1856
[52] Ortlepp WD (2000) Observation of mining-induced faults in an intact rock mass at depth. Int J Rock Mech Min Sci 37:423–426
[53] Ortlepp WD, Stacey TR (1994) Rockburst mechanisms in tunnels and shafts. Tunn Undergr Space Technol 9:59–65
[54] Ostermeijer GA, Mitchell TM, Aben FM, Dorsey MT, Browning J, Rockwell TK, Fletcher JM, Ostermeijer F (2020) Damage zone heterogeneity on seismogenic faults in crystalline rock; a field study of the Borrego Fault, Baja California. J Struct Geol 137:104016
[55] Popandopoulos GA, Lukk AA (2014) The depth variations in the b-value of frequency-magnitude distribution of the Earthquakes in the Garm Region of Tajikistan. Phys Solid Earth 50:273–288
[56] Rivera L, Kanamori H (2005) Representations of the radiated energy in earthquakes. Geophys J Int 162:148–155
[57] Sainoki A, Mitri HS (2014a) Dynamic behaviour of minig-induced fault slip. Int J Rock Mech Min Sci 66:19–29
[58] Sainoki A, Mitri HS (2014b) Dynamic modelling of fault slip with Barton’s shear strength model. Int J Rock Mech Min Sci 67:155–163
[59] Sainoki A, Mitri HS, Yao M, Chinnasane D (2016) Discontinuum modelling approach for stress analysis at a seismic source: case study. Rock Mech Rock Eng 49:4749–4765
[60] Sainoki A, Schwartzkopff AK, Jiang L, Mitri H (2021) Numerical modeling of complex stress state in a fault damage zone and its implication on near-fault seismic activity. JGR Solid Earth. https://doi.org/10.1029/2021JB021784
[61] Santis F, Contrucci I, Kinscher J, Bernard P, Renaud V, Gunzburger Y (2019) Impact of geological heterogeneities on induced seismicity in a deep sublevel stoping mine. Pure Appl Geophys 176:697–717
[62] Scholz CH (2015) On the stress dependence of the earthquake b value. Geophys Res Lett 42:1399–1402
[63] Schorlemmer D, Wiemer S (2004) Earthquake statistics at Parkfield. J Geophys Res 109:B12307
[64] Smith SVL (2018) Retrospective analysis of mine seismicity: Glencore, Kidd mine. Master, Laurentian University
[65] Urban P, Lasocki S, Blascheck P, Nascimento AF, van Giang N, Kwiatek G (2016) Violations of Gutenberg–Richter relation in anthropogenic seismicity. Pure Appl Geophys 173:1517–1537
[66] Venkatarangan A, Kanamori H (2004) Observational constrains on the fracture energy of subduction zone earthquakes. J Geophys Res. https://doi.org/10.1029/2003JB002549
[67] Venkatarangan A, Beroza G, Ide S, Imanishi K, Ito H, Iio Y (2006) Measurements of spectral similarity for microearthquakes in western Nagano, Japan. J Geophys Res. https://doi.org/10.1029/2005JB003834
[68] Wang Z, Li X, Shang X (2019) Distribution characteristics of mining-induced seismicity revealed by 3-D ray-tracing relocation and the FCM clustering method. Rock Mech Rock Eng 53:183–197
[69] Wei CH, Zhang C, Canbulat I (2020) Numerical analysis of fault-slip behaviour in longwall mining using linear slip weakening law. Tunn Undergr Space Technol 104:103541
[70] White BG, Whyatt JK (1999) Role of fault slip on mechanisms of rock burst damage, Lucky Friday Mine, Idaho, USA. In: 2nd Southern African rock engineering symposium. implementing rock engineering knowledge. Johannesburg, S. Africa
[71] Wiemer S, Wyss M (2000) Minimum magnitude of completeness in earthquake catalogs: examples from alaska, the western United States, and Japan. Bull Seismol Soc Am 90:859–869
[72] Woessner J, Wiemer S (2005) Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty. Bull Seismol Soc Am 95:684–698
[73] Wu G, Gao L, Zhang Y, Ning C, Xie E (2019) Fracture attributes in reservoir-scale carbonate fault damage zones and implications for damage zone width and growth in the deep subsurface. J Struct Geol 118:181–193
[74] Yamada T, Mori J, Ide S, Abercrombie RE, Kawakata H, Nakatani M, Iio Y, Ogasawara H (2007) Stress drops and radiated seismic energies of microearthquakes in s South African gold mine. J Geophys Res. https://doi.org/10.1029/2006JB004553
[75] Zhang P, Yang T, Yu Q, Zhu W, Liu H, Zhou J, Zhao Y (2014) Microseismicity induced by fault activation during the fracture process of a crown pillar. Rock Mech Rock Eng 48:1673–1682

About this article

Cite this article

Sainoki, A., Schwartzkopff, A.K., Jiang, L. et al. Numerical modelling of spatially and temporally distributed on-fault induced seismicity: implication for seismic hazards.Int J Coal Sci Technol 10, 4 (2023).
  • Received

    29 March 2022

  • Revised

    29 March 2022

  • Accepted

    30 December 2022

  • DOI

    https://doi.org/10.1007/s40789-022-00560-7

  • Share this article

    Copy to clipboard

For Authors

Explore