The behavior of rock mass damage and failure processes are affected significantly by the stress path the rock is experienced (Bai and Tu 2019; Cai 2008b; Diederichs et al. 2004; Kaiser et al. 2001; Wang et al. 2019; Zhang et al. 2018). Consequently, for a given condition, the evolution of damage and failure around an excavation depends on the stress path. The monotonically increased loading path, which usually applies in standard laboratory tests (such as uniaxial or triaxial compression), may not represent in situ conditions. Underground excavations not only produce stress loading in the rock mass but also unloading near the excavation boundaries, where rock fracturing and failure usually occur (Cong et al. 2020; Zheng et al. 2017). Also, we cannot relate the laboratory results to the in situ conditions because laboratory results are obtained from a different stress state and history from the in situ environment.

Stress path tests or simulations have previously been used to help understand excavation-induced rock damage. In numerical simulations, Cai (2008b) reported the deformation and yielding zone near the excavation surface, which is closely related to the unloading rate that was used to simulate the excavation. Using Discrete Element Method (DEM) models, Li et al. (2014) investigated how the unloading rate and initial stress release path affect the deformation, damage, and kinetic energy evolution during a circular tunnel excavation. In a detailed three-dimensional finite-element study, Eberhardt (2001) demonstrated the evolution of near-field stress paths during the progressive advancement of a tunnel face. In the study, the emphasis was placed on the rotation of the principal stress and its effect on microfracture initiation and propagation, brittle fracture damage, and rock strength degradation. In field simulations, Diederichs et al. (2004) employed an excavation-induced stress path to explain how the near-face stress rotation affects the in situ strength degradation and the damage development. Field measured and predicted stress paths were used by Kaiser et al. (2001) to explain the nature of different failure modes in an underground mine.

To capture the realistic rock mass responses, the stress path should be correctly represented in a test or model. Several studies have conducted stress path tests to simulate in situ problems. Ivars et al. (2011) illustrated an example of predicting block caving fragmentation by applying the mining-induced stress path to SRM (synthetic rock mass) models. Zhang et al. (2016) presented the results of mechanical and permeation properties on coal samples under the mining-induced stress path and found the stress-dependent permeability agrees well with field measurements. Still, the results are quite different from those of the conventional triaxial tests under monotonic loading regimes. Recently, a series of experiments were conducted using an in situ 3D stress path (Nasseri et al. 2016) to understand the strength, deformation, and seismic response during the Posiva’s Spalling Experiment (POSE). The experiment makes it possible to accurately obtain the rock responses (e.g., strength, deformation, damage process, etc.) using intensive monitoring systems in a laboratory environment but under an in situ stress state. The experiments on the homogenous pegmatite have shown the calculated P-wave velocities during the stress path tests were in reasonable agreement with those measured in the field (Young et al. 2020). The research also highlighted σ₂ plays a vital role in the damage of the veined rocks under some circumstances, indicating the characteristics of mechanical behaviors concluded in these tests cannot be replicated with conventional triaxial testing where σ₂ = σ₃.

This paper presents a further investigation by applying the in situ 3D stress path to a sample to simulate the roof damage process of the Mine-by tunnel. With this approach, the damage process during the excavation of the tunnel can be examined directly in the laboratory experiment and numerical simulation, something generally difficult to obtain in conventional triaxial experiments and numerical simulations. The comparison with the field observations on the damage development, seismicity, and fracture mechanism verifies the validity of this method.

2.1 Stress path at the crown of the Mine-by tunnel

The Mine-by tunnel was designed to investigate damage development around a circular tunnel in a highly stressed rock mass. It has been analyzed and well documented by numerous researchers (Martin et al. 1997; Martin 1997; Read et al. 1998). The stress path for a point of 1.0 cm from the tunnel crown as the tunnel develops was calculated using a FLAC3D model (Bai et al. 2019). Eight tunnel face locations in relation to the monitoring point, A (7 m), B (0.26 m), C (0.013 m), D (− 0.013 m), E (− 0.038 m), F (− 0.063 m), G (− 0.138 m) and H (− 10 m), are used to mark the character of the stress path. The conceptual model is shown in Fig. 1. The simulated stress was also verified with field measurements (Martin 1997) at radial distances of 2.65 m and 3.34 m (Bai et al. 2019). Velocity measurements during the experiment were equated with the in situ results at similar stress states (Tibbo 2018) and showed the stress path experiment was simulating this parameter during the tunnel excavation process.

