Zirui Yin, Hongwei Huang, Fengshou Zhang,*, Lianyang Zhang, Shawn Maxwell

a Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai, 200092, China

b Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai, 200092, China

c Department of Civil and Architectural Engineering and Mechanics, University of Arizona, Tucson, AZ, 85721, USA

d MaxSeis LLC., Calgary, Alberta, T2P 3N9, Canada

ABSTRACT This paper presents a three-dimensional fully hydro-mechanical coupled distinct element study on fault reactivation and induced seismicity due to hydraulic fracturing injection and subsequent backflow process, based on the geological data in Horn River Basin, Northeast British Columbia, Canada. The modeling results indicate that the maximum magnitude of seismic events appears at the fracturing stage.The increment of fluid volume in the fault determines the cumulative moment and maximum fault slippage, both of which are essentially proportional to the fluid volume. After backflow starts, the fluid near the joint intersection keeps flowing into the critically stressed fault, rather than backflows to the wellbore. Although fault slippage is affected by the changes of both pore pressure and ambient rock stress, their contributions are different at fracturing and backflow stages. At fracturing stage, pore pressure change shows a dominant effect on induced fault slippage.While at backflow stage,because the fault plane is under a critical stress state, any minor disturbance would trigger a fault slippage. The energy analysis indicates that aseismic deformation takes up a majority of the total deformation energy during hydraulic fracturing. A common regularity is found in both fracturing- and backflow-induced seismicity that the cumulative moment and maximum fault slippage are nearly proportional to the injected fluid volume. This study shows some novel insights into interpreting fracturing- and backflowinduced seismicity,and provides useful information for controlling and mitigating seismic hazards due to hydraulic fracturing.

Keywords:Induced seismicity Fault reactivation Hydraulic fracturing Backflow Geomechanical modeling Distinct element method

1. Introduction

As an emerging issue worldwide, fluid injection-induced seismicity has greatly raised public’s concern over the past few years(e.g. Davies et al., 2013; Ellsworth, 2013; Bao and Eaton, 2016;Elsworth et al., 2016; Lei et al., 2017; Foulger et al., 2018; Grigoli et al., 2018; Kim et al., 2018; Rubinstein et al., 2018). Typical industrial injection activities include hydraulic fracturing in lowpermeability hydrocarbon reservoir, hydraulic stimulation in enhanced geothermal system (EGS), fluid injection for carbon dioxide sequestration,and disposal of wastewater(Evans et al.,2012;Le, 2018). Although the number of induced seismicity events has significantly increased in recent years,the mechanisms of induced seismicity have not been well understood. What is the role of injected fluid volume or injection rate(Davies et al.,2013;McGarr,2014;Farmahini-Farahani and Ghassemi,2016)and how to predict the seismic energy (van der Elst et al., 2016; Galis et al., 2017;Stanˇek and Eisner, 2017; McGarr and Barbour, 2018; De Barros et al., 2019) are the two major concerns. The injected fluid is not the only factor that affects fault instability. For example, Scuderi and Collettini (2018) pointed out that the stability of a fault may be controlled by the coupling effect of fluid pressure and the rateand-state friction parameters. It is also proved that fluid injection can induce seismicity and aseismic creep (Guglielmi et al., 2015;Wei et al., 2015; Cappa et al., 2018, 2019; McGarr and Barbour,2018; Bhattacharya and Viesca, 2019). Laboratory experiments showed that fluid injection may induce aseismic creep before seismic slippage (Cappa et al., 2019). Goodfellow et al. (2015)concluded that the majority of deformation energy is aseismic during hydraulic fracturing,and De Barros et al.(2019)pointed out that the seismic energy takes up less than 1%of the total hydraulic energy during an injection activity, suggesting that most of the energy released in shear deformation is aseismic.

Due to the complexity of hydro-mechanical interactions in fractured rock masses (Rutqvist and Stephansson, 2003;Detournay, 2016), numerical simulation has become a major tool in investigating the mechanism of induced seismicity.Cappa et al.(2006) and Guglielmi et al. (2008) studied the hydro-mechanical behavior of fractured rock masses and the corresponding seismic response during hydraulic fracturing by distinct element and finite element modeling. Distinct element modeling results showed an explicit correlation between microseismicity and hydraulic fracturing treatment (Maxwell, 2013). Maxwell et al.(2015) also indicated that slowing down injection process is likely to control the total cumulative seismic energy and mitigate seismic hazards. Moreover,the distinct element modeling results by Grob et al. (2016) pointed out that a larger area undergoing high pore pressure on the fault plane indicates a larger average magnitude of seismic activities, which is consistent with the results in Shapiro and Dinske (2010). Some other studies using coupled TOUGH-FLAC model also proved that hydraulic fracturing only gives rise to microseismicity instead of larger earthquake(e.g. Rutqvist et al., 2013, 2015; Rinaldi and Rutqvist, 2019).However, the reason that some cases show a direct correlation with large seismic activities (ML≥4.0) induced by hydraulic fracturing has not been convincingly explained yet (Bao and Eaton, 2016; Lei et al., 2017). McClure and Horne (2010, 2011)developed a coupled fluid flow and rate/state model solved by finite volume method and two-dimensional (2D) displacement discontinuity method, and proposed a feasible method for reduction of the maximum event magnitude by decreasing injection pressure in an EGS project. Similar conclusion was obtained by Farmahini-Farahani and Ghassemi (2016) who used a fast multipole displacement discontinuity method (FMDDM) for simulation of geothermal and petroleum reservoirs. Baisch et al.(2010) simulated the seismic events at the geothermal site of Soultz-sous-Forêts using finite element method, and compared the magnitude-frequency characteristics of induced seismicity with field observations. Although numerical analysis has some limitations such as presence of geological uncertainties(Fairhurst,2013), it is still a vital tool for better understanding of induced seismicity.

