Numerical studies of interactions between hydraulic and natural fractures by Smooth Joint Model

Accepted Manuscript Numerical studies of interactions between hydraulic and natural fractures by Smooth Joint Model Jian Zhou, Luqing Zhang, Zhejun Pan, Zhenhua Han PII:

S1875-5100(17)30335-9

DOI:

10.1016/j.jngse.2017.07.030

Reference:

JNGSE 2272

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 29 March 2017 Revised Date:

28 June 2017

Accepted Date: 30 July 2017

Please cite this article as: Zhou, J., Zhang, L., Pan, Z., Han, Z., Numerical studies of interactions between hydraulic and natural fractures by Smooth Joint Model, Journal of Natural Gas Science & Engineering (2017), doi: 10.1016/j.jngse.2017.07.030. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

1

Numerical studies of interactions between hydraulic and natural

2

fractures by Smooth Joint Model

3

Jian Zhou 1, Luqing Zhang 1,*, Zhejun Pan 2, Zhenhua Han 1, 3,

4

1.

Academy of Science, Beijing 100029, China

RI PT

5

Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and Geophysics, Chinese

2.

CSIRO Energy, Private Bag 10, Clayton South VIC 3169, Australia

7

3.

College of Earth Science, University of Chinese Academy of Sciences, Beijing 100049, China

SC

6

8

10 11

* Correspondence: [email protected]; Tel.: +86-10-8299-8642; Fax:

M AN U

9

+86-10-6201-0846 Abstract

The hydraulic fracturing technology is a key for enhancing the permeability of

13

reservoirs during the production of natural gas. Significant progress has been made in

14

hydraulic fracturing theory and practices, however, precise prediction of the fracture

15

patterns in naturally fracture reservoirs still requires further study. The mechanisms

16

influencing the interaction between hydraulic fracture (HF) and natural fracture (NF)

17

should be well understood to achieve a wider application of hydraulic fracturing. In

18

this paper, it is the first time to simulate the details of the interactions between

19

hydraulic fractures (HFs) and natural fractures (NFs) using two-dimensional particle

20

flow code (PFC2D) with embedded Smooth Joint Model (SJM). Firstly, the ability of

21

SJM to emulate the natural rock joint was validated, and the fluid-mechanically

22

coupled mechanism was introduced. Secondly, the interactions between a driven HF

23

and a pre-existing NF were simulated considering fracturing fluid injection in a

24

borehole. Lastly, the influence of approach angles and in-situ differential stress were

25

studied. It is found that the model is capable of simulating the variety of interactions

26

between HFs and NFs such as direct crossing, crossing with an offset and HFs

AC C

EP

TE D

12

ACCEPTED MANUSCRIPT arrested by NFs. Under high approach angles and high differential stresses, the HF tend

28

to cross pre-existing NFs; on the contrary, a HF is more favorable to propagate along

29

the NF and re-initiate at the weak point or the tip of NF. Moreover, our numerical

30

results agree well with the analytical results based on the modified Renshaw and

31

Pollards’ criterion. Therefore, this numerical method may become a useful and

32

powerful tool for optimizing fracture design in naturally fractured shale gas reservoirs.

33

RI PT

27

Keywords: Hydraulic fracturing; Pre-existing natural fracture; Numerical simulation;

35

Particle flow method; Smooth Joint model 1. Introduction

M AN U

36

SC

34

Shale gas, an unconventional gas resource, is widely explored and developed in the

38

North America and other places in the world in recent years. The production of shale

39

gas increased rapidly due to the successful application of hydraulic fracturing

40

stimulation which is designed to generate complex fractures network in low

41

permeability reservoir (Warpinski et al., 2005; Maxwell et al., 2002). It is found that

42

NFs can affect the propagation direction of induced HFs, resulting in crossing,

43

branching and offset at the NF and consequently lead to complex fractures network

44

(Lamont and Jessen, 1963; Blanton, 1982; Warpinski and teufel, 1987; Renshaw and

45

Pollard, 1995; Zhou et al., 2008; Gu and Weng, 2010; Liu et al., 2014). It is widely

46

accepted that pre-existing NFs are widely distributed in most shale reservoirs, such as

47

the Barnett shale (Gale et al., 2007). Due to the variety of in-situ stress states,

48

properties and attitude of NFs, the HF propagation path is undecided. Thus, more

49

details about fluid-driven fracturing processes, including the interaction between NFs

50

and HFs, fracture patterns are essential to be understood for achieving a wider

51

application.

AC C

EP

TE D

37

52

In the past few years, a large amount of experimental, theoretical and numerical

53

works have been carried out on the interaction between HFs and NFs to understand

ACCEPTED MANUSCRIPT the mechanics of hydraulic fracturing in naturally fractured reservoirs (Lamont and

55

Jessen, 1963; Blanton, 1982, 1986; Warpinski and teufel, 1987; Renshaw and Pollard,

56

1995; Zhou et al., 2008; Gu et al., 2011; Liu et al., 2014; Zhang and Jeffery, 2006,

57

2008; Xu et al., 2009; Zhao and Young, 2011; Weng et al., 2011; Weng, 2015; Zhou et

58

al., 2016a, 2016b; Yushi et al., 2016; Dehghan et al., 2017; Zeng and Yao, 2016; Zhang

59

et al., 2017). Lamont and Jessen (1963) conducted a series experiments to study the

60

fracture crossing phenomenon, and found that all pre-existing fractures with the angle

61

of inclination from 0º to 45ºcan be crossed by HFs. Daneshy (1974) stated that

62

existing flaws whose dimensions are small compared with the driven HF have no

63

effect on the HF propagation, while the interactions between larger-scale NFs and

64

NFs are complex. In addition, he argued that the weakest planes encountered in oil

65

reservoir are of the closed type, which can be crossed by HFs. In the laboratory

66

experiments on the fractured Devonian shale and hydrostone, Blanton (1986)

67

observed that the induced fracture either crossed the NF or opened it, partly

68

depending on the approach angle and differential stress. These experiments suggest

69

that high approach angles and high differential stress favor crossing, otherwise, favor

70

opening. Warpinski and teufel (1987) studied the effect of single artificial fracture on

71

extending HF and suggested that HF crossed the artificial fracture at angles of 60ºor

72

higher under the high differential stress of 10MPa. Zhou et al. (2008) performed a

73

number of hydraulic fracturing experiments in which three different surface roughness

74

papers emulated the pre-existing fractures with different frictional coefficients. Three

75

scenarios such as arrest, dilate and cross were observed, depending on the approach

76

angles, in-situ differential stress and shear strength of NFs. Liu et al. (2014) based on

77

experiments concluded that HF propagates following the most preferential

78

propagation, the least resistance and the shortest propagation path. Recently, some

79

sophisticated technologies, such as X-ray computer tomography technology (Zuo

80

yushi, 2016; Guo tiangui, 2016) and acoustic emission technology (Stanchits et al.,

81

2015) have successfully been used to investigate the interactions between HF and NF

AC C

EP

TE D

M AN U

SC

RI PT

54

ACCEPTED MANUSCRIPT 82

in the laboratory. Based on the experimental results aforementioned, various criteria have been

84

