Inversion of seepage channels based on mining-induced microseismic data

International Journal of Rock Mechanics & Mining Sciences 126 (2020) 104180

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences journal homepage: http://www.elsevier.com/locate/ijrmms

Inversion of seepage channels based on mining-induced microseismic data Yong Zhao a, b, Tianhong Yang a, b, *, Penghai Zhang a, b, **, Haiyan Xu c, Shuhong Wang b a

Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang, Liaoning, 110819, China Center of Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang, Liaoning, 110819, China c School of Information Science and Engineering, Northeastern University, Shenyang, Liaoning, 110819, China b

A R T I C L E I N F O

A B S T R A C T

Keywords: Moment tensor Fracture network Seepage channel Fluid flow

Water inrush in mines has always been a major technical problem restricting the safe production of mines with abundant water. Seepage channels formed by fractures in rock masses provide the necessary conditions for water inrush to occur. Effectively controlling the formation of fractures in rock masses is essential for preventing water inrush disasters. In response to the shortcomings of conventional water inrush monitoring methods, a set of methods for the inversion of seepage channels in rock masses based on mining-induced microseismic (MS) data was proposed in this paper. A discrete fracture network (DFN) was built based on moment tensor (MT) theory, and a method for extracting reasonable fractures associated with different failure mechanisms was proposed. To enhance the applicability of the fractures derived from MS data (MS-derived fractures), the MS-derived fractures were transformed into a graph structure by graph theory. Subsequently, the hydraulic properties of each node and edge in the graph structure were obtained based on Darcy’s law. In addition, the fluid flow time in the MSderived fractures was weighted to detect the shortest seepage channel. Finally, the method was applied to a mine that experienced a water inrush disaster. The seepage channel formation process, water pressure distribution and main seepage channel in the study area of the mine were analyzed. The results may provide references for decision-makers regarding when and where to implement water shutoff measures in mines.

1. Introduction Water inrush disasters often occur in both coal mines1 and metal mines2,3 in China. According to statistics, more than 1,000 water inrush accidents, resulting in more than 4,000 deaths, have occurred in the past twenty years in China.4,5 At present, with the increase of mining depth, water inrush disasters are increasingly seriously influenced by high in-situ stresses and high-water pressures. To effectively avoid and con­ trol water inrush disasters, seepage channels must be accurately deter­ mined, and their evolution must be monitored in real time for early warning purposes. At present, many studies on water inrush in mines are mainly based on electromagnetic methods,6 infrared detection methods,7 compre­ hensive geophysical exploration,8,9 and many useful results have been obtained. However, these methods mainly focus on the detection of high-pressure and water-rich areas, water-conducting faults and karst collapse columns, but they lack dynamic monitoring of the development

process of water inrush channels. Generally, the developing process of water inrush channels is accompanied by rock mass failure, which leads to slip, dilation or contraction along pre-existing fractures and results in the release of energy in the form of seismic waves. These seismic waves are the mechanical response of water inrush, which provide useful in­ formation for analyzing seepage and water inrush. Therefore, the microseismic (MS) monitoring can achieve the purpose of dynamically monitoring the development process of water inrush channels. Over the past ten years, MS monitoring has become a potentially useful approach for understanding water inrush behavior.3,10–12 These studies have mainly focused on the MS events location, waveform characteristics and focal parameters, which can be used to macroscopically describe chan­ nels and hydraulic mechanisms of water inrush and therefore produce an early warning. In practice applications, accurately detecting fractures in rock masses is a feasible method for understanding seepage and water inrush channels. Currently, the methods of detecting fractures are usually

* Corresponding author. Center of Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang, Liaoning, 110819, China. ** Corresponding author. Center of Rock Instability and Seismicity Research, School of Resources and Civil Engineering, Northeastern University, Shenyang, Liaoning, 110819, China. E-mail addresses: [email protected] (Y. Zhao), [email protected] (T. Yang), [email protected] (P. Zhang). https://doi.org/10.1016/j.ijrmms.2019.104180 Received 16 November 2018; Received in revised form 9 December 2019; Accepted 15 December 2019 Available online 27 December 2019 1365-1609/© 2019 Elsevier Ltd. All rights reserved.

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

divided into rock mass surface measurements and rock mass interior measurements. The seepage and water inrush channels that form due to rock mass failure often exist in the interior of the rock mass, and the method of rock mass surface measurement fails to accurately assess these channels. Although rock mass interior measurements (including borehole video and borehole radar methods) can measure internal fractures in the rock mass, they cannot measure the formation process of water seepage channels in real time. The continuous monitoring of seismicity can help illuminate unmapped fractures that are seismically reactivated by stress redistribution.13,14 We can determine the reac­ tivated fractures through the processing of MS signals and moment tensor (MT) theory. MT inversion can be used to not only obtain the failure mechanisms of the seismic source but also determine the normal and slip vectors of the fractures. In summary, this paper aims to propose a set of methods for inversing seepage channels based on mining-induced MS data. In Section 2, the set of methods is described step by step. In Section 2.1, the inversion method combined with the principle component analysis (PCA) method with a hybrid MT method is described and the detailed implementation steps are given to establish the relationship between the MS data and fractures. In Section 2.2, the method of determining a reasonable frac­ ture from MS data under the conditions of known and unknown prin­ cipal stress is proposed. Under unknown principal stress, a stress inversion method considering various failure mechanisms is given, and a MS-derived fracture as output at the same time of the stress inversion by the stress inversion method, and detailed implementation steps are given. Sections 2.1 and 2.2 provide the foundation for Section 2.3. In Section 2.3, based on the graph theory method, the method of inversion seepage channels by MS-derived DFN is proposed. With the help of Darcy’s law, the hydraulic conductivity weight of the MS-derived frac­ ture is defined, and then the method of determining the main seepage channel is given. Finally, in Section 3, the set of methods for inversing seepage channels based on mining-induced MS data is applied to the Shirengou Iron Mine, which has experienced water inrush. The MS signals in the study areas of Shirengou Iron Mine are processed, and the MS-derived DFN is constructed. Subsequently, the formation process of seepage channels, water pressure distribution and the main seepage channels are analyzed.

where u is the n-by-1 matrix representing the far-field displacement; M is the 1-by-6 matrix containing six independent MT components (m11, m22, m33, m12, m13, m23); and the Green’s function G is the n-by-6 matrix representing the properties of the medium. One of the main applications of MT inversion in mines is the analysis of rock mass failure mechanism, that is, to determine the failure mode of seismic source. Feignier and Young developed Eq. (2) for assessing the failure mechanisms of the source in geotechnical engineering by the combination of the observation of internal fractures, the principal stresses, and the MT results.27–29 R¼

(2)

where m*i represents the eigenvalues of the deviatoric MT; trðMÞ is the trace of the MT, trðMÞ ¼ M1 þ M2 þ M3 , represents the volume change of the source. M1, M2 and M3 are the eigenvalues of the MT, where the relationship is M1 � M2 � M3. The value of R ranges from 100% to 100%. The failure mechanism of MS events is to be compressive failure if R< 30%, shear failure if 30% � R � 30%, tensile failure if R > 30%. R > 0% indicates fracture opening, R ¼ 0% indicates fracture shearing slip, and R < 0% indicates fracture closure. The failure mech­ anisms can be interpreted intuitively by Fig. A1. A detailed description of the relationship between the rock failure mechanism and MT can be found in Refs. 30 and 31, which will not be discussed here. 2.1.2. Method of MT inversion Although the MT inversion has been increasingly applied in mines in recent years, the results are sometimes not reliable. The main reasons for the inaccurate results are as follows30,31: (1) the waveforms are complex and high frequency, reflecting small-scale inhomogeneities in rock masses; (2) the moment magnitude of MS events in mines is usually low with noisy waveforms; and (3) because the scope of MS monitoring is usually tens or hundreds of meters in mining engineering, the presence of unfavorable geological bodies in the monitoring area can make the estimated Green’s function more inaccurate. To handle the issues mentioned above, we use the PCA method proposed by Vavry�cuk et al. to estimate the MTs for handling the issues with complex and relatively low signal-to-noise ratio waveforms.31 And we use the hybrid MT method proposed by Linzer30 for weakening the influences of the propagation medium on the solution of the Green’s function. The key step of PCA method is separate the common wavelet from other interference signals and reduce the errors of the observed amplitudes at each station using the PCA of the original waveforms. The mathematical description of this method is as follows:

