- Email: [email protected]
Unloading-induced rock fracture activation and maximum seismic moment prediction
Engineering Geology 262 (2019) 105352
Contents lists available at ScienceDirect
Engineering Geology journal homepage: www.elsevier.com/locate/enggeo
Unloading-induced rock fracture activation and maximum seismic moment prediction
T
Yinlin Jia, Wei Wua,*, Zhihong Zhaob a b
School of Civil and Environmental Engineering, Nanyang Technological University, 639798, Singapore Department of Civil Engineering, Tsinghua University, Beijing, 100084, China
ARTICLE INFO
ABSTRACT
Keywords: Rock fractures Frictional instability Induced seismicity Seismic moment
Hundreds of anthropogenic earthquakes have recently occurred worldwide due to underground space creation and energy extraction. The mechanism behind the human-induced geohazards is most likely associated with the reduction of normal stress on pre-existing fractures and faults. This study reports a series of laboratory experiments to investigate the mechanism of unloading-induced fracture activation, and proposes a simple approach to predict the maximum seismic moment for a critically stressed fracture. The unloading-driven shear test results exhibit that the unloading process induces the stress states of the sawcut and natural fractures to approach the Mohr-Coulomb failure envelope, and the normal stress unloading rate influences the peak slip rate. The fracture instability is dependent on the relationship between the stiffness of the system and the slip weakening rate of the fracture, and the shear dilation mainly occurs after the fracture activation. The test results also show that the critical shear stress of the sawcut fracture during the unloading-driven shear test is approximately equal to the residual shear strength after the displacement-driven fracture slip. This relationship inspires us to develop a new approach to estimate the maximum seismic moment. Our data demonstrate that the maximum seismic moments for both the fractures obtained from the unloading-driven shear tests are all below the upper limit lines, indicating that the proposed approach is reasonable. The uncertainty analysis shows that the accurate estimation of fault size can improve the maximum seismic moment prediction.
1. Introduction
principle, fracture activation occurs when the shear stress on a rock fracture exceeds the friction coefficient multiplied by the effective normal stress, which is defined as the difference between the normal stress and pore pressure. The fracture instability is accompanied by friction loss, which can be described by the slip weakening friction law based on a spring-slider system with one degree of freedom (Dieterich, 1978; Rice, 1983). The fracture stably slides when the stiffness of the spring-slider system is larger than the slip weakening rate of the fracture, otherwise the fracture exhibits an unstable motion. Besides the classic models, in-situ and ex-situ investigations have been carried out to explore the mysteries behind the induced earthquakes related to underground mass extraction. Microseismic monitoring is commonly used to observe the occurrence of fracture activation based on the spatial and temporal evolution of microseismic events (Rutlegde et al., 1998; Xu et al., 2016). A series of unloading-induced direct-shear tests is conducted to study the unloading-induced fracture instability (Wu et al., 2014), and demonstrates that the fracture instability is accompanied by the reduction of both normal and shear stresses on the fracture (Wu et al., 2017). The fracture instability is also associated
Underground resource development is expected to meet our everincreasing demand for economic development without disturbing the environment. However, several anthropogenic earthquakes have recently occurred owing to underground space creation and energy extraction. The extraction of mass (e.g., rock and fluid) from the subsurface perturbs the frictional equilibrium of pre-existing fractures and faults, resulting in the human-induced geohazards (Foulger et al., 2018). Excavation-induced seismicity has been reported during the constructions of the Gotthard Base Tunnel in Switzerland (Husen et al., 2013) and the Jinping II hydropower station in China (Zhang et al., 2013). Induced earthquakes associated with oil and gas extraction also occurred at the Wilmington oil field in the United States (Segall, 1989), the Gazli gas reservoir in Uzbekistan (Simpson and Leith, 1985), and the Groningen gas field in Netherlands (Vlek, 2018). These geohazards have caused incalculable economic losses and even led to project abandonment. According to the Mohr-Coulomb failure criterion and effective stress
⁎
Corresponding author. E-mail addresses: [email protected] (Y. Ji), [email protected] (W. Wu), [email protected] (Z. Zhao).
https://doi.org/10.1016/j.enggeo.2019.105352 Received 16 June 2019; Received in revised form 11 October 2019; Accepted 18 October 2019 Available online 21 October 2019 0013-7952/ © 2019 Elsevier B.V. All rights reserved.
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
with the change of principal stress orientation with respect to fracture orientation (Duan et al., 2019b). The volumetric contraction of reservoir rock due to gas extraction can change the local stress on the fracture (Segall et al., 1994) and induce the readjustment of crustal stress (Grasso, 1992). The occurrence of the seismic events due to gas extraction can also be correlated with the stiffness between seal- and reservoir rocks, pressure drop, and fracture density (van Eijs et al., 2006). Understanding the mechanism of unloading-induced fracture activation allows us to manage the risk of the anthropogenic geohazards. Numerical models combined with microseismic data have been used to evaluate the fracture instability during underground excavation (Xu et al., 2015). A seismic risk analysis using a risk matrix approach indicates that the Groningen gas field has a much larger seismic risk than the other onshore gas fields in the Netherlands (van Thienen-Visser et al., 2018). The total seismic moment for a fractured rock mass can be estimated by the modulus of rigidity and the volume of extracted rock in mining engineering (McGarr, 1976), and by the modulus of rigidity, the mass of removed fluid, the density of upper crust, and the fraction of seismogenic upper crust in petroleum engineering (McGarr, 1991). The seismic moment for a rectangular fracture can also be assessed based on the stress drop (Stein and Wysession, 2009). This study aims to investigate the mechanism of unloading-induced fracture activation and to forecast the maximum seismic moment for a critically stressed fracture. We perform the unloading-driven shear tests on sawcut and natural fractures to study the unloading-induced fracture instability, and discuss the mechanism of the fracture activation, in terms of slip weakening friction, fracture shear dilation, and peak slip rate. A simple approach is proposed to predict the maximum seismic moment and verified using the test results. We also conduct the uncertainty analysis to quantify the variability of the maximum seismic moment, and discuss the possible application of this approach in engineering projects.
porosity, and water content are 2660 kg/m3, 0.26%, and 0.07%, respectively. The Young’s modulus, Poisson’s ratio, cohesion, and friction angle are 74 GPa, 0.15, 38 MPa, and 54°, respectively. The physical and mechanical properties of the granite were measured based on the test standards suggested by International Society for Rock Mechanics (ISRM, 2015). The intrinsic permeability of the granite matrix is 1.3 μD measured based on the pressure pulse decay method (Brace et al., 1968; Zoback and Byerlee, 1975). Sawcut and natural fractures were prepared from 50-mm-diameter granite cores (Fig. 1a). The sawcut fracture was obtained by cutting a 100-mm-long core sample at 30° to the core axis using a diamond saw. The fracture surface roughness was controlled by fine sandpaper with 25.6 μm particle size. The natural fracture was produced by compressing a 124-mm-long core sample containing a naturally formed but unopened fracture at 27° to the core axis until the sample purely failed along the fracture. The fracture surfaces were imaged before and after the shear tests using a structured-light 3D scanner. Two Teflon layers were used to secure the sample on the coreholders and to isolate the sample from the confining oil. Two 2-mm-diameter boreholes were drilled at the sample ends to facilitate fluid flow from the coreholders to the fracture. Two 130-μm-thickness and 0.45-μm-pore size filter papers were inserted between the coreholders and sample ends to prevent sheared-off particles from entering the pore system. The sample assembly was installed in the MTS 815 rock mechanics test system (Fig. 1b). Silicon oil and distilled water were used as the confining and pore fluids, respectively. The normal and shear stresses on the inclined fracture were servo-controlled by the axial stress (σ1) and confining pressure (σ3), and corrected by considering the reduction of fracture area and the deformation of Teflon layers (Kohli and Zoback, 2013). The shear stress (τ) and normal stress (σn) can be calculated as:
2. Experimental method
where Ψ is the fracture inclination angle with respect to the core axis, σ1-σ3 is the deviator stress directly measured by an in-vessel load cell.
=(
1
=
3
n
2.1. Experimental setup
3)sin
+(
1
cos
(1)
2 3)cos
(2)
2.2. Experimental procedure
Bukit Timah granite, underlying the central region of Singapore Island, is of extreme high-quality and suitable for cavern development (Zhao et al., 1999). The granite is composed of 62% feldspar, 32% quartz, 5% biotite, and 1% hornblende. The average bulk density,
We carried out a series of displacement-driven and unloadingdriven shear tests. The displacement-driven shear tests were used to obtain the shear strength of the fracture, and the unloading-driven Fig. 1. (a) Sawcut (left) and natural (right) fractures, and (b) schematic drawing of sample assembly. Normal (σn) and shear (τ) stresses on the fracture with inclination angle (Ψ) with respect to the core axis are servo-controlled by axial stress (σ1) and confining pressure (σ3). Pore pressure is applied on the fracture through two boreholes.
2
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
shear test SDD1 and the succeeding unloading-driven shear test SDU1 as examples. First, under 5 MPa initial normal stress, the displacementdriven shear test was carried out at an axial displacement rate of 1 μm/s until the fracture reached its shear strength. Second, we fixed the shear stress as 85% of the shear strength, and assumed the fracture close to a critical stress state. Finally, the normal stress was reduced at an unloading rate of 0.01 MPa/s until the fracture was activated. For the dry sawcut fracture, seven unloading-driven shear tests were performed under different combinations of initial normal stress (e.g., 5, 7.5, and 10 MPa), shear stress level (e.g., 85%, 90%, and 95% of the shear strength), and normal stress unloading rate (e.g., 0.005, 0.01, and 0.05 MPa/s). Under 11 MPa initial normal stress, the drained and undrained conditions were considered during the unloading-driven shear tests on the saturated sawcut fracture at 0.01 MPa/s unloading rate. The pore pressure was maintained at 1 MPa under the drained condition, while the pore volume was constant with 1 MPa initial pore pressure under the undrained condition. For the natural fracture, the displacement-driven shear tests were conducted under the low normal stresses of 1, 2, and 3 MPa to minimize asperity damage and to construct the failure envelope, which was subsequently used to calculate the shear strength under high normal stresses. A critical stress state of the fracture was then simulated with the shear stress equal to 85% of the calculated shear strength. Under 10 MPa initial normal stress, three unloading-driven shear tests on the dry natural fracture were carried out at different unloading rates (0.005, 0.01, and 0.05 MPa/s). Under 11 MPa initial normal stress and 1 MPa pore pressure, the drained and undrained conditions were also studied during the unloading-driven shear tests on the saturated natural fracture at 0.01 MPa/s unloading rate. 3. Results and discussion 3.1. Experimental results The Mohr-Coulomb failure envelopes for the sawcut and natural fractures under the dry and saturated conditions are constructed based on the displacement-driven shear strengths (Fig. 3). The friction coefficients of the dry and saturated sawcut fractures are 0.83 and 0.68, respectively. The cohesion of the sawcut fractures is negligible. The friction coefficient of the natural fracture also reduces from 0.86 under the dry condition to 0.77 under the saturated condition. The cohesions of the dry and saturated natural fractures are 0.07 and 0.12 MPa, respectively. The reduction of friction coefficient is due to the lubrication effect of pore water (Cornelio et al., 2019). Fig. 3 also shows that the stress states of both the fractures are below the failure envelopes at critically stressed states and close to the failure envelopes at the onset of fracture activation. The initial normal and shear stresses significantly influence the change of stress state, whereas the effect of normal stress unloading rate is relatively minor. The drained and undrained conditions also slightly affect the unloading-driven instability of both the fractures. The sawcut and natural fractures exhibit dynamic slip and quasidynamic sliding, respectively. The dynamic slip of the sawcut fracture induced by normal stress reduction is characterized by a sudden increase in shear displacement and an abrupt drop of shear stress (Fig. 4a). After that, the fracture arrests, and the shear stress reverts to its initial value. The slip rate and acceleration are the first- and secondtime derivatives of shear displacement, respectively, which instantaneously jump during the dynamic slip. The unloading-induced dynamic slip of the sawcut fracture resembles the coseismic slip during the stick-slip cycles of natural faults, in terms of stress drop and
Fig. 2. (a) Experimental principle, and (b) typical stress paths of displacementdriven and unloading-driven shear tests. In the unloading-driven shear test, the fracture is initially subjected to a critical shear stress (τc) determined by the shear strength (τs) obtained from the displacement-driven shear test. Δτ is equal to τs minus τc. c and μ are the cohesion and friction coefficient of the fracture, respectively.
shear tests were conducted to simulate the unloading-induced activation of the fracture initially approaching a critical stress state. During the unloading-driven shear tests, the Mohr circle moves towards the Mohr-Coulomb failure envelope due to the reduction of the normal stress on the fracture, which finally results in the fracture activation (Fig. 2a). During the unloading-driven shear tests on the sawcut fracture (Table 1), the effects of initial normal stress, critical shear stress, normal stress unloading rate, and fluid drainage were considered to investigate the mechanism of unloading-induced fracture activation. Fig. 2b demonstrates the experimental procedure, taking the displacement-driven
3
4
σn = 11, P = 1
SSD3
Stable sliding Stable sliding Stable sliding
Stable sliding Stable sliding Stable sliding
Saturated natural fracture NSD1 σn = 1, P = 0 NSD2 σn = 2, P = 0 NSD3 σn = 3, P = 0
5.62
–
–
–
9.91
5.73
3.60
– – –
– – –
Residual shear strength (MPa)
Dry natural fracture NDD1 σn = 1 NDD2 σn = 2 NDD3 σn = 3
Stick-slip
Stick-slip
Stick-slip
σn = 12.5
SDD4
σn = 8.5, P = 1
Stick-slip
σn = 10
SDD3
SSD2
Stick-slip
σn = 7.5
SDD2
Stick-slip
Stick-slip
Dry sawcut fracture SDD1 σn = 5
Saturated sawcut fracture SSD1 σn = 6, P = 1
Failure mode
Displacement-driven shear test
NPD NPU
5.87
5.87
11.19
10.60
10.01
10.01
10.01
5.69
3.67
Critical shear stress (MPa)
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Failure mode
σni = 11, P = 1, τc = 85%τs, Ur = 0.01, drain σni = 11, P = 1, τc = 85%τs, Ur = 0.01, undrain
σni = 10, τc = 85%τs, Ur = 0.005 σni = 10, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.05
σni = 11, P = 1, τc = 85%τs, Ur = 0.01, drain σni = 11, P = 1, τc = 85%τs, Ur = 0.01, undrain
σni = 5, τc = 85%τs, Ur = 0.01 σni = 7.5, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.005 σni = 10, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.05 σni = 10, τc = 90%τs, Ur = 0.01 σni = 10, τc = 95%τs, Ur = 0.01
NDU1 NDU2 NDU3
SPU
SPD
SDU7
SDU6
SDU5
SDU4
SDU3
SDU2
SDU1
Unloading-driven shear test
7.00 7.00
6.70 6.70 6.70
7.98
8.38
8.89
8.13
7.40
7.59
7.54
5.99
4.06
Normal stress at failure (MPa)
0.15
0.13
0.16
0.15
0.45
0.16
0.16
0.06
0.07
Quasi-dynamic sliding Quasi-dynamic sliding
8.63 8.65
8.17 7.96 7.66
Shear displacement (mm)
Quasi-dynamic sliding Quasi-dynamic sliding Quasi-dynamic sliding
1.61
1.47
1.65
1.86
2.48
2.14
2.01
0.83
0.79
Shear stress drop (MPa)
– –
– – –
2.25
2.11
2.62
2.47
3.68
2.56
2.32
1.14
1.05
– –
– – –
Peak slip rate (mm/s)
0.22 0.23
0.24 0.25 0.29
110.51
156.72
146.97
196.97
178.25
143.84
167.49
58.73
82.14
Peak acceleration (mm/s2)
39.70 34.53
40.92 57.16 35.50
Table 1 Summary of experimental conditions and results. A rock fracture is subjected to normal stress (σn in MPa) and pore pressure (P in MPa) in the displacement-driven shear test, and initial normal stress (σni in MPa), critical shear stress (τc in MPa) equivalent to a certain percentage of displacement-driven shear strength (τs in MPa), pore pressure (P in MPa) and normal stress unloading rate (Ur in MPa/s) in the unloading-driven shear test.
Y. Ji, et al.
Engineering Geology 262 (2019) 105352
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 3. Mohr-Coulomb failure envelopes and stress states before and at the onset of unloading-induced fracture activation for (a) dry sawcut fracture, (b) saturated sawcut fracture, (c) dry natural fracture, and (d) saturated natural fracture.
displacement jump (Karner and Marone, 2000). For the natural fracture, the shear stress is constant during the quasi-dynamic sliding (Fig. 4b). The shear displacement exponentially increases, accompanied by the slight changes of slip rate and acceleration. The peak slip rate of the natural fracture is similar to that of the slow slip events of natural faults (Tinti et al., 2016). Fig. 5 presents the topographic contours of the sawcut and natural fractures before and after the displacement-driven and unloadingdriven shear tests. The lowest point in the borehole is unaffected during the shear process and considered as the zero datum, which makes the contours comparable. We calculate the root mean square (RMS) asperity height to quantify the evolution of surface roughness (Gadelmawla et al., 2002). The RMS asperity height of the sawcut fracture reduces from 1.31 mm before the shear tests to 0.65 mm after
the tests. The RMS asperity height of the natural fracture also decreases from 2.32 mm to 2.06 mm. The reductions of RMS asperity height influence the stiffness of the system, which controls the fracture behaviors. 3.2. Mechanism of unloading-induced fracture activation Our results show that the unloading-induced fracture activation is dependent on fracture surface roughness, and can be explained by a spring-slider model with one degree of freedom (Dieterich, 1978). The slip weakening friction law describes that the friction coefficient decreases with larger shear displacement and becomes nearly constant after the fracture experiences a slip weakening distance (Rice, 1983). The fracture activation is associated with the relationship between the
5
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 4. (a) Unloading-induced dynamic slip of sawcut fracture (test SDU4), and (b) unloading-induced quasi-dynamic sliding of natural fracture (test NDU3).
stiffness of the spring-slider system and the slip weakening rate of the fracture. In this study, the stiffness of the system can be considered as the combined rigidity of the test system and the fracture, which is estimated as the slope of the linear elastic portion of the shear stress
versus axial displacement curve obtained from the displacement-driven shear test. The slip weakening rate of the fracture is estimated as the slope of the linear portion of the shear stress versus shear displacement curve obtained from the unloading-driven shear test. Additionally, we
6
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 5. Topographic contours of (a) sawcut fracture and (b) natural fracture before and after displacement-driven and unloading-driven shear tests.
use the fracture strength to describe the change in shear resistance of the fracture, which is determined as the product of the friction coefficient and the effective normal stress. When the stiffness of the system is lower than the slip weakening rate of the sawcut fracture (Fig. 6a), the
fracture experiences the dynamic slip. After that, the fracture strength recovers due to the reformation and strengthening of asperity contacts (Marone, 1998). The natural fracture slides in a quasi-static manner due to the stiffness of the system greater than the slip weakening rate of the
7
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
fluid volume indicates the normal dilation of the fracture. The hydraulic apertures of both the sawcut and natural fractures slightly increase before the fracture activation, but dramatically rise with larger shear displacement (Fig. 7). The results indicate that both the fractures experience the remarkable dilation during the shear process, which is due to one asperity overriding another and pore space expansion among sheared-off particles (Morgan and Boettcher, 1999; Li et al., 2019). The slight decline in shear displacement before the onset of fracture activation is presumably due to the elastic unloading of rock blocks (Duan et al., 2019a). Fig. 8 exhibits that the peak slip rates of both the sawcut and natural fractures as a function of normal stress unloading rate. The peak slip rate of the sawcut fracture is the maximum rate value of the unloadinginduced dynamic slip, and that of the natural fracture is the rate value at 0.25 mm shear displacement, at which the stable sliding is achieved. The peak slip rate under the same stress state increases with higher unloading rate owing to a larger amount of strain energy released at the onset of fracture activation (Wu et al., 2017). The peak slip rate under the undrained condition exceeds that under the drained condition, because the pore pressure at fracture slip is reduced under the undrained condition and the amount of released energy increases with larger effective normal stress (Lei et al., 2016). 4. Prediction of maximum seismic moment The similarity between the unloading-induced dynamic slip of the sawcut fracture and the coseismic slip of natural faults allows us to predict the maximum seismic moment due to normal stress reduction. The seismic moment describes earthquake size and is determined as (Aki, 1966): (4)
M0 = GAd
where A is the fracture area, d is the shear displacement, and G is the combined rigidity of the test system and the fracture. To estimate the combined rigidity of the test system and the sawcut fracture with 3.9 × 10−9 km2 surface area, we derive the combined rigidity in each dynamic slip case based on the empirical relationships among the moment magnitude (Mw), fracture area, and seismic moment (Wells and Coppersmith, 1994; Kanamori and Brodsky, 2001):
Mw = 4.07+0.98log(A)
Mw =
(6)
We then constrain the combined rigidity of the seven dynamic slip cases using the least squares method and obtain a best-fit value of 1.2 GPa (Fig. 9a). The derived rigidity is validated based on the relationships among the earthquake scaling parameters. For a circular fracture with a radius (R), the shear displacement (d) and stress drop (Δτ) can be expressed as (Stein and Wysession, 2009):
Fig. 6. Analysis of unloading-induced instability of (a) sawcut fracture (test SDU4) and (b) natural fracture (test NDU3).
natural fracture (Leeman et al., 2016). The stress imbalance between the maintained shear stress and decreasing fracture strength leads to the increase in slip rate (Fig. 6b). The shear dilation of rock fractures also influences the unloadinginduced fracture activation. During the unloading-driven shear tests on the sawcut and natural fractures, the undrained condition with constant pore volume provides an avenue to estimate the hydraulic aperture of the fracture, which can represent the mechanical aperture (Fang et al., 2017). The change in pore fluid volume (ΔV) can be related to the change in hydraulic aperture (Δe):
V=A e
log M0 -6.06 1.5
(5)
d
cM0/(G (2R)2) =
= cM0/(2R)3 =
7M0 8GR2
7 M0 16 R3
(7) (8)
where c is the shape factor of a circular fracture and equal to 3.5. The seven dynamic slip cases are plotted in Fig. 9b. The source dimension is considered as a fracture radius of 0.075 m, which is the average length of the major and minor axes of the oval-shaped fracture. The seven cases fall within a stress drop range of 0.5–3 MPa and a shear displacement range of 50–500 μm, which are consistent with the results directly measured from the unloading-driven shear tests (Table 1). Hence, the combined rigidity is reasonable to determine the seismic
(3)
where A is the fracture area. The water compressibility is neglected at constant temperature. For a given fracture area, the increase in pore
8
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 7. Shear stress, shear displacement, pore pressure, and hydraulic aperture change during the unloading-driven shear tests on (a) sawcut fracture (test SPD) and (b) natural fracture (test NPD).
moments for the sawcut and natural fractures. Table 1 also shows that the critical shear stress of the sawcut fracture during the unloading-driven shear test (τc) is approximately equal to the residual shear strength after the displacement-driven fracture slip, which is equal to the shear strength (τs) minus stress drop (Δτ) (see
Fig. 2b). Thus, the normal stress reduction required to activate the fracture can be calculated as: n
9
=
µ
=
s
c
µ
(9)
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
1978). The standard deviations of both the parameters are considered as one-tenth of the average of the maximum and minimum values. The linear correlation between the two parameters, known as Pearson correlation coefficient, is 0.0002. The first-order Sobol indices, which quantify the variance of the increment rate from one of the parameters, are 0.92 for the fracture radius and 0.07 for the friction coefficient, respectively. The second-order Sobol index, which measures the joint effect of the two parameters, is 0.01. The results show that the firstorder index of fracture radius is much larger than the other indices, and suggest that the successful prediction of the maximum seismic moment is highly dependent on the accurate evaluation of fault size. In engineering applications, the success of predicting the maximum seismic moment for unloading-induced seismicity using Eq. (10) relies heavily on the determination of input parameters, i.e., the friction coefficient and size of pre-existing faults, and the normal stress reduction due to underground activities. The friction coefficient can be assumed as 0.6, which is a typical value obtained from laboratory experiments and has been used to forecast the maximum seismic moment for injection-induced seismicity (Zoback, 2007; McGarr, 2014). The fault size can be estimated from the moment magnitude of microseismic events by means of the seismological scaling factors (Kuang et al., 2017; Stein and Wysession, 2009). The change in normal stress is related to the change in volume (e.g., removed rock and extracted fluid) (McGarr, 1976). Therefore, this study provides a possible solution to predict the maximum seismic moment for field-scale seismic events associated with underground space creation and energy extraction. The proposed approach can be further verified by field investigations with reasonable input parameters (e.g., Calò et al., 2014; Guglielmi et al., 2015) and numerical simulations with efficient computational methods (e.g., Areias et al., 2018; Rabczuk and Belytschko, 2004). Additionally, the proposed approach is established based on the uniform reduction of normal stress on pre-existing faults. Underground activities may trigger the instability of large-scale faults through local disturbances, e.g., rock burst (Manouchehrian and Cai, 2018; Duan et al., 2019a) and stress rotation (Duan et al., 2019b). Our approach may not be suitable for this case, in which the maximum seismic moment is determined not only by the coefficient of friction, fault size, and normal stress reduction, but also by fault structure and rupture propagation.
Fig. 8. Peak slip rates of sawcut and natural fractures as a function of normal stress unloading rate.
where μ is the friction coefficient of the fracture and obtained from Fig. 3. According to Eqs. (8) and (9), the maximum seismic moment can be estimated based on the normal stress reduction at fracture failure:
M0 =
16µR3 7
n
(10)
The unloading-induced moment release of the sawcut fracture is mostly seismic (Fig. 4a). The combined rigidity of the test system and the fracture, the fracture area, and the shear displacement at fracture failure are thus used to calculate the maximum seismic moment. The friction coefficients of the dry and saturated sawcut fractures are 0.83 and 0.68, respectively, and used to plot the two upper limit lines (Eq. (10)). Based on the unloading-induced shear test results, the maximum seismic moment generally increases with larger normal stress reduction (Fig. 10a). The seismic moments for the dry and saturated fracture cases are all below the corresponding upper limit lines, which verifies the prediction of the maximum seismic moment. For the natural fracture, the unloading-induced moment release is progressively aseismic (Fig. 4b). The seismic moment is thus dependent on shear displacement and calculated as the product of the combined rigidity, the fracture area, and the cumulative shear displacement. The friction coefficient used for the upper limit line is 0.82, which is the average values of the dry and saturated natural fractures. Fig. 10b demonstrates that the unloading-induced cumulative seismic moment continuously increases with larger normal stress reduction until it reaches the upper limit line. This figure also indicates that Eq. (10) is valid to forecast the seismic moment. The variability of the maximum seismic moment due to normal stress reduction is determined by the fracture radius and friction coefficient (Eq. (10)). We carry out the uncertainty analysis to quantify the variability based on the Sobol’s method (Sobol, 1993; Vu-Bac et al., 2016; Hamdia et al., 2017). The increment rate of the maximum seismic moment as a function of normal stress reduction (i.e., 16μR3/7) is considered in this analysis. The statistical properties of the fracture radius and friction coefficient are listed in Table 2. The fault radius is assumed to follow a normal distribution with the minimum and maximum values of 10−6 m and 106 m, respectively, which represent natural faults in a size range from micrometers to thousands of kilometers. The friction coefficient also follows a normal distribution with the minimum and maximum values of 0.6 and 0.85, respectively (Byerlee,
5. Conclusions We investigate the unloading-induced instability of sawcut and natural fractures based on the unloading-driven shear tests. The test data show that the unloading-induced fracture behaviors obey the Mohr-Coulomb failure criterion, and are strongly affected by fracture surface roughness. The fracture instability is controlled by the relationship between the stiffness of the system and the slip weakening rate of the fracture, and the shear dilation mainly occurs after the fracture activation. The peak slip rates of both the fractures increase with higher normal stress unloading rate. The drained and undrained conditions play a minor role in the fracture instability. The maximum seismic moment is estimated based on the fact that the residual shear strength after the displacement-driven fracture slip is comparable to the critical shear stress during the unloading-driven fracture slip. The maximum seismic moment for the sawcut fracture obtained from the unloading-driven shear test can be bounded by the upper limit line. The time-dependent seismic moment for the natural fracture increases with larger normal stress reduction, and finally approaches the upper limit line. The proposed approach for the maximum seismic moment prediction can be used in engineering projects by reasonably determining the input parameters, and the prediction accuracy can be improved by accurately estimating fault size.
10
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 9. (a) Recovery of the combined rigidity of the test system and the sawcut fracture based on the empirical relationship between the moment magnitude and source area of natural earthquakes. (b) Comparison between the unloading-induced dynamic slips of the sawcut fracture and natural earthquakes in terms of stress drop and slip displacement. s denotes standard deviation.
11
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 10. Seismic moment as a function of normal stress reduction for (a) sawcut fracture and (b) natural fracture.
12
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Table 2 Statistical properties of fracture radius and friction coefficient. Parameter
Distribution type
Maximum value
Minimum value
Standard deviation
Friction coefficient Fault size (m)
Normal Normal
0.85 106
0.60 10−6
0.07 5 × 104
fracturing. Geophysics 82 (1), KS1–KS11. Scuderi, M.M., Leeman, J.R., Saffer, D.M., Marone, C., 2016. Laboratory observations of slow earthquakes and the spectrum of tectonic fault slip modes. Nature Comm 7, 11104. Lei, X., Funatsu, T., Ma, S., Liu, L., 2016. A laboratory coustic emission experiment and numerical simulation of rock fracture driven by a high-pressure fluid source. J. Rock Mech. Geotech. Eng. 8, 27–34. Li, Y.C., Wu, W., Tang, C.A., Liu, B., 2019. Predicting the shear characteristics of rock joints with asperity degradation and debris backfilling under cyclic loading conditions. Int. J. Rock Mech. Min. Sci. 120, 108–118. Marone, C., 1998. Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Planet. Sci. 26 (1), 643–696. McGarr, A., 1976. Seismic moment and volume changes. J. Geophys. Res. 81 (8), 1487–1494. McGarr, A., 1991. On a possible connection between three major earthquakes in California and oil production. Bull. Seismol. Soc. Am. 81 (3), 948–970. McGarr, A., 2014. Maximum magnitude earthquakes induced by fluid injection. J. Geophys. Res. Solid Earth 119 (2), 1008–1019. Morgan, J.K., Boettcher, M.S., 1999. Numerical simulations of granular shear zones using the distinct element method, 1. Shear zone kinematics and the micromechanics of localization. J. Geophys. Res. 104 (B2), 2703–2719. Rabczuk, T., Belytschko, T., 2004. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int. J. Numer. Methods Eng. 61 (13), 2316–2343. Rice, J.R., 1983. Constitutive relations for fault slip and earthquake instabilities. Instabilities in Continuous Media. Springer. Rutlegde, J.T., Phillips, W.S., Schuessler, B., 1998. Reservoir characterisation using oilproduction-induced microseismicity, Clinton County, Kentucky. Tectnophysics 289 (1–3), 129–152. Segall, P., 1989. Earthquakes triggered by fluid extraction. Geology 17, 942–946. Segall, P., Grasso, J.R., Mossop, A., 1994. Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France. J. Geophys. Res. 99 (B8), 15423–15438. Simpson, D.W., Leith, W., 1985. The 1976 and 1984 Gazli, USSR, earthquakes – were they induced? Bull. Seismol. Soc. Am. 75, 1465–1468. Sobol, I., 1993. Sensitivity analysis for non-linear mathematical models. Reliab. Eng. Syst. Safe. 52, 1–17. Stein, S., Wysession, M., 2009. An Introduction to Seismology, Earthquakes, and Earth Structure. John Wiley & Sons. Tinti, E., Scuderi, M.M., Scognamiglio, L., Di Stefano, G., Marone, C., Collettini, C., 2016. On the evolution of elastic properties during laboratory stick-slip experiments spanning the transition from slow slip to dynamic rupture. J. Geophys. Res. Solid Earth 121 (8569), 8594. van Eijs, R.M.H.E., Mulders, F.M.M., Nepveu, M., Kenter, C.J., Scheffers, B.C., 2006. Correlation between hydrocarbon reservoir properties and induced seismicity in the Netherlands. Eng. Geol. 84, 99–111. van Thienen-Visser, K., Roholl, J.A., van Kempen, B.M.M., Muntendan-Bos, A.G., 2018. Categorizing seismic risk for the onshore gas fields in the Netherlands. Eng. Geol. 237, 198–207. Vu-Bac, N., Lahmer, T., Zhuang, X., Nguyen-Thoi, T., Rabczuk, T., 2016. A software framework for probabilistic sensitivity analysis for computationally expensive models. Adv. Eng. Softw. 100, 19–31. Vlek, C., 2018. Induced earthquakes from long-term gas extraction in Groningen, the Netherland: statistical analysis and prognosis for acceptable-risk regulation. Risk Anal. 38 (7), 1455–1473. Wells, D.L., Coppersmith, K.J., 1994. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bull. Seismol. Soc. Am. 84 (4), 974–1002. Wu, W., Zhao, Z., Duan, K., 2017. Unloading-induced instability of a simulated granular fault and implications for excavation-induced seismicity. Tunn. Undergr. Space Technol. 63, 154–161. Wu, W., Zou, Y., Li, X., Zhao, J., 2014. An unload-induced direct-shear model for granular gouge friction in rock discontinuities. Rev. Sci. Instrum. 85 (9), 093902. Xu, N.W., Li, T.B., Dai, F., Li, B., Zhu, Y.G., Yang, D.S., 2015. Microseismic Monitoring and stability evaluation for the large scale underground caverns at the Houziyan hydropower station in Southwest China. Eng. Geol. 188, 48–67. Xu, N.W., Li, T.B., Dai, F., Zhang, R., Tang, C.A., Tang, L.X., 2016. Microseismic monitoring of strainburst activities in deep tunnels at the Jinping II hydropower station, China. Rock Mech. Rock Eng. 49 (3), 981–1000. Zhang, C., Feng, X.T., Zhou, H., Qiu, S., Wu, W., 2013. Rockmass damage development following two extremely intense rockbursts in deep tunnels at Jinping II hydropower station, southwestern China. Bull. Eng. Geol. Environ. 72 (2), 237–247. Zhao, J., Zhou, Y.X., Hefny, A.M., Cai, J.G., Chen, S.G., Li, H.B., Liu, J.F., Jain, M., Foo, S.T., Seah, C.C., 1999. Rock dynamics research related to cavern development for ammunition storage. Tunn. Undergr. Space Technol. 14 (4), 513–526. Zoback, M.D., 2007. Reservoir Geomechanics. University Press, Cambridge. Zoback, M.D., Byerlee, J.D., 1975. The effect of microcrack dilatancy on the permeability of Westerly granite. J. Geophys. Res. 80 (5), 752–755.
Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgement Wei Wu gratefully acknowledges the support of Start-Up Grant from Nanyang Technological University, Singapore. References Aki, K., 1966. Generation and propagation of G waves from the Niigata earthquake of June 16, 1964, 2, estimation of earthquake moment, release energy, and stress-strain drop from the G wave spectrum. Bull. Earthquake Res. Inst. Tokyo Univ. 44, 73–78. Areias, P., Reinoso, J., Camanho, P.P., Cesar de Sa, J., Rabczuk, T., 2018. Effective 2D and 3D crack propagation with local mesh refinement and the screened Poisson equation. Eng. Fract. Mech. 189, 339–360. Brace, W., Walsh, J., Frangos, W., 1968. Permeability of granite under high pressure. J. Geophys. Res. 73 (6), 2225–2236. Byerlee, J.D., 1978. Friction of rocks. Pure Appl. Geophys. 116 (4), 615–626. Calò, M., Dorbath, C., Frogneux, M., 2014. Injection tests at the EGS reservoir of Soultzsous-Forêts. Seismic response of the GPK4 stimulations. Geothermics 52, 50–58. Cornelio, C., Spagnuolo, E., Di Toro, G., Nielsen, S., Violay, M., 2019. Mechanical behaviour of fluid-lubracated faults. Nat. Commun. 10, 1274. Dieterich, J.H., 1978. Time-dependent friction and the mechanics of stick-slip. Pure Appl. Geophys. 116 (4–5), 790–806. Duan, K., Ji, Y., Wu, W., Kwok, C.Y., 2019a. Unloading-induced failure of brittle rock and implications for excavation-induced strain burst. Tunn. Undergr. Space Technol. 84, 495–506. Duan, K., Ji, Y., Xu, N., Wan, Z., Wu, W., 2019b. Excavation-induced fault instability: possible causes and implications for seismicity. Tunn. Undergr. Space Technol. 92, 103041. Fang, Y., Elsworth, D., Wang, C., Ishibashi, T., Fitts, J.P., 2017. Frictional stability-permeability relationships for fractures in shales. J. Geophys. Res. Solid Earth 122, 1760–1776. Foulger, G.R., Wilson, M.P., Gluyas, J.G., Julian, B.R., Davies, R.J., 2018. Global review of human-induced earthquakes. Earth-Sci. Rev. 178, 438–514. Gadelmawla, E.S., Koura, M.M., Maksoud, T.M.A., Elewa, I.M., Soliman, H.H., 2002. Roughness parameters. J. Mater. Process. Technol. 123 (1), 133–145. Grasso, J.R., 1992. Mechanics of seismic instabilities induced by the recovery of hydrocarbons. Pure Appl. Geophys. 139 (3-4), 507–534. Guglielmi, Y., Cappa, F., Avouac, J.P., Henry, P., Elsworth, D., 2015. Seismicity triggered by fluid injection-induced aseismic slip. Science 348 (6240), 1224–1226. Hamdia, K.M., Silani, M., Zhuang, X., He, P., Rabczuk, T., 2017. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. Int. J. Fract. Mech. 206 (2), 215–227. Husen, S., Kissling, E., von Deschwanden, A., 2013. Induced seismicity during the construction of the Gotthard Base Tunnel, Switzerland: hypocenter locations and source dimensions. J. Seismol. 17 (1), 63–81. ISRM, 2015. ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007–2014. Springer, Cham, Switzerland. Manouchehrian, A., Cai, M., 2018. Numerical modeling of rockburst near fault zones in deep tunnels. Tunnel. Undergr. Space Technol. 80, 164–180. Kanamori, H., Brodsky, E.E., 2001. The physics of earthquakes. Phys. Today 54 (6), 34–40. Karner, S.L., Marone, C., 2000. Effects of loading rate and normal stress on stress drop and stick-slip recurrence interval. GeoComplexity and the Physics of Earthquakes. pp. 187–198. Kohli, A.H., Zoback, M.D., 2013. Frictional properties of shale reservoir rocks. J. Geophys. Res. Solid Earth 118 (9), 5109–5125. Kuang, W.H., Zoback, M.D., Zhang, J., 2017. Estimating geomechanical parameters from microseismic plane focal mechanisms recorded during multistage hydraulic
13
Contents lists available at ScienceDirect
Engineering Geology journal homepage: www.elsevier.com/locate/enggeo
Unloading-induced rock fracture activation and maximum seismic moment prediction
T
Yinlin Jia, Wei Wua,*, Zhihong Zhaob a b
School of Civil and Environmental Engineering, Nanyang Technological University, 639798, Singapore Department of Civil Engineering, Tsinghua University, Beijing, 100084, China
ARTICLE INFO
ABSTRACT
Keywords: Rock fractures Frictional instability Induced seismicity Seismic moment
Hundreds of anthropogenic earthquakes have recently occurred worldwide due to underground space creation and energy extraction. The mechanism behind the human-induced geohazards is most likely associated with the reduction of normal stress on pre-existing fractures and faults. This study reports a series of laboratory experiments to investigate the mechanism of unloading-induced fracture activation, and proposes a simple approach to predict the maximum seismic moment for a critically stressed fracture. The unloading-driven shear test results exhibit that the unloading process induces the stress states of the sawcut and natural fractures to approach the Mohr-Coulomb failure envelope, and the normal stress unloading rate influences the peak slip rate. The fracture instability is dependent on the relationship between the stiffness of the system and the slip weakening rate of the fracture, and the shear dilation mainly occurs after the fracture activation. The test results also show that the critical shear stress of the sawcut fracture during the unloading-driven shear test is approximately equal to the residual shear strength after the displacement-driven fracture slip. This relationship inspires us to develop a new approach to estimate the maximum seismic moment. Our data demonstrate that the maximum seismic moments for both the fractures obtained from the unloading-driven shear tests are all below the upper limit lines, indicating that the proposed approach is reasonable. The uncertainty analysis shows that the accurate estimation of fault size can improve the maximum seismic moment prediction.
1. Introduction
principle, fracture activation occurs when the shear stress on a rock fracture exceeds the friction coefficient multiplied by the effective normal stress, which is defined as the difference between the normal stress and pore pressure. The fracture instability is accompanied by friction loss, which can be described by the slip weakening friction law based on a spring-slider system with one degree of freedom (Dieterich, 1978; Rice, 1983). The fracture stably slides when the stiffness of the spring-slider system is larger than the slip weakening rate of the fracture, otherwise the fracture exhibits an unstable motion. Besides the classic models, in-situ and ex-situ investigations have been carried out to explore the mysteries behind the induced earthquakes related to underground mass extraction. Microseismic monitoring is commonly used to observe the occurrence of fracture activation based on the spatial and temporal evolution of microseismic events (Rutlegde et al., 1998; Xu et al., 2016). A series of unloading-induced direct-shear tests is conducted to study the unloading-induced fracture instability (Wu et al., 2014), and demonstrates that the fracture instability is accompanied by the reduction of both normal and shear stresses on the fracture (Wu et al., 2017). The fracture instability is also associated
Underground resource development is expected to meet our everincreasing demand for economic development without disturbing the environment. However, several anthropogenic earthquakes have recently occurred owing to underground space creation and energy extraction. The extraction of mass (e.g., rock and fluid) from the subsurface perturbs the frictional equilibrium of pre-existing fractures and faults, resulting in the human-induced geohazards (Foulger et al., 2018). Excavation-induced seismicity has been reported during the constructions of the Gotthard Base Tunnel in Switzerland (Husen et al., 2013) and the Jinping II hydropower station in China (Zhang et al., 2013). Induced earthquakes associated with oil and gas extraction also occurred at the Wilmington oil field in the United States (Segall, 1989), the Gazli gas reservoir in Uzbekistan (Simpson and Leith, 1985), and the Groningen gas field in Netherlands (Vlek, 2018). These geohazards have caused incalculable economic losses and even led to project abandonment. According to the Mohr-Coulomb failure criterion and effective stress
⁎
Corresponding author. E-mail addresses: [email protected] (Y. Ji), [email protected] (W. Wu), [email protected] (Z. Zhao).
https://doi.org/10.1016/j.enggeo.2019.105352 Received 16 June 2019; Received in revised form 11 October 2019; Accepted 18 October 2019 Available online 21 October 2019 0013-7952/ © 2019 Elsevier B.V. All rights reserved.
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
with the change of principal stress orientation with respect to fracture orientation (Duan et al., 2019b). The volumetric contraction of reservoir rock due to gas extraction can change the local stress on the fracture (Segall et al., 1994) and induce the readjustment of crustal stress (Grasso, 1992). The occurrence of the seismic events due to gas extraction can also be correlated with the stiffness between seal- and reservoir rocks, pressure drop, and fracture density (van Eijs et al., 2006). Understanding the mechanism of unloading-induced fracture activation allows us to manage the risk of the anthropogenic geohazards. Numerical models combined with microseismic data have been used to evaluate the fracture instability during underground excavation (Xu et al., 2015). A seismic risk analysis using a risk matrix approach indicates that the Groningen gas field has a much larger seismic risk than the other onshore gas fields in the Netherlands (van Thienen-Visser et al., 2018). The total seismic moment for a fractured rock mass can be estimated by the modulus of rigidity and the volume of extracted rock in mining engineering (McGarr, 1976), and by the modulus of rigidity, the mass of removed fluid, the density of upper crust, and the fraction of seismogenic upper crust in petroleum engineering (McGarr, 1991). The seismic moment for a rectangular fracture can also be assessed based on the stress drop (Stein and Wysession, 2009). This study aims to investigate the mechanism of unloading-induced fracture activation and to forecast the maximum seismic moment for a critically stressed fracture. We perform the unloading-driven shear tests on sawcut and natural fractures to study the unloading-induced fracture instability, and discuss the mechanism of the fracture activation, in terms of slip weakening friction, fracture shear dilation, and peak slip rate. A simple approach is proposed to predict the maximum seismic moment and verified using the test results. We also conduct the uncertainty analysis to quantify the variability of the maximum seismic moment, and discuss the possible application of this approach in engineering projects.
porosity, and water content are 2660 kg/m3, 0.26%, and 0.07%, respectively. The Young’s modulus, Poisson’s ratio, cohesion, and friction angle are 74 GPa, 0.15, 38 MPa, and 54°, respectively. The physical and mechanical properties of the granite were measured based on the test standards suggested by International Society for Rock Mechanics (ISRM, 2015). The intrinsic permeability of the granite matrix is 1.3 μD measured based on the pressure pulse decay method (Brace et al., 1968; Zoback and Byerlee, 1975). Sawcut and natural fractures were prepared from 50-mm-diameter granite cores (Fig. 1a). The sawcut fracture was obtained by cutting a 100-mm-long core sample at 30° to the core axis using a diamond saw. The fracture surface roughness was controlled by fine sandpaper with 25.6 μm particle size. The natural fracture was produced by compressing a 124-mm-long core sample containing a naturally formed but unopened fracture at 27° to the core axis until the sample purely failed along the fracture. The fracture surfaces were imaged before and after the shear tests using a structured-light 3D scanner. Two Teflon layers were used to secure the sample on the coreholders and to isolate the sample from the confining oil. Two 2-mm-diameter boreholes were drilled at the sample ends to facilitate fluid flow from the coreholders to the fracture. Two 130-μm-thickness and 0.45-μm-pore size filter papers were inserted between the coreholders and sample ends to prevent sheared-off particles from entering the pore system. The sample assembly was installed in the MTS 815 rock mechanics test system (Fig. 1b). Silicon oil and distilled water were used as the confining and pore fluids, respectively. The normal and shear stresses on the inclined fracture were servo-controlled by the axial stress (σ1) and confining pressure (σ3), and corrected by considering the reduction of fracture area and the deformation of Teflon layers (Kohli and Zoback, 2013). The shear stress (τ) and normal stress (σn) can be calculated as:
2. Experimental method
where Ψ is the fracture inclination angle with respect to the core axis, σ1-σ3 is the deviator stress directly measured by an in-vessel load cell.
=(
1
=
3
n
2.1. Experimental setup
3)sin
+(
1
cos
(1)
2 3)cos
(2)
2.2. Experimental procedure
Bukit Timah granite, underlying the central region of Singapore Island, is of extreme high-quality and suitable for cavern development (Zhao et al., 1999). The granite is composed of 62% feldspar, 32% quartz, 5% biotite, and 1% hornblende. The average bulk density,
We carried out a series of displacement-driven and unloadingdriven shear tests. The displacement-driven shear tests were used to obtain the shear strength of the fracture, and the unloading-driven Fig. 1. (a) Sawcut (left) and natural (right) fractures, and (b) schematic drawing of sample assembly. Normal (σn) and shear (τ) stresses on the fracture with inclination angle (Ψ) with respect to the core axis are servo-controlled by axial stress (σ1) and confining pressure (σ3). Pore pressure is applied on the fracture through two boreholes.
2
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
shear test SDD1 and the succeeding unloading-driven shear test SDU1 as examples. First, under 5 MPa initial normal stress, the displacementdriven shear test was carried out at an axial displacement rate of 1 μm/s until the fracture reached its shear strength. Second, we fixed the shear stress as 85% of the shear strength, and assumed the fracture close to a critical stress state. Finally, the normal stress was reduced at an unloading rate of 0.01 MPa/s until the fracture was activated. For the dry sawcut fracture, seven unloading-driven shear tests were performed under different combinations of initial normal stress (e.g., 5, 7.5, and 10 MPa), shear stress level (e.g., 85%, 90%, and 95% of the shear strength), and normal stress unloading rate (e.g., 0.005, 0.01, and 0.05 MPa/s). Under 11 MPa initial normal stress, the drained and undrained conditions were considered during the unloading-driven shear tests on the saturated sawcut fracture at 0.01 MPa/s unloading rate. The pore pressure was maintained at 1 MPa under the drained condition, while the pore volume was constant with 1 MPa initial pore pressure under the undrained condition. For the natural fracture, the displacement-driven shear tests were conducted under the low normal stresses of 1, 2, and 3 MPa to minimize asperity damage and to construct the failure envelope, which was subsequently used to calculate the shear strength under high normal stresses. A critical stress state of the fracture was then simulated with the shear stress equal to 85% of the calculated shear strength. Under 10 MPa initial normal stress, three unloading-driven shear tests on the dry natural fracture were carried out at different unloading rates (0.005, 0.01, and 0.05 MPa/s). Under 11 MPa initial normal stress and 1 MPa pore pressure, the drained and undrained conditions were also studied during the unloading-driven shear tests on the saturated natural fracture at 0.01 MPa/s unloading rate. 3. Results and discussion 3.1. Experimental results The Mohr-Coulomb failure envelopes for the sawcut and natural fractures under the dry and saturated conditions are constructed based on the displacement-driven shear strengths (Fig. 3). The friction coefficients of the dry and saturated sawcut fractures are 0.83 and 0.68, respectively. The cohesion of the sawcut fractures is negligible. The friction coefficient of the natural fracture also reduces from 0.86 under the dry condition to 0.77 under the saturated condition. The cohesions of the dry and saturated natural fractures are 0.07 and 0.12 MPa, respectively. The reduction of friction coefficient is due to the lubrication effect of pore water (Cornelio et al., 2019). Fig. 3 also shows that the stress states of both the fractures are below the failure envelopes at critically stressed states and close to the failure envelopes at the onset of fracture activation. The initial normal and shear stresses significantly influence the change of stress state, whereas the effect of normal stress unloading rate is relatively minor. The drained and undrained conditions also slightly affect the unloading-driven instability of both the fractures. The sawcut and natural fractures exhibit dynamic slip and quasidynamic sliding, respectively. The dynamic slip of the sawcut fracture induced by normal stress reduction is characterized by a sudden increase in shear displacement and an abrupt drop of shear stress (Fig. 4a). After that, the fracture arrests, and the shear stress reverts to its initial value. The slip rate and acceleration are the first- and secondtime derivatives of shear displacement, respectively, which instantaneously jump during the dynamic slip. The unloading-induced dynamic slip of the sawcut fracture resembles the coseismic slip during the stick-slip cycles of natural faults, in terms of stress drop and
Fig. 2. (a) Experimental principle, and (b) typical stress paths of displacementdriven and unloading-driven shear tests. In the unloading-driven shear test, the fracture is initially subjected to a critical shear stress (τc) determined by the shear strength (τs) obtained from the displacement-driven shear test. Δτ is equal to τs minus τc. c and μ are the cohesion and friction coefficient of the fracture, respectively.
shear tests were conducted to simulate the unloading-induced activation of the fracture initially approaching a critical stress state. During the unloading-driven shear tests, the Mohr circle moves towards the Mohr-Coulomb failure envelope due to the reduction of the normal stress on the fracture, which finally results in the fracture activation (Fig. 2a). During the unloading-driven shear tests on the sawcut fracture (Table 1), the effects of initial normal stress, critical shear stress, normal stress unloading rate, and fluid drainage were considered to investigate the mechanism of unloading-induced fracture activation. Fig. 2b demonstrates the experimental procedure, taking the displacement-driven
3
4
σn = 11, P = 1
SSD3
Stable sliding Stable sliding Stable sliding
Stable sliding Stable sliding Stable sliding
Saturated natural fracture NSD1 σn = 1, P = 0 NSD2 σn = 2, P = 0 NSD3 σn = 3, P = 0
5.62
–
–
–
9.91
5.73
3.60
– – –
– – –
Residual shear strength (MPa)
Dry natural fracture NDD1 σn = 1 NDD2 σn = 2 NDD3 σn = 3
Stick-slip
Stick-slip
Stick-slip
σn = 12.5
SDD4
σn = 8.5, P = 1
Stick-slip
σn = 10
SDD3
SSD2
Stick-slip
σn = 7.5
SDD2
Stick-slip
Stick-slip
Dry sawcut fracture SDD1 σn = 5
Saturated sawcut fracture SSD1 σn = 6, P = 1
Failure mode
Displacement-driven shear test
NPD NPU
5.87
5.87
11.19
10.60
10.01
10.01
10.01
5.69
3.67
Critical shear stress (MPa)
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Dynamic slip
Failure mode
σni = 11, P = 1, τc = 85%τs, Ur = 0.01, drain σni = 11, P = 1, τc = 85%τs, Ur = 0.01, undrain
σni = 10, τc = 85%τs, Ur = 0.005 σni = 10, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.05
σni = 11, P = 1, τc = 85%τs, Ur = 0.01, drain σni = 11, P = 1, τc = 85%τs, Ur = 0.01, undrain
σni = 5, τc = 85%τs, Ur = 0.01 σni = 7.5, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.005 σni = 10, τc = 85%τs, Ur = 0.01 σni = 10, τc = 85%τs, Ur = 0.05 σni = 10, τc = 90%τs, Ur = 0.01 σni = 10, τc = 95%τs, Ur = 0.01
NDU1 NDU2 NDU3
SPU
SPD
SDU7
SDU6
SDU5
SDU4
SDU3
SDU2
SDU1
Unloading-driven shear test
7.00 7.00
6.70 6.70 6.70
7.98
8.38
8.89
8.13
7.40
7.59
7.54
5.99
4.06
Normal stress at failure (MPa)
0.15
0.13
0.16
0.15
0.45
0.16
0.16
0.06
0.07
Quasi-dynamic sliding Quasi-dynamic sliding
8.63 8.65
8.17 7.96 7.66
Shear displacement (mm)
Quasi-dynamic sliding Quasi-dynamic sliding Quasi-dynamic sliding
1.61
1.47
1.65
1.86
2.48
2.14
2.01
0.83
0.79
Shear stress drop (MPa)
– –
– – –
2.25
2.11
2.62
2.47
3.68
2.56
2.32
1.14
1.05
– –
– – –
Peak slip rate (mm/s)
0.22 0.23
0.24 0.25 0.29
110.51
156.72
146.97
196.97
178.25
143.84
167.49
58.73
82.14
Peak acceleration (mm/s2)
39.70 34.53
40.92 57.16 35.50
Table 1 Summary of experimental conditions and results. A rock fracture is subjected to normal stress (σn in MPa) and pore pressure (P in MPa) in the displacement-driven shear test, and initial normal stress (σni in MPa), critical shear stress (τc in MPa) equivalent to a certain percentage of displacement-driven shear strength (τs in MPa), pore pressure (P in MPa) and normal stress unloading rate (Ur in MPa/s) in the unloading-driven shear test.
Y. Ji, et al.
Engineering Geology 262 (2019) 105352
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 3. Mohr-Coulomb failure envelopes and stress states before and at the onset of unloading-induced fracture activation for (a) dry sawcut fracture, (b) saturated sawcut fracture, (c) dry natural fracture, and (d) saturated natural fracture.
displacement jump (Karner and Marone, 2000). For the natural fracture, the shear stress is constant during the quasi-dynamic sliding (Fig. 4b). The shear displacement exponentially increases, accompanied by the slight changes of slip rate and acceleration. The peak slip rate of the natural fracture is similar to that of the slow slip events of natural faults (Tinti et al., 2016). Fig. 5 presents the topographic contours of the sawcut and natural fractures before and after the displacement-driven and unloadingdriven shear tests. The lowest point in the borehole is unaffected during the shear process and considered as the zero datum, which makes the contours comparable. We calculate the root mean square (RMS) asperity height to quantify the evolution of surface roughness (Gadelmawla et al., 2002). The RMS asperity height of the sawcut fracture reduces from 1.31 mm before the shear tests to 0.65 mm after
the tests. The RMS asperity height of the natural fracture also decreases from 2.32 mm to 2.06 mm. The reductions of RMS asperity height influence the stiffness of the system, which controls the fracture behaviors. 3.2. Mechanism of unloading-induced fracture activation Our results show that the unloading-induced fracture activation is dependent on fracture surface roughness, and can be explained by a spring-slider model with one degree of freedom (Dieterich, 1978). The slip weakening friction law describes that the friction coefficient decreases with larger shear displacement and becomes nearly constant after the fracture experiences a slip weakening distance (Rice, 1983). The fracture activation is associated with the relationship between the
5
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 4. (a) Unloading-induced dynamic slip of sawcut fracture (test SDU4), and (b) unloading-induced quasi-dynamic sliding of natural fracture (test NDU3).
stiffness of the spring-slider system and the slip weakening rate of the fracture. In this study, the stiffness of the system can be considered as the combined rigidity of the test system and the fracture, which is estimated as the slope of the linear elastic portion of the shear stress
versus axial displacement curve obtained from the displacement-driven shear test. The slip weakening rate of the fracture is estimated as the slope of the linear portion of the shear stress versus shear displacement curve obtained from the unloading-driven shear test. Additionally, we
6
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 5. Topographic contours of (a) sawcut fracture and (b) natural fracture before and after displacement-driven and unloading-driven shear tests.
use the fracture strength to describe the change in shear resistance of the fracture, which is determined as the product of the friction coefficient and the effective normal stress. When the stiffness of the system is lower than the slip weakening rate of the sawcut fracture (Fig. 6a), the
fracture experiences the dynamic slip. After that, the fracture strength recovers due to the reformation and strengthening of asperity contacts (Marone, 1998). The natural fracture slides in a quasi-static manner due to the stiffness of the system greater than the slip weakening rate of the
7
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
fluid volume indicates the normal dilation of the fracture. The hydraulic apertures of both the sawcut and natural fractures slightly increase before the fracture activation, but dramatically rise with larger shear displacement (Fig. 7). The results indicate that both the fractures experience the remarkable dilation during the shear process, which is due to one asperity overriding another and pore space expansion among sheared-off particles (Morgan and Boettcher, 1999; Li et al., 2019). The slight decline in shear displacement before the onset of fracture activation is presumably due to the elastic unloading of rock blocks (Duan et al., 2019a). Fig. 8 exhibits that the peak slip rates of both the sawcut and natural fractures as a function of normal stress unloading rate. The peak slip rate of the sawcut fracture is the maximum rate value of the unloadinginduced dynamic slip, and that of the natural fracture is the rate value at 0.25 mm shear displacement, at which the stable sliding is achieved. The peak slip rate under the same stress state increases with higher unloading rate owing to a larger amount of strain energy released at the onset of fracture activation (Wu et al., 2017). The peak slip rate under the undrained condition exceeds that under the drained condition, because the pore pressure at fracture slip is reduced under the undrained condition and the amount of released energy increases with larger effective normal stress (Lei et al., 2016). 4. Prediction of maximum seismic moment The similarity between the unloading-induced dynamic slip of the sawcut fracture and the coseismic slip of natural faults allows us to predict the maximum seismic moment due to normal stress reduction. The seismic moment describes earthquake size and is determined as (Aki, 1966): (4)
M0 = GAd
where A is the fracture area, d is the shear displacement, and G is the combined rigidity of the test system and the fracture. To estimate the combined rigidity of the test system and the sawcut fracture with 3.9 × 10−9 km2 surface area, we derive the combined rigidity in each dynamic slip case based on the empirical relationships among the moment magnitude (Mw), fracture area, and seismic moment (Wells and Coppersmith, 1994; Kanamori and Brodsky, 2001):
Mw = 4.07+0.98log(A)
Mw =
(6)
We then constrain the combined rigidity of the seven dynamic slip cases using the least squares method and obtain a best-fit value of 1.2 GPa (Fig. 9a). The derived rigidity is validated based on the relationships among the earthquake scaling parameters. For a circular fracture with a radius (R), the shear displacement (d) and stress drop (Δτ) can be expressed as (Stein and Wysession, 2009):
Fig. 6. Analysis of unloading-induced instability of (a) sawcut fracture (test SDU4) and (b) natural fracture (test NDU3).
natural fracture (Leeman et al., 2016). The stress imbalance between the maintained shear stress and decreasing fracture strength leads to the increase in slip rate (Fig. 6b). The shear dilation of rock fractures also influences the unloadinginduced fracture activation. During the unloading-driven shear tests on the sawcut and natural fractures, the undrained condition with constant pore volume provides an avenue to estimate the hydraulic aperture of the fracture, which can represent the mechanical aperture (Fang et al., 2017). The change in pore fluid volume (ΔV) can be related to the change in hydraulic aperture (Δe):
V=A e
log M0 -6.06 1.5
(5)
d
cM0/(G (2R)2) =
= cM0/(2R)3 =
7M0 8GR2
7 M0 16 R3
(7) (8)
where c is the shape factor of a circular fracture and equal to 3.5. The seven dynamic slip cases are plotted in Fig. 9b. The source dimension is considered as a fracture radius of 0.075 m, which is the average length of the major and minor axes of the oval-shaped fracture. The seven cases fall within a stress drop range of 0.5–3 MPa and a shear displacement range of 50–500 μm, which are consistent with the results directly measured from the unloading-driven shear tests (Table 1). Hence, the combined rigidity is reasonable to determine the seismic
(3)
where A is the fracture area. The water compressibility is neglected at constant temperature. For a given fracture area, the increase in pore
8
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 7. Shear stress, shear displacement, pore pressure, and hydraulic aperture change during the unloading-driven shear tests on (a) sawcut fracture (test SPD) and (b) natural fracture (test NPD).
moments for the sawcut and natural fractures. Table 1 also shows that the critical shear stress of the sawcut fracture during the unloading-driven shear test (τc) is approximately equal to the residual shear strength after the displacement-driven fracture slip, which is equal to the shear strength (τs) minus stress drop (Δτ) (see
Fig. 2b). Thus, the normal stress reduction required to activate the fracture can be calculated as: n
9
=
µ
=
s
c
µ
(9)
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
1978). The standard deviations of both the parameters are considered as one-tenth of the average of the maximum and minimum values. The linear correlation between the two parameters, known as Pearson correlation coefficient, is 0.0002. The first-order Sobol indices, which quantify the variance of the increment rate from one of the parameters, are 0.92 for the fracture radius and 0.07 for the friction coefficient, respectively. The second-order Sobol index, which measures the joint effect of the two parameters, is 0.01. The results show that the firstorder index of fracture radius is much larger than the other indices, and suggest that the successful prediction of the maximum seismic moment is highly dependent on the accurate evaluation of fault size. In engineering applications, the success of predicting the maximum seismic moment for unloading-induced seismicity using Eq. (10) relies heavily on the determination of input parameters, i.e., the friction coefficient and size of pre-existing faults, and the normal stress reduction due to underground activities. The friction coefficient can be assumed as 0.6, which is a typical value obtained from laboratory experiments and has been used to forecast the maximum seismic moment for injection-induced seismicity (Zoback, 2007; McGarr, 2014). The fault size can be estimated from the moment magnitude of microseismic events by means of the seismological scaling factors (Kuang et al., 2017; Stein and Wysession, 2009). The change in normal stress is related to the change in volume (e.g., removed rock and extracted fluid) (McGarr, 1976). Therefore, this study provides a possible solution to predict the maximum seismic moment for field-scale seismic events associated with underground space creation and energy extraction. The proposed approach can be further verified by field investigations with reasonable input parameters (e.g., Calò et al., 2014; Guglielmi et al., 2015) and numerical simulations with efficient computational methods (e.g., Areias et al., 2018; Rabczuk and Belytschko, 2004). Additionally, the proposed approach is established based on the uniform reduction of normal stress on pre-existing faults. Underground activities may trigger the instability of large-scale faults through local disturbances, e.g., rock burst (Manouchehrian and Cai, 2018; Duan et al., 2019a) and stress rotation (Duan et al., 2019b). Our approach may not be suitable for this case, in which the maximum seismic moment is determined not only by the coefficient of friction, fault size, and normal stress reduction, but also by fault structure and rupture propagation.
Fig. 8. Peak slip rates of sawcut and natural fractures as a function of normal stress unloading rate.
where μ is the friction coefficient of the fracture and obtained from Fig. 3. According to Eqs. (8) and (9), the maximum seismic moment can be estimated based on the normal stress reduction at fracture failure:
M0 =
16µR3 7
n
(10)
The unloading-induced moment release of the sawcut fracture is mostly seismic (Fig. 4a). The combined rigidity of the test system and the fracture, the fracture area, and the shear displacement at fracture failure are thus used to calculate the maximum seismic moment. The friction coefficients of the dry and saturated sawcut fractures are 0.83 and 0.68, respectively, and used to plot the two upper limit lines (Eq. (10)). Based on the unloading-induced shear test results, the maximum seismic moment generally increases with larger normal stress reduction (Fig. 10a). The seismic moments for the dry and saturated fracture cases are all below the corresponding upper limit lines, which verifies the prediction of the maximum seismic moment. For the natural fracture, the unloading-induced moment release is progressively aseismic (Fig. 4b). The seismic moment is thus dependent on shear displacement and calculated as the product of the combined rigidity, the fracture area, and the cumulative shear displacement. The friction coefficient used for the upper limit line is 0.82, which is the average values of the dry and saturated natural fractures. Fig. 10b demonstrates that the unloading-induced cumulative seismic moment continuously increases with larger normal stress reduction until it reaches the upper limit line. This figure also indicates that Eq. (10) is valid to forecast the seismic moment. The variability of the maximum seismic moment due to normal stress reduction is determined by the fracture radius and friction coefficient (Eq. (10)). We carry out the uncertainty analysis to quantify the variability based on the Sobol’s method (Sobol, 1993; Vu-Bac et al., 2016; Hamdia et al., 2017). The increment rate of the maximum seismic moment as a function of normal stress reduction (i.e., 16μR3/7) is considered in this analysis. The statistical properties of the fracture radius and friction coefficient are listed in Table 2. The fault radius is assumed to follow a normal distribution with the minimum and maximum values of 10−6 m and 106 m, respectively, which represent natural faults in a size range from micrometers to thousands of kilometers. The friction coefficient also follows a normal distribution with the minimum and maximum values of 0.6 and 0.85, respectively (Byerlee,
5. Conclusions We investigate the unloading-induced instability of sawcut and natural fractures based on the unloading-driven shear tests. The test data show that the unloading-induced fracture behaviors obey the Mohr-Coulomb failure criterion, and are strongly affected by fracture surface roughness. The fracture instability is controlled by the relationship between the stiffness of the system and the slip weakening rate of the fracture, and the shear dilation mainly occurs after the fracture activation. The peak slip rates of both the fractures increase with higher normal stress unloading rate. The drained and undrained conditions play a minor role in the fracture instability. The maximum seismic moment is estimated based on the fact that the residual shear strength after the displacement-driven fracture slip is comparable to the critical shear stress during the unloading-driven fracture slip. The maximum seismic moment for the sawcut fracture obtained from the unloading-driven shear test can be bounded by the upper limit line. The time-dependent seismic moment for the natural fracture increases with larger normal stress reduction, and finally approaches the upper limit line. The proposed approach for the maximum seismic moment prediction can be used in engineering projects by reasonably determining the input parameters, and the prediction accuracy can be improved by accurately estimating fault size.
10
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 9. (a) Recovery of the combined rigidity of the test system and the sawcut fracture based on the empirical relationship between the moment magnitude and source area of natural earthquakes. (b) Comparison between the unloading-induced dynamic slips of the sawcut fracture and natural earthquakes in terms of stress drop and slip displacement. s denotes standard deviation.
11
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Fig. 10. Seismic moment as a function of normal stress reduction for (a) sawcut fracture and (b) natural fracture.
12
Engineering Geology 262 (2019) 105352
Y. Ji, et al.
Table 2 Statistical properties of fracture radius and friction coefficient. Parameter
Distribution type
Maximum value
Minimum value
Standard deviation
Friction coefficient Fault size (m)
Normal Normal
0.85 106
0.60 10−6
0.07 5 × 104
fracturing. Geophysics 82 (1), KS1–KS11. Scuderi, M.M., Leeman, J.R., Saffer, D.M., Marone, C., 2016. Laboratory observations of slow earthquakes and the spectrum of tectonic fault slip modes. Nature Comm 7, 11104. Lei, X., Funatsu, T., Ma, S., Liu, L., 2016. A laboratory coustic emission experiment and numerical simulation of rock fracture driven by a high-pressure fluid source. J. Rock Mech. Geotech. Eng. 8, 27–34. Li, Y.C., Wu, W., Tang, C.A., Liu, B., 2019. Predicting the shear characteristics of rock joints with asperity degradation and debris backfilling under cyclic loading conditions. Int. J. Rock Mech. Min. Sci. 120, 108–118. Marone, C., 1998. Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Planet. Sci. 26 (1), 643–696. McGarr, A., 1976. Seismic moment and volume changes. J. Geophys. Res. 81 (8), 1487–1494. McGarr, A., 1991. On a possible connection between three major earthquakes in California and oil production. Bull. Seismol. Soc. Am. 81 (3), 948–970. McGarr, A., 2014. Maximum magnitude earthquakes induced by fluid injection. J. Geophys. Res. Solid Earth 119 (2), 1008–1019. Morgan, J.K., Boettcher, M.S., 1999. Numerical simulations of granular shear zones using the distinct element method, 1. Shear zone kinematics and the micromechanics of localization. J. Geophys. Res. 104 (B2), 2703–2719. Rabczuk, T., Belytschko, T., 2004. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int. J. Numer. Methods Eng. 61 (13), 2316–2343. Rice, J.R., 1983. Constitutive relations for fault slip and earthquake instabilities. Instabilities in Continuous Media. Springer. Rutlegde, J.T., Phillips, W.S., Schuessler, B., 1998. Reservoir characterisation using oilproduction-induced microseismicity, Clinton County, Kentucky. Tectnophysics 289 (1–3), 129–152. Segall, P., 1989. Earthquakes triggered by fluid extraction. Geology 17, 942–946. Segall, P., Grasso, J.R., Mossop, A., 1994. Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France. J. Geophys. Res. 99 (B8), 15423–15438. Simpson, D.W., Leith, W., 1985. The 1976 and 1984 Gazli, USSR, earthquakes – were they induced? Bull. Seismol. Soc. Am. 75, 1465–1468. Sobol, I., 1993. Sensitivity analysis for non-linear mathematical models. Reliab. Eng. Syst. Safe. 52, 1–17. Stein, S., Wysession, M., 2009. An Introduction to Seismology, Earthquakes, and Earth Structure. John Wiley & Sons. Tinti, E., Scuderi, M.M., Scognamiglio, L., Di Stefano, G., Marone, C., Collettini, C., 2016. On the evolution of elastic properties during laboratory stick-slip experiments spanning the transition from slow slip to dynamic rupture. J. Geophys. Res. Solid Earth 121 (8569), 8594. van Eijs, R.M.H.E., Mulders, F.M.M., Nepveu, M., Kenter, C.J., Scheffers, B.C., 2006. Correlation between hydrocarbon reservoir properties and induced seismicity in the Netherlands. Eng. Geol. 84, 99–111. van Thienen-Visser, K., Roholl, J.A., van Kempen, B.M.M., Muntendan-Bos, A.G., 2018. Categorizing seismic risk for the onshore gas fields in the Netherlands. Eng. Geol. 237, 198–207. Vu-Bac, N., Lahmer, T., Zhuang, X., Nguyen-Thoi, T., Rabczuk, T., 2016. A software framework for probabilistic sensitivity analysis for computationally expensive models. Adv. Eng. Softw. 100, 19–31. Vlek, C., 2018. Induced earthquakes from long-term gas extraction in Groningen, the Netherland: statistical analysis and prognosis for acceptable-risk regulation. Risk Anal. 38 (7), 1455–1473. Wells, D.L., Coppersmith, K.J., 1994. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bull. Seismol. Soc. Am. 84 (4), 974–1002. Wu, W., Zhao, Z., Duan, K., 2017. Unloading-induced instability of a simulated granular fault and implications for excavation-induced seismicity. Tunn. Undergr. Space Technol. 63, 154–161. Wu, W., Zou, Y., Li, X., Zhao, J., 2014. An unload-induced direct-shear model for granular gouge friction in rock discontinuities. Rev. Sci. Instrum. 85 (9), 093902. Xu, N.W., Li, T.B., Dai, F., Li, B., Zhu, Y.G., Yang, D.S., 2015. Microseismic Monitoring and stability evaluation for the large scale underground caverns at the Houziyan hydropower station in Southwest China. Eng. Geol. 188, 48–67. Xu, N.W., Li, T.B., Dai, F., Zhang, R., Tang, C.A., Tang, L.X., 2016. Microseismic monitoring of strainburst activities in deep tunnels at the Jinping II hydropower station, China. Rock Mech. Rock Eng. 49 (3), 981–1000. Zhang, C., Feng, X.T., Zhou, H., Qiu, S., Wu, W., 2013. Rockmass damage development following two extremely intense rockbursts in deep tunnels at Jinping II hydropower station, southwestern China. Bull. Eng. Geol. Environ. 72 (2), 237–247. Zhao, J., Zhou, Y.X., Hefny, A.M., Cai, J.G., Chen, S.G., Li, H.B., Liu, J.F., Jain, M., Foo, S.T., Seah, C.C., 1999. Rock dynamics research related to cavern development for ammunition storage. Tunn. Undergr. Space Technol. 14 (4), 513–526. Zoback, M.D., 2007. Reservoir Geomechanics. University Press, Cambridge. Zoback, M.D., Byerlee, J.D., 1975. The effect of microcrack dilatancy on the permeability of Westerly granite. J. Geophys. Res. 80 (5), 752–755.
Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgement Wei Wu gratefully acknowledges the support of Start-Up Grant from Nanyang Technological University, Singapore. References Aki, K., 1966. Generation and propagation of G waves from the Niigata earthquake of June 16, 1964, 2, estimation of earthquake moment, release energy, and stress-strain drop from the G wave spectrum. Bull. Earthquake Res. Inst. Tokyo Univ. 44, 73–78. Areias, P., Reinoso, J., Camanho, P.P., Cesar de Sa, J., Rabczuk, T., 2018. Effective 2D and 3D crack propagation with local mesh refinement and the screened Poisson equation. Eng. Fract. Mech. 189, 339–360. Brace, W., Walsh, J., Frangos, W., 1968. Permeability of granite under high pressure. J. Geophys. Res. 73 (6), 2225–2236. Byerlee, J.D., 1978. Friction of rocks. Pure Appl. Geophys. 116 (4), 615–626. Calò, M., Dorbath, C., Frogneux, M., 2014. Injection tests at the EGS reservoir of Soultzsous-Forêts. Seismic response of the GPK4 stimulations. Geothermics 52, 50–58. Cornelio, C., Spagnuolo, E., Di Toro, G., Nielsen, S., Violay, M., 2019. Mechanical behaviour of fluid-lubracated faults. Nat. Commun. 10, 1274. Dieterich, J.H., 1978. Time-dependent friction and the mechanics of stick-slip. Pure Appl. Geophys. 116 (4–5), 790–806. Duan, K., Ji, Y., Wu, W., Kwok, C.Y., 2019a. Unloading-induced failure of brittle rock and implications for excavation-induced strain burst. Tunn. Undergr. Space Technol. 84, 495–506. Duan, K., Ji, Y., Xu, N., Wan, Z., Wu, W., 2019b. Excavation-induced fault instability: possible causes and implications for seismicity. Tunn. Undergr. Space Technol. 92, 103041. Fang, Y., Elsworth, D., Wang, C., Ishibashi, T., Fitts, J.P., 2017. Frictional stability-permeability relationships for fractures in shales. J. Geophys. Res. Solid Earth 122, 1760–1776. Foulger, G.R., Wilson, M.P., Gluyas, J.G., Julian, B.R., Davies, R.J., 2018. Global review of human-induced earthquakes. Earth-Sci. Rev. 178, 438–514. Gadelmawla, E.S., Koura, M.M., Maksoud, T.M.A., Elewa, I.M., Soliman, H.H., 2002. Roughness parameters. J. Mater. Process. Technol. 123 (1), 133–145. Grasso, J.R., 1992. Mechanics of seismic instabilities induced by the recovery of hydrocarbons. Pure Appl. Geophys. 139 (3-4), 507–534. Guglielmi, Y., Cappa, F., Avouac, J.P., Henry, P., Elsworth, D., 2015. Seismicity triggered by fluid injection-induced aseismic slip. Science 348 (6240), 1224–1226. Hamdia, K.M., Silani, M., Zhuang, X., He, P., Rabczuk, T., 2017. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. Int. J. Fract. Mech. 206 (2), 215–227. Husen, S., Kissling, E., von Deschwanden, A., 2013. Induced seismicity during the construction of the Gotthard Base Tunnel, Switzerland: hypocenter locations and source dimensions. J. Seismol. 17 (1), 63–81. ISRM, 2015. ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007–2014. Springer, Cham, Switzerland. Manouchehrian, A., Cai, M., 2018. Numerical modeling of rockburst near fault zones in deep tunnels. Tunnel. Undergr. Space Technol. 80, 164–180. Kanamori, H., Brodsky, E.E., 2001. The physics of earthquakes. Phys. Today 54 (6), 34–40. Karner, S.L., Marone, C., 2000. Effects of loading rate and normal stress on stress drop and stick-slip recurrence interval. GeoComplexity and the Physics of Earthquakes. pp. 187–198. Kohli, A.H., Zoback, M.D., 2013. Frictional properties of shale reservoir rocks. J. Geophys. Res. Solid Earth 118 (9), 5109–5125. Kuang, W.H., Zoback, M.D., Zhang, J., 2017. Estimating geomechanical parameters from microseismic plane focal mechanisms recorded during multistage hydraulic
13