Stress paths (σ₁–σ₃ and σ₂–σ₃) at the location of 1.0 cm from the tunnel crown with the eight key tunnel face locations. The crack initiation strength (σ_ci) (Martin 1997) and crack damage strength (σ_cd) (Martin 1997) using fitted Hoek–Brown (H–B) failure criterion are also given in the figure. A conceptual model to obtain the stress path is shown. The notched failure profile of the Mine-by tunnel is illustrated

The stress path (σ₁–σ_3, and σ₂–σ₃ plots) is shown in Fig. 1 according to the eight reference locations. There are some points of interest in this plot. σ₁ exceeds the crack initiation strength (σ_ci) at point B, which means cracks initiate in front of the tunnel face. However, at the area immediately in front of the tunnel face (point C), both σ₁ and σ₂ are below σ_ci. At the location directly behind the tunnel face (point D), stress levels for both σ₁ and σ₂ are above the σ_ci but much smaller than the crack damage strength (σ_cd) proposed by Martin (1997). After that, the minimum principal stress approaches tensile (point F) and reaches the level of the tensile strength of the LdB granite (9.3 MPa) (Martin 1997). Then σ₃ increases and becomes compressive at a location 0.16 m behind the face (point H). Meanwhile, σ₁ increases and reaches σ_cd.

Since tensile stress cannot be applied to the sample in the true triaxial test (TTT), the minimum magnitude of σ₃ was set to 2 MPa (e.g., σ₃ at points B, F, G, and H is 2 MPa) to maintain a functional coupling between AE transducers and loading platens (Nasseri et al. 2014). The applied stress path in the TTT is summarized in Table 1. The same stress path was used in the numerical modeling.

2.2 Numerical modeling

Compared with continuum methods, the Discrete Element Method (DEM) explicitly simulates crack initiating and propagating from microscale to macroscale without applying complex constitutive laws. This advantage makes it widely used to observe cracking and failure in brittle rocks. The Parallel Bond Model (PBM) (Potyondy and Cundall 2004) has been successfully used to simulate failure processing of brittle rocks both in laboratory tests and field practice (Potyondy 2015; Wang and Tonon 2009; Xu et al. 2020) and the effects of heterogeneity and anisotropy on rock behaviors (Duan and Kwok 2016; Peng et al. 2017; Tian and Yang 2017). However, three inherent problems are encountered when using this model (Cho et al. 2007; Potyondy and Cundall 2004): (1) the unrealistic low ratio of uniaxial compressive strength (UCS) to tensile strength, (2) the excessively low friction angle, and (3) the linear strength envelope. The unrealistic microstructure features used in the PBM model are recognized as one of the primary reasons causing these intrinsic problems (Cho et al. 2007; Wu and Xu 2016). The irregular bonded grains in the real rock are simplified mimicked as spheres (or circular disks in 2D) in the BPM, which may not adequately reflect the geometry-dependent interlocking that does not fully mobilize the frictional strength between particles (Bahrani et al. 2011; Potyondy 2010). In BPM, contacts between particles are bonded across their entire length. After bond breaks, they are removed and lose the ability to resist relative rotation (Potyondy 2012), which underestimates the internal frictional strength (Wu and Xu 2016).

