**2. Model description**

#### **2.1. Background: Synthetic rock mass approach**

Over past years, the SRM has been developed [2] as a more realistic representation of me‐ chanical behavior of the fractured rock mass compared to conventional numerical models. The SRM consists of two components: (1) the bonded particle model (BPM) of deformation and fracturing of intact rock, and (2) the smooth joint model (SJM) of mechanical behavior of discontinuities.

The BPM, originally implemented in *PFC*, is created when the contacts between the particles (disks in 2D and spheres in 3D) are assigned certain bond strength (both in tension and shear). It was found that BPM quite well approximates mechanical behavior of the brittle rocks [1]. The elastic properties of the contacts (i.e., contact shear and normal stiffness) can be calibrated to match the desired elastic properties (e.g., Young's modulus and Poisson's ratio) of the assembly of the particles. Similarly, the tensile and shear contact strengths can be adjusted to match the macroscopic strengths under different loading conditions (e.g., direct tension, unconfined and confined compression).

**1. Introduction**

820 Effective and Sustainable Hydraulic Fracturing

naturally fractured reservoirs (NFRs).

isotropic and homogeneous media.

A new generation tool that uses the bonded particle model (BPM) [1] and the synthetic rock mass (SRM) concept [2] has been developed to model hydraulic fracture (HF) propagation in

Most rock mass formations, and shale gas reservoirs in particular, consist of a large volume of fractured rock in which propagation of an HF involves both fracturing of intact rock and opening or slip of pre-existing discontinuities (joints). The pre-existing joints can significantly affect the HF trajectory, the pressure required to propagate the fracture, but also the leak-off from the fracture into the surrounding formation. None of these effects can be simulated using conventional hydraulic fracturing simulation methods, based on assumptions of continuous,

To address this challenge, a numerical approach called SRM method [2] has been developed recently based on the distinct element method. SRM method usually is realized as a bondedparticle assembly representing brittle rock containing multiple joints, each one consisting of a planar array of bonds that obey a special model, namely the smooth joint model (SJM). The SJM allows slip and separation at particle contacts, while respecting the given joint orientation rather than local contact orientations. Overall fracture of a synthetic rock mass depends on

Previous SRM models have used the general-purpose codes *PFC2D* and *PFC3D* [3,4], which employ assemblies of circular/spherical particles bonded together. Much greater efficiency can be realized if a "lattice," consisting of point masses (nodes) connected by springs, replaces the balls and contacts (respectively) of *PFC3D*. The lattice model still allows fracture through the breakage of springs along with joint slip, using a modified version of the SJM. The new 3D program, *HF Simulator* described in this paper, is based on such a lattice representation of brittle rock. *HF Simulator* overcomes all main limitations of the conventional methods for simulation of hydraulic fracturing in jointed rock masses and is computationally more efficient

The formulation of the code is described in this paper. The examples of code verification and

Over past years, the SRM has been developed [2] as a more realistic representation of me‐ chanical behavior of the fractured rock mass compared to conventional numerical models. The SRM consists of two components: (1) the bonded particle model (BPM) of deformation and fracturing of intact rock, and (2) the smooth joint model (SJM) of mechanical behavior of

both fracture of intact material (bond breaks), as well as yield of joint segments.

than *PFC*-based implementations of the SRM method.

**2.1. Background: Synthetic rock mass approach**

application are also presented.

**2. Model description**

discontinuities.

In the BPM, the contact behavior is perfectly brittle. Breakage of the bond, a function of the forces in the contact and the bond strength, corresponds to formation of a microcrack. An example of unconfined compression test conducted using *PFC2D* is illustrated in Figure 1, which shows recorded axial stress-strain response and the model configuration with generated microcracks. The shear microcracks are black; the tensile microcracks are red. Shown is the state when the sample is loaded beyond its peak strength. The stress-strain curve exhibits characteristics typical of brittle rock response. For the load levels less than ~80% of the peak strength, the stress-strain response is linearly elastic, with the slope of the line equal to the Young's modulus. Some microcracks, randomly distributed within the sample, start develop‐ ing at the load levels greater than ~40% of the peak strength. Significant non-linearity develops as the load exceeds 80% of the peak strength. In this phase, the microcracks begin to coalesce, forming fractures on the scale of the sample. After the peak strength is reached, the material starts to soften (i.e., to lose the strength). At this stage, as shown in Figure 1, the failure mechanism and the "shear bands" are well developed. It is interesting that in the unconfined compression test, the majority of cracks are tensile (red lines in Figure 1). The "shear bands" on the scale of the sample are formed by coalescence of a large number of tensile microcracks.

