## Abstract

A numerical modeling framework is described that is able to calculate the coupled processes of fluid flow, geomechanics, and rock failure for application to general engineering problems related to reservoir stimulation, including hydraulic fracturing and shear stimulation. The numerical formulation employs the use of an embedded fracture modeling approach, which provides several advantages over more traditional methods in terms of computational complexity and efficiency. Specifically, the embedded fracture modeling strategy avoids the usual requirement that the discretization of the fracture domain conforms to the discretization of the rock volume surrounding the fractures.

#### Authors

##### Jack H. Norbeck^{1} · Mark W. McClure^{3} · Jonathan W. Lo^{2} · Roland N. Horne^{1}

^{1}Department of Energy Resources Engineering, Stanford University, Stanford, CA 94305, USA. ^{2}Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA. ^{3}McClure Geomechanics LLC, Palo Alto, CA 94301, USA.

Received: 26 March 2015 / Accepted: 20 October 2015 / Published online: 13 November 2015

© The Author(s) 2015.

As fluid is exchanged between the two domains, conservation of mass is guaranteed through a coupling term that appears as a simple source term in the governing mass balance equations. In this manner, as new tensile fractures nucleate and propagate subject to mechanical effects, numerical complexities associated with the introduction of new fracture control volumes are largely negated. In addition, the ability to discretize the fractures and surrounding rock volume independently provides the freedom to choose an acceptable level of discretization for each domain separately. Three numerical examples were performed to demonstrate the utility of the embedded fracture model for application to problems involving fluid flow, mechanical deformation, and rock failure. The results of the numerical examples confirm that the embedded fracture model was able to capture accurately the complex and nonlinear evolution of reservoir permeability as new fractures propagate through the reservoir and as fractures fail in shear.

## 1 Introduction

In many reservoir engineering applications, it is imperative to incorporate a realistic description of the geologic structure of the reservoir into conceptual models and numerical models in order to establish appropriate interpretations of reservoir behavior. Several examples include hydraulic fracture treatment design, interpretation of microseismic monitoring data, and development of reservoir management strategies related to induced seismicity. In each of these cases, the interaction between fluid flow and the geomechanical response of fractured and faulted rock will have a direct influence on the reservoir behavior, and therefore also on the engineering decisions that must be made.

In geologic settings where fractures and faults are expected to have first-order impacts in terms of flow behavior, it is important to recognize that the reservoir systems are mechanically active. During development and operation of a resource, local-scale and reservoir-scale permeability and storativity can evolve as fractures deform and fail in shear, or as intact rock fails in tension. The local state of stress throughout the reservoir controls the manner in which the permeability and storativity changes manifest. It is necessary to make use of numerical modeling to investigate these types of reservoir processes for practical applications, but many traditional reservoir models neglect geomechanical processes or are based upon a set of limiting assumptions that obviate the influence of significant physical mechanisms.

In this paper, we introduce an efficient numerical modeling framework that is able to model the coupled physical processes of fluid flow and mechanical deformation of fractures, faults, and surrounding matrix rock. The framework is able to incorporate an explicit representation of the geologic structure of the reservoir by using an embedded fracture modeling (EFM) strategy [21]. The EFM approach is an extension of more traditional finite-volume-based discrete fracture models (DFM), suggesting relatively straightforward integration with industry-standard reservoir simulators. A fracture mechanics-based approach to mechanical modeling allows for accurate calculation of the complex stress distributions that arise near fracture tips, so fracture propagation problems are approached in a rigorous manner. Detailed models of friction evolution along fracture and fault surfaces are included in order to model shear failure and seismicity. To accommodate different types of rock, the model is flexible enough to incorporate a range of constitutive relationships necessary to describe permeability and storativity evolution of fracture networks due to changes in effective stress and shear failure.

The integration of the EFM approach into a geomechanical and fracture propagation model is the principal achievement of the present work. The embedded fracture approach is a numerical method that provides the critical translation necessary to attack problems that would other- wise be intractable from a computational standpoint. The key element of the EFM formulation is the treatment of the fracture system and the surrounding matrix rock volume as two separate computational domains [19, 21]. This allows for the two domains to be discretized completely independently, negating the cumbersome requirement of a matrix discretization that must conform to the fracture discretization associated with traditional discrete fracture models.