To overcome the inherent problems, the Flat-Joint Model (FJM) was proposed (Potyondy 2012). In the model, well-connected intergranular interlocks between rigid grains are produced by assigning an appropriate initial surface gap. The FJM provides the behavior of a finite-size, linear elastic, and either bonded or frictional interface. The interface is discretized into elements and can be deformable, breakable, and partially damaged. The interface evolves from a fully bonded state to a fully unbonded and frictional state. Partially damaged interface, even a fully broken interface, continues to resist relative rotation because the flat joint is not removed (Potyondy 2012). Besides, pre-existing cracks and open-pores can also be presented in the microstructure of a flat-jointed material by introducing a unbonded slit and unbonded gapped contacts, respectively (Potyondy 2017). Details of the model formulation, behaviors, force–displacement laws, and failure criterion can be found in Itasca Consulting Group Inc (2017) and Potyondy (2012, 2017).

In PFC, each bond contact breakage forms a microcrack and releases part of the stored strain energy in the form of a seismic wave (Al-Busaidi et al. 2005; Hazzard and Young 2002, 2004). The released strain causes the change of the contact forces around the breakage, which can be used to calculate the moment tensor of the acoustic emission (AE) (Hazzard and Young 2004). In reality, slipping on existing joints or faults would also produce seismicity (e.g., natural earthquakes); however, slipping on failed contacts (i.e., crack slipping) is not regarded as an AE event in this numerical study. Individual AEs should be clustered into large ones to obtain realistic magnitude distributions. AEs that fall within a specific space window and their duration (time window, i.e., computational time step in the PFC code) overlap are considered part of a larger event. The space window (SW) and time window (TW) for clustering affect the size of events and the distribution of event magnitudes (Bai et al. 2021; Khazaei et al. 2015b). Therefore, SW and TW should be properly calibrated to obtain realistic results.

The 3D stress path (listed in Table 1) was applied to a cubic LdB granite sample, which was drilled from the nearby location of the Mine-by tunnel, to simulate the roof damage development. The experiment has been previously described in Bai et al. (2019). In this study, the seismic data were further processed to obtain the complete AE time-series by revisiting the continuously recording waveforms (Young et al. 2013), which had not been reported in the previous study (Bai et al. 2019).

In the former analysis (Bai et al. 2019), 308 AE were recorded based on the triggered system. In this study, the continuous waveforms identified 4027 AEs. Due to the complex velocity structure during the experiment, 904 AEs were successfully located within the specimen. A plot of the magnitude distribution for the located AEs is shown in Fig. 2. The b-value was estimated using the maximum likelihood method (Aki 1965) and considering the binning correction (Amorèse et al. 2010):

where m_c is the magnitude of completeness, which is defined as the lowest magnitude above which all seismic events within a certain region are reliably detected (Woessner and Wiemer 2005); 〈m〉 is the mean magnitude with m ≥ m_c. Δm = 0.1 is the magnitude bin size. A correct estimate of m_c is crucial to calculate the b-value. We assess m_c by using the maximum curvature method (Wiemer and Wyss 2000). In this experimental test, the slope of the frequency-magnitude plot (b-value) is 1.38.

Frequency-magnitude distribution (FMD) plot for AE events showing a b-value of 1.38. Filled circles are the cumulative number of AE, and hollow squares are the non-cumulative number of AE. The dashed vertical line and straight blue line indicate the magnitude of completeness and the estimated Gutenberg–Richter relation (Gutenberg and Richter 1944), respectively. Note the low estimation quality due to the gradually curved FMD

In the first 3 stages (see Table 1), few AE events occurred. 423 AEs occurred in the top left part of the sample in Stage 4, as shown in Fig. 3a. Additional 209 AEs formed at Stage 5 (see Fig. 3b). In Stages 6, 7, and 8, a total of 266 AEs were recorded within the sample. Most of these events are located in the same region as that formed in Stages 4 and 5, indicating microcracks propagated along the existed cracks that formed in the former stages. The AE cloud determined a damage plane using the best fitting method described in Ingraham et al. (2013). The damage plane orientates 34° from the σ₁ direction (approximately paralleling to σ₂, 4°), which is consistent with the previous analysis (Bai et al. 2019).