developed to predict the propagation direction of induced fractures (Blanton, 1986;

85

Warpinski and teufel, 1987; Renshaw and Pollard, 1995; Gu and Weng, 2010).

86

Sarmadivaleh and Rasouli (2014) presented a summarization of the analytical criteria

87

for interaction between hydraulic and natural fractures, including criteria of Blanton,

88

Warpinski and teufel, Renshaw and Pollard and modified Renshaw and Pollard. There

89

are some basic assumptions among these analytical models, such as the flat geometry

90

of a single NF and uniformed roughness, and impermeable sides of the NF. However,

91

they are still useful to obtain an initial knowledge of the interaction mechanism.

SC

RI PT

83

Limited to the facts that laboratory work are expensive and work intensive and

93

analytical models apply simple assumptions, numerical simulations are often

94

employed.

M AN U

92

In the past decades, lots of numerical models have been developed to simulate

96

fluid-driven complex fracture networks in subsurface formation. However, this kind of

97

simulation is still full of challenge because the complexity of fluid-mechanically

98

coupled mechanisms, heterogeneity of reservoirs, and changing boundary conditions

99

have to be taken into account. To deal with the simple problems, the

100

Perkins-Kern-Nordgren model (PKN) (Perkins and Kern, 1961; Nordgren, 1972) and

101

the Khristianovic-Geertsma-de-Klerk model (KGD) (Khristianovic and Zheltov, 1955;

102

Geertsma, 1969) have been developed. With the development of computing power of

103

computers, the interactions between NFs and HFs in geomechanics-reservoir models

104

were studied. McClure (2012) investigated the fracture propagation in a pre-existing

105

discrete fracture network based on two-dimensional displacement discontinuous

106

method (DDM). Wu and Olson (2013), using the modified 2D-DDM approach, studied

107

the influence of stress shadow on the fracture extending pattern when considering

108

several perforation clusters generated in a horizontal well. Chen (2012) used the

109

commercial finite element code, ABAQUS, where the stress intensity model is replaced

AC C

EP

TE D

95

ACCEPTED MANUSCRIPT by a cohesive-zone fracture tip model, to simulate fracture propagation driven by fluid,

111

and their modeling results agreed with the 2D PKN and KGD solutions. Wang (2015)

112

also used an ABAQUS embedded extended finite element method and cohesive zone

113

method to model fracture initiation and propagation in different brittle rocks, and the

114

results show that fluid pressure field and fracture geometry are significantly affected by

115

the rock in-elastic deformations.

RI PT

110

In recent years, the discrete element method (DEM) (Cundall, 1971) has been

117

developed to investigate the process of hydraulic fracturing. In DEM the model

118

represents an assemblage of deformable blocks or rigid particles, which are in contact

119

and connected by bonds to represent continuum. Compared with continuum methods,

120

such as FEM, the principle characteristic of DEM is that the blocks or particles can

121

break if the bond reaches its strength, and that they can contact each other again. Thus,

122

DEM is more natural to simulate fracture evolution processes (initiation, propagation,

123

intersection and coalescence). The two dimensional particle flow code (PFC2D) is a

124

typical DEM (Potyondy and Cundall 2004), where the model is composed of an

125

assemblage of bonded circular particles, whereas the complex empirical constitutive

126

model can be substituted by a simple contact logic. In recent years it has been widely

127

applied in geomechanics to model failure processes of geomaterials (Zhang and Wong,

128

2012; Li et. al., 2012), and hydraulic fracturing (Hazzard et al., 2002; Al-Busaidi et al.,

129

2005 ;Zhao and Young, 2011; Shimizu et al., 2011; Eshiet et al., 2013; Yoon et al.,

130

2015a, 2015b; Zhou et al., 2016a, 2016b). Yoon et al. used the PFC2D and studied the

131

seismicity during developing an Enhanced Geothermal System by fluid injection into

132

deep underground (Yoon et al., 2015a), and simulated multistage hydraulic fracturing

133

processes in a fractured reservoir, then investigated the optimization of a HF network

134

according the stress shadow (Yoon et al., 2015b). Zhou et al. (2016a, 2016b) based on a

135

modified fluid-mechanically coupled mechanism studied hydraulic fracturing

136

processes in isotropic and laminated reservoir, and this model has been validated

137

comparing numerical value of breakdown pressures against analytical solutions. In

AC C

EP

TE D

M AN U

SC

116

ACCEPTED MANUSCRIPT 138

addition, some influence factors on the driven fracture patterns and geometries were

139

analyzed such as fluid viscosity, fluid injection rate and in-situ stress state. Although the bonded particle method has been used in hydraulic fracturing studies as

141

aforementioned, the modeling of the details of interactions between HFs and NFs by

142

PFC2D is less reported, which is a basic and important study for wider application of

143

bonded particle method in hydraulic fracturing. In this article, the Smooth Joint Model

144

(SJM) is firstly used to investigate the details of interactions between HFs and NFs in

145

PFC2D based on our fluid-mechanically coupled code (Zhou et al., 2016a, 2016b). First

146

of all, the SJM was introduced and validated for its ability to mimic shear behavior of

147

NF; secondly, the processes of Direct Across, Across with an Offset, and Arrest were

148

modeled and analyzed in details; thirdly, a series of numerical simulations was carried

149

out to analyze the influence of approach angle and differential stress on the interaction

150

between hydraulic and NF, and the results were compared to the analytical solutions of

151

modified Renshaw and Pollard’s Criteria.

152

2. Numerical modeling methodology

153

2.1. Particle flow code

TE D

M AN U

SC

RI PT

140

PFC2D is one kind of distinct element modeling method in which the solid materials

155

represent as an assembly of circular particles. Although the PFC2D is based on the

156

discontinuum method, after adding some bond models at the contacts between round

157

particles, it can also be used to model the deformation behaviors of the continuum

158

exactly. The bond with the properties of normal and shear stiffness, and shear and

159

tensile strength can simulate deformation and micro-cracks developing based on the

160

relationship presented by Potyondy and Cundall (2004). In rock mechanics research the

161

parallel bond model is one of most widely used, for which microscale properties and

162

deformation and failure behaviors are presented in Fig. 1. Normally, the Young’s

163

modulus E of emulated rock sample is related to the specified contact micro-stiffness;

164

the Poisson’s ratio ν is affected by the ratio of normal-to-shear stiffness. The microscale

AC C

EP

154

ACCEPTED MANUSCRIPT parameters in this method are different from these macroscale parameters such as E and

166

ν, which can be directly measured in laboratory. So, before modeling these parameters

167

have to be calibrated according to the laboratory’s results from confined/unconfined

168

compressive test and Brazilian indirect tensile test. Generally, by adjusting the

169

micro-stiffness and micro-strength of particles the realistic rock can be reproduced.

170

Under the applied load, when the maximum shear stress or tensile stress acting on the

171

bond exceeds the shear or tensile strength, the bond will breaks in the way of shear or

172

tensile, resulting in shear or tensile micro-cracks, respectively. By further generation of

173

micro-cracks, a macro-fracture can be formed by linking of these individual

174

micro-cracks.

AC C

EP

TE D

M AN U

SC

RI PT

165

175 176