In the EFM discretization, conservation of mass is enforced as fluid is exchanged between the two domains through the application of physics-derived coupling terms that appear as source terms in the continuity equations. In this manner, tensile fractures that nucleate and propagate subject to geomechanical considerations can be incorporated into the numerical model during a simulation with negligible amount of computational overhead.

Previous authors have applied EFM to investigate geome chanical effects in fractured reservoirs, but this work has been limited in scope to simplified models that embody the geomechanics into empirical relationships [16, 26]. In this work, the reservoir model introduced by McClure [23] and McClure and Horne [24] was extended to incorporate the EFM strategy in order to combine the effects of matrix-fracture mass exchange and a rigorous treatment of geomechanics under a unified framework. The fluid flow and geomechanical calculations are performed in a coupled sequential-implicit manner. Fracture propagation is permit- ted, and is based upon evaluating the mode-I stress intensity factor near fracture tips. Shear failure of preexisting fractures is permitted subject to a modified Mohr-Coulomb criterion [13]. The mechanical interaction between fractures as they deform is included in the model.

Reservoir-scale permeability evolution emerges as a result of deformation of individual fracture and fault planes, shear failure on preexisting fractures and faults, and propagation of tensile fractures. Limitations of the model arise from the following important assumptions: the mechanical properties of the rock are homogeneous and constant, elastic deformation is quasistatic, poroelastic and thermoelastic deformations are neglected, and the domain is two-dimensional (2-D) so all fractures in the model have the same out-of-plane dimension.

The primary purpose of this work was to verify the accuracy of the EFM approach for applications in which geomechanical effects have a first-order impact on reservoir behavior. The embedded fracture-based model developed in this work was compared to a traditional DFM for three different numerical examples of practical interest. When appropriate, the numerical results were also compared against analytical or semianalytical solutions. In addition, the accuracy of two extremely efficient approximate models, a one-dimensional (1-D) leakoff model and a zero leakoff model, were also compared in order to help quantify their range of practical applicability for reservoir stimulation modeling.

The remainder of this paper is organized as follows. In Section 2, the numerical formulation for the EFM framework is presented. The traditional DFM, 1-D leakoff approximation model, and zero leakoff approximation model that were used for comparison are also discussed. In Section 3, the major components in the geomechanical module are described. In Section 4, we demonstrate that the models achieved good agreement with an analytical solution for a problem involving injection into an infinite-conductivity fracture. Next, we modeled a shear stimulation treatment of a relatively complex fracture network, and the results are presented in Section 5. Significant nonlinear effects due to shear slip-enhanced fracture permeability

made this a challenging numerical problem. Finally, the utility of the EFM framework was highlighted by modeling mode-I fracture propagation of a vertical, two-wing fracture. The results of the fracture propagation study are shown in Section 6. In Section 7, several practical applications for the present model are illustrated, and concluding remarks are discussed.

## 2 Fluid flow module

The reservoir model introduced by McClure [23] and McClure and Horne [24] was developed originally to model reservoir stimulation treatments in low-permeability settings, such as hydraulic fracturing in shale gas reservoirs or shear stimulation in geothermal reservoirs. The model assumed that intrinsic matrix permeability in these settings is usually low enough to justify neglecting mass transfer between the fracture systems and surrounding matrix rock. In the present work, this model was extended to include the effects of matrix-fracture mass transfer. An EFM strategy was adopted in order to overcome several severe numerical and practical limitations of more traditional DFM approaches for application to reservoir stimulation problems.

In this section, the numerical formulation for the embedded fracture model is described in detail. In the sections that follow, we describe studies where we compared the EFM approach with a DFM approach that had been implemented previously in the McClure [23] model. Here, the details of the DFM model are described for the reader’s reference. In addition, the EFM and DFM models were compared with two relatively computationally efficient approximate models, namely the 1-D leakoff approximation model and the zero leakoff approximation model. The details of these models are described briefly.

### 2.1 Embedded fracture model description

In traditional reservoir simulation, fractured reservoirs are commonly modeled using the double porosity model [17, 42]. This model is applicable if the fracture orientations and lengths are relatively randomly distributed and the fracture system is connected extensively. More importantly, application of the double porosity model typically assumes that the properties and geometry of the fracture network remain relatively constant.

