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

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


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

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

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 u ˙i (t ) and u i (t ) are the velocity and position (respectively) of component i (i = 1, 3) at time t, ΣF i 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: 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 F N > F Nmax , then F N = 0 , F i S = 0, and a "fracture flag" is set.

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 a is hydraulic aperture, μ is viscosity of the fluid, p A and p B are fluid pressures at nodes "A" and "B", respectively, z A and z B are elevations of nodes "A" and "B", respectively, and ρ w is fluid density.The relative permeability, k r , is a function of saturation, s: 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.

Hydro-mechanical coupling
In HF Simulator, the mechanical and flow models are fully coupled.

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

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

3.
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 K R a / R, where K R 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.

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- Effective and Sustainable Hydraulic Fracturing 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 (10 −5 m).The crack has zero normal strength, and the in-situ stresses are also 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 m 3 /s; the dynamic viscosity is 0.001 Pa×s.The mechanical properties of the rock are characterized by Young's modulus of 7 × 10 10 Pa 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 annulus of the flow disk.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).

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 σ xx = 1 MPa, σ yy = 12 MPa and σ zz = 10 MPa.The least principal stress is aligned with the horizontal section of the borehole.This stress state favors crack propagation in the direction normal to the horizontal section of the borehole.In order to initiate the fluid calculation, fluid-filled joints have been placed at the center of each cluster; these joints are slightly larger than the cluster size.The initial apertures in these joints have been set to 0.1 mm.Both example problems use this model configuration.The example shown on the left in Figure 6 simulates the response of an unfractured medium to fluid injection.Three discrete joints that interact with the induced fractures are introduced in the example on the right in Figure 6.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.

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.Effective and Sustainable Hydraulic Fracturing 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.

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

Figure 4 .
Figure 4. Aperture profiles for three times

Figure 7 .Figure 8 .
Figure 7. Hydraulic fractures generated in a homogeneous medium (dark blue disks are microcracks)