2. Methodology Before the occurrence of a seepage and water inrush disaster, there must be a process from emergence to development. In this process, the rock mass breaks and forms fracture channels, which provides condi­ tions for the occurrence of water inrush disaster, and this process is key to effectively grasping the evolution process and formation mechanism of the fractures under the dynamic disturbance of the mine to prevent and control the water inrush disaster, which has practical value in en­ gineering. The MT provides a feasible approach for illuminating un­ mapped fractures,13,14 analyzing the failure mechanism of the source,15–17 imaging damage zones,18 describing fault activation,19 etc.

XPCA ¼ XW

(3)

where X is the data matrix, and each column represents a waveform received by one station. Single component sensor is widely used in mines. For three component sensor, only Z component is taken for calculation in this paper. The data matrix X can be mapped into an orthogonal space and generate the matrix XPCA by multiplying the ma­ trix W of orthonormal eigenvectors of the covariance matrix XTX. X PCA contains several common wavelets of the original waveforms. Seismic events with the similar source location and failure mecha­ nism, usually have similar waveforms,32 namely, matrix X is usually

2.1. Theory and method of MT inversion 2.1.1. Theory of MT inversion In recent years, MT inversion, an effective tool for analyzing failure mechanisms and processes, has been widely used in laboratory acoustic emission tests,20,21 induced seismicity,22,23 mine MS,24,25 etc. In the process of the rock mass failure, the information of failure mechanism and fracture geometry can be determined by MT inversion. A seismic source or microcrack can be induced by nine equivalent force couple in a matrix, called the MT. Any MS source can be mathematically repre­ sented by a MT, which can be calculated by the MT inversion method. An MT can be expressed as follows26: GM ¼ u

trðMÞ P ⋅100% ðjtrðMÞj þ jm*i jÞ

formed by a set of similar traces, and the first principal component X PCA is dominant and corresponds to the common wavelet (see Fig. 1). Then, Eq. (3) can be written as follows: ð1Þ

ð1Þ

X ¼ XPCA W Tð1Þ

(4)

where the matrix W(1) indicates the weights of the dominant common wavelet (the first principal component) of the individual traces, which are the coefficients of the PCA decomposition and can be set as effective amplitudes for solving MTs containing the polarities.

(1) 2

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 1. (a) Filtered waveforms of one MS event (bandpass 50–600 Hz). (b) the first principal component together with original waveforms. (c) The first, second, third and fourth principal components.

After calculating amplitudes and polarities of the waveforms by PCA method, we can determine the MT. The hybrid MT method proposed by Linzer30 could weak the influences of the propagation path and sensor. This method assumes that MS events in a cluster share similar ray paths, thus, the Green’s function solutions between a station and clustered MS events are the same. The MS events were clustered by the K-means clustering technique in this paper.33 The reason for clustering analysis is that MS events in a certain time and a certain space often have similar formation mechanisms.34 For more details on the application of K-means clustering to MS events, please refer to Ref. 35. In this paper, the MS

event location (x, y, z) and occurrence time t were used as the cluster input parameters. In addition, the widely used Krzanowski-Lai index36 as another input parameter was used for determining an adequate number of MS clusters in this paper. The solution process of the hybrid MT method is listed in Section 2.1.3. 2.1.3. Implementation steps To clearly understand the content of the MT inversion based on PCA method and hybrid MT inversion method used in this paper, the detailed implementation steps are summarized as follows, and the flowchart is

Fig. 2. Flowchart of the inversion of seepage channels based on MS data. 3

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

shown in Fig. 2a and b:

source associated with different failure mechanisms under the condition of known principal stresses according to the relationship between failure mechanisms and principal stresses. The criteria proposed in Ref. 40 can be summarized as follows. 8 < Ts ¼ τ=σ n ; Shear failure Tt ¼ ðσ1 σ n Þ=ðσ1 σ 3 Þ; Tensile failure (6) : Tc ¼ ðσ n σ3 Þ=ðσ1 σ 3 Þ; compressive failure

(1) Data preprocessing. The input is MS waveforms of the rock mass failure received by sensors. The signal-to-noise ratio of the waveform is enhanced by the bandpass filter method. The P-wave onset is obtained with an automatic picking algorithm, the Sus­ pension Bridge Picking (SBPx) algorithm of Meier.37 The MS event location method we used is an iterative method, the sim­ plex algorithm.38 (2) Calculation of the PCA amplitudes and the MT inversion. As described in Section 2.1.2, the PCA method is used to decompose the waveforms and extract the first principal component of multiple waveforms of MS events. Then, the PCA coefficients of the first principal component can be solved by Eq. (4), which serve as the effective amplitude (including their polarity) used for the amplitude MT inversion. The MTs of all MS events calculated from this step is the output of the algorithm described in Fig. 2a. (3) MT correction based on hybrid MT inversion method

where

τ ¼ ½ðσ 1

th obs uobs ij (rij ¼ uij =uij ) is used to verify the validity of the result of MT for

�2

�

σ 3 m2 n2 þ σ 1

�2

σ3 l2 n2

i1=2

(8)

2.2.2. Under the condition of unknown principal stresses It should be emphasized that the above method (Eqs. (6)–(8)) as­ sumes that the stress field distribution is known. If the principal stress distribution is unknown, the MS-derived fracture can be obtained by a stress inversion method based on focal mechanism. Many scholars have proposed many useful methods for stress inversion using the focal mechanism41,42 based on the Wallace-Bott assumption43,44; namely, the slip vector of the seismic source is parallel to the resolved shear stress on the fracture. Four parameters of the stress tensor, the orientations of the three principal stresses (σ1 , σ 2 and σ3 ) and the stress ratio Rs ¼ ðσ 1 σ2 Þ=ðσ 1 σ 3 Þ, can be obtained by the stress inversion method. The calculation of this method is described as follows: The traction T related to stress tensor σij acts on the fracture with normal n, which can be written as follows:

sensor i and event j. Finally, dynamically update the input amplitudes at each sensor i by the following Eq. (5). Repeat the above steps until the ratio correction factor ~ri 1 satisfies a certain end condition such as ð~ri 1Þ < ε(ε takes 0.02 in this paper). The final MTs of all MS events are corrected, which correspond to the output of the hybrid MT inversion method described in Fig. 2b. 1Þ

σ2 Þ2 l2 m2 þ σ2

where l, m and n are the cosines of the normal direction of the fracture with respect to the principal stress axes, σ 1 , σ2 and σ 3 , respectively; Ts , Tt and Tc are called the failure tendency indexes of seismic source under shear, tensile and compressive failure mechanisms, respectively. When the principal stresses in the source are known, the failure tendency indexes of two focal planes for the same source can be obtained under the corresponding failure mechanism from Eq. (6) and the focal plane with a larger failure tendency index can be regarded as a reasonable MS-derived fracture. For more detailed solution process, please refer to Ref. 40.

The solution process of the hybrid MT inversion method is summa­ rized as follows, and the flowchart is shown in Fig. 2b. First, the output of step (2) is clustered by K-means clustering method and used as the input of hybrid MT inversion method. Secondly, the results of all the MS events are used to back-calculate the amplitudes at individual sensor by Eq. (1). The ratio rij of theoretical amplitude uth ij and observed amplitude

old old unew ri ij ¼ uij þ wi uij ð~

(7)

σ n ¼ σ 1 � l2 þ σ2 � m2 þ σ3 � n2

(5)

where wi is the weighting factor resulting in the increase and decrease of amplitude; ~ri is the median ratio of sensor i. 2.2. Determine the reasonable fracture from MS data

