# AAPG Bulletin

- Copyright ©2012. The American Association of Petroleum Geologists. All rights reserved.

## Abstract

In this series of studies, we develop a numerical tool for modeling finite deformation of reservoir rocks. We present an attempt to eliminate the main limitations of idealized methods, for example, elastic or kinematic, that cannot account for the complexity of rock deformation. Our approach is to use rock mechanics experimental data and finite element models (Abaqus). To generate realistic simulations, the present numerical rheology incorporates the dominant documented deformation modes of rocks: (1) rock mechanics experimental observations, including finite strength, inelastic strain hardening, strength dependence on confining pressure, strain-induced dilation, pervasive and localized damage, and local tensile or shear failure without macroscopic disintegration; and (2) field observations, including large deformation, distributed damage, complex fracture networks, and multiple zones of failure.

Our analysis starts with an elastic–plastic damage rheology that includes pressure-dependent yield criteria, stiffness degradation, and fracturing via a continuum damage approach, using the Abaqus materials library. We then use experimental results for Berea Sandstone in two configurations, four-point beam and dog-bone triaxial, to refine and calibrate the rheology. We find that damage and fracturing patterns generated in the numerical models match the experimental features well, and based on these observations, we define damage fracturing, the fracturing process by damage propagation in a rock with elastic–plastic damage rheology. In part 2, we apply this rheology to investigate fracture propagation at the tip of a hydrofracture.

## INTRODUCTION

Fracturing in the crust strongly depends on the rheology of host rocks. In addition to the original lithology, the rheology is affected by nonlinear time-dependent diagenesis and deformation processes that lead to nonuniform damage-prone rock properties. Classical models of linear elastic fracture mechanics (Irwin, 1958) and crack-tip shielding models (Barenblatt, 1962) may not be applicable to in-situ growth of natural fractures because they are based on linear elasticity, require idealized flaws (e.g., elliptical cracks), and are based on arbitrary selection of interaction configurations (e.g., en echelon segments).

Our main objective is to analyze the processes of fracture initiation and fracture propagation in rocks with complex and more realistic rheology that can undergo finite strain and inelastic deformation. We explore the spontaneous generation of complex patterns, including fracture segmentation (Delaney et al., 1986), damage bridging (Weinberger et al., 2000), branching (Sagy et al., 2001), and arrest-and-rupture features, without assuming preexisting seeds of geometric complexity. We refer to the fracturing process by damage propagation in a rock with elastic–plastic damage rheology as damage fracturing, which differs from the concepts of linear elastic fracture mechanics and shielded fracture mechanics as discussed in following sections.

This analysis is for nonlinear rheology with progressive deformation that is sensitive to in-situ loading conditions and is capable of brittle failure by fracture propagation. To meet these requirements, we used the finite element (FE) method (Fish and Belytschko, 2007) that allows solving nonlinear transient problems. We used the Abaqus FE software that has a library of material models suitable for rocks and is flexible in implementing geologic boundary conditions. The present models focus on rock behavior under upper crustal conditions with elastic-plastic behavior, progressive damage, and failure in an attempt to realistically simulate rock deformation under in-situ conditions. We develop a numerical model with material parameters that are derived from tensile and compressive stress-strain curves of the Berea Sandstone, and then benchmark simulations of two experimental configurations are performed to calibrate the model rheology.

In related studies (Busetti, 2009; part 2, Busetti et al., 2012), we apply this damage rheology to hydraulic fracturing of a rock layer and find that tendency to fracture and fracture geometry are determined by damage distribution, stress path, and dynamical rupture processes. We also applied the numerical simulations to examine fracture networks formed during hydraulic fracturing and their relations to subsurface field data (borehole image logs, pressure-time logs, and microseismic distribution) (Busetti, 2009, 2010). These studies provide insight on the activation of preexisting fractures in dilation, shear, or closure and their effect on fracture propagation by rock damage. The applicability of the rheology to model fold-scale deformation was demonstrated in a study of intralayer deformation above a ramp structure (Heesakkers et al., 2009).

In this article, we first outline the concept of damage fracturing and then describe damage patterns observed in rock mechanics experiments along with schemes for damage quantification. Next, we discuss the FE damage model based on a continuum damage theory and describe the implementation of this damage rheology for the rock experiments. A product is a constitutive model and failure criteria for Berea Sandstone that incorporates elastic, plastic, and damage parameters and can be implemented numerically to simulate complex deformation.

## ROCK FRACTURING

### Fracturing of an Elastic Solid

A central concept of fracture mechanics is that a large fracture grows from a small initial flaw that amplifies the local stresses, and the large fracture propagates when these local stresses exceed the strength of the solid (Griffith, 1921; Lawn, 1993). The above concept was formulated by linear elastic fracture mechanics (LEFM), in which a mode I model consists of “crack-tip, near-field” solutions for stress and displacement fields around a slitlike plane crack (Irwin, 1958). The crack tip is assumed to be sharp, and the material is idealized as a linear elastic solid with unspecified strength (Lawn, 1993). The stresses at the proximity of a fracture are scaled by the stress intensity factor, *K*, which is a function of the remote stress (e.g., tectonic) and/or fracture internal loads, fracture geometry, and the loading configuration. According to LEFM, mode I propagation occurs when the stress intensity factor *K*_{I} reaches the critical value of the fracture toughness, *K*_{IC}, a material property; subcritical crack growth may occur under lower *K*_{I} (Lawn, 1993). Equivalently, one can consider the energy associated with fracture propagation. When the fracture energy from loading, *G*, reaches the intrinsic limit of the material, *G*_{C}, the fracture extends by length, *d*_{c}, where *E* is the isotropic Young's modulus and v is Poisson's ratio. The term *U* is the mechanical energy of the system and includes the elastic potential within the material and the potential energy from external loading.

### Fracturing of an Elastic Solid with Crack-Tip Yielding

