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.
An outline of the code hydro-mechanical formulation is presented and examples are provided to illustrate the code capabilities.
A new generation tool that uses the bonded particle model (BPM)  and the synthetic rock mass (SRM) concept  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  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
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  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
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
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,
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
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 and are the velocity and position (respectively) of component at time
where “N” denotes “normal,” “S” denotes shear, is spring stiffness and is the spring force. If the force exceeds the calibrated spring strength, the spring breaks and the microcrack is formed. In other words, if , then and a “fracture flag” is set.
2.3. Fluid flow
Fluid-flow model and hydro-mechanical coupling are essential parts of
where is hydraulic aperture, is viscosity of the fluid, and are fluid pressures at nodes “A” and “B”, respectively, and are elevations of nodes “A” and “B”, respectively, and is fluid density. The relative permeability, is a function of saturation, :
Clearly, when the pipe is saturated, 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 . The calibrated relation between and the resolution is built into the code.
2.4. Hydro-mechanical coupling
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 , where is rock bulk modulus and is the lattice resolution, is implemented in
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 . This condition is used for verification of
In the simulated example, fluid is injected at a constant rate into a penny-shaped crack of low initial aperture (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 . The injection rate is 0.01 m3/s; the dynamic viscosity is 0.001 Pa×s. The mechanical properties of the rock are characterized by Young’s modulus of 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 ). 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 m3/s. The assumed stress state is anisotropic with and 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.
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.
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.