PETROLEUM EXPLORATION AND DEVELOPMENT

Volume 47, Issue 5, October 2020

http://www.sciencedirect.com/journal/petroleum-exploration-and-development

Cite this article as: PETROL. EXPLOR. DEVELOP., 2020, 47(5): 1117–1130. RESEARCH PAPER

Numerical simulation of hydraulic fracture propagation in

laminated shale reservoirs
ZHOU Tong1, WANG Haibo1, LI Fengxia1, LI Yuanzhao2, ZOU Yushi3,*, ZHANG Chi2
1. Research Institute of Petroleum Exploration and Development, Sinopec, Beijing 100083, China;
2. Shale Gas Exploration and Development Co., Ltd., Sinopec, Chongqing 408014, China;
3. China University of Petroleum, Beijing 102249, China

Abstract: The main area of the Jiaoshiba anticline of the Fuling shale gas field was taken as the research object, laboratory rock me-
chanical experiments and direct shear experiments were conducted to clarify the mechanical anisotropy characteristics and parameters of
rock samples with rich beddings. Based on the experimental results, a 3D fracture propagation model of the target reservoir taking me-
chanical anisotropy, weak bedding plane and vertical stress difference into account was established by the discrete element method to
analyze distribution patterns of hydraulic fractures under different bedding densities, mechanical properties, and fracturing engineering
parameters (including perforation clusters, injection rates and fracturing fluid viscosity). The research results show that considering the
influence of the weak bedding plane and longitudinal stress difference, the interlayer stress difference 3–4 MPa in the study area can con-
trol the fracture height within the zone of stress barrier, and the fracture height is less than 40 m. If the influence of the weak bedding
plane is not considered, the simulation result of fracture height is obviously higher. Although the opening of high-density bedding frac-
tures increases the complexity of hydraulic fractures, it significantly limited the propagation of fracture height. By reducing the number of
clusters, increasing the injection rate, and increasing the volume and proportion of high-viscosity fracturing fluid in the pad stage, the re-
striction on fracture height due to the bedding plane and vertical stress difference can be reduced, and the longitudinal propagation of
fractures can be promoted. The fracture propagation model was used to simulate one stage of Well A in Fuling shale gas field, and the
simulation results were consistent with the micro-seismic monitoring results.

Key words: shale; lamina; hydraulic fracturing; fracture propagation law; longitudinal stress difference; Jiaoshiba anticline; shale gas reservoir

Introduction that can simulate propagation of multiple fractures, including

the line network model[1], discrete fracture network model[2–3],
Shale reservoirs are typically characterized by low porosity
unconventional fracture model[4], and other models based on
and ultra-low permeability. The key to effectively developing
the finite element method[5–6], boundary element method[7–10],
shale oil and gas is to form a large-scale complex fracture
network by multi-stage and multi-cluster fracturing in hori- extended finite element method[11–12], discrete element
zontal wells. Different from conventional reservoirs, shale, as method[13–16] and phase-field method[17–18]. The above models
a type of layered sedimentary rock, has strong heterogeneity. simplify the rock as an isotropic medium and primarily focus
Meanwhile, shale reservoirs have a large number of local dis- on the fracture distribution under the influence of high-angle
continuities, such as faults, bedding planes and natural frac- natural fractures. However, these models generally fail to
tures, which result in significant structural anisotropy that consider the influences of typical characteristics of shale such
exerts a certain impact on the hydraulic fracture propagation. as bedding and mechanical anisotropy on the morphology of
The research on the expansion mechanism of complex fracture propagation. As demonstrated in the laboratory frac-
fracture networks in shale reservoirs has become an important turing simulation and field fracture monitoring, bedding frac-
topic in the field of unconventional oil and gas development. tures have an important impact on the vertical extension and
In recent years, a lot of theoretical research on the propagation propagation morphology of hydraulic fractures[19–23]. When
of hydraulic fractures in shale reservoirs has been carried encountering bedding surfaces, hydraulic fractures may result
out, and a series of fracturing models have been established in penetration, diversion, termination, or stepped extension,

Received date: 28 Nov. 2019; Revised date: 18 Aug. 2020.

* Corresponding author. E-mail: [email protected]
Foundation item: Supported by the China National Science and Technology Major Project (2016ZX05060001-032).
https://doi.org/10.1016/S1876-3804(20)60122-7
Copyright © 2020, Research Institute of Petroleum Exploration & Development, PetroChina. Publishing Services provided by Elsevier B.V. on behalf of KeAi Com-
munications Co., Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