The stress singularity at the crack tip under LEFM (e.g., Kanninen and Popelar, 1985) was replaced by a small zone with a nonlinear cohesive stress function (Barenblatt, 1962). This cohesive zone accounts for the fact that materials have a finite strength governed by interatomic stresses (Barenblatt, 1962) or plastic yield strength (Dugdale, 1960). The existence of a cohesive zone changes the crack profile from parabolic to a cuspate form, where the fracture narrows at the crack tip. According to this model, propagation occurs when the stress intensity factor, *K*, reaches the critical value of the fracture toughness of the cohesive zone, or the intrinsic toughness, given for a straight crack (Lawn, 1993). The J-integral technique (Rice, 1968) effectively solves for the energy balance around the nonlinear cohesive zone as long as the stress-strain characteristics are reversible. If dissipative processes are considered, additional steps are required, for example, limiting plastic deformation by enforcing monotonic loading or by introducing a shielding zone (Lawn, 1993). The shielding concept (Thompson, 1986) accounts for the fact that many materials exhibit toughening at the crack tip that is independent of external loading conditions and results in an effective “blunting” of the crack tip. Shielding mechanisms include a zone of dislocations, microcracks, and phase transformations in front of the crack tip.

The model of a cohesive crack-tip and a shielding zone is relevant to fracture propagation in rocks that exhibit macroscopic propagation via coalescence and linkage of microcracks within a damage front. For example, experimental observations indicate that the physical location of the crack tip may be ambiguous because of the development of a network of microcracks in the crack frontal zone. The process follows several stages: (1) initial state, (2) initiation, (3) small growth and linkage, and finally, (4) large growth and linkage (Takahashi and Abe, 1987), and can generate complex fracture patterns (Figure 1).

### Damage Fracturing of an Inelastic Solid

The present analysis examines fracture propagation in rocks with damage rheology. Using the context of continuum damage mechanics eliminates the need to a priori define multiple zones (elastic, cohesive, and shielding) and allows damage to initiate and propagate according to the intrinsic constitutive rheology of the rock. This rheology is based on the experimental observation that many rocks undergo pervasive damage from the initial stages of loading and fractures form within the zones of damage. The nonlinear damage evolution degrades the rock stiffness primarily by microcrack interaction (e.g., Kachanov et al., 1990; Lyakhovsky et al., 1993; Katz and Reches, 2004). Unlike the insertion of cohesive or shielding zones as discussed above, damage fracture propagation is localized within weakening zones that are spontaneously determined by the material constitutive behavior. Furthermore, rock mechanics experiments indicate that, during progressive deformation, the critically stressed regions are not necessarily restricted to the tip of a distinct initial crack and may occur randomly within the host because of innumerable randomly distributed microflaws. In brittle rocks, yielding is characterized by nonlinear inelasticity associated with stress-induced damage accumulation (e.g., Ashby and Hallam, 1986; Lockner et al., 1992; Reches and Lockner, 1994). A field application of damage fracturing was presented by Weinberger et al. (2000) who compared damage simulations with field observations of fracture patterns and damage distribution around dike segments in sandstones. They used a damage model (Lyakhovsky et al., 1997) that considers a linear elastic matrix populated with penny-shaped cracks (Kachanov, 1992). In this scheme, the energy expression considers the effective moduli, *λ*^{e} (*λ*^{d}, *μ*^{d}, *λ*, strain) and *μ*^{e} (*λ*^{d}, *μ*^{d}, *γ*, strain), for a medium damaged by a mode I crack, where *λ* and *μ* are the first and second Lame parameters, *γ* is an additional modulus, and the superscripts e and d reflect the moduli in their original and damaged states. The crack energy release rate is given by the J-integral: *E* is Young's modulus, *ν* is Poisson's ratio, and *K*_{I} is the mode I fracture toughness. The energy of the damaged rock is given by *ρ* is the microcrack density, *I*_{1} and *I*_{2} are strain invariants, and *β* reflects the state of strain (Lyakhovsky et al., 1997).

Rock fracturing is analyzed here by simulating deformation in solids with damage rheology that may realistically represent the in-situ rock behavior. We refer to this process as damage fracturing to distinguish it from the LEFM and shielded-tip fracture mechanics outlined above. The analysis is based, in part, on previous damage models (below); however, our approach differs in one central point: we focus on geometry and pattern of damage fractures, whereas previous analyses focused on the general rheologic character attributed to damage evolution.

This approach has several advantage points. First, it was demonstrated in many works (e.g., Lockner et al., 1992; Weinberger et al., 1994; Katz and Reches, 2004) that rocks undergo pervasive and localized damage starting under small strain (<1%). Second, field and experimental studies display complex fractures and networks of fractures that cannot be explained by elastic analysis, and as demonstrated here and in part 2 (Busetti et al., 2012), the damage fracturing analysis spontaneously predicts such modes. Third, damage fracturing does not require any special assumption, such as initial perturbations or nonrealistic high stresses. Fourth, damage fracturing does not suffer from the present computational limitations of local element enrichment formulations (i.e., the extended finite element method [XFEM]) such as limits on the number of fractures that can cross an enriched element or the difficulty in synthesizing pressure-dependent yielding with traction-separation definitions for crack propagation.

## DAMAGE RHEOLOGY IN ROCKS

### Rock Inelasticity and Damage