Fig. 1. Bonded particle model in PFC2D. (a) Schematic of bonded particle model. (b)

177

The micro-scale parameters of the bond at contacts, and (c) the deformation and

178

failure mechanisms of the bond under tensile and shear stress. (Modified from

179

Bahaaddini et al., 2015)

180

ACCEPTED MANUSCRIPT 181

2.2. Smooth Joint Model The general method to emulate natural joints in PFC is the bond removal method,

183

in which the particle contacts lying on a joint track are left unbonded. This approach

184

has been used to investigate the shear behavior of rock joints in a number of studies

185

(Asadi et al., 2012; Park and Song, 2013; Jing and Stephansson, 2007). Bahaaddini et

186

al. (2015) argued that the bond removal method is limited to mimic the realistic shear

187

deformation characteristics of rock joints because of the circular shape, unequal size

188

of the particles and non-uniform distribution. In order to overcome the shortcomings

189

of this approach Pierce et al. (2007) introduced the Smooth Joint Model into PFC.

190

Recently, the SJM was used in the modeling of hydraulic fracturing by Yoon et al.

191

(2015a, 2015b). However, its ability of modeling of interactions between HFs and NFs

192

never be validated.

M AN U

SC

RI PT

182

As shown in Fig. 2(a), the SJM emulates the behavior of macroscale rock joint

194

(black line) by a series of micro-scale surfaces (short red lines) at the contacts

195

between circular particles. The parallel bonds are removed at these contacts and new

196

bonding model is assigned with pre-defined orientation and each pair of contacting

197

particles can overlap and pass through each other (Piece et al., 2007). Since the details

198

of fundamental algorithm of Smooth joint model can be acquired form the software

199

manuals (Itasca Consulting Group, 2008), only a summary of this model is given here.

200

In SJM, the relative displacement U and force F at contact point can be given as

202 203

EP

AC C

201

TE D

193

follows (Itasca Consulting Group, 2008): U = U n nˆ j + Us

(1)

F=Fn nˆ j + Fs

(2)

204

where nˆ j is the normal unit vector at contact as shown in Fig. 2(a). The positive

205

values of Un and Fn represent the overlapping of particles and compression force

206

respectively; Us and Fs represent the shear displacement vector and shear force vector

207

respectively. As shown in Fig. 2(b), the force-displacement relationship of each SJM

ACCEPTED MANUSCRIPT contact point is assumed following the Coulomb sliding model with dilation.

209

Micro-scale parameters such as SJM normal stiffness, k nj , SJM shear stiffness, k sj ,

210

SJM friction coefficient, µ j and SJM dilation angle, ϕ j , reflect the mechanical

211

behaviors of SJM contact with area of cross section A.

RI PT

208

212

The increments of joint force are calculated by the elastic components of the

213

displacement increment multiplying the SJM normal and shear stiffness. The normal

214

force and shear force are updated as:

216

Fs′ := Fs − k sj A∆U s

217 218

(3)

SC

Fn := Fn + k nj A∆U n

(4)

M AN U

215

If Fs′ ≤ ( Fs* = µ j Fn ) then Fs = Fs′ . If Fs′ > Fs* , then sliding occurs at the SJM contact and normal force and shear force are updated as: Fn := Fn + [∆U s* tan ϕ j ]knj A = Fn + (

220

Fs = Fs*

Fs′ − Fs* k sj

)knj tan ϕ j

(5)

(6)

AC C

EP

TE D

219

221 222

Fig. 2. Smooth Joint Model: (a) Sketches of SJM, (b) force-displacement relationship at

223

SJM contacts (Itasca Consulting Group Inc., 2008)

ACCEPTED MANUSCRIPT 224 In order to investigate the ability of the smooth joint model to reproduce the shear

226

behavior of rock joints, a series of direct shear tests were simulated. Fig. 3 shows the

227

model of which the length and width are 0.4m and 0.4m respectively. A joint modeled

228

by SJM was generated in the model center, which divided the model to upper and

229

lower parts. The lower part was fixed by three non-movable walls, and a constant

230

normal stress was applied on the up wall of the upper part, while a constant horizontal

231

speed of 0.001m/s was applied to two side walls of upper part (Fig. 3). The shear

232

displacement was recorded at the center of left wall of upper part, and the shear force

233

was recorded according the unbalance force of right wall of lower part. The

234

coefficient of friction of joint was set to 0.38, and its cohesion was 0. The direct shear

235

test results were shown in Fig. 4. Under four different normal stresses, the shear stress

236

increase until to a stable value with shear displacement increasing. The stable shear

237

stress are 1.14MPa, 1.90MPa, 2.66MPa and 3.80MPa corresponding to the normal

238

stress of 3MPa, 5MPa, 7MPa and 10MPa, respectively (Fig. 4(a)). The linear

239

regulation relationship between shear stress and normal stress shows that the friction

240

coefficient is 0.38 (Fig. 4(b)), which is the same as the input. According to the test

241

results, it can be concluded that the SJM can well reproduce the shear behaviors of

242

rock joint.

AC C

EP

TE D

M AN U

SC

RI PT

225

243 244

Fig. 3. Rock model with a joint in the center for direct shear test to validate the

ACCEPTED MANUSCRIPT 245

ability of SJM

(10MPa)

2.66

(7MPa)

1.90

(5MPa)

1.14

0

(3MPa)

0.5 1 1.5 2 Shear displacement (mm)

2.5

3.8

4 3.5 3 2.5 2 1.5 1 0.5 0

2.66 1.9

1.14

2

RI PT

3.80

y = 0.38 x + 0.00 R² = 1.00

4 6 8 Normal stress (MPa)

SC

4 3.5 3 2.5 2 1.5 1 0.5 0

Shear stress (MPa)

Shear stress (MPa)

246

(b)

M AN U

(a) 247

Fig. 4. The results from direct shear test. (a) Shear stress vs shear displacement under

248

normal stress of 3MPa, 5MPa, 7MPa and 10MPa respectively; (b) The linear

249

regression relationship between shear stress and normal stress.

250

252

2.3. Fluid flow in bonded granular matrix

TE D

251

The original concept of fluid flow algorithm in bonded granular matrix was first proposed by Cundall (Hazzard, et al., 2002), which is based on the network flow

254

model. There are two assumptions in Cundall’s fluid flow algorithm. One is that fluid

255

domain is surrounded by contacted round particles, and the other is that these fluid

256

domains are linked by flow channels at the contacts. This algorithm has been modified

257

and applied in a number of studies (Al-Busaidi et al., 2005; Hazzard et al., 2002;

258

Shimizu et al., 2011; Zhao and Young, 2011; Eshiet et al., 2013; Yoon et al. 2014,

259

2015a, 2015b, 2016; Chang M and Huang R C, 2015, Zhou et al., 2016a, 2016b,

260

2017). The details of fluid-mechanically coupled mechanism can be found in the

261

literatures (Al-Busaidi et al., 2005; Zhou et al., 2016a, 2016b, and 2017), so it is not

262

described in this section. However, two assumptions need to be mentioned. One is

263

that the pre-existing NF and rock matrix have the same permeability if the pre-exiting

264