which cause uncertainty of the shape of the final fracture Formation, below which is a set of light-gray nodular lime-
network. As a result, most numerical models can not reflect stone of the Upper Ordovician Jiancaogou Formation. Ac-
the essential difference between hydraulic fracture propaga- cording to the reservoir physical properties, shale layers in the
tion in laminated shale reservoirs and that in naturally frac- Longmaxi and Wufeng Formations are divided into 9 sublay-
tured reservoirs. ers from bottom to top. Gray-black carbonaceous and sili-
Field coring data shows that shale reservoirs have highly ceous shale take dominance in Wufeng Formation-Longmaxi
dense bedding fractures (several to hundreds per meter) Formation. The main gas-bearing layers have rich beddings,
commonly. The strong anisotropy brings great challenges to and the numbers of beddings decrease from the bottom to the
the simulation of fracture propagation in shale reservoirs. It is top[28]. Among them, the sublayer ① is developed strongly in
accurate to analyze the influence of the weak surface of the bedding fractures, with hundreds of bedding fractures per
natural fracture on fracture propagation using the cohesive meter. In contrast, the sublayers ⑧ and ⑨ at the top have
element of the finite element[24]. Nevertheless, when the natu- fewer bedding fractures. The core observation results are
ral fracture density is set too high, the calculation stability and shown in Table 1. It can be seen that the high-angle natural
fractures in the reservoir are generally small in scale and not
convergence are significantly poor. For the shale with abun-
connected with each other.
dant weak bedding planes, the rock mass medium is more like
The vertical stress profile of Well Jiaoye A was interpreted
a discrete body of discontinuous media. As a numerical simu-
based on the corrected logging curve obtained during the
lation for discontinuous medium, the discrete element method
Kaiser in-situ stress test (Fig. 1). The minimum horizontal
has great advantages in dealing with large deformation prob-
principal stress varies in the longitudinal direction. The in-situ
lems of discontinuous medium such as rock and soil fractures
stresses of sublayers ①–④ fluctuate within a small range
(like highly complex fracture network)[25]. Zhao et al.[13] used
(only 1–2 MPa), with an average of 49 MPa; the in-situ stress
the two-dimensional particle discrete element method to sim-
of sublayer ⑤ is significantly higher (up to 53 MPa), and
ulate the behavior of hydraulic and natural fractures, in which there is a stress difference of 3–4 MPa between sublayer ⑤
they simulated natural fractures by weakening the bond and sublayers ①–④; and the in-situ stress of sublayer ⑥ in-
strength between particles. Zangeneh et al.[14] used the creases further and reaching 56 MPa in maximum. In addition,
two-dimensional discrete element method to simulate the hy- the Jiancaogou Formation that underlies the target layer is a
draulic fracture network. In their model, the stratum was di- good stress-shielding layer, resulting in a stress difference of
vided by multiple sets of joints, bounded by deformable rock more than 12 MPa.
blocks. Nagel et al.[15] studied the types and influencing fac-
tors of hydraulic fractures in naturally fractured reservoirs
2. Anisotropy of shale
using the 3D discrete element method. However, they didnot Taking the Fuling shale gas field as the research object,
consider the influence of bedding fracture. Zou et al.[26–27] triaxial rock mechanics tests of shale samples taken from dif-
established a three-dimensional discrete element fracture ferent sublayers and different coring directions of the same
network model to explore the influence of bedding fractures
on the propagation of hydraulic fractures in shale reservoirs at Table 1. Development of beddings in different sublayers.
engineering scale, but they did not take into account the in- Density of bed-
Sub- Bond
fluence of the longitudinal reservoir stress heterogeneity on Core photo ding fracture/
layer strength
the hydraulic fracture propagation pattern. (numberm1)
Taking the Longmaxi shale reservoir in the main area of the
Jiaoshiba anticline in the Fuling shale gas field as the research ⑨ 10 Strong
object, we carried out a series of experimental tests to inves-
tigate the mechanical anisotropy of the shale reservoir in this
study. Then, a three-dimensional complex fracture propaga-
⑧ 32 Strong
tion model was built and solved based on the discrete element
method. Moreover, the propagation model was combined with
the indoor experimental test to investigate the law of hydrau-
lic fracture propagation in shale reservoirs under the influ- ⑥–⑦ 65 Medium
ences of longitudinal stress difference and bedding.
1. Geological background Weak-
⑤ 90
In this study, the main area of the Jiaoshiba anticline is medium
taken as the research object. The organic-rich, gas-bearing
shale intervals of the Upper Ordovician Wufeng Formation
and the Lower Silurian Longmaxi Formation are mainly con- ①–④ 165 Weak
centrated in the bottom of Wufeng Formation-Longmaxi
 1118 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 1. Comprehensive evaluation of the target interval in Well Jiaoye A.

well were carried out with a dynamic rock triaxial test system bedding fractures, the fracturing patterns of core samples from
under high temperature and high pressure conditions, and the different bedding directions are different from each other.
results are shown in Fig. 2. It can be seen that elastic modulus Generally, longitudinal tensile fractures occur more frequently
Eh measured in the direction parallel to the bedding (0°) is in core samples with parallel bedding, while conjugate shear
generally higher than Ev measured in the direction perpen- failures with tensile and shear fractures are predominant in
dicular to the bedding (90°). The Eh/Ev value decreases from core samples with vertical bedding.
1.26 in sublayer ① with many bedding fractures to 1.06 in To study the mechanical properties of bedding and their in-
sublayer ⑧ with few bedding fractures (Fig. 2a). When the fluences on the shear strength anisotropy of shale, the shear
loading stress acts vertically on the bedding surface, the com- strengths of standard samples parallel to and perpendicular to
paction effect of the external force will lead to the closure of the bedding direction were tested using a direct shear testing
the bedding microfractures. A larger axial strain generally machine. The internal friction angle and cohesion of the sam-
corresponds to smaller elastic and deformation moduli. The ples were approximately obtained by linear regression (Fig. 3).
more developed the bedding fractures,the more obvious the As shown in the Mohr-Coulomb strength envelope curves
phenomenon will be. As the confining pressure increases, and obtained by the direct shear test, the shear stress required for
the compaction effect of the interlayer fractures enhances, the shear failure of the rock mass in the case that the direction of
elastic moduli measured in different directions tend to be the shear stress is parallel to that of the bedding direction is
more similar (Fig. 2b). In addition, due to the influence of lower than that in the case that the shear stress direction is

 1119 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

splitting and peeling.