Cumulative AE locations in a Stage 4, b Stage 5, c Stage 8, and d the projection in σ₁–σ₃ plane. Events are colored according to time (cold color for early occurrence and hot color for later occurrence)

The stress path (Table 1) used in the TTT produced a potential damage plane, but the sample did not fail, which is not consistent with the field observations. Monitoring in the field showed that microseismic events initiated 0.5–1.0 m ahead of the tunnel face (stress path B and C or Stage 2 and 3, see Table 1). Microcracks propagated and coalesced into small-scale flaking and then resulted in thin slabs in the area approximately 0.5–1.0 m behind the tunnel face. The slabs progressed in an unstable manner with larger slabs and finally produced stable notches in the roof and floor (Martin et al. 1997; Martin 1997), as shown in Fig. 1. It is likely that the modified σ₃ used in the TTT contributed to the disagreement between field observation and the experiment (Bai et al. 2019). Field observations show that damages around the tunnel were susceptible to confining stress (Read and Martin 1996). However, tensile stress and low compressive stress were replaced by a confining pressure of 2 MPa (see Table 1), which should be provided to make efficient contact between AE sensors and the sample. Therefore, the modified stress path listed in Table 1 cannot fully reflect the real manner of the field conditions, although similar stress paths (not include the tensile part) were widely used to explain the notch formation of the Mine-by tunnel, e.g., Diederichs et al. (2004), Li et al. (2017) and Martin (1997).

The numerical model is calibrated to present the macro properties of the LdB granite. Numerical calibrations were conducted on cylinder samples with a size of 50 mm in diameter and 125 mm in length. The sample consists of approx. 9400 particles. Table 2 summarizes the calibrated microparameters for the LdB granite. The calibration results are illustrated in Fig. 4. The simulated strength envelope is nonlinear and is consistent with the H–B failure criterion (Read et al. 1998). The numerical results are also summarized in Table 3 and compared with laboratory results and the H–B model. The results reveal that the errors between the laboratory properties and the simulated values are within ± 10%.

Comparison of the simulated results with the H-B failure criterion given by Read et al. (1998)

The UCS test produced shear failure with the primary fracture traveling through the sample with an angle of approx. 60°, as shown in Fig. 5a. Most of the larger AE events concentrated along the failure plane. The crack initiation (σ_ci) agrees well with the onset of microcracks and AE events, as shown in Fig. 5b. This phenomenon is consistent with laboratory results (Eberhardt et al. 1998; Martin 1993). The axial stress sharply decreased after the peak value with a quick increase of AE events (Fig. 5c). Therefore, the FJM used in the study successfully reproduced the brittle failure of the LdB granite.

Summary of numerical results of the UCS test. a Particle displacement vector and AE distribution and magnitude, b Stress–strain curve, volumetric-strain curve, and AE number-strain curve accompanying the σ_ci and σ_cd. σ_ci is defined as the onset of dilatancy − increase in volume relative to elastic changes, and σ_cd is defined as the point of the volumetric strain reversal

The space window (SW) of 2 times averaged particle diameter is suggested (Hazzard and Young 2004; Yoon et al. 2015) to calculate the AE moment tensor in the PFC model. The time window (TW) was determined by calculating the b-value using different TWs, as shown in Fig. 6. It can be seen the 30 and 40 TW produced a huge event with a magnitude of − 4.5 and a source dimension of 5.1 cm, which clustered most of the microcracks near the failure plane, so there are few events along the failure plane. There is also a large magnitude gap between the largest and the second-largest events (see Fig. 6). Such a large event is unrealistic and is quite different from laboratory experiments. Laboratory tests showed a large number of AEs always spatially aligned on primary fractures (Lockner et al. 1991; Sellers et al. 2003). Events in 10 TW were smaller than − 5.7 and therefore generated a larger b-value of 1.56. The 20 TW produced more realistic AE behaviors (see Fig. 5), and the b-value agrees well with laboratory results which are near unit for hard rocks (Lei et al. 2016; Lockner et al. 1991; Sellers et al. 2003). This time window was employed in the simulation.