(9)

Ti ¼ σ ij nj

The fracture formed by mining disturbance can be obtained by MT analysis, which facilitates the establishment of the relationship between MS data and rock mass fracture. However, a single seismic source can yield two interchangeable focal planes, which is caused by the basic assumption of the seismic source model.39 In seismology, one focal plane, called the fault plane, is the fracture plane needed for research, and the other is called the auxiliary plane. In Appendix A, the rela­ tionship between focal planes and failure mechanisms was discussed. It is stated that the focal plane is parallel only when pure tensile or pure compressive failure occurs, and the focal plane forms a certain angle, not parallel or vertical, under other failure mechanisms. Therefore, we need a method to extract a reasonable fracture from two interchangeable focal planes. Determining the fracture from focal planes of seismic sources is of great significance for analyzing the propagation of fractures and the formation of seepage channels. The methods involved are described in detail below, and the detailed flowchart is shown in Fig. 2c. Hereafter, the term “fracture” refers to the “fault plane”, the term “MS-derived fracture” refers to the “reasonable fracture” derived from an individual MS event, and the MS-derived fracture is assumed to be planar and penny-shaped.

with its normal and shear components σn and τ: (10)

σ n ¼ Ti ni ¼ σ ij ni nj τNi ¼ Ti

2.2.1. Under the condition of known principal stresses Because different failure mechanisms have different stress relations (Fig. 3), In Ref. 40, the authors proposed criteria for extracting the MS-derived fracture from two interchangeable focal planes of a single

σn ni ¼ σ kj nk ðδik

ni nk Þ

(11)

Fig. 3. Griffith-Coulomb diagrams illustrating different failure modes: pure tensile failure, tensile-shear failure, pure shear failure and compressiveshear failure.40 4

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

where δ is the Kronecker delta and N is the vector direction of τ. The basis of this method is that the vector direction N is parallel to the slip vector v. Because this method cannot determine absolute stress values, τ is normalized to be 1 in Eq. (11). Therefore, the equation of stress inversion can be written as following: v ¼ σkj nj ðδik

(2). The geometric parameters (size, normal, slip vector and aperture) of focal planes of every MS event are determined by Appendixes A and B. (2) Determine the MS-derived fracture from MS data. If the principal stress is known, the MS-derived fracture can be determined by calculating the failure tendency indexes of the focal planes of each seismic source by Eqs. (6)–(8). If the principal stress is un­ known, then the MS-derived fracture should be determined by the stress inversion method described in Section 2.2.2.

(12)

ni nk Þ

Subsequently, Eq. (12) can be expressed in matrix form: (13)

Aσ ¼ v

According to Eq. (18), if the fracture normal of a single source is unknown, each focal plane has a 50% chance of being selected in the calculation process. Obviously, the incorrect selection of the focal plane will affect the inversion result of the principal stresses. To solve this problem, we developed a new package for obtaining the fractures associated with different failure mechanisms based on the “STRESSIN­ VERSE” package.42 In the new package, we replaced the principal stress inversion equation (Eq. (13)) (used in “STRESSINVERSE” package) with Eq. (17). In addition, the " STRESSINVERSE 00 package uses an index similar to shear failure tendency Ts in Eq. (6) to determine all fractures, assuming that all fractures are caused by shear failure. In the new package, we determine the fracture according to corresponding failure mechanism by Eq. (6). The detailed procedure of the new package is described below. Initially, one of the two focal planes of each MS event are randomly selected for the inversion of the principal stresses by Eq. (17). According to the inverted principal stresses in this step and the criteria listed in Eqs. (6)–(8), MS-derived fractures with different failure mechanisms are selected. Then, the MS-derived fractures are used for the next step of stress inversion. The above operations are repeated until the userdefined termination condition (user-defined iterations) is satisfied. Finally, the MS-derived fractures are determined. Thus, an MS-derived fracture, including the orientation (strike and dip), can be determined by the above method.

where σ is the reduced stress tensor: (14)

σ ¼ ½σ 11 σ 22 σ 33 σ 23 σ 13 σ 12 �T

where A is a 3 � 6 matrix calculated from the fracture normal n: 2 3 � � � 2 2 2 2 7 6 n1 1 n2 n1 n3 n1 n2 2n1 n2 n3 n3 1 2n1 n2 1 2n1 7 1 6 � � �7 6 A ¼ 6 n2 n21 n2 1 n22 n2 n23 n3 1 2n22 2n1 n2 n3 n1 1 2n22 7 7 6 � � � 4 n3 n2 n3 n2 n3 1 n2 2n1 n2 n3 5 n2 1 2n2 n1 1 2n2 1

2

3

3

3

(15) For m MS-derived fractures determined by m MS events, 3 � m linear equations can be solved through the linear least square method for obtaining principal stress by Eq. (13). Vavry�cuk42 developed a MATLAB package “STRESSINVERSE” based on this method for solving the three principal stresses. The stress inversion method mentioned above assumes that the failure mechanism is only shear failure. However, the slip vectors of the sources are not always parallel to the fracture surface especially in mines; thus, this assumption is not always valid. Therefore, the stress inversion method described above needs to be modified. Referring to a modification method developed in Ref. 45, the traction T is no longer decomposed but expressed as a resultant force; in other words, the di­ rection of the traction T is considered to be parallel to the slip vector determined by MT, which is not parallel to the fracture surface under tensile and compressive mechanisms. Therefore, Eq. (9) can be changed into the following form45: � T ei ¼ σ eij nj ¼ σij δij P nj ¼ v (16)

2.2.4. Numerical test We used the MS data and the stress field data of Shirengou Iron Mine in Section 3 to verify the availability of the above stress inversion method. By the measurement of in situ stress, the orientation of the maximum horizontal stress in Shirengou Iron Mine is NE 60� 70� . Therefore, the orientations of σ 1 , σ 2 and σ 3 can be assumed to be plunge 0� , azimuth 65� ; plunge 0� , azimuth 155� ; and plunge 90� , azimuth 90� , respectively. The principal stress ratio is Rs ¼ 0.75. According to the data of principal stresses in Shirengou Iron Mine, the MS-derived fractures can be determined from focal planes using three indexes, Tt , Ts and Tc . Then, the stress inversion method described in Section 2.2.2 and Fig. 2c is used to determine the principal stresses (σ1 , σ 2 , σ3 and Rs). As described in Section 2.2.3, we can determine the MS-derived fractures from focal planes while inversing the principal stress by the new package. Next, a comparative analysis is performed between the fractures determined directly by the measured principal stress and the fractures determined by the inverted principal stress through the new package, and the results are shown in Fig. 4. In this figure, the blue dot shows that the MS-derived fracture determined by the measured principal stress and the inverted principal stress is consistent for the same seismic source, and the red dot indicates that the MS-derived fracture determined by the measured principal stress and the inverted principal stress is opposite, that is, the other focal plane is selected as the MS-derived fracture by the stress inversion method for the same source. The proportion of red dots is 10.36%, which indicates that there is an error of 10.36% of MS-derived fractures obtained from stress inversion. Meanwhile, Fig. 4 also shows the orientations of measured and inversed principal stresses. The yellow dot is the measured principal stress orientation, and the green dot represents the principal stress orientation inverted from the MS data by the stress inversion method. The errors in the principal stress orientations between

where σ eij is the effective stress tensor, and P is the fluid pressure.

Subsequently, Eq. (16) can be transformed into a matrix form, as follows: A σe ¼ v

(17)

where 2

3 0 0 n3 n2 n1 0 A ¼ 4 0 n2 0 n3 0 n1 5 0 0 n3 n2 n1 0

(18)