Following the boom development of unconventional resources(Zhao et al., 2018), increasing number of seismic cases was reported in association with hydraulic fracturing. Some of these cases have brought serious damage to the public properties(Holland,2013;Clarke et al.,2014;Bao and Eaton,2016;Lei et al.,2017). The Sichuan Basin (China) and the Horn River Basin (Canada) are the two typical project sites that have experienced growing induced seismicity during hydraulic fracturing.Statistics of recorded earthquakes in the Changning-Weiyuan shale gas production region in the Sichuan Basin shows the number of earthquakes with magnitude (ML) larger than zero has increased dramatically from an average number of 1525 times per year(2009—2013) to 9539 times in 2018 (Fig. 1). Similarly, the Horn River Basin,the largest shale gas reservoir site in Canada,has also experienced a soaring number of earthquakes. The recorded earthquakes with magnitudeML≥1.0 within 100 km of the Fort Nelson station ascended from an average value of 24 times in 2002 and 2003(before the shale gas exploitation)to 131 times in 2011.In addition, the maximum magnitude rises fromML= 2.9 toML= 3.6, showing the correlation between increasing seismic hazards and fracturing operation(Farahbod et al.,2015).However,given the significant variations of field parameters such as fluid injection rate and volume, the mechanism of large induced seismicity during hydraulic fracturing needs further investigation(Clarke et al., 2014; Bao and Eaton, 2016; Lei et al., 2017).

Fig.1. The number of earthquakes annually with a magnitude ML ≥0.0 from 2009 to 2018 in the Changning-Weiyuan shale gas production region,Sichuan Basin,Southwest China (102.25°E to 105.43°E, and 27.83°N to 30.03°N). Data source:China Earthquake Networks Center (CENC).

Nevertheless, it has been reported that induced seismicity may occur during injection period, and some large events may take place in post-injection period (H?ring et al., 2008; Baisch et al.,2010; Ellsworth, 2013; McClure, 2015), or even months after shut-in (Grigoli et al., 2018; Keranen and Weingarten, 2018; Kim et al., 2018). The triggering mechanism of induced seismicity during post-injection is still not mature (De Simone et al., 2017). For hydraulic fracturing treatment,fluid would backflow for the recycle and protection of environment, with proppant remained in the stimulated fractures. In this situation, fluid pressure drops rapidly near the wellbore,and stress state in the fault may be redistributed,which is possible to trigger more seismic events. McClure (2015)pointed out that the backflow fluid can enter other larger natural fractures and trigger additional seismic events. However, most of the previous studies focused on fracturing period, and few have been reported on the induced seismicity during backflow. Therefore, it is of significance to investigate different triggering mechanisms during fracturing and backflow periods so that a comprehensive understanding of induced seismicity during hydraulic fracturing can be achieved.

This paper focuses on fault reactivation and induced seismicity due to hydraulic fracturing and backflow. For this, a threedimensional (3D) fully hydro-mechanical coupled distinct element model is employed to quantitatively evaluate the synthetic microseismicity for the two different periods (Damjanac and Cundall, 2016). As one of the two typical fracturing-induced seismicity sites, the Horn River Basin is selected with associated geological parameters (Farahbod et al., 2015; Woo et al., 2017;Zhang et al., 2019). Several indicators in the modeling are analyzed, including fault slippage and seismic energy, to show the spatio-temporal evolution of fault reactivation. Then the energy budget of synthetic seismic events is estimated. In addition, the mechanism of induced seismicity during fluid backflows and different seismicity characteristics between injection and backflow periods are discussed.

2. Numerical scheme

2.1. Numerical method

The 3D distinct element code (3DEC) is employed to construct the model of fault reactivation and induced seismicity due to hydraulic fracturing (Cundall, 1988; Hart et al., 1988). 3DEC is suitable for analysis of discontinuous medium, which is represented as an assemblage of discrete blocks.Each block is subdivided into a mesh of finite difference elements consisting of tetrahedral zones and nodes.Newton’s law of motion is applied at all nodes to update the velocity,displacement and nodal force at each time step.In addition, discontinuities are treated as boundary conditions between blocks. A contact detection algorithm is performed, and the integration of the law of motion is used to obtain new positions of blocks and contacts at the next time step.3DEC has the capability of analyzing fluid flow in fractures in hydro-mechanically fully coupled manner (Cappa et al., 2006; Guglielmi et al., 2008; Zhang et al., 2013; Maxwell et al., 2015; Konietzky et al., 2017; Zhang and Mack, 2017).

2.2. Hydro-mechanical coupling model

Transient flow of a compressible fluid is used as a basic algorithm in 3DEC (Itasca, 2015). Within a domain network structure,each domain is a triangular element with a certain fluid pressure.Fluid flow is controlled by the pressure difference between adjacent domains ΔP, which is given by

where ρwis the fluid density;P1andP2are the domain pressures at two adjacent domains, respectively;gis the gravitational acceleration; andz1andz2are the depths of the two adjacent domains,which are the centroid depths of the triangular elements.The cubic law for fluid flow in a planar fracture is applied(Witherspoon et al.,1980), where the flow rateqis expressed by

whereais the domain aperture;μ is the fluid viscosity;andlis the domain length, which is defined as the distance between the centroids of two adjacent domains.The domain apertureais updated in each mechanical step using

wherea0is the initial domain aperture;Δunand Δusare the normal and shear displacement increments, respectively; and ψ is the dilation angle. In each mechanical step, the domain pressurePis updated according to the domain fluid volume change:

whereP0is the initial domain pressure;Kwis the fluid bulk modulus;Qis the sum of flow rate flowing into the domain from its adjacent domains; Δtis the current time step;VandV0are the domain volumes at the new step and the old step,respectively;and ΔV=V-V0is the domain volume change.For each time step,the newly added fluid volume is first stored in the injected domains,and then the domain pressure and the domain aperture are updated, forcing the fluid to flow into adjacent domains and changing the domain volumes, which subsequently influences the domain pressure.

2.3. Method validation

The method of modeling hydraulic fracturing with 3DEC is first validated by comparing the modeling results with theoretical solution of a radial hydraulic fracture(HF)(Detournay,2016;Dontsov,2016; Dontsov and Zhang, 2018; Zhang and Dontsov, 2018). Radial HF is selected herein because there is a rigid theoretical solution for such a problem.In distinct element modeling,tensile strength and model I fracture toughness are related by mesh size(Potyondy and Cundall, 2004):

where σtis the tensile strength,KIcis the model I fracture toughness,dis the average edge size of tetrahedral zones, and σt0is a dimensionless parameter(σt0= 1 by default).For the convenience of mathematical expression, four parameter groups are defined as

where μ is the fluid viscosity;Eand ν are the Young’s modulus and Poisson’s ratio of the rock mass,respectively;andCLis the Carter’s leak-off coefficient (Economides and Nolte,2000).

For a penny-shaped HF,fracture propagation can be divided into four regimes: toughness-dominated (K), viscosity-dominated (M),leak-off viscosityand leak-off toughnessThese regimes represent different dominant mechanisms of energy dissipation.Viscosity-dominated regime means fluid viscosity is the main factor controlling the solution of fracture propagation,and toughnessdominated regime means fracture toughness plays a dominant role.Two variables, dimensionless variable τ and leak-off parameter φ(Detournay, 2016; Dontsov, 2016), are written as

wheretmkis the dimensionless time, by which the regime of fracture propagation is determined in the absence of leak-off, andAs noted by Dontsov and Zhang(2018),a uniform tensile strength used in a 3DEC model may lead to nonuniform apparent toughness at the fracture tip,but the effect is trivial in the viscosity-dominated regime. Hence, the tensile strength of the HFs is set to be zero,i.e.the fracture toughnessKIc=0 according to Eq. (5). Since fluid leak-off is ignored in the validation model, the HF propagates in the viscosity-dominated regime which can be described by the M-vertex solution (viscosity-dominated regime):

wherewmis the fracture aperture as a function of scaled spatial coordinate from the injection point ρ (ρ =r/R, whereRdenotes the location of the fracture tip andrrefers to the coordinate of a point inside the fracture) and injection timet;Q0is the injection rate;pmis the fluid pressure at the pointr;Rmis the fracture radius;and the function F is defined as

wheresis also a scaled spatial coordinate (s= (R-r)/R), which denotes the relative distance from a pointrto the fracture tip;λ is a variable for the convergence of iteration procedure, which is initially set to be 1/2 to capture the elliptic solution for a uniformly pressurized crack; andM(ρ，s) refers to a piecewise function as

Note that the functionsK(·) andE(·) are the complete elliptic integrals of the first and the second types,respectively.is a slowly varying function from the tip asymptotic solution which is expressed as

A radial HF model is constructed using mechanical parameters in the original model(Fig.2).This model is a cuboid with the size of 120 m × 120 m × 160 m, and only contains one horizontal plane crossing the central point, which represents the HF. A uniform compressive stress acts perpendicularly to the plane and is equal to the initial minimum stress (Shmin= 53.75 MPa) at the depth of 2500 m. The injection point is at the centroid of this model. Note that the initial minimum stress is consistent with the one in the subsequent fault reactivation model.

Besides the four parameter in Eq.(6),the M-vertex solution also depends on injection timetand injection rateQ0. Therefore, we consider injection duration of 200 s with a rate ofQ0= 0.027m3/s for the validation model. Young’s modulus (E), fluid viscosity (μ)and Poisson’s ratio (ν) of the rock mass are derived from the subsequent fault reactivation model. These parameters are:E=39.6 GPa，μ = 0.02 Pa s， ν = 0.221，CL= 0，σt= 0，Shmin=53.75 MPa，Q0= 0.027 m3/s， andt= 200 s.

Note that there is a singularity in φ andtmkifK′= 0 (Eq. (7)).Thus,when applying numerical solution of Eq.(8),the scaled model I fracture toughness is set to be a small value ofK′= 1 × 105MN/m3/2for approximation. In the absence of fluid leak-off, the validation model is in viscosity-dominated regime where φ = 0 and τ = 1.43×10-18(Fig. 3).

Modeling results are shown in Fig. 4. A circular fracture propagation centered in the injection point can be identified. Given the discretization nature of the tetrahedral mesh grid, there are some minor notches at the fracture tip.However,the fracture is overall in a radial shape.

Comparison of 3DEC solution and M-vertex solution(viscositydominated regime) is shown in Fig. 5. The red dot lines are footprints of theoretical solutions while the black solid lines are 3DEC solutions.Results in Fig.5 demonstrate that 3DEC can well predict the radial fracture size, fracture aperture and fluid pore pressure.

2.4. Model set-up