In order to honor more realistic representations of fractured reservoir geology, discrete fracture approaches have been developed. For example, Karimi-Fard et al. [14] presented a DFM in which the geometry of the fractures and faults are captured by discretizing them explicitly in lower-dimensional space, and creating a matrix discretization that conforms to the fractures. In general, DFM approaches are useful in settings where production is dominated by flow through fractures (e.g., formations with low matrix permeability) or if fractures tend to have preferred orientations. However, traditional DFM techniques are subject to several well-documented drawbacks. Most notable is that a matrix discretization that conforms to the fractures inevitably results in a large number of “small” matrix control volumes in areas where there are many fractures or where fractures intersect at low angles. In some cases, the geologic structure of the reservoir must be sacrificed for numerical convenience.

Moreover, DFM techniques are not well suited for fracture propagation problems in which the fracture networks are growing over time. Previous work has been done in the area of developing models that apply adaptive grid refinement as fractures propagate [12, 33]. This requires a significant level of computational overhead, and the numerical results have been observed to be grid-dependent. Alter-natively, it is possible to define planes where hydraulic fractures may potentially propagate in advance of a simulation, and then prediscretize the system around these potentially forming planes. Naturally, this approach will require an unnecessarily high number of additional degrees of freedom, and, perhaps worse, involves making implicit assumptions about the mechanics of fracture propagation.

#### 2.1.1 Embedded fracture model overview

In the present work, the use of traditional DFM techniques was avoided, and instead the EFM approach was adopted. In the EFM approach, the fracture and matrix domains are treated as separate computational domains. The two systems are discretized completely independently (i.e., a conforming mesh is not required; see Fig. 1), and mass conservation is strictly enforced through physics-derived coupling terms. In fact, EFM is conceptually very similar to dual porosity or dual permeability models, but is able to maintain a more realistic representation of complex geologic features. As will be demonstrated in the numerical examples that follow, the ability to define realistic representations of the geologic structure of a reservoir and the use of a nonconforming grid are the critical features that make the EFM approach an attractive modeling strategy to perform rigorous geomechanical analyses.

**Fig. 1** Schematic of the embedded fracture discretization strategy. The solid blue lines are natural fractures, and the dashed red lines are hydraulic fractures. The circles represent the centers of fracture con- trol volumes, and the diamonds represent the centers of matrix control volumes. The matrix control volumes that will pick up EFM coupling terms are shaded gray

The EFM approach was introduced originally by Lee et al. [19], and later expanded upon by Li and Lee [21]. Karvounis [15] developed a heat and mass transfer geothermal model based on EFM, and demonstrated that EFM can obtain a suitable degree of accuracy with improved computational performance compared to traditional models. Hajibeygi et al. [9] incorporated EFM into an iterative multiscale finite volume scheme. Several studies have compared DFM to EFM for multiphase flow problems and demonstrated that EFM was able to capture a high degree of accuracy at a reduced computational expense [25, 32]. Ding et al. [4] drew upon EFM fundamentals and calculated the matrix-fracture transmissibility numerically to be able to capture pressure transients in the near-fracture region more accurately.

Recently, EFM has also been used in geomechanics applications. Moinfar et al. [26] incorporated a treatment for calculating fracture permeability evolution due to changes in effective stress within an EFM framework, but did not include a formal treatment for geomechanics. Karvounis et al. [16] extended their model to include a proxy geomechanical model based on changes in pore pressure to investigate injection-induced seismicity. Norbeck et al. [27] and Nor beck and Horne [28] introduced an EFM-based model that integrated a rigorous treatment of fluid flow, fractured reservoir mechanics, and fracture propagation. Norbeck and Horne [29] performed a study of porothermoelastic effects on injection-induced seismicity using a rate-and-state earth-quake model.

It should be noted that in some of the works cited previously, the concept of EFM was applied in a context related to upscaling techniques [9, 21]. In those applications, it was assumed that fractures that existed at a relatively small scale could be homogenized in order to obtain effective “damaged matrix rock” properties, and the geometries of the larger fracture systems expected to contribute to flow at a reservoir scale were maintained explicitly. For the purposes of this paper, it is sufficient to recognize this distinction purely at the conceptual level. In the remainder of this paper, it is assumed that any reference to matrix permeability may imply an effective upscaled permeability, and any fracture domain is representative of a scale of practical engineering interest.