**Figure 1.** Example of unconfined compressive test using bonded particle model (BPM).

In order to model a typical rock mass in the BPM, it is also necessary to represent pre-existing joints (discontinuities). A straightforward approach is to simply break or weaken the bonds (in the contacts between the particles) intersected by the pre-existing joints. The created discontinuity will have roughness with the amplitude and wavelength related to the resolu‐ tion, or the particle size of the BPM. The mechanical behavior of discontinuities is very much affected by their roughness. The problem is that the selected particle size (or resolution) typically is not related to actual roughness of the pre-existing joints. The SJM overcomes this limitation. The contacts in the BPM model are oriented in the direction of the line connecting the centers of the particles involved in the contact. The SJM contacts are oriented perpendicular to the fracture plane irrespective of the relative position of the particles. Consequently, the particles can slide relative to each other in the plane of the fracture as if it is perfectly smooth.

cal relations to estimate the rock mass properties and to account for the size effect considering

Three-Dimensional Numerical Model of Hydraulic Fracturing in Fractured Rock Masses

http://dx.doi.org/10.5772/56313

823

The new code, *HF Simulator*, is based on implementation of the SRM in the lattice, which is a simplified, but also a computationally more efficient version of particle flow code (*PFC*). Despite simplifications, the lattice approach represents all physics important for simulation of

The lattice is a quasi-random array of nodes (with given masses) in 3D connected by springs. It is formulated in small strain. The lattice nodes are connected by two springs, one represent‐ ing the normal and the other shear contact stiffness. The springs represent elasticity of the rock mass. In *HF Simulator*, the calibration factors for spring stiffness are built-in and the user may specify typical macroscopic elastic properties as it is done for other conventional numerical models. The tensile and shear strengths of the springs control the macroscopic strength of the lattice. As for elastic constants, calibration factors are built-in for the strength parameters.

The model simulation is carried out by solving an equation of motions (three translations and three rotations) for all nodes in the model using an explicit numerical method. The following

is the central difference equation for the translational degrees of freedom:

+D -D +D +D

& &

( /2) ( /2) ( ) ( ) ( ) ( /2)

*u uu t*

/ *tt tt t ii i tt t tt i ii*

> NN N S S SS

*F F uk t F F uk t* ¬+ D ¬+ D &

*i ii*

formed. In other words, if *<sup>F</sup>* <sup>N</sup> <sup>&</sup>gt; *<sup>F</sup>* Nmax, then *<sup>F</sup>* <sup>N</sup> =0 , *Fi*

*N*

where "N" denotes "normal," "S" denotes shear, *k* is spring stiffness and *F* is the spring force. If the force exceeds the calibrated spring strength, the spring breaks and the microcrack is