This model is based on a published case which provides both geological parameters and recorded seismicity information in the Horn River Basin, Canada (Sneiling et al., 2013; Sneiling and de Groot, 2014; Zhang et al., 2019). Fig. 6 shows a schematic view of the numerical model made of a cubic block with the centroid coordinate of (0, 0, -2500) m. The length, width and depth of the cuboid are 1920 m,1200 m and 512 m,respectively.In this model,the directions of all HF planes are predefined.The height of the HFs is set to be 160 m. A fault plane is constructed with the dip of 90°and the dip direction of 45°. The size of the fault plane is 700 m × 320 m (length × width), intersecting withy-axis at the point of(0,150)m.The wellbore alongx-axis contains five clusters with spacing of 20 m. A strike-slip stress regime is applied to the model. The maximum and minimum horizontal stresses, with stress gradients of 35 kPa/m and 21.5 kPa/m,are used alongx-andy-axis, respectively. The vertical principal stress and the pore pressure gradients are set to be 25 kPa/m and 12 kPa/m, respectively. This model focuses on the interaction between the HFs and the fault.Thus,the discrete fracture network(DFN)is not included.Note that the HFs only account for these parts with tensile failure due to fluid injection.All planes are initially glued together with a given tensile strength before fluid injection. To sustain uniform toughness at the fracture tip,the given tensile strength is set to be zero in this model,which is valid for the propagation of HFs in the viscosity-dominated regime. Parameters of model geometry and stress are listed in Table 1.

Fig. 3. Variation diagram of limiting fracture propagation regimes versus dimensionless time τ and leak-off parameter φ (Detournay, 2016; Dontsov, 2016). The black dot represents the viscosity-dominated propagation regime based on the selected injection parameters for the hydraulic fracturing validation model.

Fig. 4. Modeling results of the hydraulic fracturing validation model. Upper row: fracture aperture after 100 s injection (left) and 200 s injection (right); and Bottom row: pore pressure after 100 s injection (left) and 200 s injection (right).

Fig.5. Simulation results of a radial hydraulic fracture in 3DEC and the comparison with M-vertex theoretical solutions.Upper row:fracture aperture(left)and fluid pore pressure(right) after 100 s injection; and Bottom row: fracture aperture (left) and fluid pore pressure (right) after 200 s injection.

Fig. 6. The 3D model geometry (left) and the joints configuration (right). The heights of fault plane and hydraulic fracture planes are 320 m and 160 m, respectively.

Table 1Parameters of model geometry and stress.

Rock matrix is modeled as an elastic body(Zhao et al.,2016).The density of the block material is 2550 kg/m3, the bulk modulus is 23.66 GPa, and the Poisson’s ratio is 0.221. The fluid viscosity is 0.02 Pa s and the fluid density is 1200 kg/m3(Table 2). A Carter’s leak-off model is applied with leak-off coefficient of 1×10-6m/s1/2(Economides and Nolte, 2000). The simulation scheme is 1-h injection with an injection rate of 0.17 m3/s, followed by a 3-h backflow with a fixed borehole pressure of 36 MPa. The backflow borehole pressure is slightly higher than the initial pore pressure(30 MPa)at the depth of injection point (-2520 m).

The initial joint aperture of the fault plane and the HFs are 0.2 mm and 0.05 mm, respectively. To avoid shear failure on the HFs, the friction angle of the HFs is set to be 40°, larger than the value on the fault plane(30°).Besides,a cohesion of 50 MPa is also assumed on the HFs to ensure the occurrence of only tensile failure.

2.5. Fault constitutive model

The Coulomb slip model has been the most widely applied joint constitutive model for hydro-mechanical coupled problems.However, the standard Coulomb slip model is not suitable for modeling seismic slip behavior. Therefore, this study uses an improved nonlinear softening-healing Coulomb slip model(Itasca,2015;Zhang et al.,2020).This model is simplified from the‘slip-weakening model’proposed by Rice(1983)for shear faulting.It overcomes two drawbacks of the standard Coulomb slip model.First, when the joint starts to slip, the shear strength will gradually descend from the peak value to the residual one, following a specified exponential curve. This improvement avoids instantaneous stress drop and raises the threshold of the seismic slip.Second, when the slipping stops, the residual shear strength will recover back to the peak value,equivalent to the case that a static friction coefficient is re-assigned. With the assumption that the joint stiffness is constant,the increments of the normal and shear stresses are proportional to the corresponding displacements.Thus,the changes of the normal(Δσn)and shear(Δτ)stresses are given by

whereknandksare the normal and shear stiffnesses of the joint,respectively; and Δunand Δusare the normal and shear displacements,respectively.Note that the tensile stress change is limited by the specified tensile strength.If the tensile strength is exceeded,the tensile stress change (Δτ) will drop to zero. The shear strength is given by the Coulomb model function (Raleigh et al.,1976; Jaeger et al., 2007):

where τsis the shear strength,cis the cohesion,and φ is the friction angle. Note that the peak cohesion and friction angle values are used when the shear stress is smaller than the shear strength;once the slipping distance exceeds a critical value,the residual cohesion and friction angle values are used instead. At the beginning of slipping, the shear strength is proportional to the shear slip at a slopeksbefore shear failure occurs. In the transition section after the shear resistance reaches the peak value given by Eq. (14), the shear strength evolution follows a nonlinear curve:

where the superscripts ‘peak’ and ‘resid’ represent the peak and residual values, respectively;denotes the plastic slip in the transition section, beginning when the shear strength reaches the peak value;Dcis a critical slipping distance from the peak to theresidual shear strength; and α is a coefficient controlling the amount of the strength drop, and α> 1. Note that the curve is not smooth at the pointbut the difference can be ignored if α is large enough (>2). Thus, the coefficient α is set to be 3 in this model. The schematic of the nonlinear softening-healing Coulomb slip model is shown in Fig.7.Wong(1982)inferred that the abscissa value of centroid of the triangle-like area enclosed by the power function curve, the slope line and the horizontal residual shear strength line is about the order of 0.5 mm. Therefore, the critical slipping distanceDcis set to be 1 ×10-3m in this model.