### 2.1.2 Numerical formulation of the embedded fracture model

The key insight introduced by Li and Lee [21], that the fracture and matrix domains can be discretized independently, is leveraged by expressing the mass conservation equations for the matrix and fracture domains separately. For a porous medium saturated with single-phase quations can be written, for flow in the matrix domain, as:-

and, for flow in the fracture domain, as:

Here, p^{m} is fluid pressure in the matrix domain, p^{f} is fluid pressure in the fracture domain, ρ is fluid density, λ is inverse of fluid viscosity, **k**_{m} is the diagonal matrix permeability tensor, **k**_{f }is the diagonal fracture permeability tensor, φ is matrix porosity, e is fracture hydraulic aperture, E is fracture void aperture, m^{wm} is a normalized mass source term related to wells in the matrix domain, and is a normalized mass source term related to wells in the fracture domain. In addition to the usual terms related to flux, wells, and storage, the terms Ψ^{fm} and Ψ^{mf} are introduced to account for mass transfer between the two domains. To ensure continuity upon integration over the respective control volumes, these mass transfer terms take the following form [9]:

where the parameter ϒ is a transmissibility called the fracture index and is analogous to the Peaceman well index [31]. The normalization parameters in Eqs. 3 and 4 are V, the bulk volume of the matrix control volume, and A, the surface area of the fracture control volume. Comparing terms in Eqs. 1–4, the units indicate that the fracture mass balance equations are expressed on a lower-dimensional manifold.

Similar to the treatment of wells in traditional reservoir simulators, the fracture index serves to capture subgrid behavior of the pressure gradient near fractures. In this work, the derivation provided by Li and Lee [21] was followed to calculate the fracture index. The assumptions in the derivation are as follows: (a) flow in the vicinity of the fractures is linear, (b) the fractures fully penetrate the matrix control volume in the out-of-plane direction, and (c) the matrix and fracture pressures represent average pressures over their respective control volumes.

The rate of mass transfer from a fracture control volume into a matrix control volume is defined as:

This term has units of mass per time. Assuming that flow is linear (i.e., 1-D) in the local region near the fracture, the mass exchange rate can be described alternatively by integrating the Darcy flux over the surface area of the fracture:

where A^{f }is the total fracture surface area (i.e., both faces of the fracture), k^{∗} is an effective fracture-normal permeability, n is the unit normal vector to the fracture face, and the pressure gradient term is:

Here, L represents the average normal distance from the fracture surface in the matrix control volume, which can be calculated numerically for complex fracture and matrix control volume geometries [9]. The effective permeability, k^{∗}, represents the ability for fluid to flow in the direction perpendicular to the fracture. In reservoir settings with low matrix permeability and highly conductive fractures, k^{∗} will likely be dominated by the matrix rock permeability. In geologic settings where relatively thick faults with sealing properties exist, the fracture (fault) permeability can affect k^{∗} significantly. In general, k^{∗} can be considered a harmonic mean of the matrix and fracture permeabilities:

where k_{⊥}^{m} and k_{⊥}^{f }are the matrix and fracture permeabilities, respectively, resolved in the direction normal to the fracture, and W is the physical half-width of the fracture. In Eq. 8, it was assumed that L >> W. Equating the right-hand side expressions in Eqs. 5 and 6 allows for the determination of the fracture index:

(9) where I is a grid dependent property with units of length that can be calculated as:

With the matrix-fracture mass flux terms that appear in Eqs. 1 and 2 now fully defined, the utility of the EFM approach for problems that involve fracture propagation is revealed. The coupling between the fracture and matrix domains has been reduced to a collection of simple source terms. Numerical complexities associated with conforming mesh approaches that would tend to make fracture propagation problems become intractable for problems with many fractures are now avoided with EFM. Standard finite volume schemes can be used to discretize Eqs. 1 and 2 [1, 14].

### 2.2 Discrete fracture model description

The DFM was implemented using the finite volume method and a conforming mesh of the rock volume around the fractures. The volume around the fractures was discretized with triangular control volumes, aided by the program Triangle [35]. Triangle is an algorithm designed to create Delaunay triangularizations of 2-D regions. The finite volume method was implemented according to the method described by Karimi-Fard et al. [14].