Fluid-flow model and hydro-mechanical coupling are essential parts of *HF Simulator*, as a code for simulation of hydraulic fracturing. The fluid flow occurs through the network of pipes that connect fluid elements, located at the centers of either broken springs or springs that represent pre-existing joints (i.e., springs intersected by the surfaces of pre-existing joints). (The code

*u u F t m*

= +D S =+ D

(*t*) are the velocity and position (respectively) of component *i* (*i* =1, 3) at time

is the sum of all force-components *i,* acting on the node of mass *m,* with time step *Δ t*. The relative displacements of the nodes are used to calculate the force change in the springs:

& (1)

& (2)

<sup>S</sup> =0, and a "fracture flag" is set.

the size of the samples tested in the laboratory and the scale of interest in the model.

hydraulic fracturing.

**2.2. Lattice**

where *u*˙*<sup>i</sup>*

*t*, *ΣFi*

(*t*) and *ui*

**2.3. Fluid flow**

The SRM and its components are shown in Figure 2. The BPM represents the intact rock, its deformation and damage. The pre-existing joints are represented explicitly, using the SJM. They can be treated deterministically, by specifying each discontinuity by its position and orientation as mapped in the field. However, typically, for practical reasons, it is not possible to treat the DFN deterministically. Instead, fracturing in the rock mass is characterized statistically. The synthetic DFNs that are statistically equivalent (i.e., fracture spacing, orien‐ tation and size) to fracturing of the rock mass are generated and imported into the SRM using SJM (Figure 2). Very often a reasonable compromise is to represent few dominant structures (faults) with their deterministic position and orientation and the rest of the fracturing in the rock mass (smaller structures) using a synthetic DFN.

**Figure 2.** Synthetic rock mass (SRM).

One of the advantages of the SRM is that the components, the intact rock and the joints, can be mechanically characterized by standard laboratory tests. The mechanical response of the rock mass and the size effect are the model results, functions of the model size, DFN characteris‐ tics and mechanical properties of the components. Thus, it is not necessary to rely on empiri‐ cal relations to estimate the rock mass properties and to account for the size effect considering the size of the samples tested in the laboratory and the scale of interest in the model.

The new code, *HF Simulator*, is based on implementation of the SRM in the lattice, which is a simplified, but also a computationally more efficient version of particle flow code (*PFC*). Despite simplifications, the lattice approach represents all physics important for simulation of hydraulic fracturing.

#### **2.2. Lattice**

In order to model a typical rock mass in the BPM, it is also necessary to represent pre-existing joints (discontinuities). A straightforward approach is to simply break or weaken the bonds (in the contacts between the particles) intersected by the pre-existing joints. The created discontinuity will have roughness with the amplitude and wavelength related to the resolu‐ tion, or the particle size of the BPM. The mechanical behavior of discontinuities is very much affected by their roughness. The problem is that the selected particle size (or resolution) typically is not related to actual roughness of the pre-existing joints. The SJM overcomes this limitation. The contacts in the BPM model are oriented in the direction of the line connecting the centers of the particles involved in the contact. The SJM contacts are oriented perpendicular to the fracture plane irrespective of the relative position of the particles. Consequently, the particles can slide relative to each other in the plane of the fracture as if it is perfectly smooth.

The SRM and its components are shown in Figure 2. The BPM represents the intact rock, its deformation and damage. The pre-existing joints are represented explicitly, using the SJM. They can be treated deterministically, by specifying each discontinuity by its position and orientation as mapped in the field. However, typically, for practical reasons, it is not possible to treat the DFN deterministically. Instead, fracturing in the rock mass is characterized statistically. The synthetic DFNs that are statistically equivalent (i.e., fracture spacing, orien‐ tation and size) to fracturing of the rock mass are generated and imported into the SRM using SJM (Figure 2). Very often a reasonable compromise is to represent few dominant structures (faults) with their deterministic position and orientation and the rest of the fracturing in the

One of the advantages of the SRM is that the components, the intact rock and the joints, can be mechanically characterized by standard laboratory tests. The mechanical response of the rock mass and the size effect are the model results, functions of the model size, DFN characteris‐ tics and mechanical properties of the components. Thus, it is not necessary to rely on empiri‐

rock mass (smaller structures) using a synthetic DFN.

822 Effective and Sustainable Hydraulic Fracturing

**Figure 2.** Synthetic rock mass (SRM).

The lattice is a quasi-random array of nodes (with given masses) in 3D connected by springs. It is formulated in small strain. The lattice nodes are connected by two springs, one represent‐ ing the normal and the other shear contact stiffness. The springs represent elasticity of the rock mass. In *HF Simulator*, the calibration factors for spring stiffness are built-in and the user may specify typical macroscopic elastic properties as it is done for other conventional numerical models. The tensile and shear strengths of the springs control the macroscopic strength of the lattice. As for elastic constants, calibration factors are built-in for the strength parameters.

The model simulation is carried out by solving an equation of motions (three translations and three rotations) for all nodes in the model using an explicit numerical method. The following is the central difference equation for the translational degrees of freedom:

$$\begin{cases} \dot{u}\_{l}^{(t \circ \Lambda t/2)} = \dot{u}\_{l}^{(t \circ \Lambda t/2)} + \Sigma F\_{l}^{(t)} \,\Delta t \,/\, m \\ u\_{l}^{(t \circ \Lambda t)} = u\_{l}^{(t)} + \dot{u}\_{l}^{(t \circ \Lambda t/2)} \Delta t \end{cases} \tag{1}$$

where *u*˙*<sup>i</sup>* (*t*) and *ui* (*t*) are the velocity and position (respectively) of component *i* (*i* =1, 3) at time *t*, *ΣFi* is the sum of all force-components *i,* acting on the node of mass *m,* with time step *Δ t*. The relative displacements of the nodes are used to calculate the force change in the springs:

$$\begin{cases} F^N \leftarrow F^N + \dot{\boldsymbol{u}}^N \boldsymbol{k}^N \boldsymbol{\Delta t} \\ F\_i^s \leftarrow F\_i^s + \dot{\boldsymbol{u}}\_i^S \boldsymbol{k}^S \boldsymbol{\Delta t} \end{cases} \tag{2}$$

where "N" denotes "normal," "S" denotes shear, *k* is spring stiffness and *F* is the spring force. If the force exceeds the calibrated spring strength, the spring breaks and the microcrack is formed. In other words, if *<sup>F</sup>* <sup>N</sup> <sup>&</sup>gt; *<sup>F</sup>* Nmax, then *<sup>F</sup>* <sup>N</sup> =0 , *Fi* <sup>S</sup> =0, and a "fracture flag" is set.

#### **2.3. Fluid flow**

Fluid-flow model and hydro-mechanical coupling are essential parts of *HF Simulator*, as a code for simulation of hydraulic fracturing. The fluid flow occurs through the network of pipes that connect fluid elements, located at the centers of either broken springs or springs that represent pre-existing joints (i.e., springs intersected by the surfaces of pre-existing joints). (The code also can simulate the porous medium flow through unfractured blocks as a way to represent the leakoff. This capability is not discussed further in this paper.) The flow pipe network is dynamic and automatically updated by connecting newly formed microcracks to the existing flow network. The model uses the lubrication equation to approximate the flow within a fracture as a function of aperture. The flow rate along a pipe, from fluid node "A" to node "B," is calculated based on the following relation:

$$q = \rho k\_r \frac{a^3}{12\,\mu} \left[ p^A - p^B + \rho\_w g \left( z^A - z^B \right) \right] \tag{3}$$

off, the response is viscosity-dominated, which corresponds to the "M-asymptote" identified

In the simulated example, fluid is injected at a constant rate into a penny-shaped crack of low

zero. Thus, the test conditions approximate those of the analytical solution for the no-lag case (i.e., no fluid pressure tension cut-off) provided by [5]. The injection rate is 0.01 m3

dynamic viscosity is 0.001 Pa×s. The mechanical properties of the rock are characterized by Young's modulus of 7×1010Pa and Poisson's ratio of 0.22. Figure 3 provides a visualization of the state of the model at 10 s of elapsed time. Note that pressures are negative in the outer

**Figure 3.** View of pressure (Pa) field (icons, colored according to magnitude) and cross-section of displacement (m)

Figure 4 shows the aperture profiles at three times during the simulation — averaged numer‐ ical results (for 30 radial distances), together with asymptotic solutions (derived from the equations of [5]). Figure 5 shows the pressure profile at 10 s, together with the asymptotic solution. Note that there is a lack of match at small and large radial distances: at small distances, the numerical source is a finite volume, rather than a point source (which is assumed in the exact solution); at large distances, the finite initial aperture allows seepage (compared to zero

seepage in the exact solution, which assumes zero initial aperture).

m). The crack has zero normal strength, and the in-situ stresses are also