Effect of TW on b-values in the UCS test. a 10 timestep, b 20 timestep, c 30 timestep and d 40 timestep. Note the low estimation quality due to the gradually curved FMD

5.1 Microcrack development

Overall, numerical simulation generated the same AE performance as the laboratory results, as shown in Fig. 7. Most AE events occurred in Stages 4 and 5. Events started to increase in Stage 4 but had a sharp increase between Stage 4 and 5. This agrees well with the crack initiation (σ_ci) and the crack damage (σ_cd) thresholds of the LdB granite. At the middle of Stage 4, both σ₁ and σ₂ exceeded σ_ci, but much smaller than σ_cd because of the high confining stress, although σ₁ reached its maximum value, as shown in Fig. 7b. The stress state produced hundreds of normally distributed microcracks (see Fig. 8) with moment magnitude below − 7.0, consistent with many laboratory experiments that stable microcracks development occurs when σ₁ exceeds σ_ci (Brace et al. 1966; Martin and Chandler 1993). In Stage 5, σ₁ exceeded σ_cd because of the quickly decreasing of σ₃. At this stage, microcracks grew unstable (Turichshev and Hadjigeorgiou 2016), causing a sharp increase of AEs, as shown in Fig. 7a. Some cracks clustered into large AE events, as can be seen by the difference between the AE number and the microcrack number (Fig. 7a). In this stage, microcracks close to each other started to coalesce but did not form macrocracks. Therefore, AE events remained normally distribution in the sample, as illustrated in Fig. 8b. In the next three stages, few AEs occurred since σ₁ decreased to a low level.

Comparison of numerical simulation and the laboratory experiment. a Numerical AE and microcracks at different stages, b Principal stresses used in the test accompanying σ_ci (Martin 1997) and σ_cd (Read et al. 1998), c AE rate and numbers at different stages in the laboratory experiment

Numerical simulation of cumulative AE locations at a Stage 4, b Stage 5, and c Stage 8 and d the corresponding source mechanisms

The scattered distributed AEs in Fig. 8c disagreed with the laboratory result, which formed a potential damage plane initiating from an edge of the sample (see Fig. 3). This difference was considered to result from the end effect. In the laboratory experiment, a 4 mm chamfer was machined onto all edges of the sample to meet the edge-sealing purposes (Young et al. 2013). Therefore, loading platens were slightly smaller than the rock specimen’s ends, so stress concentration occurred at the edges. The elastic mismatch between the sample and plates and the friction between them also leads to a non-uniform distribution of stress near the ends (Pan et al. 2012; Xu et al. 2017). These shortcomings caused a shear damage plane initiated from the edge and associated with the concentrated distribution of the AE events in the sample (Fig. 3c). While in the modeling, the perfect loading condition (i.e., no end and corner effects) was considered, so microcracks scattered distributed in the homogeneous rock specimen when the stress exceeds the strength of the particle contacts.

Compared to the laboratory experiment, AE events in the PFC3D modeling have higher moment magnitudes. Instant breakage of bonds in the model, as opposed to gradual bond, decays in nature are known as the primary factor overestimating the moment magnitudes (Khazaei et al. 2016). Khazaei et al. (2015a) suggested that the numerical AE energies are 6.6 orders of magnitude larger than the real specimens in UCS tests when using the Gutenberg–Richter formula to calculate event magnitudes.

During the laboratory test, it is not possible to record all the events due to practical limitations. It should record many more events with magnitudes smaller than − 7.0 based on the frequency-magnitude relationship shown in Fig. 2. Events with magnitudes lower than − 7.0 have a weak signal to background noise and were filtered by the monitoring system (Lockner 1993). While PFC model can record every single event, these two facts caused the numerical modeling to record many more events than the laboratory experiment.