NF has no failure, or otherwise, its permeability increases two orders of magnitude of

AC C

EP

253

10

ACCEPTED MANUSCRIPT 265 266

its original value. 3. Validation of SJM modeling HF-NF interactions To form complex patterns of HF network, which results from the interaction

268

between HF and NF, is one of the critical factors during shale gas development.

269

Numerous NFs exist in shale formation, which are closed, open or sealed by

270

precipitated materials such as quartz and calcite. When the HF encounters NF, the

271

direction of its propagation is very difficult to predict owing to the factors including

272

in-situ stress, approach angle, cohesion and interfacial friction coefficient, fracturing

273

fluid viscosity and fluid injection rate. HF-NF interactions have been investigated in

274

the field and laboratory, and using numerical simulations (Blanton, 1982; Renshaw

275

and Pollard, 1995; Zhang and Jeffery, 2006; Zhou et al., 2008, 2010; Gu and Weng,

276

2010; Liu et al., 2014). The results indicate that three possibilities will take place

277

when HFs intersect with NFs (Fig. 5): (a) the HF propagates and crosses the NF

278

directly without change of direction; (b) the HF crosses the NF with an offset, where

279

HF will deflect into NF for a certain distance and resume its propagation along the

280

maximum stress direction at some weak points of NF; (c) the HF is arrested and

281

propagates along the path of NF until to its end and then makes NF propagating into

282

the reservoir.

AC C

EP

TE D

M AN U

SC

RI PT

267

283 284

Fig.5. Possible sceneries when HFs intersect with NFs. (a) The HF cross the NF; (b)

285

the HF cross the NF with an offset; (c) the HF is arrested by the NF and propagates

286

off at the end of the NF.

287

ACCEPTED MANUSCRIPT In this section these three interaction processes when a HF approaching a NF were

289

simulated by PFC2D with SJM emulating NF. As presented in Fig. 6, a rock model was

290

generated with injection borehole in the center for viscous fluid injection, which is an

291

assembly of 22600 circular particles with their radii ranging from 2 cm to 3 cm. The

292

height of rock specimen is 6 meters and its width is 8 meters. A NF with length of 2

293

meters is simulated by SJM, of which the cohesion is close 0 and the friction angle is

294

30.96°. The dip angle of the NF is θ, which is the same as approach angle of HF while

295

the direction of maximum principle stress is in horizontal direction. In order to mimic

296

the in-situ stress, a kind of servo-mechanism was applied by using four movable walls

297

surrounding the rock model (Itasca Consulting Group Inc., 2008). The maximum

298

principle stress σ1 is in horizontal direction, and the minimum principle stress σ3 is in

299

vertical direction. The microscale parameters are listed in Table 2, which are

300

optimized according to macroscale mechanical properties (Table 1) of cement mortar

301

tested by Zhou et al. (2008). The stress σ1 is kept constant at 10.0 MPa, while σ3 is

302

5.0 MPa. In the simulation, the fluid viscosity is 1.0×10−3 Pa·s, and the rate of

303

fracturing fluid injection is 2.0×10−4 m3/s. With continuous fracturing fluid injection,

304

the HF initiates in the surrounding rock of the injection borehole, and propagates

305

along the direction of the maximum principle stress. Under different dip angle of NF,

306

the HF will cross NF directly, cross NF with an offset or be arrested by NF.

SC

M AN U

TE D

EP

AC C

307

RI PT

288

SC

RI PT

ACCEPTED MANUSCRIPT

308

Fig. 6. Schematic of a HF approaching NF and numerical model

310 311

Table 1 The macroscale mechanical properties of cement mortar (Zhou et al., 2008) Cement mortar parameters Young’s modulus

Permeability Porosity

Laboratory values

Numerical values

E

GPa

8.40

8.34

0.23

0.23

σc

MPa

28.34

28.7

σt

MPa

—

3.73

K

mD

0.1

—

1.85%

—

φ

EP

Table 2 The input microscale parameters for PFC model Microscopic input parameters

Symbols

Ratio of particle radius

Rmin/ Rmax

AC C

313

Units

TE D

Unconfined compressive strength Tensile strength

Symbols

ν

Poisson’s ratio

312

M AN U

309

Lower bound of particle radius

Rmin

Units

Values 2.6

mm

4.0 3

Particle density

ρ

kg/m

2350.0

Young’s modulus of the particle

Ec

GPa

6.7

Ratio of stiffness of the particle

kn/ks

2.6

Friction coefficient of particle

µ

0.50

Young’s modulus of the parallel bond

Ec

Ratio of stiffness of the parallel bond

kn / ks

Tensile strength of the parallel bond

σt

MPa

8.0

Cohesion of the parallel bond

c

MPa

18.0

Friction angle of the parallel bond

φ

°

30.0

Radius multiplier

λ

1.0

Moment contribution factor

β

0.2

GPa

6.7.0 2.6

ACCEPTED MANUSCRIPT Normal stiffness of SJM

kn

GPa/m

3000.0

Shear stiffness of SJM

ks

GPa/m

500

Tensile strength of SJM

σ t′

MPa

0.05

Friction angle of SJM

φ′

°

30.96

Cohesion of SJM

c′

MPa

0.05

314 3.1. HF arrested by NF

RI PT

315

An example that a produced HF was arrested by a pre-existing NF was simulated.

317

In this case, the HF approaches the NF with an incidence angle of 30°. The simulated

318

results are presented in Fig. 7 including a borehole pressure history curve and fracture

319

patterns and pore pressure distributions at 5 different times (0.6s, 1.2s, 1.8s, 3.0s and

320

3.6s) in fracturing process. The borehole pressure history curve presented in Fig. 7

321

shows that with continuous injection of fracturing fluid, the borehole pressure

322

increased firstly, and then slowly dropped down to a certain stress level after a macro

323

HF was generated in surrounding rock. The fracture patterns in Fig. 7 were formed by

324

newly induced microscale HFs and reinitiated the pre-exiting NFs. In models, almost

325

all newly induced fractures were tensile fractures shown by magenta short lines;

326

pre-existing HF reinitiated in slippage presented by black short lines. The scenario of

327

HF arrested by NF is described clearly by these fracture patterns at different times.

328

Initially, an induced bi-wing fracture symmetrically propagated in both directions from

329

the borehole in maximum principle direction; fracturing fluid was penetrating into

330

surrounding rock from the wall of borehole and induced fracture. Although HF did not

331

encounter NF, most microscale parts of the pre-existing HF had slipped due to the

332

deformation caused by near-tip stress field of the HF. Under continuous fracturing fluid

333

driving, the bi-wing fracture propagated along its original direction and pore pressure

334

distributed along the HF. When the HF encountered the NF, fracturing fluid leak off

335

into the damaged NF more easily, which results in NF opening. The opening NF

336

became a part of HF. The fluid diversion causes HF propagating along NF until one of

337

its tip. The fracture patterns at 3.0s and 3.6s show that the HF restarts at one tip of NF

338

with its orientation rerouted from the dip direction of NF to the direction of maximum

AC C

EP

TE D

M AN U

SC

316

ACCEPTED MANUSCRIPT 339

principle stress.

340

SC

1.8s

M AN U

0.6s