Table 2Parameters of rock and fluid materials in the numerical model.

In this study,the standard Coulomb slip model is implemented for the HFs, while the nonlinear softening-healing Coulomb slip model is applied on the fault plane. Parameters of the two joint materials are summarized in Table 3.

2.6. Model of synthetic microseismicity

Numerical method has the advantage of quantitative evaluation for the simulation results.In this study,microseismic events can be calculated according to the shear displacement of slipping nodes and the attributes of the joint material(Zhang and Mack,2017).For dissociation of kinematic equation,each contact face is subdivided into multiple triangle zones, where the nodes denote the triangle vertices.Each node has an area in charge,which is defined as onethird the sum of the areas of all the surrounding triangles. Therefore,the seismic moment can be calculated by

Fig. 7. Nonlinear softening-healing Coulomb slip model for fault plane. The shear strength increases linearly before reaching the peak value, and then drops within a specified distance Dc. When slipping stops, the shear strength resumes to the peak value instantaneously.

whereGis the shear modulus of the rock mass,Adenotes the total area of slipping nodes in one seismic event, andis the average plastic slip of all slipping nodes in this event and is weighted by node areas. Note that Eq. (16) only considers shear failure,i.e.seismic events triggered by tensile failure are neglected.This way will possibly overestimate the proportion of large seismic events in the model.Some slipping behavior does not cause seismic activities;therefore,only seismic slipping nodes can be counted as seismic events. A specified criterion is applied for combining seismic slipping nodes to larger seismic events.In the spatial scale,each slipping node as a single event has a radiusRfor recognition.The subsequent slipping nodes in this range can be counted in the same event; otherwise they belong to different events (Fig. 8). In the time scale, the slipping time of a large event is determined by the beginning of the first slipping node,and the end of the last one,as long as these nodes have overlaps in slipping duration.After the combination, the microseismic moment is estimated by Eq. (16).The slipping area of one event is represented by the total area of slipping nodes in it.The corresponding moment magnitude can be calculated by Hanks and Kanamori (1979):

In this study, a slipping node is regarded as a seismic slip only when two criteria are met.One condition is that the normal stress is larger than 5 MPa, and the other is that the slipping velocity is above 0.5 mm/s. A positive normal stress guarantees a slipping event under compression, and a critical velocity wipes off the low speed creep. Any event which does not meet one of these two requirements is treated as aseismic slip (Zhang et al., 2020).

3. Simulation results of fault reactivation

The modeling is carried out at two stages: a 1-h hydraulic fracturing stage and a 3-h backflow stage. Fig. 9 illustrates the borehole pressure history. The borehole pressure at the fracturing stage reaches 55 MPa,about 1 MPa higher than the initial minimum principal horizontal stress at depth of -2520 m(Shmin=54.18 MPa),which is consistent with the field observations during hydraulic fracturing(Haimson and Fairhurst,1967;Hickman and Zoback, 1983; Economides and Nolte, 2000). At the entire backflow stage, the borehole pressure value is fixed at 36 MPa.

Fig.10 plots the joint aperture and the pore pressure at the end of each stage.After 60-min injection,the maximum aperture on the HFs reaches about 3.5 mm, and the maximum pore pressure is about 56 MPa.The elevated pore pressure does not cross the fault,because shear failure along the fault occurs before the propagation of the HFs on the back of it. At the end of backflow, the HFs are nearly closed near the injection point.However,at places far away from the wellbore, pore pressure inside the HFs is much higher than the initial value, indicating that the HFs still propagate away from the wellbore at the backflow stage. Moreover, the pore pressure difference near the intersection is relatively higher at the beginning of backflow, forcing the fluid nearby to flow into the fault.

The pore pressure and the shear displacement on the fault plane at different stages are shown in Fig.11.The pink spheres represent the location of the fracturing clusters. Note that at 180 min of injection, only the slippage appearing at the backflow stage is counted. The result shows that the fault slippage at the fracturing stage is about 2.5 times that at the backflow stage (20 mm vs.8 mm). Besides, at the backflow stage, slippage only occurs on thetwo ends of the fault plane, in accordance with the pore pressure diffusion. This finding indicates that fluid diffusion plays an important role in fault slippage during backflow.

Table 3Parameters of joint material in the numerical model.

Fig. 8. Spatio-temporal combination of seismic slipping nodes for a large event: (a) Spatial scale combination. Note that the radius for recognition (R) is set to be two times the average length of the modeling zone edges(Edge_Size);(b)Schematic of the node area(An).The area of each node is one-third of the total area of the surrounding triangular zones(At); and (c) Time scale combination. Nodes in the range of recognition radius and with overlapped slipping time will be regarded as one large event.

Fig. 9. Borehole pressure history in the whole simulation process.

The evolution of fluid volume in the HFs and the fault is shown in Fig.12. After resetting the borehole pressure to a relatively low value at the beginning of backflow, the fluid volume drops rapidly in the HFs, showing that the fluid is extracted through wellbore.However, the fluid volume in the fault keeps increasing in the backflow stage.

Synthetic microseismic events in the fracturing and backflow stages are shown in Fig. 13. Note that the semitransparent background is the shear slip result as shown in Fig. 11. The size of spheres represents the moment magnitude of events,and the color represents the time of occurrence.Note that the time is reset to be zero when plotting microseismicity at the backflow stage. It is obvious that the microseismicity spreads along the trace of the pore pressure diffusion, and some large events appear at the rim of the slipping contour,showing a sudden release of the seismic energy at nearly 40 min of injection. At the backflow stage, most of the microseismic events occur at the first 1 h. It can be estimated that the moment magnitude of events at the fracturing stage is larger than that at the backflow stage.