The mechanical parameters of shale, such as elastic pa-
rameters and shear strength, are anisotropic due to the exis-
tence of weak planes or bedding. For the shale gas reservoir
with developed bedding, accurately describing the develop-
ment degree of bedding and the anisotropic characteristics
under its influence are the key to reflecting the difference of
hydraulic fracture propagation in shale reservoirs and natu-
rally fractured reservoirs.
3. The three-dimensional fracture propagation
model
The governing equation of the model is mainly composed
of the fracturing fluid flow equation, the rock mass deforma-
tion equation, and the fracture failure criterion. This govern-
ing equation is solved by combining the finite element and
discrete element methods[26–27]. According to the discrete ele-
ment method[29], the solution domain of the stratum model (f)
is discretized into several triangular-prism block elements that
are linked by virtual springs to transfer the interaction force.
The break of spring represents the rock fracture. The joint
elements between all the contact block elements constitute the
connected fracture network for fracturing fluid flow. The dis-
tribution of the fluid pressure in the contact block is calculated
using the finite element method. The pressure is taken as ex-
Fig. 2. Mechanical test results of shale samples in different cor- ternal load acting on the fracture surface (i.e. the contact sur-
ing directions.
face of the blocks), and the deformation of the block and the
stress state of the spring are calculated. The spring break
(fracture propagation) is determined by the maximum tensile
stress criterion and the Mohr-Coulomb criterion. At the same
time, the constitutive equation of transversely isotropic linear
elastic material is used to replace the isotropic constitutive
equation to study the influence of the mechanical anisotropy
of shale rock on the morphology of fracture propagation[26–27].
3.1. Flow equation
The fluid flow in the fracture is regarded as the plane plate
flow of incompressible Newtonian fluid, which satisfies the
continuity equation without considering the effect of grav-
Fig. 3. Mohr-Coulomb strength envelope curves in different
ity[28]:
shear directions.  w3 p
 qx  
 12 x
perpendicular to the bedding direction, under the same normal  (1)
q   w p
3
stress. At the same time, more developed bedding corresponds
 y 12  y
to lower shear strength in the direction parallel to the bedding.
When sheared in the direction parallel to the bedding, the w qx q y
cohesion of sublayer ① with extremely developed bedding is
   ql  0 (2)
t x y
5.24 MPa, while that of sublayer ⑧ with few bedding is 15.18
Due to extremely low permeability of the shale matrix, the
MPa. Besides, affected by the degree of bedding development,
matrix block is regarded impermeable and the filtration of the
the fracture surface is flat and smooth when the sample is
fracturing fluid is ignored, i.e., ql=0. Then, the global mass
sheared in the direction parallel to the bedding. In contrast,
balance equation in the fracture network is as follows:
when the shearing direction is perpendicular to the bedding
direction, composite failure modes such as matrix tensile  f
wds  tQ (3)
splitting and bedding plane shear slip occur, the fracture sur- When N cracks initiate simultaneously, and the total fluid
face is unsmooth and some fragments are generated by local rate entering all the fractures is Q, the fluid flow into each
 1120 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

fracture depends on the width and pressure of the fracture in action force disappears, namely:
the process of its extension. The model assumes the fractures  Fn( n )  0
are completely filled with fluid, and there is no flow at the end  (n) (9)
 Fs  0
of the hydraulic fractures.
When Fs    Fs,max  AS0  Fn  n  tan , shear failure oc-
n

3.2. Deformation equation of the rock mass curs, accompanied by shear displacement between adjacent
blocks (Δun≤0, Δus>0). At this moment, different from the
The linear elastic dynamic equilibrium equation is as fol-
tensile failure, the shear failure is associated with compressive
lows[30]:
 ij , j  bi   ai ,t   ui ,t  0 (4) stress and friction resistance between blocks:
 Fn  Fn
( n 1)
( n)
 k n Δun( n )
The boundary of the model block is fixed, and ui=0; in the  ( n) (10)
 Fs  Fs,re  Fn tan 
(n)
condition that a contact force is applied on the contact surface
of the block, the fluid pressure pi ( pi   ij  n j ) is applied on
3.3. Iterative coupling of equations
the fracture wall when the hydraulic fracture occurs. Since the
stress-strain relationship conforms to the linear elastic consti- Under the action of the fluid pressure in the fracture, after
tutive equation[30], then: the matrix block deforms, the width of the fracture changes.
 ij  Dijkl  kl (5) Meanwhile, the fluid rate in the fracture also changes with the
As a type of laminated sedimentary rock, shale can be re- width, which affects the fluid pressure in the fracture. In other
garded as a type of transversely isotropic material, namely, its words, the fluid pressure in the fracture and the fracture width
elastic characteristics are the same in the bedding planes, but affect each other. The weak coupling method is used to realize
different in the direction perpendicular to the bedding plane[31]. the iterative process of the fracturing fluid flow and solid de-
Five elastic constants, Eh, Ev, υh, υv, and Gv, are used to charac- formation in fracture, while the flow equation and rock mass
terize the linear elastic characteristics of transversely isotropic deformation equation are discretized respectively, and solved
rock. If the shale layer is horizontal, its flexibility coefficient sequentially and iteratively[26–27]. Before the tensile or shear
matrix is as follows: failure occurs in the matrix block, the permeability of the ini-
 1 h v  tial joint unit is equal to that of the matrix block. The perme-
 E   0 0 0  ability of the matrix block is equivalent to the initial width of
 h Eh Ev  the fracture, that is, w0=(12K0)1/2, and the original fracture
 h 1 v 
  0 0 0  width is designed for the initial flow of the fluid[33]. Then, in
 Eh Eh Ev  the current time step, it is necessary to select appropriate ex-
  v 1  perimental solutions pm and wm to solve the pressure pm+1 in
 v  0 0 0 
 Ev Ev Ev  the next time step, then wm+1 is calculated with equation (11).
A  (6)
2 1   h  Generally, 0< ≤0.5 is taken, and the time step is small enough
 0 0 0 0 0 
  to facilitate the iterative convergence of the equation[34].
Eh
     wm 
 0 1  pm  12  f
0 0 0 0  
  wm   tm
 Gv    (11)
   p  (1   ) p   p
 0 1   m 1 m m  12
 0 0 0 0
 Gv 
