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.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).

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.

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).

where, a is a scaling exponent dictating the ratio between smaller and larger fractures; l_max and l_min 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).

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, l_max = 10 m, l_min = 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°.

Discrete fracture network in a fault damage zone

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 l_min 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.

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.

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; K_n and K_s 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.

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 (E_intact) 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/m³ is applied to the model.

As for the fracture stiffnesses K_n and K_s 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, K_n and K_s 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.

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.

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.

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.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.

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.

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).

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.

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.

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 × 10⁷ N m and 3.5 × 10⁸ N m in 20th and 25th stages, respectively. This discrepancy corroborates the result shown in Fig. 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).

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.

where, E_R denotes seismically radiated energy; τ, τ₁, and τ₂ 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.

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.

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 M_c. For the determination, the present study employ the maximum curvature method (Wiemer and Wyss 2000). Although this method is known to underestimate M_c 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 M_c is interpreted as the minimum earthquake that can be adequately modelled in the model. That is, M_c 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 M_c 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.

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.

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.

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.

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 m³. 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 m³. 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 m³. 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.

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.

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.

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.