3.0s

RI PT

1.2s

3.6s

Fig. 7. The borehole pressure history curve, fracture patterns and pore pressure

342

distributions at different times during the process of HF arrested by NF when the dip

344 345

angle of NF is 30°.

3.2. HF crossing NF with an offset

EP

343

TE D

341

Fig. 8 presents an example where a HF crossed a NF with an offset. In this case,

347

except that the dip angle of NF is 60°, the other parameters are the same as those in

348

section 3.1. The shape of the borehole pressure history curve presented in Fig. 8 is

349

similar to the one in Fig. 7. It also shows that the borehole pressure increased firstly,

350

and then slowly dropped down to a limit level. According to the simulated results in

351

Fig. 8, including fracture patterns and pore-pressure distributions in model at different

352

times, the process of a HF crossing pre-existing fracture with an offset can be well

353

discussed. Initially, the fracture pattern at 0.6s shows the fracturing fluid drove a

354

bi-wing fracture in horizontal direction on both sides of borehole, and also shows

355

some local parts of the HF have slipped (shown by black short lines). Under fluid

AC C

346

ACCEPTED MANUSCRIPT continuous injection condition, the bi-wing fracture propagated resume along the

357

horizontal direction until it encountered the NF. At the intersection fracturing fluid

358

diverted to the pore of NF, which decreased the effective normal stress and opened it

359

consequently. However, the fracture pattern and pore pressure distribution at 1.8s

360

show fracturing fluid did not completely divert into one branch of NF as in Fig. 7

361

although the NF has deboned. The reason is that normal constraint force acting on NF

362

limits its dilated deformation and the real permeability of contact part of NF is still low.

363

From 1.8s to 3.0s, the borehole pressure increased which drives the NF going to dilate

364

until near its tip. Because of complexity of stress state the HF re-initiates not at the NF

365

tip but at a weak point along the NF, then keeps propagating close to the original

366

direction.

M AN U

SC

RI PT

356

3.0s

EP

TE D

1.2s

367 368

AC C

0.6s

1.8s

3.6s

Fig.8. The borehole pressure history curve, fracture patterns and pore pressure

369

distributions at different times during the process of HF crossing NF with an offset

370

when the dip angle of NF is 60°.

371

ACCEPTED MANUSCRIPT 372

3.3. HF direct crossing NF The scenario of HF direct crossing NF is presented in Fig. 9. In this model the dip

374

angle of NF is 80°. Comparing borehole pressure history curves, it can be found that

375

the breakdown pressure in this scenario is a little higher than other twos’ in Fig. 7 and

376

Fig. 8. The fracture pattern and pore pressure distribution at 0.6s in Fig. 9 also show

377

that HF propagated in horizontal direction as a form of bi-wing fracture. At that time

378

the NF was not disturbed by the near-tip stress filed of HF. The fracture pattern at 1.2s

379

in Fig. 9 suggests the right branch of HF approached NF and crossed it directly.

380

Although the local part of NF deboned in the intersection, the fracturing pressure is

381

not high enough to overcome the constraint force perpendicular to the NF. In this

382

scenario, the HF propagates without interruption, maintaining its orientation in the

383

direction of maximum principle stress.

TE D

M AN U

SC

RI PT

373

1.2s

AC C

EP

3.0s

0.6s

1.8s

3.6s

384 385

Fig. 9. The borehole pressure history curve, fracture patterns and pore pressure

386

distributions at different times during the process of HF crossing NF directly when the

387

dip angle of NF is 80°.

ACCEPTED MANUSCRIPT 388

4. Influence of approach angle and in-situ differential stress The interaction between HF and NF is a complex process, which is related to lots

390

of influence factors, such as formation rock mechanical properties and strength of NF,

391

approach angle, in-situ stress, and permeability of NF. According to the scenarios

392

described in the preceding section, the influence of approach angle and in-situ stress

393

are studied in this section, which tend to provide a basis conclusion for selecting

394

stimulation methods in hydraulic fracturing process. 4.1. Influence of approach angle

SC

395

RI PT

389

Firstly, the influence of approach angle on the HF propagation is studied. The

397

model used is the same as aforementioned in section 3. The in-situ stress ratio in the

398

modeling is chosen as σ1 σ3 =2.0, and σ1 is 10MPa in the horizontal direction. The

399

viscous fluid with viscosity of 1.0×10-3Pa·s was injected from the borehole at the rate

400

of 2.0×10-4m3/s. Generally, the HF propagates along the direction of maximum

401

principle stress σ1 in homogeneous rocks. By changing the dip angle of NF, the HF

402

approaching the pre-existing NF at different angle were simulated. Fig. 10 presents

403

the simulation results that include the fracture patterns and pore-pressure distributions

404

in models where the approach angles are about 40°, 50°, 60°, 70°, 80° and 90°

405

respectively. When the approach angles were about 40° and 50°, the driven HFs were

406

arrested by the NF and continuously propagated from the NF tips. While, the HF

407

crossed the NF with an offset between the intersection and NF tip, if the approach

408

angles were about 60° and 70°, respectively. With the approach angle increasing (θ=80°

409

or 90°), the HF crossed the NF directly at the intersection. Therefore, the results suggest

410

that for non-orthogonal approach, the approach angle will greatly affect the fluid driven

411

fracture pattern. With small approach angle, it is more favorable for opening and

412

dilating the NF and generating complex fracture network; on the contrary, the HF tend

413

to cross the HF and remain bi-wing planar geometry.

AC C

EP

TE D

M AN U

396

θ=50º

θ=60º

M AN U

SC

θ=40º

RI PT

ACCEPTED MANUSCRIPT

θ=70º

θ=80º

θ=90º

414

Fig. 10. The interactions between produced HF and pre-exiting NF with different dip

415

angle θ

416

4.2. Influence of in-situ differential stress

Secondly, the effect of differential in-situ stress is investigated. In the simulation,

418

the model with a pre-exiting fracture with the dip angle of 60ºwas chose. The

419

maximum principle stress σ1 is constant as 10MPa, and the minimum principle stress

420

σ 3 varies from 2 to 8MPa, so that the conditions under various in-situ differential

421

stress were generated. The fluid viscosity and fluid injection rate are the same as

422

aforementioned, respectively. The simulation results are shown in Fig. 11. It is found

423

that the HFs initiated and propagated mainly along the direction of σ1 before they

424

intersected with NFs. While, under low in-situ differential stress, the direction of the

425

initiation of the HF propagation was a small offset from the direction of the far field

426

maximum principle stress and consistent with the parallel direction of NF; when HF

427

came close to the NF it rotated to the normal direction of NF; at last, the HF was

428

arrested. This modeling phenomenon is consistent with the experimental results of Liu

429

et al. (2014). The reason is the transformation of direction of local maximum principle

430

stress caused by relative deformation of two walls of NF. With increasing of

AC C

EP

TE D

417

ACCEPTED MANUSCRIPT differential stress up to 6MPa, the HF crossed the NF with an offset from the

432

intersection. Under larger differential stress (∆σ = 7.0MPa, 8.0MPa), the NF prevents

433

the propagation of the right wing of the generated bi-wing fracture even though

434