Finally, the direction of the three principal stresses and the stress ratio Rs can be obtained by solving Eq. (17). Obviously, when the principal stress is solved, the MS-derived fracture can be determined according to the corresponding failure tendency index (Eqs. (6)–(8)). 2.2.3. Implementation steps To clearly understand the method of determining the MS-derived fracture from MS data used in this paper, the detailed implementation steps are summarized as follows, and the flowchart is shown in Fig. 2c. (1) Determine the failure mechanism of MS events and the geometric parameter of focal planes. The corrected MT obtained from Sec­ tion 2.1 is used as the input of the algorithm described in Fig. 2c. The failure mechanism of every MS event is determined by Eq. 5

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 5. (a) Spatial relationship of fractures. (b) Fractures and intersection lines.

2.3.2. Transforming a MS-derived DFN into a graph structure The intersection and connection of multiple MS-derived fractures in the MS-derived DFN forms many complex connected channels. These complex connected channels can be determined by transforming the MSderived DFN into a graph structure through graph theory. In mathe­ matics, a graph is a structure comprising a set of nodes with pairs of nodes linked by edges.49 According to this approach, the MS-derived DFN can be represented as a graph. The geometrical information about the MS-derived fractures, such as the apertures, the length of the intersections and flow hydraulic properties, including permeability, pressure drop and porosity, are assigned to the edges connecting the nodes (see the right of Fig. 6a). For details of graph theory, please refer to Refs. 49–52. In the graph, only those connecting the inlet surface and outlet surface can become the seepage channels. These seepage channels with the significance of hydraulic conductivity can be selected from the complex connected channels by the graph search algorithm. The depth first search (DFS) algorithm was used in this paper. The MS-derived DFN generated by rock mass damage can determine the seepage channel through the above operations. In the study of fracture seepage, a three-dimensional DFN is usually transformed into a pipe model. A pipe for a MS-derived fracture can be established between its centroid and the center of the intersection line with another MS-derived fracture, as shown in the left of Fig. 6a. A schematic diagram of the transformation from a DFN to graph structure is also shown in the right of Fig. 6a. The node is shown as a red sphere, the intersection line is shown as a gray tube and the pipe is shown as a blue tube. Adjacency list is used to store the graph, which is a common method and occupies less memory space.53 The data in Fig. 6a can be the converted into an adjacent list, as shown in Fig. 6b. The ‘data’ in Fig. 6b represents the data domain and stores vertex information, likes the vertex A in the right of Fig. 6a. The ‘dest’ refers to the adjacency point domain, likes the end point B of edge AB in the right of Fig. 6a. The ‘weight’ is used to record geometric and flow hydraulic properties of edge, likes ωAB in the right of Fig. 6a.

Fig. 4. Stereographic projection figure containing the inverted principal stresses, measured principal stresses and the orientation of MS-derived fractures determined by the measured principal stress and the inverted principal stress.

the inverted and measured results are within 10� , and the error in the principal stress ratio Rs is 0.12. Therefore, we think the error is acceptable under the condition of unknown principal stresses, and the stress inversion method can effectively invert the principal stresses and determine the MS-derived fracture. 2.3. Seepage channel analysis of the MS-derived DFN Through the above method, the MS-derived fracture group can be determined from a set of MS data, and this group can be called as MSderived DFN. The size and aperture of the MS-derived fracture in the MS-derived DFN can be calculated by Appendix B. The pre-existing fractures are seismically reactivated (MS-derived fractures) due to external disturbances, which will change the fluidity of fluid in the fractures. In this study, only the MS-derived fractures were considered, and un-activated fractures were neglected. Fluid transport in the MSderived DFN is influenced by the intersection, topology and the geo­ metric properties of MS-derived fractures.46,47 2.3.1. Intersection of MS-derived fractures Accurately evaluating the connectivity of MS-derived fractures is the primary task when determining the seepage channel. The intersection of MS-derived fractures forming seepage channels play an important role in fluid transport for low-permeability rock masses.48 Fluids flow from one MS-derived fracture to another through the intersection segment of the MS-derived fractures. A way of evaluating the connectivity of the MS-derived fractures is intersection analysis. There are four types of spatial relations among two MS-derived fractures in three-dimensional (see Fig. 5a): intersection, connection, intercalation, and separation, which can be converted to mathematical problems. The data of each MS-derived fracture consist of center coordinates x, y, z, radius r, dip α and strike β. According to these geometric parameters of two MS-derived fractures, geometric analysis can be used to determine whether there is intersection line between two MS-derived fractures. If an intersection line exists, the equation of the intersection line and the coordinates of the intersection points can then be determined. A three-dimensional DFN with 100 randomly generated fractures is used to shows the intersection analysis of fractures, as shown in Fig. 5b, in which the red lines represent the intersection lines.

2.3.3. Fluid flow through MS-derived fractures Darcy’s law50 is usually used to effectively evaluate the fluid flow through a single fracture. Two basic assumptions are needed for this paper: the rock material is impermeable, inert and incompressible, namely, only seepage in MS-derived fractures is considered, and each MS-derived fracture is comprised of two parallel flat and smooth plates having a specified aperture and equivalent permeability. For the graph structure of the MS-derived DFN, the equilibrium area is formed through node i and k edges connected to node i, and the flow equation is obtained as follows:

6

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 6. (a) Transforming a DFN into a graph structure. The left is the schematic diagram for determining the pipe model by the DFN. The right is the schematic diagram of the graph structure determined by the DFN. (b) Sketch map of the adjacency list.

!

k X

qj j¼1

qj ¼ Cij ΔHj

(19)

þ Qi ¼ 0 ði ¼ 1; ⋅ 2; ⋅ ⋅ ⋅ ⋅ ⋅ kÞ

where

i

where qj represents the flow rate of inlet or outlet of node i; k represents the number of edges which connect with node i; and Qi represents the source or sink term at node i. In the graph structure of the MS-derived DFN, there are three types of nodes, inlet nodes, inner nodes and outlet nodes. The pressure head values of inlet and outlet nodes set to a certain value, such as H1 ¼ 1 for inner nodes and H2 ¼ 0 for outlet nodes. If there are N nodes in the graph structure of MS-derived DFN, Eq. (19) can be turned into N equation sets. Combining with boundary conditions, a mathematical model for seepage in the MS-derived DFN can be expressed in matrix form. 8 < Aq þ Q ¼ 0 H ¼ H1 (20) : inlet Houtlet ¼ H2

Cij ¼

(23)

2.3.4. Shortest fluid flow channel Because of anisotropy in a rock mass, the MS-derived fractures that form in response to an external disturbance usually have different geo­ metric and hydraulic properties due to diverse fluid flow properties in various seepage channels. There are usually multiple seepage channels between storage points and water inrush points in rock mass engineer­ ing. Some MS-derived fractures with larger apertures can form a shorter seepage channel, making the flow transport easier in this channel. As shown in the right of Fig. 6a, edges between different nodes have different flow weights because of the geometric and hydraulic properties of the MS-derived fractures. In this paper, the flow time Tf in the MSderived fracture is defined as the weight of the conductivity. The

T

n P

Cij Hi Hj ¼ i¼1 n P Cij

ga3 b 12νL

where Hj is the head pressure and Cij is the conductance between nodes i and j calculated based on gravitational acceleration g, the aperture of edge a, the length of intersection line b, the kinematic viscosity v, and the length of edge L. Because the aperture of different MS-derived fractures in the seepage network is inconsistent, the average of the aperture of MS-derived fractures corresponding to nodes i and j is taken as the aperture of the edge ij in this paper. According to Eqs. (20)-(23), the flow rate q of each edge and the head pressure H of each node were obtained using numerical methods.

� � where Q ¼ ½Q1 ; Q2 ; :::; Qn � ; q ¼ ½q1 ; q2 ; :::; qm � ; A ¼ aij n�m is a connectivity matrix, describing the connected relation between nodes and edges in the graph structure of MS-derived DFN. The unknown quantity in Eq. (20) is the head pressure H of inner nodes in the graph structure of the MS-derived DFN, which can be solved by the following equations from Darcy’s law: T

(22)

(21)

i¼1

7

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

derivation process of weight ω is shown as follows. The fluid velocity in an edge is as follows: V¼

Q KΔH ¼ A νL

(24)