3.4. Model validation
where the shear modulus Gv is[32]:
Eh E v The propagation patterns of symmetrical vertical fracture
Gv 
Eh 1   v   Ev and horizontal radial fracture were simulated by the model,
and compared with the analytical solutions of the classical
There are normal and tangential springs between adjacent
fracturing models. The main input parameters were E=35 GPa,
blocks, and the springs can be broken by tensile or shear slip,
ν=0.2, =5 MPaꞏs, and Q=5 m3/min. The results are shown in
i.e. tensile failure or shear failure. Whether a spring is broken
Fig. 4. In the classical PKN model and the radial model, the
is judged according to the maximum tensile stress criterion
energy consumption during fracture propagation is assumed to
and the Mohr-Coulomb criterion. When −Fn(n)<Fn,max (the
(n) be mainly used for fluid flow in the fracture, and the influ-
tensile stress is negative) or Fs <Fs,max , there is no failure
ences of mechanical fracture parameters, such as fracture
between adjacent contact points. The normal stress and shear
toughness and tensile strength of the rock, are not considered.
stress in step n are solved with equations (7) and (8).
The tensile strength is set at zero in the numerical model. The
Fn( n )  Fn( n 1)  k n Δun( n ) (7)
results of the numerical model are basically consistent with
Fs( n )  Fs( n 1)  ks Δus( n ) (8) those of the PKN and radial models, with a difference of 3.6%
(n)
When Fn Fn,max=AT0, tensile failure occurs, the spring and 4.9% respectively, which verifies the reliability of the
breaks, adjacent blocks are separated (Δun>0), and the inter- numerical model.
 1121 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 4. Comparison of results of the numerical model and clas-

sical models.

4. Simulation and analysis

Based on the reservoir parameters of the Jiaoshiba anticline
in the Fuling shale gas field, a layered formation model was
established. Considering that the formation is of horizontal Fig. 5. The minimum horizontal principal stress network model.
occurrence, the boundary of the triangular-prism grid element weak bedding planes with the interval of 2–10 m before frac-
is along the preset trace line of the interlayer interface (weak turing. At the same time, in order to study the influence of
bedding plane). In the early stage of simulation, it was found natural fractures on the morphology of the fracture growth,
that the fractures were not affected by the upper sublayers discrete high angle natural fractures 20 m long and 10 m high
⑦–⑨ in the simulation of sublayers ①–⑨. Out of the con- were set up in random distribution in the model. Given the
sideration of computational efficiency of the model, the in- fact that high-angle natural fractures in the reservoir are of
fluences of sublayers ⑦–⑨ are not considered in the follow- small scales and not connected with each other, their linear
ing modeling. The model was 200 m wide, 500 m long and 50 density was set at 0.05/m. In reference to results of the labo-
m thick in X, Y, and Z directions, respectively. According to ratory rock mechanics test, the mechanical parameters of the
the results of laboratory experiments and in-situ stress inter- shale matrix, natural fracture and weak bedding plane such as
pretation from well logging, the 9 sublayers were divided into tensile strengths, cohesion and internal friction angles were
3 intervals with different stresses vertically (Table 2). Among determined (Table 3).
them, the total thickness of sublayers ①–④ at the bottom was
Table 3. Main input parameters of fracturing simulation.
28 m; the thickness of sublayer ⑤ in the middle part was 10
m, and the stress difference with sublayer ⑥ was 3 MPa; the Category Parameter value
thickness of sublayer ⑥ in the upper part was 12 m. The hor- Density 2 600 kg/m3
Elastic modulus Eh=35 GPa and Ev=30 GPa
izontal well passed through the middle of sublayers ①–④,
Poisson's ratio υh=0.200 and υv=0.175
and the grid model is shown in Fig. 5. To simulate 3-cluster Shale
Permeability 0.001×103 μm2
fracturing in a single stage, it was assumed that only one hy- matrix
Tensile strength 10 MPa
draulic fracture was produced in each perforation cluster. The
Shear strength 30 MPa
basic fracturing simulation parameters were 3 clusters at the
Internal friction angle 40°
spacing of 25 m in one stage, displacement of 14 m3/min, and
Permeability 0.01×103 μm2
fracturing fluid viscosity of 2.5 MPaꞏs.
Bedding Tensile strength 110 MPa
In the study area, bedding fractures are developed in gen- fracture Shear strength 5–30 MPa
eral, especially in the sublayers ①–④ at the bottom of the Internal friction angle 5°–40°
Wufeng-Longmaxi Formations, bedding fractures reach the Permeability (0.01–10.00)×103 μm2
density of over one hundred per meter. During the simulation Natural Tensile strength 0–5 MPa
of the fracture propagation, it is necessary to use an equivalent fracture Shear strength 0–20 MPa
method to approximate the bedding fractures to "explicit" Friction angle 0°–35°
Maximum horizontal
Table 2. Parameters of the in-situ stress profile. 56 MPa
principal stress
Sublayer Z/m σH/MPa σh/MPa σv/MPa p0/MPa Stress Minimum horizontal
49–56 MPa
⑥ 38–50 56 56 58 40 state principal stress
⑤ 28–38 56 53 58 40 Vertical stress 58 MPa
①–④ 0–28 56 49 58 40 Pore pressure 40 MPa

 1122 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