An important problem is that Delaunay triangularization does not guarantee uniform or smoothly varying line- segment length along domain edges, which in this case are the fracture elements. This can create problems in the boundary element mechanical calculations described in Section 3, which are inaccurate unless fracture element length is uniform or gradually varying. Therefore, Triangle could not be used to generate a true Delaunay triangularization.

To guarantee uniform fracture element size, the fractures were discretized first by imposing a constant element length (with some minor and unavoidable deviation from constant length at fracture tips and intersections). Next, a uniform grid of matrix nodes was superimposed over the fracture network. Third, nodes were identified that were in close proximity to fracture elements, and they were removed. Finally, the list of fracture elements and matrix nodes was provided to Triangle, which was used to produce a constrained Delaunay triangularization.

The triangularization was “constrained” in the sense that the algorithm was required to use only the fracture elements and matrix nodes provided and was not permitted to subdivide the fracture elements. Because of the constraint, the mesh was not guaranteed to be truly Delaunay, which degraded the quality of the mesh and the accuracy of the calculations of flow in the matrix. Despite this problem, the approach was used because it was more important to avoid inaccuracy in the mechanical calculations due to unevenly sized fracture elements than inaccuracy in matrix flow calculations due to high aspect ratio triangles. It should be noted that this entire issue is avoided with the EFM approach, which does not require a conforming mesh between the matrix and fracture elements. Figure 2 is an example of the type of conforming mesh used in the calculations.

**Fig. 2** Example of a triangular discretization used for the conforming mesh discrete fracture model. The red lines are the natural fractures, and the blue lines are the edges of the triangular matrix control volumes

### 2.3 One-dimensional leakoff approximation model description

For a well connected to a highly conductive fracture within an infinite reservoir, flow near the fracture has been shown to be linear at early times [8]. This suggests that a useful approximation to model the leakoff behavior near fractures is to assume 1-D flow away from the fracture. In this work, the semianalytical method of Vinsome and Westerveld [41] was used to develop an approximate model that is relatively efficient computationally compared to the EFM and DFM approaches. The purpose of the simplified model is to avoid numerical discretization of the volume of rock surrounding the fractures, while still accounting for fluid exchange between the two domains.

The model treats fluid leakoff at each fracture control volume using a sink term that is independent from all other fracture control volumes. The key advantage of the Vinsome and Westerveld [41] method is that it gives a highly accurate and efficient solution to the diffusivity equation in 1-D, even for arbitrarily varying pressure in the fracture. In contrast, the Carter leakoff model assumes constant pressure in the fracture [11], a simplifying assumption that seriously reduces model generality.

The Vinsome and Westerveld [41] method was created originally as a model of heat loss due to conduction into cap rock. However, the equation for heat conduction is identical to the equation for single-phase fluid flow in a porous media with constant pore and fluid compressibilities, matrix permeability, and fluid viscosity. Therefore, the method can be adapted easily by changing the variables in the original equations of Vinsome and Westerveld [41] to their equivalents for flow in porous media, assuming that the aforementioned variables are considered as constants.

The assumptions of 1-D leakoff and no interference between fractures are justified if the fracture spacing is sufficiently large relative to the penetration distance of the pressure signal and if the injection duration is short enough to preclude a change in the flow regime (e.g., towards late-time radial flow). The penetration distance, l, can be estimated as [2]:

where the total compressibility, c, is the sum of the pore and fluid compressibilities, μ is fluid viscosity, and φ_{0} is the initial porosity of the matrix rock. If fractures are in pressure communication with other nearby fractures, then the 1-D leakoff approximation tends to overestimate the amount of leakoff. Relatively high leakoff suppresses pressure within the fracture domain, which discourages shear failure and fracture propagation.

### 2.4 Zero leakoff approximation model description

In geologic settings where the intrinsic permeability of the matrix rock is extremely low, a useful approximation is that the matrix rock is impermeable. In this case, fluid flow can occur only within a network of connected fractures, and no fluid is able to leakoff into the surrounding rock. Under this assumption, there is no flow in the matrix rock, and the volume surrounding the fractures does not require discretization. The improvement in computational efficiency that can be achieved by avoiding discretization of the matrix rock volume can be tremendous.