Fig.14 plots the moment magnitude and the cumulative seismic moment at the fracturing and the backflow stages. The maximum magnitude of microseismic events occurs at the fracturing stage and most of the events fall into the range of magnitude-0.5 to 1.A dramatic rise of the cumulative moment occurs at about 40 min of injection, showing a sudden release of the strain energy, which is also illustrated by the large spheres at the rim of slipping contour(Fig. 13). While at the backflow stage, the cumulative seismic moment ascends slowly,and is consistent with the tendency of the fluid volume change in the fault(Fig.12).The low number of events in the range of-1.5 to 0 is probably due to the small total number of events(348)at a relatively coarser discretization of the fault plane.

To explore the effect of mesh size on the simulated seismicity,two additional simulations with a mesh size of 7 m and 11 m,respectively, are carried out. For clarification, ‘zone size’ means‘mesh size’in this context.Fig.15 shows the magnitude distribution of seismic events at different mesh sizes. It can be seen that the magnitude distribution is sensitive to the mesh size.The number of medium-sized events(-1 to 0)increases significantly as the mesh size decreases. An explanation is that a finer mesh brings more events and leads to more normally distributed magnitudes. However, many other factors such as velocity threshold and critical normal stress also affect the magnitude distribution. Since this study mainly compares the total moment in the fracturing and backflow processes, further investigations on the magnitude distribution are needed.

The details of shear displacement,cumulative moment and fluid volume from the numerical model are presented in Table 4.Regardless of the number of events, the maximum shear displacement and the cumulative seismic moment are almost proportional to the fluid volume change in the fault,indicating that the fluid volume increment is a crucial factor in inducing fault slippage and the release of seismic energy.

Fig.10. Simulated aperture and pore pressure at the end of fracturing and backflow stages.

4. Energy budget analysis

Seismic energy is analyzed in this section to evaluate the seismicity at the fracturing stage. Theoretical predictions include McGarr and Barbour (2018) model for cumulative moment, Galis et al. (2017) model for maximum moment, and the hydraulic energy mentioned by De Barros et al.(2019).To convert energy scale into moment scale,the following relation is used(Aki and Richards,2002):

whereEsis the seismic energy, and Δσ is the background stress drop. Considering the slipping area and the maximum magnitude,Δσ is estimated to be 100 kPa(Zoback and Gorelick,2012).In Galis et al.(2017)model,the dynamic friction coefficient is set to be 0.4 and the characteristic length is defined as the fault height(320 m).The hydraulic energy is expressed by the integral of fluid pressure and flow rate:

whereEhis the hydraulic energy (J), which can be converted to hydraulic moment through Eq. (18); andPis the fluid pressure which is estimated to be 55 MPa from Fig. 9.

The theoretical predictions and the modeling results are shown in Fig.16. Only 0.0114% of the hydraulic energy is converted to the seismic energy.Both the Galis et al.(2017)and McGarr and Barbour(2018) modeling results are at least one order of magnitude larger than the simulated seismic moment,which means that the seismic deformation only takes up a small fraction of the total deformation energy (De Barros et al., 2019). However, selection of the two thresholds in Section 2.6 is a critical factor affecting the calculation of the seismic energy. The threshold velocity of 0.5 mm/s is a conservative estimation to ensure that enough seismic events are detected. The normal stress threshold of 5 MPa may also be underestimated for the pressurized slipping area, leading to overestimation of the seismic events.

5. Mechanisms of fault reactivation during injection and backflow

To investigate the mechanism of fault reactivation during injection and backflow, they-direction block stress and the displacement contour are plotted onx—yplane at depth of -2500 m (Fig.17). Five parallel lines represent the traces of the HFs,and the oblique line represents the fault plane(Fig.17a).Note that the compressive stress is positive. The result shows that the upper side of the fault slides upward while the lower side goes downward during fracturing. At the backflow stage, fault slippage occurs at the two ends of the fault plane(Fig.11),possibly impacted by the closure of the HFs.Nevertheless,at the end of the backflow stage (180 min), the block stress iny-direction at the upper left of the fault plane reduces by about 2 MPa compared with the one at the end of fracturing (60 min), indicating the decrease of normal stress on the fault plane. This stress reduction relaxes the upper section of the fault plane, and promotes the fault slippage at the two ends at the backflow stage.

Fig.18 shows the measurement points on the fault plane and the HF planes. The red dots (F1—F5) represent the points outside the fracturing region and are used to monitor fault slippage and ambient block stress change. The other ten points at the intersection of the fault and the HFs are set to study the evolution of the aperture and the normal stress on the fault plane (FP1—FP5) and the HF planes (HF1—HF5), respectively. The depth of the crosssection is -2500 m. To distinguish the intersection points on the HF planes and those on the fault plane,the measurement points on the HF planes are 1.5 times the average zone size away from the HF/fault intersection points.

Fig.11. Pore pressure and fault slippage at the fracturing and backflow stages. Note that at 180 min of injection, only the fault slippage at the backflow stage is counted.

Fig.12. Evolution of fluid volume in the hydraulic fractures(left)and the fault(right).At the beginning of backflow stage,the fluid volume drops rapidly in the hydraulic fractures,while the fluid volume inside the fault keeps increasing.

Fig.13. Synthetic microseismic events at the fracturing and backflow stages. Note that the time is reset to be zero when plotting microseismicity at the backflow stage.

