Conventional methods for simulation of hydraulic fracturing are based on assumptions of continuous, isotropic and homogeneous media. These assumptions are not valid for most rock mass formations, particularly shale gas reservoirs, as these typically consist of a large volume of naturally fractured rock in which propagation of a hydraulic fracture (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 and also the leak-off from the fracture into the surrounding formation. None of these effects can be simulated using conventional methods.

*HF Simulator* is a new three-dimensional numerical code that can simulate propagation of hydraulic fracture in naturally fractured reservoirs, accounting for the interaction between the hydraulic fracture and pre-existing joints. In *HF Simulator*, fracture propagation occurs as a combination of intact-rock failure in tension, and slip and opening of joints. The code uses a lattice representation of brittle rock consisting of point masses (nodes) connected by springs. The pre-existing joints are derived from a user-specified discrete fracture network (DFN).

*HF Simulator* can model fluid injection or production from one or multiple boreholes each with one or multiple clusters. Non-steady, hydro-mechanically coupled fluid flow and pressure within the network of joint segments and the rock matrix are considered.

An outline of the code hydro-mechanical formulation is presented and examples are provided to illustrate the code capabilities.

**Keywords:** Numerical model, naturally fractured, rock mass

## 1. Introduction

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 naturally fractured reservoirs (NFRs).

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, isotropic and homogeneous media.

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 bonded-particle 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 both fracture of intact material (bond breaks), as well as yield of joint segments.

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 than *PFC*-based implementations of the SRM method.

The formulation of the code is described in this paper. The examples of code verification and application are also presented.

## 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 mechanical 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).

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 developing 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.

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 resolution, 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, orientation 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.

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 characteristics and mechanical properties of the components. Thus, it is not necessary to rely on empirical 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

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 representing 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:

where *t*, *i,* acting on the node of mass *m,* with time step

where “N” denotes “normal,” “S” denotes shear,

### 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:

where

Clearly, when the pipe is saturated,

### 2.4. Hydro-mechanical coupling

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

Fracture permeability depends on aperture, or on the deformation of the solid model.

Fluid pressure affects both deformation and the strength of the solid model. The effective stress calculations are carried out.

The deformation of the solid model affects the fluid pressures. In particular, the code can predict changes in fluid pressure under undrained conditions.

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

## 3. Verification test: Penny-shaped crack propagation in medium with zero toughness

The non-steady response of rock to injection of fluid depends on fracture toughness, the viscosity of the fluid and the rate of leak-off. In the case of zero fracture toughness and no leak-off, the response is viscosity-dominated, which corresponds to the “M-asymptote” identified by [5]. This condition is used for verification of *HF Simulator*.

In the simulated example, fluid is injected at a constant rate into a penny-shaped crack of low initial aperture (^{3}/s; the dynamic viscosity is 0.001 Pa×s. The mechanical properties of the rock are characterized by Young’s modulus of

Figure 4 shows the aperture profiles at three times during the simulation — averaged numerical 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).

## 4. Example application

Two example problems are discussed in this section. Fracture propagation in a homogeneous (unfractured) and fractured media is analyzed. These two problems involve a horizontal borehole segment with two injection clusters with centers at 4.8 m distance (Figure 6). The model domain is 18 m × 18 m × 18 m, and the lattice resolution was set to 0.5 m. Fluid is injected into the clusters at rate of 0.01 m^{3}/s. The assumed stress state is anisotropic with

The induced microcracks in the homogeneous model after 15 s of injection are shown in Figure 7. The microcracks form two roughly circular (penny-shape) hydraulic fractures. In this example, the fractures are not parallel. There is a slight trend of fractures curving away from each other as a result of stress interaction.

In the second example, the HF propagation is clearly affected by the pre-existing joints, as shown in Figure 8. When the HF intersects the pre-existing joint, the fluid is diverted into the pre-existing joints. (In general case, the HF can cross or be diverted into the pre-existing joint, depending on a number of parameters, including stress state, strength and permeability of the pre-existing joint.) The propagation continues by reinitiation along the edges of pre-existing joints.

## 5. Conclusion

*HF Simulator* is a powerful 3D simulator for hydraulic fracturing in jointed rock mass that allows the main mechanisms (nonlinear mechanical response, fluid flow in joints and coupled fluid-mechanical interaction) to be reproduced. The formulation of *HF Simulator* is based on a quasi-random lattice of nodes and springs.

The springs between the nodes break when their strength (in tension) is exceeded. Breaking of the springs corresponds to the formation of microcracks, and microcracks may link to form macrofractures. The SJM (smooth joint model) is used to represent pre-existing joints in the model. Thus, the SJM allows simulation of sliding of a pre-existing joint in the model, unaffected by the apparent surface roughness resulting from lattice resolution and random arrangement of lattice nodes.

The model is fully coupled hydro-mechanically. There are several ways in which fluid interacts with the rock matrix. First, fluid pressures may induce opening or sliding of the fractures. Second, mechanical deformation of fractures causes changes in joint pressures. Third, the mechanical deformation changes the permeability of the rock mass as the joint apertures change.

The new code is a promising tool for simulation and understanding of complex processes, including propagation of HF and its interaction with DFN, during stimulation of unconventional reservoirs.

## Acknowledgments

The development of the numerical code described in this paper was funded by BP America. The authors would like to thank BP America for their support. Matt Purvance, Jim Hazzard and Maurilio Torres of Itasca Consulting Group, Inc. are thanked for their valuable work on HF Simulator.