where Q represents the flow through edge, ΔH represents the head dif­ ference of edge, v represents the kinematic viscosity, and A, L and K represent the cross-sectional area, the length and permeability of the edge, respectively. The time Tf indicates the water flow time from one MS-derived fracture to another, that is, the flow time of fluid in the edge:

ω ¼ Tf ¼

L νL2 ¼ V KΔH

(25)

Based on Eq. (25), the water flow time is used as the weight for the seepage channels between the storage points and water inrush points. Namely, the lowest water flow time with the lowest weight suggests the best seepage channel in the grapy structure of the MS-derived DFN. The shortest weighted channel in the graph structure can be detected using path-finder algorithms, such as the Dijkstra algorithm.51 To demonstrate the above calculation process of Section 2.3, we built a DFN containing 400 random fractures as shown in Fig. 7a. By calcu­ lating the intersection lines between fractures, the pipes were built ac­ cording to the method illustrated in Fig. 6a. Subsequently, the pipe model was transformed into a graph structure, and the result is shown in Fig. 7b. Let the top surface of the fracture seepage network cube be the inlet surface H1 ¼ 1, the bottom surface be the outlet surface H2 ¼ 0, and the other surface be the impermeable boundary Hother ¼ 0. All nodes in the graph structure were calculated according to Eqs. (20)-(23), and the results containing the pressure head values at all nodes and the interpolated pressure values for all edges are presented in Fig. 7c, which shows the distribution of the fluid pressure head. Sections 2.3.1-2.3.4 are described step by step in accordance with the process of seepage channel analysis of the MS-derived DFN, as shown in Fig. 2d; therefore, the detailed implementation steps will not be described separately as in Sections 2.1.3 and 2.2.3. In addition, the procedure of fluid calculation in this section was adapted from the openaccess MATLAB code ADFNE package developed by Alghalandis,54 which is available for download. 3. Engineering application of water inrush in the fractured zone of the Shirengou Iron Mine 3.1. Engineering background and MS data The Shirengou Iron Mine is a typical mine transition from open pit mining to underground mining in China. Fig. 8 describes the geological engineering conditions of the Shirengou Iron Mine.55 In Fig. 8a, many stopes have been left at the - 60 m level due to decades of mining ac­ tivities. Fig. 8b describes the surface situation with abundant accumu­ lated water and a local landslide in the pit bottom in 2009. It should be noted that the height of the accumulated water at the pit bottom reached tens of meters with a large pressure on the crown pillar. Serious defor­ mation and many local collapses occurred at 0 m level, and there were several seepage zones on the roof of the roadway. In addition, there were some large fractures and seepage zones in the surrounding rock at the 60 m level. To sum up, the rock mass quality in blue line zone in Fig. 8a was very poor. This area was selected as a study area and deserves monitoring and our attention. To effectively monitor the seepage and water inrush disaster in the fractured zone of the Shirengou Iron Mine, ESG MS monitoring system56 was established to ensure the safe pro­ duction of the mine. Fig. 8c shows the sensors distribution around the study area. Eleven sensors (thirteen channels) arranged around the study area were centrally distributed, covering the study area well and ensuring a high locational accuracy. The MS data recorded in the study area are closely related to the development of seepage channels, thus

Fig. 7. Modeling of fluid flow in fractures. (a) A DFN. (b) Graph structure (yellow sphere indicates node and blue tube indicates edge). (c) Distribution of the pressure head. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

8

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 8. (a) 3D model of the Shirengou Iron Mine. (b) Site condition of the pit bottom in 2009. (c) Sensors arrangement in the study area (the red cylinder represents the triaxial geophone sensor and the gray cylinder represents the uniaxial geophone sensor). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 9. Distribution of the MS events (upper two panels) and the MT beachballs (lower two panels). 9

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

providing high value for our study (see Fig. 9). Fig. 9 shows the distribution of MS events and the corresponding MT beachballs recorded from October 27, 2009, to January 10, 2010. The source radius of MS event, i.e. the radius of MS-derived fracture, is 4.5–14.4 m, which was obtained from Eq. (B1). The seismic moment M0 is from 1.37E7 N.m to 7.49E9 N.m, and the moment magnitude Mw of MS events is from 1.41 to 0.51, which can be obtained from Eq. (26). Readers can refer to Ref. 55 for a detailed description about the source parameters of the MS events in the study area. Through the frequency analysis of the MS waveforms, the bandpass frequency is taken as 50–600 Hz when filtering in this paper. The MTs of all MS events were calculated based on PCA method and corrected based on hybrid MT inversion method. The accuracy of MT solutions for all MS events in the study area is shown in Fig. 10. The RMS of the MTs in Fig. 10 defined as the normalized differences between the synthetic amplitude Ath i and

3.2. Formation process of fracture channels It should be emphasized that before the arrangement of the MS monitoring system, the initial fracture formed due to the mining activity is like a black box, and we have no data to determine the role it plays. Therefore, the continuity monitoring of MS monitoring is very important to evaluate the seepage channel. According to the actual observation in the study area, we can only clarify that the previous mining activities did not form a macro seepage channel. In this section, the influence of external disturbance on the formation of seepage channel is ignored before the arrangement of MS monitoring. To study the formation process of fracture channels, the occurrence period of MS events is divided into four parts: Stage 1 (from October 27, 2009, to November 15, 2009), Stage 2 (from November 16, 2009, to November 29, 2009), Stage 3 (from November 30, 2009, to December 15, 2009) and Stage 4 (December 16, 2009, to January 1, 2010). For discussion of the temporal and spatial variations of MS events at various stages, please refer to Ref. 55. The formation process of the fracture channels obtained from MS data is shown in Fig. 11. In Fig. 11, the fracture channels are represented in the form of pipes, which is the result of the intersection of MS-derived fractures. During Stage 1, the fracture channels were mainly distributed in two concentrated areas as follows: one concentrated area formed between stopes No. I16 and No. I17, and the other formed between stopes No. 6 and No. 18. However, the frac­ ture channels in these two concentrated areas were not linked together. The reason for the formation of the first accumulation zone is that tailings were dumped into the pit bottom in the study area, resulting in the rock mass failure in the west side of the stope No. I16 due to the large overburden load. The reason for the second accumulation zone is that the residual ore recovery in stopes No. 6 and No. 18 increased the reactive degree of the pre-existing fractures. During Stage 2, the connectivity of the fracture channels in Stage 1 was enhanced by the newly derived fractures in this stage. The fracture channels of these two concentrated areas became intensively connected. The fracture channels from the pit bottom (water storage area) to the 60 m roadway initially formed. In this stage, the accumulated water seeped into stope No. 18, which was consistent with field records. Many large fractures were observed around stopes No. I16 and No. 18 increasing the damage degree of rock masses, and abundant water was present in stope No. 18. During Stage 3, the fractures became further concentrated in the pit bottom and in stopes No. I16 and No.18, corresponding to the second damage development period.55 As a result, the connectivity degree of fractures between stopes No. I16 and No. 18 increased, the number of fracture channels increased, the water conductivity of the seepage network increased, and the accumulated water volume in stope No. 18 increased. During Stage 4, the final MS-derived DFN was formed as shown in the upper two panels of Fig. 12. Because a MS-derived DFN with good water conductivity formed in Stages 1–3, the connectivity of the final MSderived DFN was greatly enhanced due to the new MS-derived frac­ tures in this stage. It can be seen from the upper two panels of Fig. 12 that MS-derived fractures with high connectivity are mainly distributed in between stopes No. I16, No. 18 and 60 m level. Thus, the final fracture seepage channel formed as shown in the lower two panels of Fig. 12. The hydraulic calculation of the MS-derived DFN in the study area was carried out using the method proposed in Section 2.3.3. During this calculation progress, the upper surface of the cuboid frame corre­ sponding to the pit bottom was set as the inlet (H ¼ 1) and the bottom surface corresponding to the 60 m level was set as the outlet (H ¼ 0). The pressure head at all nodes are demonstrated in the lower two panels Fig. 12. The network in the figure is the effective seepage network connecting the inlet surface and outlet surface. The decrease in water pressure between stope No. I16 and the pit bottom was low, namely, stope No. I16 and its surrounding roadway were prone to water inrush