The recoverable elastic deformation stage for brittle rocks corresponds to small strain magnitudes (<1%). However, in most geologic structures, deformation exceeds the elastic limit, leading to the development of permanent heterogeneous strain patterns. The deformation involves many inelastic processes: continuous distortion (e.g., twinning and dislocations), microfracturing and macrofracturing, cataclasis, brecciation, faulting, and formation of shear zones, compaction bands, and tight folds. Consequently, analysis of most geologic structures requires the consideration of finite strain and nonlinear inelasticity. A typical stress-strain curve for rock (Figure 2) can be divided into five main stages interpreted as (1) initial nonlinear stress change associated with crack and pore dilation and closure; (2) elastic stage (linear or nonlinear); (3) nonlinear strain hardening associated with the onset of brittle microcracking and plasticity; (4) continued hardening characterized by progressive crack coalescence in a fracture process zone; and finally, (5) ultimate failure, strain softening, and macroscopic crack propagation (Katz and Reches, 2004). After an initial elastic stage, the rock enters a strain-hardening stage that is dominated by the accumulation of damage, primarily by microcracking (e.g., Ashby and Hallam, 1986; Lockner et al., 1992; Reches and Lockner, 1994). The rock progressively weakens because of evolving cracks and microcracks (Walsh and Brace, 1964; Walsh, 1965), plastic deformation (Handin and Hager, 1957; Mogi, 1973), stress shielding (Thompson, 1986), and dilatancy (Brace et al., 1966; Nur, 1975). Macroscopically, this degraded stiffness is linked to the evolution of stress-induced damage (Lyakhovsky et al., 1997; Katz and Reches, 2004) that leads to local fracturing (e.g., Rice, 1975) and, eventually, to failure. The locations, timing, and amount of damage associated with the strain-hardening stage have been quantified using thin-section maps, scanning electron microscopy, and acoustic emissions (AE) logs. The experiments showed that the total amount of damage increases nonlinearly beginning at the onset of inelasticity (e.g., Lockner et al., 1992; Eberhardt, 1998; Stanchits and Dresen, 2003; Backers et al., 2005; Bobich, 2005; Chen et al., 2006). The damage is nonuniformly distributed, as shown for example by Backers et al. (2005) for sandstone under three-point beam bending and by Reches and Lockner (1994) for triaxial loading of granite.

### Deformation Modulus