Martin (1997) pointed out that a deviatoric stress threshold of 70–75 MPa is a reliable indicator for predicting the onset of damage near the surface of the tunnel where the confining pressure is low. However, our experiment and modeling results do not match this criterion. Deviatoric stress in Stages 6–8 is more than 100 MPa (see Fig. 1 and Table 1), but few microcracks (or AEs) occurred in these stages. In fact, damages in the tunnel crown are progressive. Damage produced in Stages 4 and 5 degrades the rock strength [deemed as cohesion loss (Martin 1997)], which promotes the subsequent damage development under the deviatoric stress criterion. Damage development will further cause strength degradation and thus damage propagation. This processing repeats and ultimately forms a notch in the roof. While in our experiment and modeling, because of the inconsistent stress path, the amount of damage in Stages 4 and 5 did not reach the level that could produce sufficient strength degradation to promote damage development. Since the LdB granite strength is sensitive to the amount of induced damage (Martin 1997), not enough damage may impose damage development in the following stages. If this progressive strength reduction processing was introduced to the PFC model, the notch formation might successfully reproduce (Potyondy and Cundall 1998).

5.2 Damage mechanism

Source mechanisms in Fig. 8d show that most of the events are distributed between − 0.6 < K < 0.6 in the Hudson diagram. These events were regarded as tensile cracks but containing a significant component of compensated linear vector dipole (CLVD). The interpretation of microseismic events in the field also indicated that the cracking mechanism was rich in tension (Cai et al. 1998).

In the simulation, 17608 microcracks were generated. The associated stereonet is illustrated in Fig. 9, where the two dominant orientations are visible. Most microcracks are roughly parallel with σ₁ and have dip direction within 15° with respect to σ₂. This result agrees with the major tensile failure of the microcracks. Tensile cracks usually open perpendicular to σ₃ and propagate in the σ₁–σ₂ plane.

The lower hemisphere of 17608 microcrack orientations (equal area projection, Fisher concentrations). The directions of the principal stresses are indicated

FMD plot suggested a b-value of 1.63, as shown in Fig. 10. This result is slightly higher than the laboratory experiment result, which produced a b-value of 1.38 (see Fig. 2). The potential damage plane in the experiment should contribute to the disagreement since events clustering near the damaged plane usually produced larger events, i.e., lower b-value.

Frequency-magnitude relationship plot for AE events in simulation showing a b-value of 1.63. Note the low estimation quality due to the gradually curved FMD

5.3 Comparison with field monitoring

AE events obtained from laboratory and numerical tests were converted to locations in relation to the tunnel face based on the relationship between stress states and distance to the tunnel face (see Table 1). Field monitoring of MSs occurring in the top-notch region in rounds 6–11 (Collins 1997) was also provided to compare with the laboratory and numerical results, as shown in Fig. 11. Overall, laboratory and numerical tests produced the same seismicity process as that monitoring in the field. In the experiment, moment magnitudes for all the AEs are not available since only a sub-set of events are well located. Herein, the logarithm of the averaged maximum magnitude of the waveforms from an event is used to show the relative magnitude. 84.8% AE events with magnitudes > 0.0 occurred near the tunnel face. 58 events were recorded in the laboratory experiment with a distance > 1.0 m behind the tunnel face. There are 203 events (5%) occurring 8.0–10.0 m behind the tunnel face. Field monitoring and numerical simulation also obtained a few events concentrated in the region around 10 m behind the face. It is believed that the recovery of σ₁ and decreasing of σ₂ contribute to the reactivity of the events (see Fig. 7).

The distance of events from the active tunnel face versus magnitudes for a Laboratory experiment, b Numerical model, and c Field monitoring (Collins 1997). The field events in c occurred in the upper notch region of rounds 6–11 volumes along the Mine-by tunnel, where the lithology is predominantly granite (modified from Collins (1997)). The sub-plots in a and b show the close-up near the location of 0.0 m