4.1. Results of simulation not considering the effect of 4.2. Simulation considering the effect of bedding
bedding
4.2.1. Influence of the density of weak bedding plane
To illustrate the influence of the weak bedding plane on the
According to the physical modeling of fracturing in labora-
distribution of artificial fractures in the heterogeneous interval,
tory, a "fence-shaped" complex fracture network interwoven
the simulated morphology of fracture propagation in the ho-
by bedding and hydraulic fractures is generated in the shale
mogeneous reservoir interval of a single stage of horizontal
after fracturing[16], and the density of the weak bedding plane
well (without considering bedding and natural fractures) is
would significantly affect the morphology of the propagation
taken as a comparison (Fig. 6a). The simulation with the
of artificial fractures. Under the basic fracturing operation
above fracturing parameters showed three single-wing main
parameters (NC=3, Q=14 m3/min, μ=2.5 mPaꞏs) and medium
fractures 450–480 m long were formed after fracturing in the
strength (TBP=4 MPa, SBP=15 MPa, and φBP=25°) and 2 m of
homogeneous interval. Under the influence of the stress
spacing between weak bedding planes, the fracture propaga-
shadow, the fractures in the outer clusters deflected to the
tion morphologies at different injection times are shown in Fig.
outside; the fractures in the middle clusters extended along the
fracture height direction more than the outer cluster fractures 7. In the early stage of expansion (5 min into injection in the
under the compression of the outer fractures. The fracture near simulation), the hydraulic fracture extended in the sublayers
the wellbore passed through the interface of sublayers ⑤ and ①–④, accompanied with bedding opening (Fig. 7a); after 15
⑥, and was up to 50 m high. With the increase of the distance min of injection, due to the influence of the interlayer stress
to the wellbore, the height of the artificial fracture decreased difference, the fracture height propagation was obstructed and
gradually. Restricted by the stress difference of the interlayer, stabilized at 28 m. The hydraulic fracture extend along the
it can be seen that the fracture height at the fracture tip was direction of the fracture length and the bedding opened con-
limited to the top of sublayer ④, and couldn't enter sublayer tinuously (Fig. 7b); after injection for 45 min, the hydraulic
⑤, with the height of 28 m. In the case considering the influ- fractures in the middle cluster broke through sublayer ④ and
ence of natural fractures (without considering bedding), with expanded to sublayer ⑤. Finally, they were cut off by bed-
the opening of natural fractures, the complexity of the hydrau- ding fractures, and the fracture height was stable at 38 m
lic fractures increased, and the main fractures decreased in (Fig. 7c).
length to about 330 m. Due to the existence of local Fig. 8 shows the numerical simulation results of the artifi-
high-angle cross-layer natural fractures, the fracture height cial fracture propagation at different weak plane densities. As
was very large at local positions (Fig. 6b). In general, the a large-area, continuous weak plane, beddings in shale can
near-wellbore fractures passed through two stress-shielding increase the fracture density during the hydraulic fracturing,
layers in the simulation not considering the influence of bed- and thus improve the adequacy and effect of the reservoir
ding. stimulation. The larger the bedding density, the higher the

Fig. 6. Simulated morphology of fracture propagation in the single stage of horizontal well without considering the bedding effect (with
simulated injection of 15 min).

 1123 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 7. Fracture propagation morphologies at different simulation time during the single stage fracturing in the horizontal well in the
case considering bedding effect (at weak plane spacing of 2 m).

Fig. 8. Simulated fracture growth morphologies at different bedding densities (45 min into injection).

hydraulic fracture density in unit volume will be. However, tionships between the length and height of artificial fractures
the opening of too many beddings will severely restrict the and the stimulation time are shown in Fig. 9. When dBP was 8
expansion of hydraulic fractures in the direction of length and m (weakly developed), the fracture length and height ex-
height, consequently, the stimulated reservoir volume will panded rapidly, and the fracture near the well went through
greatly reduce. Under different bedding densities, the rela- the interface of sublayers ⑤ and ⑥ at 15 min. Then, the overall

 1124 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 10. Fracture propagation dynamics under different beddin

Fig. 9. Fracture propagation dynamics under different bedding g strength conditions.
densities.
high (Fig. 11a). In contrast, when the mechanical strength of
artificial fracture height was stabilized at 40 m, and the final the bedding is moderate to high, the fractures height expanded
simulated fracture length was 400 m. When dBP was 2 m, the smoothly. When the hydraulic fractures passed through sub-
fracture length was less than 300 m, the hydraulic fractures layer ④, they stopped expanding in height for a period of time
ended at the interface of sublayers ⑤ and ⑥, and the fracture due to the stress difference between sublayer ④ and ⑤. At
near the wellbore reached a maximum height of 37 m. In gen- the later stage of fracturing simulation (40 min), the middle
eral, the fractures near the well mostly ended at the interface cluster fractures penetrated through the stress shielding layer
of sublayers ⑤ and ⑥, with a fracture height of about 38 m, under the condition of medium-high to high bedding strength
which is obviously smaller than that in the model without (Fig. 10). As the bedding strength increases, it becomes more
considering the bedding effect. In contrast, the artificial frac- difficult to open the bedding, while the fracturing fluid is
tures far from the wellbore terminated at the interface of more efficient in creating main fractures, thereby the fractures
sublayers ④ and ⑤, with a fracture height of about 28 m. increase in length (Fig. 11d). When the bedding strength is
weak, the bedding is easy to be opened by hydraulic fracture,
4.2.2. Influence of the strength of the weak bedding
then the hydraulic fracture would divert along the bedding,
plane
ending with much smaller fracture height.
The tensile strength and shear strength (including cohesion
4.2.3. Influences of fracturing engineering parameters
and internal friction angle) of beddings are used to describe
whether the bedding is easy to open or not. To study the in- To investigate the correlation between fracturing parame-
fluence of the bedding strength on the fracture propagation ters and fracture propagation morphology, under constant
morphology, four groups of experiment schemes were set up, geological parameters, medium bedding strength (TBP=4 MPa,
in which different values of TBP, SBP, and φBP were input to SBP=15 MPa, and φBP=25°), and bedding spacing of 2 m,
simulate weak bedding planes with different strengths. The simulations at different cluster numbers, displacements and
strength grading is shown in Table 4. fracturing fluid viscosities were conducted, and the results are
The bedding spacing was 2 m. The dynamic parameters of shown in Fig. 12.
the fracture propagation and the fracture morphology at the The same injection parameters, the less the clusters in one
end of the simulation are shown in Figs. 10 and 11, respec- stage, the better it is to increase the fluid pressure in the frac-
tively. It can be seen that the artificial fractures in the shale ture, thereby promote the expansion of the fracture height. In
reservoir with low-strength weak bedding planes were cut off the simulation under 2 clusters in one stage and the cluster
by the bedding fractures near the well and only propagated spacing of 25 m, two transverse fractures 400 m long were
horizontally along the bedding, and were less than 5 m formed; the near-well artificial fractures extended to the top of
sublayer ⑥, with a height of 48 m, while those far from the
Table 4. Mechanical parameters of the weak bedding plane. wellbore stopped growth at the interface of sublayers ⑤ and
Serial Grade of bedding ⑥, with a height of 38 m (Fig. 12a). As the cluster number
TBP/MPa SBP/MPa φBP/(°)
number strength increases, the fracture height decreases. In the simulation un-
1 High ≥8 ≥25 ≥35 der 3 (Fig. 7c) and 5 clusters (Fig. 12b) per stage, the artificial
2 Medium-high (5, 8) (15, 25) (25, 35) fractures near wellbore stopped growth at the interface of
3 Medium [3, 5] [10, 15] [20, 25] sublayers ⑤ and ⑥ longitudinally, with fracture heights of
4 Low <3 <10 <20 38–40 m, while those far from the wellbore terminated at the
 1125 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 11. Simulation of fracture propagation under different bedding strength conditions (45 min into injection).