fracturing fluid penetrated into and opened some local parts of NF. Last two figures in

435

Fig. 11 show that the end of induced fractures are dry. This phenomenon can be

436

attributed to the facts that under higher in-situ differential stress state HF favors more

437

to propagation in direction of maximum principle stress and it less likely to open the

438

NF, and the fluid cannot leak off into the NF consequently.

RI PT

431

To summarize, the influence of approach angle and in-situ differential stress on

440

HF pattern mainly due to the normal stresses on the NF. With increasing of approach

441

angle and in-situ differential stress, the normal stress on the NF increases, which are

442

beneficial for the driven fracture to cross the NF and disadvantage for the opening of

443

NF and fracturing fluid penetrating into it.

TE D

M AN U

SC

439

σ3 =6.0MPa (∆σ = 4.0MPa)

σ3 =5.0MPa (∆σ = 5.0MPa)

σ3 =3.0MPa (∆σ = 7.0MPa)

σ3 =2.0MPa (∆σ = 8.0MPa)

AC C

EP

σ3 =8.0MPa (∆σ = 2.0MPa)

σ3 =4.0MPa (∆σ = 6.0MPa)

444

Fig. 11. The interactions between produced HF and pre-exiting NF under a series of

445

differential stress (σ1 is 10MPa in horizontal direction)

446

ACCEPTED MANUSCRIPT 4.3. Numerical results compare against Modified Renshaw and Pollards’ Criteria

448

As aforementioned in previous section, the SJM can reproduce different scenarios

449

of interactions between the generated HF and NF. Tring to move forward in validation

450

of the ability of SJM in hydraulic fracturing modeling, a series of modeling were

451

performed and compared against analytical model under different approach angles

452

and in-situ stress states. Table 3 presents the numerical results comparing against the

453

analytical results according the Modified Renshaw and Pollards’ Criterion (Modified

454

R&P Criterion). The bold highlighted cells of the Table 3 show the cases at which the

455

model did not satisfy the modified Renshaw and Pollards’ criterion. It can be shown in

456

Table 3 that these numerical results agreed well with the analytical results based on

457

the modified Renshaw and Pollards’ criterion except in a few cases. Thus, the ability

458

of SJM modeling of HF-NF interaction can be validated again.

459

Table 3 Comparing Modified Renshaw and Pollards’ Criterion and our numerical results σ1 (MPa)

σ3 (MPa)

-10

-8

30

-10

-6

30

-10

-5

30

-10

-4

30

-10

30

-10

40

-10

σT (MPa)

τ(MPa)

τ0

Modified

µ

Numerical results R&P Criterion

-8.50

-9.50

0.87

0

0.6

No Crossing

Arrest

-7.00

-9.00

1.73

0

0.6

No Crossing

Arrest

-6.25

-8.75

2.17

0

0.6

No Crossing

Arrest

-5.50

-8.50

2.60

0

0.6

No Crossing

Arrest

EP

30

σn (MPa)

TE D

θ (º)

-3

-4.75

-8.25

3.03

0

0.6

No Crossing

Arrest

-2

-4.00

-8.00

3.46

0

0.6

No Crossing

Arrest

-8

-8.83

-9.17

0.98

0

0.6

No Crossing

Arrest

AC C

460

M AN U

SC

RI PT

447

40

-10

-6

-7.65

-8.35

1.97

0

0.6

No Crossing

Arrest

40

-10

-5

-7.07

-7.93

2.46

0

0.6

No Crossing

Arrest

40

-10

-4

-6.48

-7.52

2.95

0

0.6

No Crossing

Arrest

40

-10

-3

-5.89

-7.11

3.45

0

0.6

No Crossing

Arrest

40

-10

-2

-5.31

-6.69

3.94

0

0.6

Crossing

Arrest

50

-10

-8

-9.17

-8.83

0.98

0

0.6

No Crossing

Arrest

50

-10

-6

-8.35

-7.65

1.97

0

0.6

No Crossing

Arrest

50

-10

-5

-7.93

-7.07

2.46

0

0.6

No Crossing

Arrest

50

-10

-3

-7.11

-5.89

3.45

0

0.6

No Crossing

Arrest

50

-10

-2

-6.69

-5.31

3.94

0

0.6

Crossing

Arrest

50

-10

-4

-7.52

-6.48

2.95

0

0.6

No Crossing

Offset Crossing

60

-10

-8

-9.50

-8.50

0.87

0

0.6

No Crossing

Offset Crossing

ACCEPTED MANUSCRIPT -6

-9.00

-7.00

1.73

0

0.6

Crossing

Offset Crossing

60

-10

-5

-8.75

-6.25

2.17

0

0.6

Crossing

Offset Crossing

60

-10

-4

-8.50

-5.50

2.60

0

0.6

Crossing

Offset Crossing

60

-10

-3

-8.25

-4.75

3.03

0

0.6

Crossing

Offset Crossing

60

-10

-2

-8.00

-4.00

3.46

0

0.6

Crossing

Offset Crossing

70

-10

-6

-9.53

-6.47

1.29

0

0.6

Crossing

Offset Crossing

70

-10

-5

-9.42

-5.58

1.61

0

0.6

Crossing

Offset Crossing

70

-10

-4

-9.30

-4.70

1.93

0

0.6

Crossing

Offset Crossing

70

-10

-8

-9.77

-8.23

0.64

0

0.6

No Crossing

Direct Crossing

70

-10

-3

-9.18

-3.82

2.25

0

0.6

Crossing

Direct Crossing

70

-10

-2

-9.06

-2.94

2.57

0

0.6

Crossing

Direct Crossing

80

-10

-8

-9.94

-8.06

0.34

0

0.6

Crossing

Direct Crossing

80

-10

-6

-9.88

-6.12

0.68

0

0.6

Crossing

Direct Crossing

80

-10

-5

-9.85

-5.15

0.86

0

0.6

Crossing

Direct Crossing

80

-10

-4

-9.82

-4.18

1.03

0

0.6

Crossing

Direct Crossing

80

-10

-3

-9.79

-3.21

1.20

0

0.6

Crossing

Direct Crossing

80

-10

-2

-9.76

-2.24

1.37

0

0.6

Crossing

Direct Crossing

90

-10

-8

-10.00

-8.00

0.00

0

0.6

No Crossing

Direct Crossing

90

-10

-6

-10.00

-6.00

0.00

0

0.6

Crossing

Direct Crossing

90

-10

-5

-10.00

-5.00

0.00

0

0.6

Crossing

Direct Crossing

90

-10

-4

-10.00

-4.00

0.00

0

0.6

Crossing

Direct Crossing

90

-10

-3

-10.00

-3.00

0.00

0

0.6

Crossing

Direct Crossing

90

-10

-2

-10.00

-2.00

0.00

0

0.6

Crossing

Direct Crossing

5. Conclusions

M AN U

SC

RI PT

-10

TE D

461

60

The Smooth Joint Model shows value in modeling the behaviors of rock joints,

463

and it is capable of emulating the elastic deformation, shear failure and tensile failure.

464

After introducing a modified fluid-mechanically coupled model, the SJM is firstly used

465

to simulate the details of interactions between driven HFs and NFs in PFC2D, and its