observed amplitude Aobs i in the MT inversion, as shown in Eq. (27). A RMS value greater than 0.3 indicates unreliable MTs.31 Fig. 10 shows that the MTs calculated by combining PCA and hybrid MT inversion method has higher accuracy, and more MTs have RMS lower than 0.3. In the following study, MTs with RMS greater than 0.3 were deleted. 2 Mw ¼ logM0 3

6:0

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffi P synth 2 Ai Aobs i qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffi RMS ¼ P synth 2 Ai

(26)

(27)

The MS events are mainly distributed around stopes No. I16, No. I17 and No. 18. Fig. 9 shows that Shear failure MT beachball is dominant, tensile failure MT beachball is rare and scattered accounting for a low percentage, and compressive failure MT beachball related to blasting is the lowest. The reason for the generation of shear failure MS events is that during occurrence period of MS events, tailings were dumped into the pit bottom and the surrounding rock was subjected to a great load due to the increase in the overburden load. Moreover, the accumulation of water at the pit bottom led to the expansion and lubrication of the fractures. The tensile failure MS events are formed due to the highpressure water in the study area and some mining activities. Readers can refer to Ref. 55 for a detailed description about the reason of the formation of MS events and its focal mechanisms.

Fig. 10. Histogram of the RMS residuals of the 328 MTs calculated by the standard-amplitude MT inversion method and PCA combined with Hybrid MT inversion method. 10

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 11. Distribution of fracture channels during different stage in the form of pipe. Yellow sphere represents node and red pipe represents fracture channel. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 12. Final MS-derived DFN (upper two panels) and the seepage channels colored by pressure head (lower two panels). Upper surface corresponds to inlet, and bottom surface corresponds to outlet.

due to the high-water pressure. In addition, the water pressure near stope No. 18 was at a medium level and well connected with the above MS-derived fractures; therefore, water inrush was also likely to occur. In the areas with seepage channels, corresponding measures should be taken to prevent water inrush.

of Fig. 13, the roof of stope No. 18 was taken as the study object and three nodes around it were taken as the outlet nodes. First, taking one of the three nodes as an outlet node, the seepage channels between all the inlet nodes and this outlet node were traversed to find the shortest channel with the least weight. Then, the shortest channels of the other two outlet nodes were calculated in the same way. Channels 1, 2 and 3 in the left part of Fig. 13 indicate the three shortest channels from the open pit to the roof of stope No. 18. The colors of the channels in this figure represent the magnitudes of the weights ω, with red representing the largest weight. As evidenced by the results, these three paths had a common seepage channel, Ch 1.

3.3. Analysis of the main seepage channel As mentioned in Section 2.3.4, the channel with the shortest flow time deserves attention. The shortest channel in this case was studied with on the basis of three noteworthy positions. As shown in the left part 11

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 13. Shortest channels associated with the three study objects. Upper surface corresponds to inlet, and bottom surface corresponds to outlet.

Similarly, two nodes on the roof of stope No. 6 were used to identify the shortest seepage channels, and channels 4 and 5 were obtained, as shown in the middle part of Fig. 13. Nevertheless, the weights of channels 4 and 5 were larger, indicating that water inrush was more likely to occur in stope No. 18 than in the roof of stope No. 6. It is interesting that channels 4 and 5 also had a common seepage channel, Ch 2. In the right part of Fig. 13, all nodes at the 60 m level were taken as outlet nodes, and the seepage channel with the smallest weight, Ch M, was obtained. This channel was the shortest channel with the least weight in all channels of the study area. It was easy to conclude that many segments of the channels Ch M, Ch 1 and Ch 2 were coincident. Therefore, the coincident segments likely played a leading role in the seepage of the study area and were the main seepage channel in the study area. Reference 40 gives the results of the permeability distribution along the longitudinal plane with maximum permeability of the same study area and the seepage of actual observation points, as shown in Fig. 14a. Reference 55 gives the shear failure zone distribution of the same study area through the coupling analysis of MS data and numerical simulation, as shown in Fig. 14b. According to Fig. 14a and b, it can be seen that the main seepage channel is consistent with the distribution trend of high permeability and shear failure zones. Meanwhile, we can draw the conclusion that the connected shear failure zone is the reason for the formation of the main seepage channel. In addition, as can be seen from Fig. 14a, the field hydrologic condition of observation point 1 at the 0 m level of the roadway was extremely poor with serious seepage from the roadway wall. It is a pity that there was no observation condition for the pillar of stope No.18, but the field hydrologic condition of the area near the stope No.18 (observation point 2) was recorded. There were several large fractures with continuously seepage in the wall of observation point 2. By analyzing the field hydrologic condition records of these two observation points, it can be proved that the main seepage channel determined by the inversion method proposed in this paper is reasonable and has certain engineering guiding significance. Finally, the orientations of the MS-derived fractures in the coincident segments of the channels Ch M, Ch 1 and Ch 2 were statistically analyzed, and the results are shown in Fig. 14c. Accordingly, the water inrush in the fractured zone of the Shirengou Iron Mine was controlled

by two sets of MS-derived fractures with orientations of 159� ∠46� and 251� ∠27� . These two sets of MS-derived fractures are mainly deter­ mined by the shear failure seismic source. The shear failure tendency index of the two sets of MS-derived fractures was calculated by Eq. (6) with the measured stress conditions (see Fig. 4) and plotted in Fig. 14d. It can be seen from this figure that the two sets of intersecting MSderived fractures are mainly distributed in the region with high shear failure tendency index. Thus, the formation of the two sets of MS-derived fractures forming the main seepage channel is closely related to the stress conditions in the study area. 4. Conclusions In the present study, a set of methods for the inversion of seepage channels using mining-induced MS data was proposed and applied to Shirengou Iron Mine that has experienced water inrush. The methods were divided into three main parts. The first part was to solve the MT. To weaken the influence of interference signals on the polarity and ampli­ tude of P-wave and to weaken the influence of the wave propagation medium, the method combining PCA method with the hybrid MT inversion method was used to solve the MT. The second part involved the method of determining the MS-derived fracture. To select a reasonable fracture from two focal planes derived from a source, three indexes, Tt , Ts and Tc , were used for the condition of known principal stresses. If the principal stresses were unknown, the MS-derived fracture could be determined by iteratively combining the stress inversion method with the Tt , Ts and Tc indexes. The last part was the seepage channel analysis of the MS-derived DFN. The MS-derived DFN was transformed into a graph structure by graph theory, and the hydraulic properties of MS-derived fractures were solved according to Darcy’s law. Moreover, the weights were defined by the flow time of the fluid in the MS-derived fractures, and the Dijkstra algorithm was used to identify the shortest seepage channel. The above-mentioned methods were applied to the water inrush analysis in the fractured zone of the Shirengou Iron Mine. Through the analysis of seepage channels in the study area of the Shirengou Iron Mine, it was concluded that the roof of stope No. 18 was more prone to water inrush than that of stope No. 6, and that the MS-derived fractures between stope No. I16 and the pit bottom were complex and had high 12

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Fig. 14. (a) Permeability distribution and field observation of the study area40 (b) Shear yield zone distribution of the study area calculated by MS monitoring and numerical simulation55 (c) Orientation of MS-derived fractures in the main seepage channel (d) Shear failure tendency of MS-derived fractures in the main seepage channel.

water pressure; consequently, these fractures deserved more attention. At the same time, the study showed that the main seepage channel was controlled by two sets of MS-derived fractures with orientations of 159� ∠46� and 251� ∠27� in the study area of the Shirengou Iron Mine. The results were proved by coupling analysis of the observations of the field hydrologic condition and the distribution of high permeability and shear failure zones, and can provide a reference for the formulation of water plugging measures. The purpose of this paper is to make a preliminary attempt to use MS data to invert seepage channels. However, MS monitoring is often affected by the complexity of rock mass structure and the diversity of mechanical properties of the rock mass. This will lead to incomplete data record and inaccurate location in the MS monitoring system, which has a serious impact on the use of MS data to characterize fracture behavior. Therefore, the research on the inversion of water inrush seepage channel based on MS data needs further exploration and improvement.