Experiment and simulation show that events started at the tunnel face (although several AEs occurred within 0.4 m ahead of the face), while field monitoring highlighted MSs occurring 0.7–1.4 m ahead of the active tunnel face (Collins and Young 2000; Martin 1997). It is suggested that stress rotation plays one of the significant roles in cracks occurring in front of the tunnel face (Diederichs et al. 2004; Martin 1993, 1997). Both σ₁ and σ₃ were found to rotate near the tunnel face with an angle of approximately 15°–35° (Bai et al. 2019). Pre-existing microcracks in the rock could take advantage of stress rotations to extend in their favorable directions (Eberhardt 2001). However, only the change of stress magnitudes was considered in the laboratory and numerical simulation.

A number of seismic events happened in the region 0–2.0 m behind the tunnel face (see Fig. 11c), while these events were absent in the laboratory and numerical tests in the range 0.5–2.0 m. The different stress paths used in the tests attributed a large part to this difference. The elastic model shows that magnitude of 8.7 MPa tensile stress occurred at the crown immediately behind the tunnel face (see Fig. 1 and Table 1), which reaches the level of the tensile strength of the granite and should produce crack localization and spalling. That means the rock sample should fail if such tensile stress was applied in the laboratory and numerical test. Interpretation of microseismic events also showed that cracking mechanism was rich in tension (Cai et al. 1998), and thus tensile stress dominated the failure mechanism and the fracture size. However, this tensile stress was replaced by a 2 MPa compressive stress due to the limitation of the test facility. A detailed explanation of the roles of stress rotation and tensile stress, causing the disagreements with field observations can be found in Bai et al. (2019).

Notch formation in the crown is a progressive process; after the damage initiating in front of the face, microcracks dilated, developed, slabbed, and eventually formed the stable notch (Martin 1997) in the region 0–1.5 m behind the tunnel face. This progressive process generated events continually (see Fig. 11c). Some of the interspersed nature of the large and small source radius events in the field suggest that small events are occurring first in intact rock, followed by larger events that may be connecting previously formed cracks (Collins 1997). While in our laboratory and numerical tests, preferred cracks did not occur near the face due to the inconsistent stress state (tensile stress and stress rotation). Therefore, the subsequent crack development, coalescence, and the associated seismicity were absent, and thus no macro-fractures were obtained.

5.4 Influence of intermediate principal stress

Both laboratory tests and numerical simulations show that σ₂ plays a role in rock strength, damage process, and failure mechanism (Cai 2008a; Li et al. 2015; Liu et al. 2020). However, the deviatoric stress criterion for predicting the notch scale does not reflect any information of σ₂ (Martin 1997), i.e., in situ cracking starts when the maximum deviatoric stress (σ₁–σ₃) exceeds approximate 70 MPa. In the numerical simulation, σ₂ was gradually decreased to the same level of σ₃ after Stage 5 to check whether σ₂ has any effect on rock damage under the specific stress path.

As predicted, both microcracks and AEs gradually increased with σ₂ decreasing, but they presented quickly growing when σ₂ is lower than 20 MPa, as shown in Fig. 12a. Overall, the increased events were uniformly distributed in the sample, as shown in Fig. 12b, which means the increased microcracks did not coalesce into macro fractures, although many small events clustered into larger events in the sample. Note AE events result from contact breakages, and slipping on existing cracks (i.e., failed contacts) is not regarded as an event in this numerical study. The increased AEs produced a b-value of 1.55, which is lower than the original model before σ₂ decreasing (see Fig. 10), indicating more clusters of microcracks.

a Increased microcracks, AE events and their magnitudes with σ₂ decreasing, b Distribution of the increased AEs, c Lower hemisphere of orientations of the microcracks (equal area projection, Fisher concentrations)

Most of the increased microcracks failed in tension and had a dip angle within 15° with respect to the σ₁ direction, as shown in Fig. 12c, consistent with the original state (see Fig. 9) before σ₂ decreasing. However, the dip direction of the microcracks is absent within 30° in relation to σ₂. This is reasonable because most of the contacts within this direction had failed in the original state, as shown in Fig. 9.