Fig.19 shows the effective normal stress and the fracture/fault aperture versus time at the measurement points as shown in Fig.18. At the fracturing stage, all measurement points on the HF planes experience a reduction of normal stress and an increase of aperture except HF1, which is the farthest HF measurement point from the wellbore. HF5 experiences the largest aperture increase during injection due to the shortest fluid path from the wellbore to the fault plane.Thus,the largest amount of injection fluid flows into the fault plane through HF5.Meanwhile, a rapid aperture increase at HF5 and HF3 brings severe stress shadow effect in HF4, which possibly restricts the propagation of the fourth HF(Konietzky et al.,2017).At the backflow stage,HF5 experiences an increase of normal stress and a corresponding dramatic reduction of aperture, illustrating that a fast closure of fluid path occurs once the backflow starts (Fig. 19a and c). The fault aperture change is synchronous with the normal stress change,and the fault fluid pressure remains at a high level during backflow(Fig.19b and d).In addition,the pore pressure difference near the HF/fault intersection is relatively high at the beginning of the backflow (Fig.10). This indicates that once the backflow starts, the fluid near the HF/fault intersection keeps flowing into the fault from the HFs.Pore pressure dissipation near the HF/fault intersection could increase the effective normal stress and narrow the fluid path between them, resulting in fluid locked in the clamped fault.

Fig.14. Moment magnitude and cumulative seismic moment at the fracturing and backflow stages.

Fig.15. Magnitude distribution of seismic events with different mesh sizes.

Table 4Detailed simulation results.

To intuitively observe the effect of the fluid volume on the fault slippage,a real-time slippage monitoring method is applied on the measurement points (F1—F5). Fig. 20 plots the slippage history of five measurement points in Fig.18 and cumulative fluid volume in the fault(blue line).Once a point starts slipping,the corresponding segment of the dash line is marked as bold.The results show that all slippage occurs within the first 90 min of injection and most of the durations are within 2 min. The start time of slippage at different points is in accordance with the arrival time of the fluid pressure.At the backflow stage,measurement points F1 and F5 only slip at the beginning when the fluid volume increases rapidly,indicating that fault slippage depends on the pore pressure change. This result is consistent with the observation in Fig.11.At about 40 min,the fluid reaches measurement point F2 and induces a continuous slippage,which can explain why there is a sudden seismic energy release at this time (Fig. 14) and why the maximum event occurs near the limits of the slipping area (Fig.13).

During fluid injection, the pore pressure increment will decrease the effective normal stress on the fault plane and push the stress state closer to the critical condition.In addition,the change of block stress during fluid injection will increase the differential stress in an elastic rock formation(Keranen and Weingarten,2018).Therefore,combined pore pressure perturbations and the ambient stress changes determine the stress state on the fault(Segall et al.,1994; Segall and Lu, 2015; Keranen and Weingarten, 2018; Zhang et al., 2019). For further investigation of the relevance between fault slippage and the stress change as well as the effect of fluid at both fracturing and backflow stages, real-time Mohr circle diagrams for some of measurement points are monitored.

Note that measurement points F1—F5 are regarded as microunits which contain a 45°joint (Fig. 21a). According to the Mohr—Coulomb criterion,the joint shear strength can be expressed as Eq.(14)and the relationship between the compressive strength and the joint oblique angle can be derived as (Brady and Brown,2004; Zoback, 2007):

where β is the oblique angle between the joint plane and the direction of the maximum principal stress σ1.When β equals φ or 90°,the value of σ1is infinite,indicating that the shear failure along the joint plane only occurs where φ ＜ β ＜ 90°. The minimum of the compressive strength occurs at β = 45°+φ/2 and can be expressed as

Fig.16. Seismic moment and hydraulic energy versus injected fluid volume.

Fig.17. Block displacement and block stress in y-direction at the depth of -2500 m. Note that compressive stress is positive.

Fig.18. Schematic of the measurement points on the hydraulic fractures and the fault plane(the depth of the cross-section is-2500 m).Red dots represent the fault points which are outside the fracturing region,orange dots are the fault points at the intersection of these planes,and gray dots are the hydraulic fracturing points near the intersection of these planes. ‘Edge_Size’ means the average zone size.

Fig.19. (a) Effective normal stress obtained by the measurement points on the hydraulic fractures; (b) Effective normal stress obtained by the measurement points on the fault plane; (c) Fracture aperture obtained by the measurement points on the hydraulic fractures; and (d) Fracture aperture obtained by the measurement points on the fault plane.

In this case,the stress Mohr circle at failure is tangent with the strength line of the joint.At a certain angle β other than 45°+φ/2, the maximum stresses at failure are greater than the strength obtained from Eq. (21) and the stress Mohr circle crosses the strength line of the joint (Brady and Brown,2004; Zoback, 2007).

In this study, the friction angle of the fault φ = 30°and the oblique angle β = 45°. The red dot in Fig. 21b shows the corresponding maximum principal stress at the moment of shear failure along the joint. It is noted that the model in this study assumes elastic rock. Thus, the shear failure only occurs along the joint(fault) plane.

Fig.20. Slippage history of the five measurement points in Fig.18 and the cumulative fluid volume in the fault (blue line). The bold segments on the dash lines indicate the slipping periods.

The measurement points F2 and F3 are at different sides of the HF planes (Fig. 18), and their stress states under slipping at the fracturing stage are shown in Fig. 22. Note that the fault plane is parallel to thez-axis, which is the direction ofSv, while theShmaxandShminare alongx-andy-axis,respectively.Therefore,in the 3D Mohr circle diagram under the strike-slip stress regime, the fault stress state can be expressed by the Mohr circle whose two endpoints areShmaxandShminon the σ-axis. Besides, the dip direction of the fault plane is 45°, indicating that the corresponding stress state lies on the top of the Mohr circle.Therefore,the two cyan dots at the top of the Mohr circles (Fig. 22) represent the fault stress states before injection and under slipping, respectively.