466

ability has been investigated. Based on a serious simulation, we can arrive at

467

conclusions: PFC model can describe exactly the tensile failure of incipient HF in

468

surrounding rock of the borehole; Furthermore, the SJM is also able to reproduce the

469

variety of interactions between a driven HF and a NF such as direct crossing, crossing

470

with an offset and HF arrested by NF in hydraulic fracturing.

AC C

EP

462

471

The interactions between generated HFs and pre-existing NFs are further studied

472

considering different in-situ stress states and approach angles. The simulation results

ACCEPTED MANUSCRIPT indicate that the NF is more favorable for opening and diverting fracturing fluid under

474

low approach angles and low differential stresses, which results in that a HF

475

propagated along the NF and re-initiated at the weak point or the tip of NF. On the

476

contrary, under high approach angles and high differential stresses, the HF tend to cross

477

pre-existing NFs. Moreover, a series of modeling results under different approach

478

angles and in-situ stress states indicate it is a very good agreement comparing against

479

the analytical results based on the modified Renshaw and Pollards’ criterion except in

480

a few cases.

RI PT

473

To summarize, the SJM embedded in PFC2D was successfully applied to simulate

482

different scenarios of the interactions between HF and NF. This numerical method will

483

be a useful and powerful tool for investigating the hydraulic fracturing behavior and

484

optimizing fracture design in naturally fractured shale gas reservoirs.

M AN U

SC

481

485

Acknowledgments: This research is supported by funding from the National Natural

487

Science Foundation of China (Grants 41502307, 41572312, 41272353), China

488

Postdoctoral Science Foundation (Grants 2014M550838) and the Strategic Priority

489

Research Program of the Chinese Academy of Sciences (Grants XDB10030104).

490

References

491 492 493 494 495 496 497 498 499 500 501 502 503 504

Warpinski, N.R., Kramm, R.C., Heinze, J.R., Waltman, C.K., (2005). Comparison of singleand dual-array microseismic mapping techniques in the Barnett shale. In: Paper SPE 95568, SPE Annual Technical Conference and Exhibition, Dallas, Texas, 9–12 October. Lamont, N., & Jessen, F. W. (1963). The effects of existing fractures in rocks on the extension of hydraulic fractures. Journal of Petroleum Technology, 15(02), 203-209. Blanton, T.L., 1982. An experimental study of interaction between hydraulically induced and pre-existing fractures. In: SPE 10847, SPE/DOE Unconventional Gas Recovery Symposium, Pittsburgh, Pennsylvania, USA, 16-18 May. Blanton, T.L., 1986. Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. In: Paper SPE 15261, SPE Unconventional Gas Technology Symposium, Louisville, Kentucky, USA, 18-21 May. Warpinski, N.R., Teufel, L.W., 1987. Influence of geologic discontinuities on hydraulic fracture propagation (includes associated papers 17011 and 17074). SPE J. Pet. Technol. 39 (2), 209-220.

AC C

EP

TE D

486

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Renshaw, C.E., Pollard, D.D., 1995. An experimentally verified criterion for propagation across unbounded frictional interfaces in brittle, linear elastic-materials. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 32 (3), 237-249. Zhao, J., Chen, M., Jin, Y., Zhang, G., 2008. Analysis of fracture propagation behavior and fracture geometry using tri-axial fracturing system in naturally fractured reservoirs. Int. J. Rock Mech. Min. Sci. 45, 1143-1152. Gu, H., Weng, X., Lund, J.B., Mack, M., Ganguly, U., Suarez-Rivera R., 2011. Hydraulic fracture crossing natural fracture at non-orthogonal angles, a criterion, its validation and applications. In: Paper SPE 139984, Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, 24–26 January. Gu, H., Weng, X., Lund, J.B., Mack, M., Ganguly, U., Suarez-Rivera R., 2011. Hydraulic fracture crossing natural fracture at non-orthogonal angles, a criterion, its validation and applications. In: Paper SPE 139984, Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, 24–26 January. Liu, Z., Chen, M., & Zhang, G. (2014). Analysis of the influence of a natural fracture network on hydraulic fracture propagation in carbonate formations. Rock mechanics and rock engineering, 47(2), 575-587. Gale, J.F.W., Reed, R.M., Holder, J., 2007. Natural fractures in the Barnett shale and their importance for hydraulic fracture treatment. AAPG Bull. 91 (4), 603–622. Zhang, X., Jeffrey, R.G., 2006. The role of friction and secondary flaws on deflection and re-initiation of hydraulic fractures at orthogonal pre-existing fractures. Geophys. J. Int. 166 (3), 1454–1465. Zhang, X., Jeffrey, R.G., 2008. Reinitiation or termination of fluid-driven fractures at frictional bedding interfaces. J. Geophys. Res.-Sol. Ea. 113 (B8), B08416. Xu, W., Calvez, J.L., Thiercelin, M. 2009. Characterization of hydraulically-induced fracture network using treatment and microseismic data in a tight-gas formation: A geomechanical approach. In: Paper SPE 125237, SPE Tight Gas Completions Conference, San Antonio, Texas, USA, 15-17 June. Zhao, X., & Paul Young, R. (2011). Numerical modeling of seismicity induced by fluid injection in naturally fractured reservoirs. Geophysics, 76(6), WC167-WC180. Weng, X., Kresse, O., Cohen, C., Wu, R., Gu, H., 2011. Modeling of hydraulic fracture network propagation in a naturally fractured formation. In: Paper SPE 140253, SPE Hydraulic Fracturing Conference and Exhibition, Woodlands, Texas, USA, 24-26 January. Weng, X. (2015). Modeling of complex hydraulic fractures in naturally fractured formation. Journal of Unconventional Oil and Gas Resources, 9, 114-135. Yushi, Z., Xinfang, M., Shicheng, Z., Tong, Z., & Han, L. (2016). Numerical Investigation into the Influence of Bedding Plane on Hydraulic Fracture Network Propagation in Shale Formations. Rock Mechanics and Rock Engineering, 49(9), 3597-3614. Dehghan, A. N., Goshtasbi, K., Ahangari, K., Jin, Y., & Bahmani, A. 3D Numerical Modeling of the Propagation of Hydraulic Fracture at Its Intersection with Natural (Pre-existing) Fracture. Rock Mechanics and Rock Engineering, 1-20.

AC C