The accumulation of damage and the changing shape of the stress-strain curve (stage III, black line in Figure 2) can be macroscopically presented by the deformation modulus, *D*, which is the local slope of an experimental stress-strain curve (Johnson and Page, 1976; Katz and Reches, 2004) (red curve, Figure 2). It provides a phenomenological approximation of the material properties at different stages of deformation. The strain hardening can be described as macroscopic stiffness degradation associated with a reduction in the elastic load-carrying capacity of the rock. Stiffness, deformation, and load-carrying capacity are expressed by modifying the familiar equation (e.g., Hooke's law):*σ* is the Cauchy stress tensor, *E* is the stiffness matrix, and *ε* is the elastic strain tensor. A stress change, *dσ*, corresponds to each strain increment, *dε*, such that

Because the material stiffness changes during deformation, the initial elastic stiffness, *E*, is no longer applicable, and is replaced by a more general tensor, *D*, that includes permanent inelastic deformation attributed to stress-induced damage (Katz and Reches, 2004). The *D* accounts for the sum of the degradation behavior but does not distinguish individual weakening mechanisms. Rearranging (1b) and using differential stress (*σ*_{1} − *σ*_{3}), where *σ*_{1} and *σ*_{3} are the axial and confining stresses, in place of *σ* (for triaxial tests) defines the deformation modulus, *D*:

Before the onset of damage, when *D* = *D*_{max} (Figure 2), *D* is equivalent to the Young's modulus. At various incremental states, *D* = *D*′ < *D*_{max}. The *D*′ and *D*_{max} are the deformed and undeformed states. The ratio of the deformed to undeformed states provides a dimensionless approximation for stiffness degradation:*E*_{0} is the original Young's modulus and E′ is the apparent modulus at the deformed state.

Coalescence and shear-zone development indicate a strong dependence on the strain history; a damage “memory” that affects future events exists. A nondegrading or static damage scheme will give nonphysical results. Principally, the mechanisms for damage and the magnitude of damage interaction (i.e., inelastic energy dissipation) are transient phenomena. The solution is to introduce an evolving stress function that depends on the state of stress and strain history; clearly, this is a task for an advanced numerical method as presented below.

## PRESENT ANALYSIS

### Approach

The central goal of this study is to analyze the process of damage fracturing of rocks. We model fractures that initiate and propagate spontaneously to form patterns of natural fractures in a rock body with damage rheology. The selected rheology is valid for finite deformation and includes plastic yielding, as well as failure in tension and compression. In this scheme, pressure-dependent yield strength controls the onset of microcracking damage and internal dissipation processes (plasticity) that are coupled to control finite deformation. Plasticity is based on the Mohr-Coulomb and Drucker-Prager failure models as implemented in the FE code Abaqus that we use in this study. The formulation of the rheology is presented in Appendix 1.

We base the present models on rock mechanics experimental data for Berea Sandstone, a common reservoir rock analog. Our approach to adapt rock mechanics data into the Abaqus FE damage model includes the following steps:

Acquisition of experimental data of the relevant rock and deformation conditions

Conversion of the experimental data into input material parameters used in Abaqus; we used the material “concrete damage plasticity” (Abaqus manual Simulia, 2010b and below)

Creation of an FE model for the laboratory configuration

Calibrations by iterative model simulations and adjustment of material parameters to match the laboratory results (forward modeling)

Application of the calibrated material to other models and geologic analysis (part 2, Busetti et al., 2011)

The details of these steps are described in the next section.

### Experimental Data

#### Berea Sandstone: A Typical Reservoir Rock

The Berea Sandstone is commonly used as an analog for reservoir rock (e.g., Hart and Wang, 1995; Menendez et al., 1996). This rock consists of medium-grained, well-sorted, subangular quartz (∼80%), feldspar (∼5%), and clay (∼8%) and is cemented by calcite (∼6%) (Hart and Wang, 1995). The Berea Sandstone has well-documented mechanical properties (*E* = 19.3–27.5 GPa [2.8 × 10^{6}–4 × 10^{6} psi]; *ν* = 0.17–0.34; unconfined compressive strength = 71.3–74 MPa [10,341–10,733 psi]; unconfined tensile strength = 3.8–9.8 MPa [551–1421 psi]) (Table 1) that fall within the ranges of many reservoir rocks (e.g., Chang et al., 2006).

To calibrate the material parameters, we performed benchmark simulations of rock mechanics experiments with Berea Sandstone that induce tensile failure under *P*_{c}. The incorporation of both tensile and shear failure simultaneously in one rheological law is a central advantage of the present analysis. We attempted to accurately calibrate the critical features of stress-strain relations: (1) onset of inelasticity, (2) strain hardening, (3) ultimate strength, (4) strain softening and brittle failure, and (5) postfailure finite strain. The benchmark models are at a 1:1 scale to the laboratory tests, with only slight changes to the load rates to improve computational efficiency. We used the results of two different experimental configurations published by Weinberger et al. (1994), Ramsey and Chester (2004), and Bobich (2005).

#### Four-Point Beam Bending

Weinberger et al. (1994) deformed samples of the Berea Sandstone, Indiana limestone, and Tennessee sandstone in a four-point beam configuration. The beams were loaded while emplaced in a pressure vessel under a *P*_{c} of 5 to 50 MPa (725–7252 psi) (Figure 3). The four-point configuration induces near-uniform tension across the top center of the beam and was dimensioned to minimize interference between the lower two loading points. Rectangular beam samples were jacketed with polyolefin tubing and loaded inside a servo-controlled 138-MPa (20,015-psi) pressure vessel. An internal load cell monitored the axial load, *P*, applied on the top of the beam; beam strains were measured with beam-parallel strain gauges located on the top and at the bottom. The resultant maximum tension and compression stresses in the beam were calculated as done by Yokoyama (1988) (Figure 3).

#### Triaxial Extension of Dog-Bone Samples

Ramsey and Chester (2004) and Bobich (2005) deformed dog bone–shaped samples of Carrara marble and Berea Sandstone under *P*_{c} as part of their study on the transition from extension to shear fracturing. The dog-bone shape with a central cylindrical neck (Figure 4) was designed to produce uniform tensile stress at the sample center. A modeling clay layer along the neck part was constructed to transmit the *P*_{c} uniformly, and a polyolefin jacket was used around the entire sample. Pistons at the top and bottom of the sample moved simultaneously while recording axial force, axial displacement, and *P*_{c}. The experiments recorded acoustic emissions using a piezoelectric transducer attached to the top piston (Bobich, 2005).

### Finite Element Models of Experiments

We created three-dimensional (3-D) FE models in the dimensions of the laboratory tests (Figures 5, 6). The models were meshed with eight-node linear hexahedral elements, with mesh density increasing in the regions with high strain. We used mirror symmetry at the center of the beam to reduce computation time. The loading points (two rigid pins in the beam, and rigid platens in the dog bone) were modeled using 3-D rigid elements and low-friction (coefficient of friction, 0.01) master-slave penalty contacts. Loading was performed in two steps: (1) *P*_{c} was established on all of the free surfaces and was held constant throughout the simulations, and (2) piston motion was simulated by a constant-velocity boundary condition of the appropriate rigid surfaces. For the beam, the top (outer) pin was lowered at a constant velocity of 0.5 mm/s (i.e., a velocity of 0.5 mm/s for 1 min would lead to a total displacement of 30 mm) while holding the bottom (inner) pin fixed (Figure 5). For the dog-bone sample (Figure 6), axial pressure was removed and the two pistons retracted at a rate of 0.2 to 0.5 mm/s while maintaining constant *P*_{c}. These applied load rates are about 10 times faster than in the laboratory experiments and were used to reduce simulation time that took 6 to 12 h for a typical run. The ratio of kinetic energy to total system energy was found to be very low, indicating that using a higher load rate did not introduce adverse inertial effects. No chatter was observed between the rigid contacts and the rock during the axial loading step, thus, it was not necessary to use the no-separation contact option in Abaqus. We also confirmed that the material model was unaffected by the load rate by running a few tests at slower loading rates and comparing the stress-strain results.

### Material Parameters for Numerical Simulations

Table 1 lists the range of stress-strain values from experiments on the Berea Sandstone. Table 2 lists the parameters used in the FE models. Intermediate points defining the input strain-hardening and damage curves (Figure 7) were obtained by iteratively calibrating the benchmark models against laboratory results (Figures 3⇑⇑–6) to achieve stress-strain and behavior fitting within the experimental ranges reported for the Berea Sandstone while ensuring that realistic definitions for postyield behavior and other unknown plastic parameters were used. For damage-strain evolution, we use a quasisinusoidal curve (Figure 7) based on Eberhardt (1998) and Bobich (2005) (e.g., Figure 8). We selected a constant dilation angle of 15°. Other Abaqus parameters are a *K*_{c} value of 0.66 (see equation 9f), indicating a small dependence on the intermediate principal stress, as well as the default value for the ratio of initial equibiaxial to uniaxial compressive yield stress, fb_{0}/fb_{c} = 1.16 (see equation 9c), which contributes to the shape of the yield surface (Appendix 1; Figure 13). The *K*_{c} is the ratio of the length of the tensile to compressive meridians for a given pressure, that is, controlling the dependence of *σ*_{2}.

### Numerical Procedures

Details of the selected constitutive model and rheological parameters are described in Appendix 1. We used the explicit dynamic solver of Abaqus/Explicit; it is well suited for extreme nonlinearity such as strain localization or contact. When small stable time increments (10^{−6} s–10^{−8} s) are used, the simulation continues even during extensive damage propagation. We found that implicit solutions had difficulty converging during strain softening and terminated shortly after failure. Although the constitutive relationship is time independent (i.e., inferring quasistatic deformation), the analysis solves for acceleration and inertial effects, thus the modeled time represents true time, a feature necessary for analysis of transient damage propagation, which is covered in part 2 (Busetti et al., 2012). Because mesh dependence can be a difficulty in continuum damage simulations, we experimented with several different meshing schemes. Variations in the mesh density did not significantly affect the global solution (e.g., the stress-strain curve). Capturing the localized pattern of damage necessary to discern individual fractures required a finer mesh density with elements on the scale of a few millimeters. Damage localization patterns were generally consistent for similarly sized and patterned meshes, although the specific fracture geometry varied somewhat.

#### Macroscopic Failure and Fracturing

The Abaqus software allows implementation of the damage parameter by either specifying *d* as a function of plastic strain during uniaxial loading (Figure 7) on a load-displacement curve (Rice, 1968) or by directly specifying the fracture energy, *G*_{f}, in compression or tension. The FE implementation permits simulation of fracture propagation based on the equivalent crack concept, which states that there exists a length-scaled damage zone that is thermodynamically equivalent to a crack and vice versa (see Appendix 2) (Mazars and Pijaudier-Cabot, 1996). Using this approach, a discrete fracture is represented by an equivalent fracture zone: a path of elements that are completely damaged and thus have no strength.

## SIMULATIONS RESULTS

### Macroscopic Stress-Strain Relations

We ran numerous calibration simulations on the beam and dog-bone configurations. We plotted the calculated differential stress, Δ*σ*, as a function of strain across the top center of the beam (Figure 5) and in the central part of the dog bone (Figure 6) using averaged element integration point values. When a reasonable fit between simulation and experiments was achieved, we ran an additional refined simulations on the beam model for 10 MPa (1450 psi) *P*_{c} and on the dog-bone model for *P*_{c} = 10, 20, 30, 40, 50, 60, 80, 100, 150 MPa (1450, 2901, 4351, 5802, 7252, 8702, 11,603, 14,504, and 21,756 psi). A good fit in the simulations should capture the essential stages of strain hardening, failure, and softening at values lying within the ranges observed experimentally (Figures 3, 5, 8; Table 1). In comparison with experimental results, the simulated rheology underpredicts nonlinearity associated with strain hardening and ultimate strength and overpredicts failure strain in the beam while underpredicting failure strength in the dog bone (Figures 5, 6). We think that this issue is not a limitation of the model, and that by further refining the plastic parameters, a closer reproduction can be achieved; however, more detailed calibration is beyond the present scope.

### Deformation and Failure Features

#### Damage Propagation

The distribution of damage in the simulations shows good agreement with experimental observations. Figures 9 and 10 show contour plots of the simulated tensile damage in the two experimental configurations; for brevity, we highlight a description of the beam. In the beam simulation, the damage initiates in a broad region corresponding to maximum strain located off center and at the beam top. The off-center damage initiation can be predicated from the curvature plots, *d*^{2}*y*/*dx*^{2}, along the beam top (solid lines, Figure 11). The figure shows that the center of the beam top is not quite uniformly curved, and maximum curvature occurs about 2 cm (0.8 in.) from the beam center (most likely affected by the lower fixed points). Damage initiates (color contours in Figure 9, dashed lines in Figure 11) at 3.5 cm (1.4 in.) from the beam center and increases nonlinearly toward the region of maximum curvature, reaching a maximum value of 1 at about 1.5 cm (∼0.6 in.). Before the first fracture, the damage in the top center of the beam increases to *d* ≈ 0.1, which corresponds to 10% stiffness reduction. Multiple fractures form in regions of 15 to 20% stiffness reduction.

#### Onset of Plastic Yielding

The initial linear parts of the stress-strain curves reflect the Young's modulus of the undamaged rock (Figures 5, 6). From the dog-bone simulations, the elastic stage accounts for the initial 0.0006 to 0.0014 strain at a low *P*_{c} (10–30 MPa [1450–4351 psi]), which increases linearly to 0.0063 at a *P*_{c} = 150 MPa (21,756 psi). Calculating the linear dependence of plastic onset on the *P*_{c} yields*ε*_{PO} is the plastic onset strain and *P*_{c} is the confining pressure in pascals.

#### Strain Hardening

For a *P*_{c} less than 30 MPa (4351 psi), the strain-hardening stage is pervasive and accounts for up to 50% of the prefailure strain, displayed as the broad hump on the stress-strain curve in Figure 5 and for the curves at a lower *P*_{c} in Figure 6. For *P*_{c} = 30 to 60 MPa (4351–8702 psi), strain hardening is greatly reduced and the hump essentially disappears. With increasing *P*_{c} = 60 to 150 MPa (8702–21,756 psi), strain hardening gradually increases to approximately 10% of the total prefailure strain. Extensional strain at ultimate failure occurs at 0.0015 for a *P*_{c} less than 30 MPa (4351 psi), increasing to 0.007 at 150 MPa (21,756 psi). For the beam configuration, results are not shown for compressive stress-strain at the bottom of the beam because the elastic limit was not exceeded.

#### Strength

Consistent with the yield surface definition (Appendix 1, equation 9b), the ultimate strength increases linearly with *P*_{c}, that is, constant internal friction angle. Failure occurs in the regions where stress is locally least compressive: under local tension for a low *P*_{c} (<20 MPa [<2901 psi]) and under local compression for a *P*_{c} more than 20 MPa (2901 psi). The yield-strain linear relation to *P*_{c} is*ε*_{UY} is the ultimate yield strain and *P*_{c} is the confining pressure in pascals.

#### Softening and Finite Strain

The macroscopic failure at the ultimate strength is followed by an immediate stress drop (strain softening) with a transition to constant stress deformation (Figures 5, 6). Note that, in Figure 6, the stress drops were only observed if stress was plotted against monotonic axial strain (i.e., in the figure the strain curve was extrapolated to monotonic because, on failure of individual elements, an immediate jump in strain occurs that effectively canceled out the stress drops on a stress-strain plot. Alternatively, strain could have been measured from displacement of the rigid platens). The plastic parameters (Table 2; Appendix 1) determine the character of the stress drop, where a more rounded peak and shallower postfailure slope reflects a more ductile behavior and a sharper peak, and a steep slope indicates brittle behavior. The simulations show that increasing strain softening as well as the rate and magnitude of damage accumulation contributed to the more brittle behavior observed experimentally. Brittleness also increased using nonassociated plastic flow, that is, decreasing the dilation angle, *ψ*, to 15° (e.g., Alejano and Alonso, 2005), resulting in an additional stress drop of 1 to 2 MPa (145–290 psi).

#### Stress Path

We plot the stress paths for the dog-bone simulations to examine the variations of the stress tensor during progressive deformation. Following Harrison and Hudson (2003), Figure 12 shows the incremental change in pressure stress *P*_{c} = 10–150 MPa (1450–21,756 psi). All simulations begin at a hydrostatic stress state (*q* = 0). Retracting the pistons decreases *p*′ and increases *q* along a constant linear path during elastic deformation. The onset of strain hardening, plasticity, and damage causes the stress path to curve to the left as the deviatoric stress decreases relative to the pressure change. Ultimate yielding and macroscopic failure occur as the stress path intersects the yield surface (Figure 13) and is followed by a rapid stress drop.

## DISCUSSION

### Synthesis of Simulation Results

#### Calibrating Simulated Results to Experimental Data

A main strength of the present modeling approach is its ability to generate realistic finite deformation and to simulate stress-strain and damage curves by iterative calibration of material parameters. The final results lie within published ranges for the Berea Sandstone (Table 1); however, the match between simulated and experimental results (solid lines vs. dots in Figures 5, 6) is not perfect. Although it was possible to generate a Berea Sandstone rheology to generally fit the experimental ranges, it was not a simple task to precisely match both sets of tests results simultaneously. In addition to rheologic considerations, accurate calibrations of the numerical rheology using laboratory data depend on (1) capturing relevant global effects, such as load frame stiffness, which can affect postfailure behavior; (2) data sampling (e.g., location, size, frequency, etc.), for instance, using averaged integration point values from a few elements to calculate stress and mimic strain gauge readings versus determining stress from piston loads or measuring axial displacement at the sample ends; (3) unknown environmental factors including frictional contact parameters; and (4) scaling issues, meshing, and other numerical considerations. Despite these limitations, we find the simulations to be suitable to study the results of complex rheologic behavior as displayed in numerous field observations and rock mechanics experiments.

#### Main Simulation Results

*Deformation Stages*—The simulations capture the main stages of rock deformation, including (1) onset of inelasticity, (2) strain hardening, (3) ultimate strength, (4) strain softening and brittle failure, and (5) postfailure finite strain.

*Inelasticity Onset*—Under confining pressure, inelastic deformation initiates within the regions of lowest compressive stresses as early as 0.1% strain.

*Stiffness Reduction*—Strain-hardening initiates with the onset of inelasticity with a shown stiffness reduction of 10 to 15% before tensile fracturing (e.g., time points preceding Figure 9A). Fractured zones display a stiffness reduction of 20% or more (Figure 9B–D).

*Damage Pattern*—The simulated damage patterns resemble those of the rock mechanics tests. The four-point beam configuration displays multiple fractures that initiate at the beam's extension side and extend toward the neutral axis. In the dog-bone configuration, damage localizes in the center of the sample and the pattern transitions from simple to complex with increasing confinement.

### Damage Fracturing

#### Extension under Confining Pressure

The experimental configurations that are analyzed here, four-point beam and dog-bone triaxial, have a very important deformation feature: induced local extension under global confining pressure. This condition is likely to develop in many field situations, for example, fault-fold systems (Reches and Johnson, 1978), dike emplacement and hydrofracturing (Delaney et al., 1986), and fault tip zones (Vermilye and Scholz, 1999). The present and previous simulations (e.g., Lyakhovsky et al., 1997), as well as microstructural analyses (e.g., Lockner et al., 1992; Katz and Reches, 2004), show intense damage development in the regions of localized extension that commonly leads to macroscopic failure. We refer to this process as damage fracturing, and its product is damage fractures. The damage fracturing process is analyzed below by linking the known rheological parameters in the simulations and the simulated damage features to the observed stress-strain relations and fracture patterns in the experiments.

#### Early Damage Distribution in Fresh Rock

*Experiments*—In the beam experiments, Weinberger et al. (1994) described two modes of tensile yielding: minor early stage and major late stage fracturing. They distinguished the early-stage fracturing as the lower bound of the tensile strength and the ultimate stress as the upper bound, where large fractures were associated with a major stress drop.

*Simulations*—Early-stage damage is controlled by stress-induced yielding and the onset of inelasticity (top Figure 7) as manifested by the continuous strain hardening. Damage evolves slowly and stably (dotted lines, Figures 5, 6). We note that, in the simulated stress-strain curves plotted in Figure 5, small stress drops representing local fracture events were not distinct as in the experiments; however, some minor events do appear.

#### Late Propagation of Highly Localized Damage Zones

The coalescence of major fractures within a pervasively damaged region controls unstable late-stage deformation and failure. Thus, the active failure criterion reflects the preconditioned state of a damaged rock with degraded strength.

*Experiments*—From the beam experiments, it was interpreted (Weinberger et al., 1994) that the early fractures propagated stably from the tensile surface, whereas the late-stage fractures propagated unstably toward the beam neutral surface and vanishing tensile stresses.

*Simulations*—Late-stage propagation is controlled by an inelastic strain evolution during the early stage that preconditioned the extension zone for failure (note blue-green damage zone in Figure 9). Preconditioning damage also forms ahead of the fracture tip, where it controls the propagation direction of the fracture as it grows into less stressed zones. Figure 9 further shows that when the fractures are isolated (widely spaced), fracture propagation paths are linear and continuous (Figure 9A, B). However, with advanced fracturing, damage zones link, leading to complex fracturing and crosscutting patterns (Figure 9D). The stage is associated with a global failure of the beam, late strain softening, and follows a major stress drop.

#### Are Damage Fractures Shear or Tensile Fractures?

Traditionally, shear and tensile failures are defined by the stress state on the failure surface starting with the works of Coulomb (1773) and Griffith (1921). Thus, application of the relevant failure criterion, for example, the Coulomb criterion, analytically predicts the orientation of an idealized fracture with respect to a given stress tensor (e.g., Anderson, 1951). Damage fractures are different; their initiation and propagation depend on both stress and strain fields and are controlled by several nonlinear and partly interdependent parameters, for example, volumetric strain, tensile stresses, loading history, and strain history. Lyakhovsky et al. (1997) analytically determined the complex relations between these parameters for a series of cases with homogeneous loading and found a good agreement with the macroscopic experimental data of stress and strain. However, the shapes and orientations of damage fractures under nonhomogeneous loading can be determined only numerically. For example, Lyakhovsky et al. (2001) studied damage in a layered crust subjected to basal horizontal shear, in which the upper part has elastic damage rheology. They determined the damage and fracturing patterns in this nonhomogeneous system by numerical simulations that delineated a nontrivial fault structure, the shape of which depends strongly on the rate of healing. Similarly, in the present analysis, a nontrivial network of damage fractures developed in both the simulation of the four-point beam (Figure 9D) and dog-bone experiments (Figure 10). Most of the deformation associated with these damage fractures is likely tensile, yet they are not simple smooth tensile fractures. We will explore in detail the morphology and nature of damage fractures in part 2 (Busetti et al., 2012) of this study.

## Acknowledgments

Funding and resources for this project came from the ConocoPhillips School of Geology and Geophysics and the School of Civil Engineering and Environmental Science at the University of Oklahoma, Devon Energy, and ConocoPhillips. We thank Vincent Heesakkers, Tricia Allwardt, K. K. Muraleetharan, David Lockner, and Peter Hennings for their helpful discussions. We also thank Fred Chester and Andreas Kronenberg of Texas A&M for kindly providing details of the experimental data. The finite element simulations were conducted using the academic version of the code Abaqus (a product of Simulia, a division of Dassault Systèmes). We found the reference material in the digital Abaqus manual to be very helpful. We also are extremely thankful to the two manuscript reviewers whose comments and recommendations were most invaluable.

The AAPG Editor thanks the anonymous reviewers for their work on this paper.

## APPENDIX 1: FINITE ELEMENT DAMAGE AND CONSTITUTIVE MODEL

### Damage Formulation

We used a damage model (Lubliner et al., 1989; Lee and Fenves, 1998) included in Abaqus (Simulia, 2010a, b) that applies the concepts of deformation modulus (E) and stiffness reduction, and uses the damage parameter, *d*, as a dimensionless approximation for stiffness degradation to scale the true stress. In the initial stage, *D* = *D*_{0}, *d* = 0 (no degradation); and at failure, *d* = 1, the material is completely damaged and the effective stress drops to zero.

The incremental plastic strain includes all irreversible deformations as well as brittle microcracking damage. Decomposing strain into elastic and plastic strain components (*ε* = *ε*^{e} + *ε*^{p}) gives the stress-strain relationship

The effective-stress concept (here unrelated to pore pressure) is used to degrade the elastic stiffness, which in turn controls the shape of the yield surface. The damage parameter evolves separately as a function of plastic strain and is tracked separately for tension and compression damage (Figure 7). The result is a material that retains directional strength depending on how it is strained. For example, a region containing microcracks oriented normal to the loading direction is very strong in compression, but weak in tension. The constitutive equations couple the evolution of *d* with plastic strain evolution of the dissipation potential, *d* = *d* (*κ*_{ℵ}) (Lubliner et al. 1989):

To distinguish mode I and II damage, damage dissipation κ_{ℵ} is defined separately for compression or tension. The term *g*_{ℵ} is the normalized dissipated energy during microcracking and is set under the uniaxial stress state. For a continuum framework, *g*_{ℵ} normalizes the energy released during compressive or tensile fracturing (defined by fracture mechanics theory) *G*_{ℵ} by a localization size, *l*_{ℵ} (Lubliner et al., 1989; Lee and Fenves, 1998):

The value of *l*_{ℵ} is an “objective value” or “assumed as a material property” (Lee and Fenves, 1998), and in the Abaqus FE formulation, it appears as an element regularization parameter.

### Constitutive Model

The constitutive model incorporates a pressure-dependent yield criteria, a plastic flow rule, a hardening rule, and damage. The combination of strain softening and plastic damage permits simulation of extreme localized weakening, where fractures or fracture zones (depending on the fineness of the FE mesh) occur where material degrades to zero strength. The model is based on modifications of classical Mohr-Coulomb plasticity, as described below.

The onset of yielding that occurs at the brittle ductile transition in rocks is typically pressure dependent (Murrell, 1965; Brace et al., 1966; Jaeger and Cook, 1976). The Mohr-Coulomb criterion described this relationship by relating the shear and normal stress across a plane by the function:*τ* is the shear stress, *σ*_{n} is the normal stress, *μ* = tan(*ϕ*) gives the coefficient of friction and friction angle, and *c* is the shear strength or cohesion. Byerlee showed experimentally that a value of *μ* = 0.6 to 0.8 was a property intrinsic to most upper crustal rocks. The Mohr-Coulomb criterion (Figure 13A) assumes that failure does not depend on the intermediate principal stress, yet Mogi (1973) and Reches and Dieterich (1983) showed that this assumption is not necessarily valid in the general case. By assuming that the yield surface is fully dependent on the intermediate principal stress, sets *τ*_{oct} is constant for all rotations of θ, as stated by the Drucker-Prager criterion, where the shape of the yield surface is a cone that opens with increasing mean stress. The Drucker-Prager failure surface is a pressure-dependent variation on the von Mises criterion and is given by*α* and *k* are material constants, and the stress invariants in terms of the principal stresses *σ*_{1}, *σ*_{2}, and *σ*_{3} and the mean stress *σ*_{0} are

We use the modified yield surface of Lubliner et al. (1989), which combines positive features of both the Mohr-Coulomb and Drucker-Prager models. Figure 13B shows the Barcelona model yield surface (Alonso et al., 1990). The Drucker-Prager surface (equation 7a) fits within the form*σ*_{max} and assigning two additional parameters, *β* < *σ*_{max} > and *γ* < *σ*_{max} > to equation 9a, the Drucker-Prager circle can be modified to reduce the dependence of the intermediate principal stress and establish tensile and compressive meridians. Thus, the form presented by Lubliner et al. (1989) and Lee and Fenves (1998) is

The uniaxial compressive and tensile yield strength are given by *σ*_{c0} and *σ*_{t0} and can be taken directly from experimental data. The biaxial tensile strength *σ*_{b0} is approximately 1.3% less than the uniaxial tensile strength (Lee and Fenves, 1998), yielding typical values of 0.08 ≤ *α* ≤ 0.12 and 1.10 ≤ *β* ≤ 1.16 (Lubliner et al., 1989). The *K*_{c} is the ratio of the length of the tensile to compressive meridians for a given pressure, that is, controlling the dependence on *σ*_{2}. For *K*_{c} = 1, *β* and *γ* drop out, leaving the original DP function:*p* = I_{1}/3 = *σ*_{0}, and TM and CM are the tensile and compressive meridians.

This constitutive model has a few useful qualities for our purposes. First, it accommodates a broad range of rock types that were previously modeled by Mohr-Coulomb or Drucker-Prager failure. Second, in modeling new rock types from experimental data, particularly in reservoir applications where core is limited, it is convenient to use uniaxial failure and stress-strain data for calibrations. Third, the model is implemented in the commercial FE software Abaqus, which can be used for a wide range of geologic problems.

## APPENDIX 2: FRACTURE ENERGY EQUIVALENCE

A brief point of clarity is helpful regarding fracturing from a classical fracture mechanics theory compared with the continuum damage method used here. The fracture mechanics approach (Griffith, 1921) and subsequent refinements (e.g., Barenblatt, 1962; also see Busetti, 2009) are useful for describing the growth of discrete cracks propagating from initial flaws. Continuum damage theory is well suited when describing the global averaged behavior of fracture networks (e.g., smeared crack models) (Bazant and Oh, 1983; de Borst and Nauta, 1985; Mazars and Pijaudier-Cabot, 1989) or, locally, the evolution of microcracks (e.g., Steinmann et al., 1994; Mazars and Pijaudier-Cabot, 1996). The present work fits the latter perfectly but also expands to model discrete fracture growth (Busetti, 2009; part 2, Busetti et al., 2012). Thus, note that the two theories are thermodynamically equal according to the equivalent crack concept (Mazars and Pijaudier-Cabot, 1996), which states that there exists a damage zone that is equivalent to a discrete fracture and vice versa. In the present formulations, this manifests in the normalized dissipated energy term *g*_{ℵ}. The energy consumed by forming all of the microcracks in a volume is equivalent to decohesion of singular cracks with surface area *A*_{d}:*Y* is the damage energy release rate, *d* is the damage parameter increment, and −*G* = *G*_{f} is the fracture energy release rate. On the left side, the damaged area (or volume) is calculated directly using the energy density function. Following Mazars (1986), an equivalent crack, Ae, is attained from equation 10a:*dd*(*x*) represents evolution of the damage at a point, *x*. The total damaged area then reflects the summation of the area comprising all equivalent cracks over a volume, *V*, which is the definition of the energy density function, *ρ*. This equivalence relates continuum and discrete fracture theory, but also offers a means for comparison between numerically predicted damage and the damage observed in experiments and quantified using other stiffness reduction models.

- Manuscript receivedJanuary 27, 2011.
- Revised manuscript receivedJuly 1, 2011.
- Revised manuscript receivedSeptember 9, 2011.
- Revised manuscript receivedNovember 28, 2011.
- Revised manuscript receivedJanuary 3, 2012.
- Final acceptanceFebruary 1, 2012.

Seth Busetti received his B.S. degree in geological engineering from the Colorado School of Mines and his M.S. degree and Ph.D. in structural geology from the University of Oklahoma. Busetti is currently employed at ConocoPhillips and works in a subsurface technology position focusing on applied structural geology and geomechanics problems worldwide. He is currently involved in projects involving fracture and fault mechanics, fluid flow in fractured reservoirs, and geomechanics in nonconventional reservoirs.

Kyran Mish received his B.S. degree in mathematics, his M.S. degree in structural mechanics, and his Ph.D. in engineering, all from the University of California at Davis. Before his present position at Sandia National Laboratories, Mish was in management at Lawrence Livermore National Laboratory, served on faculty at University of California at Davis, and most recently was director of the Donald G. Fears Structural Engineering Laboratory and professor of structural engineering at the University of Oklahoma.

Ze'ev Reches received his B.S. and M.S. degrees in geology from Hebrew University, Israel, and his Ph.D. in structural geology from Stanford University, California. Reches serves as a professor of structural geology at the University of Oklahoma. His prior work includes positions at Arizona State University, Stanford University, and the U.S. Geological Survey at Menlo Park, California, and Hebrew University, Israel. His research interests include earthquake and fault processes and rock mechanics.