If only the effect of the pore pressure change is considered, the Mohr circle is supposed to move towards the fault strength line,as shown in the dashed magenta semi-circle in Fig. 22. The result demonstrates that at the fracturing stage, measurement points F2 and F3 have experienced a combined effect of pore pressure elevation and ambient block stress change, suggesting that the diameter of the Mohr circle also changes besides the lateral translation.The translation amount(Δpp)from the solid blue semicircle to the dashed magenta semi-circle is larger at measurement point F3 than that at measurement point F2, which means that measurement F3 experiences a larger pore pressure elevation due to the shorter distance from the wellbore.In addition,the distance between the two left endpoints of the magenta dashed semi-circle and the red solid semi-circle (ΔShmin) at F3 is shorter compared with the one at F2.This result indicates that reduction of minimum stress at the upper side of the HF planes(F3)is restrained due to the dilation of HFs, which is in accordance with the observation in Fig. 17. Therefore, it can be concluded that pore pressure change plays an important role in inducing fault slippage at both the upper and lower sides of the HF planes. However, the upper side experiences larger pore pressure change and smaller block stress change than the lower side does. The contributions of the pore pressure elevation and the ambient block stress change to the fault slippage vary at different fault locations.

Fig. 21. (a) Schematic of the measurement point regarded as a micro-unit containing a 45° joint; and (b) Relationship between the maximum principal stress at the moment of shear failure and the joint oblique angle.

Fig.22. Fault stress state of measurement points F2(a)and F3(b)under slipping at the fracturing stage. The black line is the Mohr—Coulomb strength line of the fault. The solid blue, dotted red and dashed magenta semi-circles represent the Mohr circle before injection, the Mohr circle at slipping, and the Mohr circle when only the effect of pore pressure change is considered, respectively. The two cyan dots at the top of Mohr circles represent the stress state before injection and under slipping,respectively.

Fig. 23 plots the stress state of measurement points F1 and F5 under slipping at backflow stage.Similar to measurement points F2 and F3,the two points have experienced a combined effect of pore pressure elevation and ambient block stress change, while the difference lies on that after shut-in,the solid blue circle is touching the fault strength line, indicating that the fault plane is under a critical state at the two points.The small translation(Δpp)from the solid blue semi-circle to the dashed magenta semi-circle shows that the impact of the fluid is not as obvious as that at the fracturing stage,but high enough for fault slippage as long as the fault is under critical stress state after fracturing stage. However, release of injection-induced compressive stress iny-direction also attributes to fault slippage as represented by ΔShminin Fig.17,making the fault stress point (cyan dot) closer to the fault strength line and more prone to slipping.

Fig. 23. Fault stress state of measurement points F1 (a) and F5 (b) under slipping at backflow stage. The black line is the Mohr—Coulomb strength line of the fault. The solid blue,dotted red and dashed magenta circles represent the Mohr circle after shutin, Mohr circle at slipping, and Mohr circle when only considering the effect of pore pressure change, respectively. The two cyan dots at the top of Mohr circles represent the stress states before backflow and under slipping, respectively.

In summary, fault slippage can be attributed to the combined effect of pore pressure elevation in the fault and ambient stress change of the rock. At the fracturing stage, pore pressure change has more effect on inducing fault slippage.While at backflow stage,the fault plane is under a critical stress state, and any minor disturbance will trigger fault slippage, including pore pressure elevation and release of injection-induced stress change in the rock mass. This finding indicates that once a backflow operation is conducted after ceasing injection,the fluid may flow into the fault,disturb the fault which is already under critical state after fracturing, and thus trigger additional microseismicity. It is hard to control the flow of the injected fluid as well as the rock stress.Therefore,the best solution is to avoid encountering a pre-existing fault when choosing the fracturing site.However,if crossing a fault is inevitable, controlling the injected fluid volume is one feasible way for mitigating seismic hazards.

6. Conclusions

A 3D discrete element modeling of fault reactivation and induced seismicity during hydraulic fracturing injection and backflow is carried out. The study analyzes the reason why microseismicity occurs during backflow and reveals the different mechanisms for induced seismicity at fracturing and backflow stages.The modeling results show that the maximum magnitude of seismic events occurs at fracturing stage. The increment of fluid volume in the fault determines the cumulative moment and maximum fault slippage,both of which are essentially proportional to the fluid volume. The energy budget analysis indicates that the aseismic deformation takes up a majority of the total deformation energy during hydraulic fracturing.

There is a difference between the mechanisms of fracturing-and backflow-induced seismicities, which is reflected in the contributions of pore pressure and ambient rock stress.At fracturing stage,pore pressure change shows a dominant effect on inducing fault slippage.While at backflow stage,because the fault plane is under a critical stress state, any minor disturbance would trigger a fault slippage.

It is noted that this paper only presents a preliminary study of the induced seismicity during fracturing and backflow. More parametric analyses can be carried out based on the established frame in this study. The synthetic seismic events may not be completely in accordance with field data due to the simplified assumptions and the limitation of model resolution. Inclusion of more complex natural fractures in the vicinity of fault may be more realistic for the simulation.

Declaration of competing interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This research is supported by the Key Innovation Team Program of Innovation Talents Promotion Plan by Ministry of Science and Technology of China(Grant No.2016RA4059),and National Natural Science Foundation of China(Grant Nos.41672268 and 41772286).The useful discussion with Jim Hazzard and Branko Damjanac in Itasca Consulting Group is also greatly appreciated.