505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Zeng, Q., & Yao, J. (2016). Numerical simulation of fracture network generation in naturally fractured reservoirs. Journal of Natural Gas Science and Engineering, 30, 430-443. Zhang, Z., Li, X., He, J., Wu, Y., & Li, G. (2017). Numerical study on the propagation of tensile and shear fracture network in naturally fractured shale reservoirs. Journal of Natural Gas Science and Engineering, 37, 1-14. Daneshy AA (1974) Hydraulic fracture propagation in the presence of planes of weakness. Paper SPE 4852. SPE European spring meeting, American Institute of Mining, Metallurgical, and Petroleum Engineers, Inc., Amsterdam, Netherlands. Jia L, Chen M, Sun L, et al. (2013). Experimental study on propagation of hydraulic fracture in volcanic rocks using industrial CT technology. Petroleum Exploration and Development, 40(3): 405-408 Zuo Y., Zhang S., et.al. (2016). Experimental investigation into hydraulic fracture network propagation in gas shales using CT scanning technology. Rock Mechanics and Rock Engineering, 49(1), 33-45. Guo, T., Zhang, S., et.al. (2014). Experimental study of hydraulic fracturing for shale by stimulated reservoir volume. Fuel, 128, 373-380. Stanchits, S., Burghardt, J., & Surdi, A. (2015). Hydraulic Fracturing of Heterogeneous Rock Monitored by Acoustic Emission. Rock Mechanics and Rock Engineering, 48(6), 2513-2527. Sarmadivaleh, M., & Rasouli, V. (2014). Modified Reinshaw and Pollard criteria for a non-orthogonal cohesive natural interface intersected by an induced fracture. Rock mechanics and rock engineering, 47(6), 2107-2115. Perkins, T.K.; Kern, L.R. Widths of hydraulic fractures. J. Pet. Technol. 1961, 13, 937–949. Nordgren, R.P. Propagation of a vertical hydraulic fracture. Soc. Pet. Eng. J. 1972, 12, 306– 314. Khristianovic, S.; Zheltov, Y. Formation of vertical fractures by means of highly viscous fluids. In Proceedings of the 4th World Petroleum Congress, Rome, Italy, 6–15 June 1955; Volume 2, pp. 579–586. Geertsma, J.; De Klerk, F. A rapid method of predicting width and extent of hydraulically induced fractures. J. Pet. Technol. 1969, 21, 1571–1581. McClure, M. W. (2012). Modeling and characterization of hydraulic stimulation and induced seismicity in geothermal and shale gas reservoirs (Doctoral dissertation, Stanford University). Wu, K., & Olson, J. E. (2013). Investigation of Critical In Situ and Injection Factors in Multi-Frac Treatments: Guidelines for Controlling Fracture Complexity. In: Paper SPE 163821, SPE Hydraulic Fracturing Conference, Woodlands, Texas, 4 - 6 February. Chen, Z. (2012). Finite element modelling of viscosity-dominated hydraulic fractures. Journal of Petroleum Science and Engineering, 88, 136-144. Wang, H. (2015). Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using XFEM with cohesive zone method. Journal of Petroleum Science and Engineering, 135, 127-140.

AC C

547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Cundall, P.A. (1971). A computer model for simulating progressive, large-scale movements in blocky rock systems. In Proceedings of the International Symposium on Rock Fracture, Nancy, France, 4–6 October. Potyondy, D. O., & Cundall, P. A. (2004). A bonded-particle model for rock. International journal of rock mechanics and mining sciences, 41(8), 1329-1364. Zhang, X. P., & Wong, L. N. Y. (2012). Cracking processes in rock-like material containing a single flaw under uniaxial compression: a numerical study based on parallel bonded-particle model approach. Rock Mechanics and Rock Engineering, 45(5), 711-737. Li, W. C., Li, H. J., Dai, F. C., & Lee, L. M. (2012). Discrete element modeling of a rainfall-induced flowslide. Engineering geology, 149, 22-34. Hazzard, J.F.; Young, R.P.; Oates, S.J (2002). Numerical modeling of seismicity induced by fluid injection in a fractured reservoir. In Proceedings of the 5th North American Rock Mechanics Symposium Mining and Tunnel Innovation and Opportunity, Toronto, ON, Canada, 7–10 July; pp. 1023–1030. Al‐Busaidi A, Hazzard J F, Young R P. Distinct element modeling of hydraulically fractured Lac du Bonnet granite[J]. Journal of Geophysical Research: Solid Earth, 2005, 110(B6). Zhao, X., & Paul Young, R. (2011). Numerical modeling of seismicity induced by fluid injection in naturally fractured reservoirs. Geophysics, 76(6), WC167-WC180. Shimizu, H., Murata, S., & Ishida, T. (2011). The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution. International Journal of Rock Mechanics and Mining Sciences, 48(5), 712-727. Eshiet, K. I., Sheng, Y., & Ye, J. (2013). Microscopic modelling of the hydraulic fracturing process. Environmental earth sciences, 68(4), 1169-1186. Yoon, J. S., Zimmermann, G., & Zang, A. (2015a). Numerical Investigation on Stress Shadowing in Fluid Injection-Induced Fracture Propagation in Naturally Fractured Geothermal Reservoirs. Rock Mechanics and Rock Engineering, 48(4), 1439-1454. Yoon, J. S., Zimmermann, G., & Zang, A. (2015b). Discrete element modeling of cyclic rate fluid injection at multiple locations in naturally fractured reservoirs. International Journal of Rock Mechanics and Mining Sciences, (74), 15-23. Zhou, J., Zhang, L., Pan, Z., & Han, Z. (2016a). Numerical investigation of fluid-driven near-borehole fracture propagation in laminated reservoir rock using PFC 2D. Journal of Natural Gas Science and Engineering, 36, 719-733. Zhou, J., Zhang, L., Braun, A., & Han, Z. (2016b). Numerical Modeling and Investigation of Fluid-Driven Fracture Propagation in Reservoirs Based on a Modified Fluid-Mechanically Coupled Model in Two-Dimensional Particle Flow Code. Energies, 9(9), 699. Zhou, J., Zhang, L., & Han, Z. (2017). Hydraulic fracturing process by using a modified two-dimensional particle flow code-method and validation. Progress in Computational Fluid Dynamics, an International Journal, 17(1), 52-62.

AC C

588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Bahaaddini, M., Hagan, P. C., Mitra, R., & Hebblewhite, B. K. (2015). Parametric study of smooth joint parameters on the shear behaviour of rock joints. Rock Mechanics and Rock Engineering, 48(3), 923-940. Jing L, Stephansson O (2007). Fundamentals of discrete element methods for rock engineering: theory and applications. Elsevier, Amsterdam. Asadi M, Rasouli V, Barla G (2012). A bonded particle model simulation of shear strength and asperity degradation for rough rock fractures. Rock Mech Rock Eng 45:649–675. Park JW, Song JJ (2013). Numerical method for the determination of contact areas of a rock joint under normal and shear loads. International Journal of Rock Mechanics and Mining Sciences, 58:8–22. Pierce M, Cundall P, Potyondy D, Mas Ivars D (2007). A synthetic rock mass model for jointed rock. Paper presented at the rock mechanics: meeting society’s challenges and demands, 1st Canada-U.S. Rock Mechanics Symposium, Vancouver. Itasca Consulting Group Inc., 2008. PFC2D – Particle Flow Code in 2 Dimensions. Version 4.0, Minneapolis.

AC C

629 630 631 632 633 634 635 636 637 638 639 640 641 642 643

ACCEPTED MANUSCRIPT

Highlights:

(1) It is first time study the details of interactions between HFs and NFs by PFC2D with SJM;

RI PT

(2) SJM in PFC2D describes well different scenarios of interactions between HFs and NFs; (3) Our numerical results agreed well with the analytical results from the modified Renshaw and Pollards’ criterion;

AC C

EP

TE D

M AN U

SC

(4) Interaction types between HFs and NFs rely on approach angles and in-situ differential stress states;