interface of sublayers ④ and ⑤, with a height of 28 m and fracture was only 168 m, as a result, the fracture-controlled
length of about 300 m. Compared with the 3 clusters fractur- reserves of the single well will drop significantly, which is
ing, the 5 clusters fracturing was weaker in maintaining the unfavorable for the long-term stable production (Fig. 12e). As
fracture height along the fracture length direction, the frac- the viscosity of the fracturing fluid increased to 2.5 mPaꞏs, the
tures far away from the wellbore dropped in height quickly to fracture height reached sublayer ⑤, and the beddings at dif-
28 m, and all fractures reduced in height. ferent positions were opened to various degrees. As the vis-
Injection displacement and fracturing fluid viscosity are cosity of fracturing fluid was further increased to 25 mPaꞏs,
also important engineering parameters to be considered in the hydraulic fractures propagated along the direction of the
shale reservoir volume fracturing. As the displacement was maximum principal stress (vertical fracture), and the fractures
increased from 12 m3/min to 16 m3/min, the propagation trend broke through to sublayer ⑥ (Fig. 12f); however, increasing
of artificial fractures was obvious, and the stimulation scope the viscosity of the fracturing fluid greatly reduced the open-
of sublayer ⑤ increased (Fig. 12c, 12d). Increasing the vis- ing degree of the bedding fractures and the complexity of the
cosity of the fracturing fluid can significantly reduce the re- fractures.
striction of the bedding and longitudinal interlayer stress dif- The aim of the hydraulic fracturing in the shale reservoir is
ference on the fracture height propagation. When the viscosity to improve the complexity of fractures and increase the vol-
was 1.0 mPaꞏs, the hydraulic fracture fully opened the bed- ume of stimulated reservoir, so as to maximize the productiv-
ding near the injection point, and the artificial fracture ity. According to the simulation results, highly developed
vertically ended at the interface of sublayers ④ and ⑤, with a bedding fractures would increase the complexity of hydraulic
fracture height of 28 m; meanwhile, due to the significant fractures, but also restrict the extension of the fracture height
diversion effect of the bedding fracture, the length of the main and reduce the volume of stimulated reservoir. For the interval
 1126 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Fig. 12. Simulation of fracture propagation at different engineering parameters (45 min into simulated injection).