Failure mechanism affected by σ₂ was observed in the laboratory test on hard rocks (Li et al. 2015). They found rock failure mode changed from slabbing to shear when σ₂ decreased to a critical value of 20 MPa. Our simulation also illustrates a similar mechanism. As the σ₂ unloaded, microcracks are found to incline at about 45° and 135° along the unloading direction, as the four dominant orientations are shown in the lower hemisphere figure (Fig. 12c). Because as σ₂ decreased, the orientation suppression by σ₂ was gradually weakening (Browning et al. 2017), and thus tensile microcracks turned near failed contacts where stress concentration usually occurs. However, these microcracks do not coalesce; therefore, macro shear failure was not observed in the sample. If the unloading stopped at the 20 MPa, and σ₁ was increased to fail the sample, shear failure would be observed as illustrated in Li et al. (2015).

Detailed analysis of the microcrack development process shows that most of the generated microcracks preferred in the two orientations when σ₂ was higher than a threshold of about 20 MPa, as shown in Fig. 13. Below this threshold, failed contacts with all dip directions quickly developed. Preferred failure orientation for contacts gradually disappeared as σ₂ decreased to the same level as σ₃, i.e., cracks generated are only parallel to the maximum principal stress and have no preferred orientation relative to σ₂ and σ₃.

a Orientation distribution of the microcracks during σ₂ unloading, b Development of microcracks for specific orientations during σ₂ unloading. The numbered orientations (e.g., #1) in plot b are also indicated in plot a

The increased AE events during σ₂ unloading highlighted that the σ₂ magnitude also affects the micro failure mechanism, which could cause different fracture mechanisms in the field. The simulation shows that lower σ₂ may produce shear failure in the roof of the Mine-by tunnel, which would differ from the observed results that tensile cracking is the dominant mechanism during the notch processing (Cai et al. 1998) at the specific stress state. Also, lower σ₂ may produce larger damage zones in the roof since σ₂ unloading promotes microcracks development. Therefore, the deviatoric stress criterion for predicting the damage region in the roof (Martin 1997) also depends on σ₂, i.e., σ₂ variance could affect the accuracy of the deviatoric stress criterion.

We have conducted a true triaxial experiment and numerical simulation on an LdB granite sample under a 3D field stress path that comes from the crown of the Mine-by tunnel. AE characteristics were investigated in the tests. In general, the numerical simulation produced the same seismicity process and AE magnitude distributions (b-value). AE events obtained from laboratory and numerical tests were converted to locations with respect to the tunnel face and were compared to the field monitoring events occurring in the top-notch region of the Mine-by tunnel. Overall, laboratory and numerical tests produced the same AE performance as monitoring in the field.

σ₂ unloading after Stage 5 reveals that microcracks continued to develop, and crack planes rotated in the direction within 30 ± 15° to σ₂ before it decreases to a level of 20 MPa. The concentration of the increased microcracks indicated shear damage planes may be forming during σ₂ unloading, although tensile cracks dominated in this process. Generated cracks had no preferred orientation when σ₂ was below 20 MPa. This mechanism indicated that σ₂ may play a role in damage development and should be considered when using the deviatoric stress criterion for predicting the damage region of underground excavations.

In the laboratory test, a potential damage plane that initiated from a sample edge was produced. In contrast, the micro-cracks produced in the numerical model do not indicate a possible damage plane. This disagreement may result from the stress concentration near the sample edges because the size of the loading platens is smaller than that of the rock specimen (Bai et al. 2019). Therefore, in future studies, we will conduct numerical simulations with the same boundary conditions to better reproduce the laboratory study. In this study, the tensile stress that appears near the tunnel face (see Table 1) was replaced by a compressional stress due to the limitation of the true triaxial device, which results in the different distributions of the seismic events in relation to the tunnel face (see Fig. 11). Although it is impossible to apply tensile load in the true triaxial tests, it is feasible to employ tensile boundary conditions in the numerical model. In future simulations, we will use the true stress path (including the tensile stress) to better reproduce the field observations and to better uncover the underlying mechanisms.