Three-Dimensional Numerical Model of Hydraulic Fracturing in Fractured Rock Masses

/s; the

825

http://dx.doi.org/10.5772/56313

by [5]. This condition is used for verification of *HF Simulator*.

initial aperture (10−<sup>5</sup>

annulus of the flow disk.

field (vectors, colored according to magnitude)

where *a* is hydraulic aperture, *μ* is viscosity of the fluid, *p <sup>A</sup>* and *p <sup>B</sup>* are fluid pressures at nodes "A" and "B", respectively, *z <sup>A</sup>* and *z <sup>B</sup>* are elevations of nodes "A" and "B", respectively, and *ρw* is fluid density. The relative permeability, *kr*, is a function of saturation, *s*:

$$k\_r = s^2 \text{ ( $3 - 2s$ )}\tag{4}$$

Clearly, when the pipe is saturated, *s* =1 and the relative permeability is 1. The dimensionless number *β* is a calibration parameter, a function of resolution, used to match conductivity of a pipe network to the conductivity of a joint represented by parallel plates with aperture *a*. The calibrated relation between *β* and the resolution is built into the code.

#### **2.4. Hydro-mechanical coupling**

In *HF Simulator*, the mechanical and flow models are fully coupled.


A new coupling scheme, in which the relaxation parameter is proportional to *KRa* / *R*, where *KR* is rock bulk modulus and *R* is the lattice resolution, is implemented in *HF Simulator*, allowing larger explicit time steps and faster simulation times compared to conventional methods that use fluid bulk modulus as a relaxation parameter.