with developed bedding and easy to open, it is necessary to with undeveloped bedding or deep shale gas reservoirs with
reduce the cluster number, increase the displacement, and large overlying stress), the target is to increase the fracture
increase the amount and proportion of high-viscosity fractur- complexity and enlarge the gas drainage fracture area. In
ing fluid in the prepad fluid stage, to reduce the restriction of this context, it is necessary to increase the cluster number,
the bedding on fracture height expansion; for the interval with increase the displacement, and reduce the viscosity of the
undeveloped bedding and difficult to open (upper gas layer fracturing fluid (using low-viscosity slick water).
 1127 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

pass through the stress shielding layer, reaching a height of

close to 50 m. In contrast, when the bedding effect is consid-
ered, there is a 3–4 MPa interlayer stress difference, and the
fractures are controlled in the stress shielding layer, and less
than 38 m high. The influences of weak bedding plane and
longitudinal stress difference should be considered when pre-
dicting fracture propagation in shale reservoirs, otherwise, the
predicted result of fracture height would be quite different
from the actual height. In fracturing of the interval with
abundant bedding, the cluster number should be reduced, the
displacement and the amount and proportion of high-viscosity
fracturing fluid in the prepad fluid stage should be increased
Fig. 13. Comparison between simulation results and micro-seis- to reduce the restriction of bedding on the fracture height ex-
mic monitoring results (red dots indicate micro-seismic signals). pansion, so as to improve the vertical extension of fractures
and increase the volume of stimulated reservoirs.
5. Comparison with micro-seismic monitoring
results Nomenclature
Taking Well Jiaoye A in the Fuling shale gas field as an
example, the horizontal section of the well mainly goes ai,t—acceleration of gravity, m/s2;
through sublayers ①–③. The well was fractured in a total of A—contact area, m2;
28 sections (16 with 2 clusters and 12 with 3 clusters of per- A—flexibility coefficient, Pa−1;
forations), at a cluster spacing of 17–25 m, average fluid con- bi—force per unit volume, N/m3;
sumption of single stage of 1965 m3 and proppant consump- dBP—weak plane spacing, m;
tion of 50.7 m3. Micro seismic monitoring was carried out in Dijkl—elastic tensor, Pa;
the 28 sections during fracturing, and the number of micro E—elastic modulus in isotropic and homogeneous reservoir, Pa;
seismic events in each stage was 5–135. The monitoring re- Eh, Ev—elasticity moduli parallel to and perpendicular to the bed-
sults show that the fractures are 229–483 m long, the fracture ding direction, Pa;
network is 34–204 m wide, and the fractures are 16–49 m f—right vector, m/s;
high (52% of the fractures are 30–40 m high). The three-di- Fn—normal force, N;
mensional numerical model of fracture propagation was used Fn, max—the normal force required for tension break of spring
to simulate each stage, and the results were compared with the between block elements, Pa;
results of micro-seismic monitoring (Fig. 13, the third section Fs—tangential force, N;
of monitoring is taken as an example for display). From the Fs, max—tangential force required for shear break of spring between
simulation, the fractures in the sections were 250–430 m long block elements, Pa;
and 25–50 m high, and fracture network was 76–167 m wide, Fs, re—residual shear resistance after spring break, N;
which are generally consistent with the results of micro-seis- Gv—shear modulus, Pa;
mic monitoring. GR—gamma ray, API;
i, j, k, l—tensor indexes, dimensionless;
6. Conclusions
kn and ks—normal and tangential spring stiffness, N/m;
Based on the discrete element method, a fracture propaga- K0—matrix permeability, m2;
tion model of shale reservoir with developed bedding which m, m+1—current iteration and next iteration;
considers the influence of the weak bedding plane and the n−1, n—current time step and next time step;
longitudinal stress difference is established. The actual pa- nj—unit vector normal to the fracture surface, dimensionless;
rameters of the horizontal shale gas well in the Longmaxi N—number of fractures;
Formation of the main area of the Jiaoshiba anticline in the Nc—number of clusters;
Fuling shale gas field were used for numerical simulation p—pressure, Pa;
analysis. It is found that high-density beddings can increase pi—fluid pressure exerted on the fracture wall, Pa;
the fracture complexity of the stimulated volume, but inhibit pm—experimental solution of pressure from the current iteration
the fracture height and fracture length obviously. When the step, Pa;
low-strength bedding is opened, artificial fractures change to pm+1/2—experimental solution of pressure from the current iteration
expand horizontally, resulting in low bottom hole pressure and step to the next iteration step, Pa;
limited fracture height expansion. When the bedding effect is p0—pore pressure, Pa;
not considered, and there is a 7 MPa interlayer stress differ- ql—filtration rate, m/s;
ence, some local hydraulic fractures in the middle cluster can qx, qy—velocity in x and y directions, m2/s;

 1128 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Q—total injection volume, m3/s; chanics, 2014, 131(2): 269–281.

RLLD—deep lateral resistivity, Ωꞏm; [6] SHIN D H. Factors controlling the simultaneous propagation
RLLS—shallow lateral resistivity, Ωꞏm; of multiple competing fractures in a horizontal well. USA:
s—fracture area, m2; Prentice Hall, 2014.
SBP—shear strength of the bedding fracture, MPa; [7] OLSON J E. Multi-fracture propagation modeling: Applica-
S0—matrix shear strength, Pa; tions to hydraulic fracturing in shales and tight gas sands.
t—time, s; ARMA 08327, 2008.
TBP—tensile strength of the bedding fracture, MPa; [8] XU Yun, CHEN Ming, WU Qi, et al. Stress interference cal-
T0—tensile strength of the matrix, Pa; culation model and its application in volume stimulation of
w—dynamic crack width, m; horizontal wells. Petroleum Exploration and Development,
w0—initial fracture width, m; 2016, 43(5): 780–786.
wm—experimental solution of fracture width of current iteration [9] ZHOU Tong, CHEN Ming, ZHANG Shicheng, et al. Simula-
step, m; tion of fracture propagation and optimization of ball-sealer
x, y—two directions of the coordinate system, m; in-stage diversion under the effect of heterogeneous stress
ui—displacement, m; field in a horizontal well. Natural Gas Industry, 2020, 40(3):
ui,t—velocity, m/s; 82–91.
3
α—damping per unit volume, (kgꞏs)/m ; [10] CHEN Ming, ZHANG Shicheng, XU Yun, et al. A numerical
β—iterative correction coefficient; method for simulating planar 3D multi-fracture propagation in
γ—experimental solution coefficient, m/(Paꞏs); multi-stage fracturing of horizontal wells. Petroleum Explora-
εkl—strain tensor, dimensionless; tion and Development, 2020, 47(1): 163–174.
μ—fluid viscosity, Paꞏs; [11] DAHI T A, OLSON J E. Numerical modeling of multi-
υ—Poisson's ratio, dimensionless; stranded hydraulic fracture propagation: Accounting for the
υh, υv—the Poisson's ratios in the directions of parallel to bedding interaction between induced and natural fractures. SPE Jour-
and vertical to bedding, dimensionless; nal, 2011, 16(3): 575–581.
3
ρ—rock density, kg/m ; [12] GORDELIY E, PEIRCE A. Implicit level set schemes for
σh—minimum horizontal principal stress, Pa; modeling hydraulic fractures using the XFEM. Computer
σH—maximum horizontal principal stress, Pa; Methods in Applied Mechanics & Engineering, 2013, 266:
σij—Cauchy tensor, Pa; 125–143.
3
σij, j—Cauchy Tensor derivative, N/m ; [13] ZHAO X, YOUNG R. Numerical simulation of seismicity
σv—vertical stress, Pa; induced by hydraulic fracturing in naturally fractured reser-
φ—internal friction angle, (°); voirs. SPE 124690, 2009.
φBP—internal friction angle of bedding fracture, (°); [14] ZANGENEH N, EBERHARDT E, BUSTIN R M. Investiga-
Ωf—solution domain; tion of the influence of natural fractures and in situ stress on
Δtm—size of an adaptive time step in Step m, s; hydraulic fracture propagation using a distinct-element ap-
Δwm—fracture width change of current iteration step, m; proach. Canadian Geotechnical Journal, 2015, 52(7): 926–946.
Δun, Δus—normal and tangential relative displacement between [15] NAGEL N B, SANCHEZ M A, ZHANG F, et al. Coupled
adjacent nodes, m. numerical evaluations of the geomechanical interactions be-
tween a hydraulic fracture stimulation and a natural fracture
References system in shale formations. Rock Mechanics and Rock Engi-
neering, 2013, 46(3): 581–609.
[1] XU W, THIERCELIN M, GANGULY U, et al. Wiremesh: A [16] ZOU Y, ZHANG S, MA X, et al. Numerical investigation of
novel shale fracture simulator. SPE 132218, 2010. hydraulic fracture network propagation in naturally fractured
[2] MEYER B R, BAZAN L W. A discrete fracture network shale formations. Journal of Structural Geology, 2016, 84: 1–13.
model for hydraulically-induced fractures: Theory, parametric [17] MIEHE C, HOFACKER M, SCHÄNZEL L M, et al. Phase
and case studies. SPE 140514, 2011. field modeling of fracture in multi-physics problems. Part II.
[3] NAGEL N, IVAN G, MARISELA S N, et al. Simulating hy- Coupled brittle-to-ductile failure criteria and crack propaga-
draulic fracturing in real fractured rock: Overcoming the limits tion in thermo-elastic-plastic solids. Computer Methods in
of pseudo 3D models. SPE 140480, 2011. Applied Mechanics and Engineering, 2015, 294(1): 486–522.
[4] WENG X, KRESSE O, COHEN C E, et al. Modeling of hy- [18] LIU Guowei, LI Qingbin, ZUO Zheng. Implementation of a
draulic fracture network propagation in a naturally fractured staggered algorithm for a phase field model in ABAQUS.
formation. SPE 140253, 2011. Chinese Journal of Rock Mechanics and Engineering, 2016,
[5] BAO J Q, FATHI E, AMERI S. A coupled finite element 35(5): 1019–1030.
method for the numerical simulation of hydraulic fracturing [19] LIU Zhanli, ZHUANG Zhuo, WANG Tao, et al. The key me-
with a condensation technique. Engineering Fracture Me- chanical problems on hydraulic fracture in shale. Chinese

 1129 
 ZHOU Tong et al. / Petroleum Exploration and Development, 2020, 47(5): 1117–1130

Journal of Solid Mechanics, 2016, 37(1): 39–49. [27] ZOU Y, MA X, ZHOU T, et al. Hydraulic fracture growth in a
[20] ZHANG Shicheng, GUO Tiankui, ZHOU Tong, et al. Fracture layered formation based on fracturing experiments and dis-
propagation mechanism experiment of hydraulic fracturing in crete element modeling. Rock Mechanics and Rock Engineer-
natural shale. ActaPetroleiSinica, 2014, 35(3): 496–503. ing, 2017, 50: 2381–2395.
[21] XU Dan, HU Ruilin, GAO Wei, et al. Effects of laminated [28] MA Yongsheng, CAI Xunyu, ZHAO Peirong. China’s shale gas
structure on hydraulic fracture propagation in shale. Petroleum exploration and development: Understanding and practice. Pe-
Exploration and Development, 2015, 42(4): 523–528. troleum Exploration and Development, 2018, 45(4): 561–574.
[22] HENG Shuai, YANG Chunhe, ZENG Yijin, et al. Experimen- [29] MUNJIZA A. The combined finite-discrete element method.
tal study on hydraulic fracture geometry of shale. Chinese Chichester: John Wiley and Sons, 2004.
Journal of Geotechnical Engineering, 2014, 36(7): 1243–1251. [30] JAEGER J C, COOK N G W, ZIMMERMAN R. Fundamental
[23] ZOU Y, ZHANG S, ZHOU T, et al. Experimental investiga- of Rock Mechanics. 4th ed. Oxford: Blackwell Publishing
tion into hydraulic fracture network propagation in gas shales Ltd., 2009.
using CT scanning technology. Rock Mechanics and Rock [31] HENG S, YANG C, ZHANGB, et al. Experimental research
Engineering, 2016, 49(1): 33–45. on anisotropic properties of shale. Rock and Soil Mechanics,
[24] CHEN Z, JEFFREY R G, ZHANG X, et al. Finite-element 2015, 36(3): 606–616.
simulation of a hydraulic fracture interacting with a natural [32] LEKHNITSKII S G. Theory of elasticity of an anisotropic
fracture. SPE 176970, 2017. body. Russian: Mir Publishers, 1981.
[25] CUNDALL P A, STRACK O D L. A discrete numerical model [33] ZIMMERMAN R W, BODVARSSON GS. Hydraulic conduc-
for granular assemblies. Geotechnique, 1979, 29(1): 47–65. tivity of rock fractures. Transport in Porous Media, 1996,
[26] ZOU Y, MA X, ZHANG S, et al. Numerical investigation into 23(1): 1–30.
the influence of bedding plane on hydraulic fracture network [34] ADACHI J, SIEBRITS E, PEIRCE A, et al. Computer simula-
propagation in shale formations. Rock Mechanics and Rock tion of hydraulic fractures. International Journal of Rock Me-
Engineering, 2016, 49: 3597–3614. chanics & Mining Sciences, 2007, 44(5): 739–757.

 1130 