Declaration of competing interest The authors (Yong Zhao, Tianhong Yang, Penghai Zhang, Haiyan Xu and Shuhong Wang) declared that they have no conflicts of interest to this work. Acknowledgments This work was supported by the National Key Research and Devel­ opment Program of China (2016YFC0801602), the National Natural Science Foundation of China (51604062, U1710253 and U1602232) and Development Program of Science and Technology in Liaoning �clav Province, China (2019JH2/10100035). The authors thank Dr. Va Vavry�cuk in the Czech Academy of Sciences for his guidance with the STRESSINVERSE package and PCA package, and Dr. Younes Fadakar Alghalandis for his considerable guidance with the ADFNE package.

Appendix A. Relationship between focal plane and failure mechanisms The eigenvectors e1 , e2 and e3 of MT give the directions of the principal strain axes T, B, and P.25 Based on their geometrical relationship with the normal n and slip vector v of focal plane, in which the T axis bisects the angle between n and v, and the T, B, and P are mutually orthogonal, vectors n and v can be derived from MTs by the following expressions39: 13

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 M1 M2 M2 M3 > > e1 þ e3 > > M1 M2 M2 M3 > :v¼ e1 e3 M1 M3 M1 M3

(A1)

To explain the characteristics of the normal n and slip vector v of focal plane with the change of failure mechanism, seven representative seismic sources were constructed according to the method proposed by Vavry�cuk,39 including two tensile failure seismic sources, three shear failure seismic sources and two compressive failure seismic sources as shown in Fig. A1b and Table A1. We used Eq. (A1) to solve the normal n and slip vector v of the seven seismic sources in Table A1, and the result is shown as the brown plane in the MT beachball in Fig. A1b. It can be seen that the normal n and slip vector v of focal planes of a seismic source are closely related to the component of MT and failure mechanism. When the seismic source is pure tensile or pure compressive failure source, the two focal planes of the seismic source are parallel. When the seismic source is pure shear failure source, the two focal planes of the seismic source are vertical. When the seismic source is shear - tensile or shear - compressive failure source, the two focal planes of the seismic source are at a certain angle, not parallel or vertical.

Fig. A1. (a) Transformation of failure mechanism with R value. (b) MT beachballs change with R value. The method for drawing full MT beachballs referred to a method proposed by Tape and Tape.57 The different colors of the MT beachballs represent different failure mechanisms, red for tensile failure, blue for shear failure and green for compressive failure. Table A1 Seven representative MTs. Event ID

m11

m22

m33

m12

m13

m23

R

Failure mechanism

S1 S2 S3 S4 S5 S6 S7

1.7500 0.2001 0.4991 0.8320 0.4991 0.2001 1.7500

1.2500 0.8791 0.7095 0.4899 0.7095 0.8791 1.2500

2.000 0.9900 0.7008 0.3420 0.7008 0.9900 2.000

0.4330 0.3181 0.2600 0.1841 0.2600 0.3181 0.4330

0.8660 0.0241 0.1688 0.3322 0.1688 0.0241 0.8660

0.5000 0.7093 0.6569 0.5754 0.6569 0.7093 0.5000

100% 43% 30% 0% 30% 43% 100%

Tensile Tensile Shear Shear Shear Compressive Compressive

Appendix B. Size and aperture of MS-derived fracture Size of MS-derived fracture plane. The source dimension is usually expressed as the radius of the MS-derived fracture and is related to the corner frequency of the P-wave or S-wave,25 which can be determined by the following equation. rs ¼

Kc β 0 2πfc

(B1)

where Kc is a constant that depends on the source model, Kc ¼ 2:01 for S-wave and Kc ¼ 1:32 for P-wave58; β0 is the S-wave velocity in the vicinity of the source; fc is the corner frequency of the P-wave or S-wave, which can be determined by displacement spectrum,26 as shown in Eq. (B2). The low frequency amplitude Ω0 and the corner frequency fc in Eq. (B1) can be determined by nonlinear fitting of the velocity spectrum of the MS waveform by 14

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180

Eq. (B2). π fr

Ωðf Þ ¼

Ω0 eVc Qc

(B2)

1 þ ðf =fc Þ2

where r is the distance between the source and the sensor, Vc is the velocity of P wave or S wave; Qc is the quality factor independent of frequency and represents the inherent attenuation effect of the medium. Aperture of MS-derived fracture. The fracture aperture change Δa can be deduced from the MT,40 as shown in Eq. (B3), and the new fracture aperture can be written as Eq. (B4) Δa ¼ u⋅cos θ ¼

Mkk ⋅4 f 2c ⋅cos θ 2 2 K c β0 ⋅ð3λ þ 2 Þn⋅v

π

(B3)

μ

(B4)

a ¼ a0 þ Δa

where θ is the angle between n and v, u is the displacement on the slip direction of the MS-derived fracture and a0 is the initial aperture of the preexisting fracture. It should be reiterated that Δa is the increase or decrease in the initial aperture a0 due to the reactivation of the fracture. Because no field measurement of the aperture of pre-existing fracture was conducted, the initial aperture a0 in Eq. (B4) was estimated to be 0.8 mm in this paper based on field measurements in the Refs. 59 and 60. Appendix C. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijrmms.2019.104180.

References

21 Kwiatek G, Goebel T, Dresen G. Seismic moment tensor and b-value variations over successive seismic cycles in laboratory stick–slip experiments. Geophys Res Lett. 2015; 41(16):5838–5846. 22 Nolen-Hoeksema RC, Ruff LJ. Moment tensor inversion of microseisms from the Bsand propped hydrofracture, M-site, Colorado. Tectonophysics. 2001;336(1):163–181. 23 Martínez-Garz� on P, Kwiatek G, Bohnhoff M, Dresen G. Impact of fluid injection on fracture reactivation at the Geysers geothermal field. J Geophys Res-Sol Ea. 2016;121: 7432–7449. 24 Trifu CI, Urbancic TI. Fracture coalescence as a mechanism for earthquakes: observations based on mining induced microseismicity. Tectonophysics. 1996;261(13):193–207. 25 Mendecki AJ. Seismic Monitoring in mines[M]. London: Chapman & Hall; 1997. 26 Aki K, Richards PG. Quantitative Seismology, Theory and methods[M]. New York: W. H. Freeman; 1980. 27 Young RP, Maxwell SC, Urbancic TI, Feignier B. Mining-induced microseismicity: monitoring and applications of imaging and source mechanism techniques. Pure Appl Geophys. 1992;139(3-4):697–719. 28 Young RP, Collins DS. Seismic studies of rock fracture at the underground research laboratory, Canada. Int J Rock Mech Min Sci. 2001;38(6):787–799. 29 Young RP, Collins DS, Reyes-Montes JM, Baker C. Quantification and interpretation of seismicity. Int J Rock Mech Min Sci. 2004;41(8):1317–1327. 30 Linzer LM. A relative moment tensor inversion technique applied to seismicity induced by mining. Rock Mech Rock Eng. 2005;38(2):81–104. 31 Vavry�cuk V. Moment tensor inversion based on the principal component analysis of waveforms: method and application to microearthquakes in West Bohemia, Czech Republic. Seismol Res Lett. 2017;88(5):1303–1315. 32 Gibbons SJ, Ringdal F. The detection of low magnitude seismic events using arraybased waveform correlation. Geophys J Int. 2010;165(1):149–166. 33 Hartigan JA, Wong MA. A K-means clustering algorithm. Appl Stat. 1979;28(1): 100–108. 34 Woodward K, Wesseloo J, Potvin Y. A spatially focused clustering methodology for mining seismicity. Eng Geol. 2017;232(1):104–113. 35 Shang XY, Li XB, Morales-Esteban A, Asencio-Cort�es G, Wang ZW. Data field-based kmeans clustering for spatio-temporal seismicity analysis and hazard assessment. Remote Sens. 2018;10(3):461. 36 Krzanowski WJ, Lai YT. A criterion for determining the number of groups in a data set using sum-of-squares clustering. Biometrics. 1988;44(1):23–34. 37 https://ww2.mathworks.cn/matlabcentral/fileexchange/51996-suspension-bri dge-picking-algorithm-sbpx?s_tid¼srchtitle. 38 Prugger AF, Gendzwill DJ. Microearthquake location: a nonlinear approach that makes use of a simplex stepping procedure. Bull Seismol Soc Am. 1988;78(2): 799–815. 39 Vavry�cuk V. Tensile earthquakes: theory, modeling, and inversion. J Geophys Res-Sol Ea. 2011;116(B12). 40 Zhao Y, Yang TH, Zhang PH, Xu HY, Zhou JR, Yu QL. Method for generating a discrete fracture network from microseismic data and its application in analyzing the permeability of rock masses: a case study. Rock Mech Rock Eng. 2019;52(9): 3133–3155. 41 Michael AJ. Determination of stress from slip data: faults and folds. J Geophys Res-Sol Ea. 1984;89(B13):11517–11526. 42 Vavry�cuk V. Iterative joint inversion for stress and fault orientations from focal mechanisms. Geophys J Int. 2014;199(1):69–77.