This approximation represents a lower bound on the amount of leakoff that would occur during a stimulation treatment. Neglecting leakoff tends to promote elevated pressure in the fracture domain, which encourages both shear stimulation and fracture propagation. Even in scenarios in which the implicit assumptions of this approximation are not strictly valid, the model can be used to provide informative constraints on reservoir behavior.

## 3 Geomechanics module

The present numerical model was built upon the framework introduced by McClure [23] for coupling fluid flow in fractures and fracture deformation. The reader is referred to McClure [23] and McClure and Horne [24], where thorough descriptions of the assumptions and the numerical formulation for the geomechanics module are explained. Here, we provide a brief overview of the major components of the model.

### 3.1 Fracture deformation

We are interested principally in modeling systems that contain a large number of fractures and faults, and therefore the displacement fields that arise as a result of fracture deformation are expected to be discontinuous. A boundary element method called the displacement discontinuity (DD) method is capable of calculating the complex displacement fields due to the mechanical interaction between fractures as they deform [3]. In addition, the DD method has been demonstrated to calculate accurate stress distributions in the vicinity of fracture tips, which is necessary to model fracture propagation [30].

The DD model assumes a 2-D faulted and fractured porous domain, saturated with a single-phase fluid. The mechanical properties of the intact matrix rock are homogeneous. Deformations occur quasistatically, and linear elasticity applies. Using the approach described by Shou and Crouch [37], a linear system of equations can be developed that relates the mode-I (opening) and mode-II (sliding) fracture deformations to the traction boundary conditions along the fractures. The primary variables that must be solved for are the mode-I and mode-II displacement discontinuities, Δe (or equivalently ΔE) and Δδ, respectively.

The traction boundary conditions are functions of the effective normal stress distributions along the fractures. The effective normal stress, σ¯ , at a specific location along a fracture surface is defined as:

where σ^{R} is the normal component of the remote tectonic stress, and Φ is the normal component of the elastic stress transfer that occurs as nearby fracture elements deform. Compressive stress has been taken as positive in this sign convention. The shear stress, τ , acting on a fracture surface is:

where τ^{R} is the shear component of the remote tectonic stress, Θ is the shear component of the elastic stress transfer that occurs as nearby fracture elements deform, η is a radiation damping coefficient, and v is sliding velocity. The radiation damping term is used to approximate inertial effects when sliding occurs very rapidly [34]. In the model, Φ and Θ are calculated using the DD method. A Mohr-Coulomb shear failure criterion is used for the fracture sliding calculations, whereby the frictional strength of a fracture, τ_{s}, is defined as [13]:

Here, f is the coefficient of friction, and S is the fracture cohesion.

Fractures can be classified into several groups depending on their local state of stress. A fracture element is classified as closed if it is bearing compression such that its walls are in physical contact. The walls of closed fractures are rough surfaces, and so these fractures have the ability to transmit and store fluid. Elastic stress transfer due to mode-I deformations of closed fractures is assumed to be negligible. A fracture element is open if the effective stress drops to zero such that the fracture walls become out of contact. Open fractures deform subject to fracture mechanics considerations. If the shear stress acting on the fracture is less than its frictional resistance to slip (i.e., τ < τ_{s}), then the fracture is classified as stuck. Mode-II deformations and the associated elastic stress transfer is assumed to be negligible for stuck fractures. If the shear stress acting on the fracture becomes equal to its frictional resistance to slip (i.e., τ = τ_{s}), the failure criterion has been met and the fracture is classified as sliding.

The boundary conditions that must be enforced for the mechanical deformation problem depend on these classifications. For closed fractures that are sliding, the following boundary condition is enforced:

For open fractures, the following boundary conditions are enforced:

For open fractures, the following boundary conditions are enforced:

Equations 1, 2, and 15–17 can be recast in residual form to develop a coupled system of nonlinear equations. The primary variables involved in these equations are p^{m}, p^{f}, Δe, and Δδ. In the model, a sequential implicit strategy is used to solve the fully coupled system of equations [18]. The equations involving fluid pressure and opening displacement are grouped together and solved simultaneously (1, 2, and 16), and the residual equations involving shear displacement are solved separately (15 and 17). The sequential strategy iterates between these two groups until all residual equations have converged to within a prescribed tolerance.