1 Xing ML, Li WP, Wang QQ, Yang DD. Risk prediction of roof bed-separation water inrush in a coal mine, China. Electron J Geotech Eng. 2015;20(1):301–312. 2 He KQ, Wang RL, Jiang WF. Groundwater inrush channel detection and curtain grouting of the gaoyang iron ore mine, China. Mine Water Environ. 2012;31(4): 297–306. 3 Zhou JR, Yang TH, Zhang PH, Xu T, Wei J. Formation process and mechanism of seepage channels around grout curtain from microseismic monitoring: a case study of Zhangmatun iron mine, China. Eng Geol. 2017;226:301–315. 4 Ma D, Miao XX, Bai HB, et al. Effect of mining on shear sidewall groundwater inrush hazard caused by seepage instability of the penetrated karst collapse pillar. Nat Hazards. 2016;82(1):73–93. 5 Sun WJ, Zhou WF, Jiao J. Hydrogeological classification and water inrush accidents in China’s coal mines. Mine Water Environ. 2016;35(2):214–220. 6 Wang Q, Wang XY, Hou QL. Geothermal water at a coal mine: from risk to resource. Mine Water Environ. 2016;35(3):294–301. 7 Zhang SH, Wang YW. Study of tunnel water detection through infrared acquisition technology based on 3D trend analysis. Adv Mater Res. 2014;838–841:1370–1380. 8 Xu ZM, Sun YJ, Dong QH, Zhang GW. Predicting the height of water-flow fractured zone during coal mining under the Xiaolangdi Reservoir. Int J Min Sci Technol. 2010; 20(3):434–438. 9 Zhai LJ. Comprehensive exploration technology of water disaster prevention and control in coal mining roof. Procedia Earth Planet Sci. 2011;3:303–310. 10 Tang SF, Tong MM, Hu JL. Study on prediction of water inrush in mine by microseismic technique. In: 2009 International Workshop on Intelligent Systems and Applications. 2009:1–4. 11 Li T, Mei TT, Sun XH, Lv YG, Sheng JQ, Cai M. A study on a water-inrush incident at Laohutai coalmine. Int J Rock Mech Min Sci. 2013;59:151–159. 12 Cheng S, Li SC, Li LP, Shi SS, Zhou ZQ, Wang J. Study on energy band characteristic of microseismic signals in water inrush channel. J Geophys Eng. 2018;15(5): 1826–1834. 13 Moriya H, Naoi M, Nakatani M, Van Aswegen G, Murakami O, et al. Delineation of large localized damage structures forming ahead of an active mining front by using advanced acoustic emission mapping techniques. Int J Rock Mech Min Sci. 2015;79: 157–165. 14 Candela T, Wassing B, Ter Heege J, Buijze L. How earthquakes are induced. Science. 2018;360(6389):598–600. 15 Mcgarr A. Some comparisons between mining-induced and laboratory earthquakes. Pure Appl Geophys. 1994;142(3-4):467–489. 16 Shigeishi M, Ohtsu M. Acoustic emission theory for moment tensor analysis. Constr Build Mater. 2001;15:311–319. � 17 Stierle E, Vavry�cuk V, Sílený J, Bohnhoff M. Resolution of non-double-couple components in the seismic moment tensor using regional networks-I: a synthetic case study. Geophys J Int. 2014;196(3):390–399. 18 Maxwell SC, Young RP. Seismic imaging of blast damage. Int J Rock Mech Min Sci. 1993;30(7):1435–1440. 19 Cheng GW, Li LC, Zhu WC, et al. Microseismic investigation of mining-induced brittle fault activation in a Chinese coal mine. Int J Rock Mech Min Sci. 2019;123:104096. 20 Ohtsu M. Acoustic emission theory for moment tensor analysis. Res Nondestruct Eval. 1995;6(3):169–184.

15

Y. Zhao et al.

International Journal of Rock Mechanics and Mining Sciences 126 (2020) 104180 53 Arifuzzaman S, Khan M. Fast parallel conversion of edge list to adjacency list for large-scale graphs. In: Proceedings of the Symposium on High Performance Computing. 2015:17–24. 54 Alghalandis YF. ADFNE: open source software for discrete fracture network engineering, two and three dimensional applications. Comput Geosci. 2017;102:1–11. 55 Zhao Y, Yang TH, Yu QL, Zhang PH. Dynamic reduction of rock mass mechanical parameters based on numerical simulation and microseismic data–A case study. Tunn Undergr Space Technol. 2019;83:437–451. 56 https://www.esgsolutions.com/. 57 Tape W, Tape C. A geometric setting for moment tensors. Geophys J Int. 2012;190(1): 476–498. 58 Madariaga R. Dynamics of an expanding circular fault. Bull Seismol Soc Am. 1976;66 (3):639–666. 59 Odonne F, L� ezin C, Massonnat G, Escadeillas G. The relationship between joint aperture, spacing distribution, vertical dimension and carbonate stratification: an example from the Kimmeridgian limestones of Pointe-du-Chay (France). J Struct Geol. 2007;29(5):746–758. 60 Larsen B, Grunnaleite I, Gudmundsson A. How fracture systems affect permeability development in shallow-water carbonate rocks: an example from the Gargano Peninsula, Italy. J Struct Geol. 2010;32(9):1212–1230.

43 Wallace RE. Geometry of shearing stress and relation to faulting. J Geol. 1951;59(2): 118–130. 44 Bott MHP. The mechanics of oblique slip faulting. Geol Mag. 1959;96(2):109–117. 45 Jia SQ, Eaton DW, Wong RCK. Stress inversion of shear-tensile focal mechanisms with application to hydraulic fracture monitoring. Geophys J Int. 2018;215(1): 546–563. 46 Huseby O, Thovert JF, Adler PM. Geometry and topology of fracture systems. J Phys A Gen Phys. 1997;30(5):1415. 47 Alghalandis YF, Dowd PA, Xu C. The RANSAC method for generating fracture networks from micro-seismic event data. Math Geosci. 2013;45(2):207–224. 48 Berkowitz B. Characterizing flow and transport in fractured geological media: a review. Adv Water Resour. 2002;25(8-12):861–884. 49 Sanderson DJ, Nixon CW. Topology, connectivity and percolation in fracture networks. J Struct Geol. 2018;115:167–177. 50 Priest SD. Discontinuity Analysis for Rock Engineering. Springer Science & Business Media; 2012. 51 Gross JL, Yellen J. Handbook of Graph Theory. London: CRC press; 2004. 52 Bondy JA, Murty USR. Graph Theory with Applications. London: Macmillan press; 1976.

16