PUBLICATIONS Journal of Geophysical Research: Solid Earth RESEARCH ARTICLE 10.1002/2014JB011805 Key Points: • We quantify the role of deformation mechanisms using new mathematical model • Formulation accounts for elastic deformation, cataclastic flow, and brittle failure • Deformation bands formation includes dilation, high angle shear, and compaction Correspondence to: V. Lyakhovsky, [email protected] Citation: Lyakhovsky, V., W. Zhu, and E. Shalev (2015), Visco-poroelastic damage model for brittle-ductile failure of porous rocks, J. Geophys. Res. Solid Earth, 120, doi:10.1002/2014JB011805. Received 3 DEC 2014 Accepted 8 MAR 2015 Accepted article online 11 MAR 2015 Visco-poroelastic damage model for brittle-ductile failure of porous rocks Vladimir Lyakhovsky1, Wenlu Zhu2, and Eyal Shalev1 1 Geological Survey of Israel, Jerusalem, Israel, 2Department of Geology, University of Maryland, College Park, Maryland, USA Abstract The coupling between damage accumulation, dilation, and compaction during loading of sandstones is responsible for different structural features such as localized deformation bands and homogeneous inelastic deformation. We distinguish and quantify the role of each deformation mechanism using new mathematical model and its numerical implementation. Formulation includes three different deformation regimes: (I) quasi-elastic deformation characterized by material strengthening and compaction; (II) cataclastic flow characterized by damage increase and compaction; and (III) brittle failure characterized by damage increase, dilation, and shear localization. Using a three-dimensional numerical model, we simulate the deformation behavior of cylindrical porous Berea sandstone samples under different confining pressures. The obtained stress, strain, porosity changes and macroscopic deformation features well reproduce the laboratory results. The model predicts different rock behavior as a function of confining pressures. The quasi-elastic and brittle regimes associated with formation of shear and/or dilatant bands occur at low effective pressures. The model also successfully reproduces cataclastic flow and homogeneous compaction under high pressures. Complex behavior with overlap of common features of all regimes is simulated under intermediate pressures, resulting with localized compaction or shear enhanced compaction bands. Numerical results elucidate three steps in the formation of compaction bands: (1) dilation and subsequent shear localization, (2) formation of shear enhanced compaction band, and (3) formation of pure compaction band. 1. Introduction Within the geological pressure and temperature range of the seismogenic zone (i.e., the upper 10–15 km), rocks often fail by shear localization or by cataclastic flow. The latter is a special semibrittle failure mode in which the distributed macroscopic deformation is associated with pervasive microcracking [e.g., Fredrich et al., 1989]. While the macroscopic manifestation of brittle faulting and cataclastic flow failure modes are drastically different (i.e., localized versus distributed deformation), the micromechanics of both regimes involves stress-induced microcracking at the grain scale [e.g., Wong et al., 1997]. Responses of low-porosity rocks to mechanical loads are porosity generation via microcracking and the resultant permeability enhancement during both the brittle faulting and cataclastic flow [e.g., Zhu and Wong, 1999]. In contrast, for porous siliciclastic rocks permeability reduction may actually be observed with an appreciable amount of dilation in the brittle faulting regime, whereas significant porosity reduction with concomitant permeability decrease is generally observed in the cataclastic flow regime [e.g., Zhu and Wong, 1997]. This is because inelastic deformation behavior of high-porosity rocks, when subjected to compressive loading, includes damage accumulation in form of microcracking coupled with reduction of preexisting porosity. Depending on the material properties, confining pressure, and loading rates, deformation behavior in porous rocks is controlled by the interplay between shear, dilation, and compaction. Under low effective pressures, damage localizes forming shear bands and the material undergoes dilatancy; under high effective pressures, the material compacts and damage are homogeneously distributed in the samples [e.g., Tembe et al., 2007]. In the intermediate effective pressures, the transition from localized shear dilation to homogeneous compaction includes localized compaction [e.g., Holcomb et al., 2007; Wong and Baud, 2012]. The style of localization from shear to pure compaction bands [Baud et al., 2004; Charalampidou et al., 2014] depends on the lithology and initial porosity of the rock [Baud et al., 2004]. Numerous field observations have shown the existence of compaction bands [Antonellini and Aydin, 1994; Eichhubl et al., 2010], shear bands [Antonellini et al., 1994; Sternlof et al., 2005; Eichhubl et al., 2010], and dilation bands [Du Bernard et al., 2002; Fossen et al., 2007]. Several types of deformation bands are also found LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 1 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 to coexist in the field [Sternlof et al., 2006; Aydin and Ahmadov, 2009; Fossen et al., 2011]. Different types of deformation bands may either act as fluid conduits or flow barriers [e.g., Zhu and Wong, 1997; Sun et al., 2011; Baud et al., 2012]. Understanding fluid flow and faulting in the upper crust requires better model on inelastic behavior and failure modes in porous rocks. Inelastic and failure behaviors of porous rocks, including brittle failure at low pressures and cataclastic flow at high pressures, are often described using an elastoplastic formulation also referred as bifurcation analysis [e.g., Rudnicki and Rice, 1975; Aydin and Johnson, 1983; Olsson, 1999; Issen and Rudnicki, 2000; Zhu and Walsh, 2006; Besuelle, 2001; Rudnicki, 2007; Chemenda, 2011; Chemenda et al., 2012]. Laboratory experiments show that the yield behavior of porous rocks evolves during deformation of porous rocks [Baud et al., 2006; Tembe et al., 2007; Wong and Baud, 2012]. Grueschow and Rudnicki [2005] discussed different models of the evolving yield cap [DiMaggio and Sandler, 1971; Carroll, 1991] and suggested new constitutive model with elliptic yield cap. The bifurcation analysis is successful in predicting the onset of deformation bands, but it does not produce temporal evolution of deformation bands, because the time dependency is not explicitly included in their formulation. Recent numerical models [Foster et al., 2005; Sun et al., 2014] included ad hoc kinetic equations for the constitutive relation that controls the temporal evolution of the yield cap. However, this approach lacks well-constrained deformation mechanisms and does not account for gradual evolution of rock elastic properties associated with microscopic damage and porosity change. Evolving rock properties are considered in the models based on the continuum breakage mechanics for granular media [Einav, 2007a, 2007b]. This approach is relevant only for cataclastic flow under relatively high pressures and does not address brittle failure [Das et al., 2013]. Following pioneering work by Kachanov [1958], many researchers have developed damage theories for brittle and ductile deformation that account for temporal evolution of rock properties. Some of these theories have been adopted for the rheology of high-porosity rocks and capture important features of different deformation regimes [e.g., Ricard and Bercovici, 2009; Cai and Bercovici, 2014; Karrech et al., 2014]. The theoretical model developed by Hamiel et al. [2004a, 2004b, 2005] is based on the thermodynamic principles [e.g., Biot, 1973; Coussy, 1995] leading to a system of coupled kinetic equations for the evolution of damage and porosity. A significant advantage of their model is the ability to simulate the entire yield curve by a single formulation. In this paper we modify the Hamiel et al. model by incorporating the evolving porosity-dependent yield conditions and simulate the deformation and failure behaviors of porous sandstones under different stress states. We demonstrate that the model is capable to quantify the deformation processes that take place during mechanical loading of porous rocks and compare the numerical results with the experimentally observed deformation behavior of the Berea sandstone samples reported by Zhu and Wong [1997] and Wong et al. [1997]. Our model establishes links between stress, microstructure, and the inelastic behavior and failure mechanisms in actively deforming porous rocks. 2. Visco-Poroelastic Damage Model 2.1. Conceptual Model Traditionally, the analysis of rock deformation and failure criteria has been formulated by the classical Coulomb-Mohr conditions that define the brittle failure and by a yield cap criterion that defines cataclastic flow [e.g., Issen and Rudnicki, 2000]. The failure/yield surface (stress-space formulation) represents the boundary between the elastic and plastic behavior of a granular rock. Loading a rock sample to stress beyond the initial yield surface will cause significant plastic strain. Consequently, if the sample is unloaded and reloaded, the new yield cap is at higher stress states. However, the stress-space formulation usually ignores connection between yield stress and amount of damage in the form of microcracks, voids, or any other type of flaws [e.g., Han and Chen, 1986]. As a result, the application of the stress-space formulation encounters difficulties in modeling the strain-softening behavior since it cannot explicitly give a clear description of the evolution of the yield cap. Several authors formulated yield conditions in terms of strains (strain-space formulation) and demonstrated significant advantages of this approach for materials with evolving yield conditions as a function of material properties and loading history [Naghdi and Trapp, 1975; Yoder and Iwan, 1981; Han and Chen, 1986; Lehane and Simpson, 2000; Puzrin and Houlsby, 2001]. Although there are some similarities between stress-space and strain-space formulations, they are not equivalent when material weakening is considered [Casey and Naghdi, 1983; Einav, 2004, 2005a, 2005b]. Porous rocks exhibit a strong strain-softening postfailure behavior, indicating a LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 2 Journal of Geophysical Research: Solid Earth Figure 1. Three deformational regimes that are defined by the yield cap and the compaction-dilation (C-D) transition. Loading path 1 represents loading in low effective pressure and occur in Regimes I and III. Loading path 2 represents loading in high effective pressure and occur in Regimes I and II. Loading in Regime II increases the yield cap as the rock is undergoing compaction (Δφc) and cataclastic flow (Δεd). The amount of the growth of the yield cap (f) is proportional to compaction rate. 10.1002/2014JB011805 significant elastoplastic coupling for the degradation of elastic modulus with increasing plastic deformation. Stress-space formulation of plasticity based on Drucker’s stability postulated for these materials encounters mathematical difficulties in modeling the softening/elastoplastic coupling behavior because the need to adjust both the elastic properties and the yield surface simultaneously. The advantage of formulating the model in the strain space is that the failure envelopes do not depend on strain softening and hardening, and thus, the mathematical formulation is robust during rock property evolution. Furthermore, numerical methods such as finite elements are formulated in terms of displacements (strains), and the description of the yield cap in the stress-space is not efficient. The continuum damage-porosity model presented here considers two thermodynamic state variables, i.e., porosity and damage intensity coupled by kinetic relations discussed in the next section. The suggested model addresses three different deformation regimes that are defined by two envelopes (Figure 1). Heavy line schematically defines the yield cap and straight-dotted line stands for the transition between compaction and dilation (C-D), i.e., onset of the inelastic porosity increase. The C-D transition coincides with Coulomb-Mohr failure criterion at very slow loading rates. These lines define three zones in the differential-volumetric strain diagram: 1. Beneath the yield cap (Regime I) the quasi-elastic deformation is accompanied by material strengthening associated with crack closure and compaction (porosity decrease). 2. Regime II is associated with stress-induced damage (microcracking) and porosity decrease. The damage accumulation starts when the loading crosses the yield cap. Distributed microcracking and grain crushing enable sliding along newly created internal surfaces, collapse of pore space, and change of grain packing. These processes are most prominent under high confining pressures usually treated as cataclastic flow associated with compaction. 3. Regime III is associated here with macroscopic brittle failure. When the loading is beyond the C-D transition, intensive damage accumulation along with significant differential strains lead to material dilation (porosity increase). Damage increase is bounded by certain critical value (damage, α equals 1) corresponding to macroscopic brittle failure, dynamic stress drop, and fast conversion of the differential elastic strain into plastic strain components. When a porous rock is loaded with low effective (confining-pore) pressure (loading path 1 in Figure 1), the strain state crosses the C-D transition, which overlaps with the yield cap at low confinement, and the rock fails by shear localization. Under typical triaxial laboratory conditions, the peak stress is reached above the onset of dilation (C-D transition). The distance of the strain state at failure from the C-D transition depends on the rock properties (damage accumulation rate) and on the loading rate. At very low loading rates or very fast damage accumulation, the rock fails just at the C-D transition, which then represents the Coulomb-Mohr failure criterion. The strain accumulated above the C-D transition increases with loading rates leading to elevated peak stress (high strength). Such connection between loading rate and measured rock strength, widely observed in rock mechanics experiments [e.g., Paterson and Wong, 2005], is not addressed by classical plasticity theory utilizing Coulomb-Mohr failure criterion. When the rock fails, a stress drop reduces the strains back to the field corresponding to Regime I. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 3 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 Under loading conditions with high confinement (loading path 2), the strain state crosses only the yield cap, where damage increases and compaction are taking place. In this path, the distance of the strain state from the cap depends on the loading rate, the rate of shear strain relaxation (effective viscosity), and compaction rate (porosity decrease). While the rock is loaded beyond the yield cap, inelastic compaction acts to reduce elastic volumetric strain and at the same time to expand the yield cap (red arrows). Cataclastic flow reduces the differential strain (Maxwell relaxation) and also shifts the loading path. Compaction and cataclastic flow continues as long as loading continues. As a result, the yield cap keeps growing, and thus, it follows the strain state which keeps shifting with response to compaction and cataclastic flow. The rates of compaction and cataclastic flow depend on the distance between the strain state and the yield cap. If initially, compaction and cataclastic flow are slower than loading rate, they both increase as the strain state distance from the yield cap grows, until they are fast enough to keep pace with loading, and their rates reach a steady state in which the yield cap grows at the same rate as the loading rate. 2.2. Mathematical Formulation In this study we extend the continuum damage-porosity model developed by Hamiel et al. [2004a, 2004b] to account for (1) non-linear elasticity that connects the effective elastic moduli to a damage variable and loading conditions; (2) evolution of the damage state variable as a function of the ongoing deformation and gradual conversion of elastic strain to permanent inelastic deformation during material degradation; (3) evolution of the material porosity, including pressure-driven compaction, as well as damage-related dilation and compaction; (4) macroscopic brittle instability at a critical level of damage and related rapid conversion of elastic strain to permanent inelastic strain and material porosity; and (5) coupling between deformation and porous fluid flow through poroelastic constitutive relationships incorporating damage rheology with Biot’s poroelasticity. Governing equations are summarized in this section (for details, see Appendix A). First we present general equations which are common for all deformation regimes and then discuss damage-porosity kinetic equation for each deformation regime separately. The free energy of a poroelastic medium, F, is a sum of the elastic energy and the poroelastic coupling term of the saturated medium: F¼ pffiffiffiffi 1 λðα; ϕ Þ 2 I1 þ μðα; ϕ ÞI2 γðα; ϕ ÞI1 I2 þ M ½αB I1 ðζ ϕ Þ2 ; 2 2 (1) where I1 = εkk and I2 = εijεij are the first and second invariants of the elastic strain tensor, defined as a t ir difference between the total and irreversible strain components, εij ¼ εij εij . The total strain tensor is expressed through the displacements, εtij ¼ 1=2 ∂ui =∂x j þ ∂uj =∂x i . The deviatoric part of the irreversible strain εirij is accumulated during the viscous, Maxwell-type, relaxation, while change of the volumetric component is equal to the porosity change (dεirii ¼ dϕ). The fluid volume content is ζ, ϕ is the material porosity, α is the damage variable, λ, and μ, are the Lame drained moduli, γ is the strain coupling modulus, M is the reciprocal of the specific storage coefficient at constant strain, and αB is the Biot coefficient. Differentiation of the poroelastic energy (1) leads to constitutive relations for the stress tensor, σ ij: ∂F γ ¼ λ I1 δij þ ð2μ γξ Þεij þ αB M½αB I1 ðζ ϕ Þδij ; σ ij ¼ ∂εij ξ (2) Following Agnon and Lyakhovsky [1995], we assume that the elastic moduli λ, μ, and γ depend linearly on an evolving damage state variable 0 ≤ α ≤ 1. λ ¼ const: μ ¼ μ0 þ γr ξ 0 α; (3) γ ¼ γr α where μ0 is the initial shear modulus and ξ 0 and γr are constant. Fluid pressure is proportional to the difference between total porosity or water content (ζ ) and material porosity (ϕ) and first invariant of the strain tensor representing the volumetric elastic deformation [e.g., Coussy, 1995]: p ¼ M½αB I1 þ ðζ ϕ Þ; LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. (4) 4 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 where M = 1/Sε is the reciprocal specific storage coefficient at constant strain. The fluid mass conservation equation combined with the Darcy law leads to pressure-diffusion equation for saturated isothermal flow [e.g., Wang, 2000]: ! k ðα; ϕ Þ ∂p ∂εkk ∇ – ð∇p þ ρgz Þ ¼ Sε þ αB ; (5) ∂t ∂t μf where k is the intrinsic permeability of the medium, μf is the fluid viscosity, ρ is the fluid density, g is the gravitational acceleration, and z is a unit vector. We assume that the permeability depends exponentially on damage and porosity: k– ¼ k– o ðϕ ÞexpðbαÞ; (6) where –k o ðϕ Þ is the porosity-dependent permeability tensor for the damage-free material and b is a constant. Since we deal with drained conditions, the permeability dependence in porosity is negligible. Equations (1)–(6) are the standard Biot [1941, 1956] equations for the porous media slightly modified to include nonlinear elasticity (2) and variable permeability (6). The core of the new formulation is coupled kinetic equations for evolving porosity and damage discussed in Appendix A. These equations are derived assuming that the free energy of poroelastic medium, F(T,εij,α,ζ ,ϕ), is a function of temperature, elastic strain, damage, water content, and porosity. dϕ ∂F ∂F ¼ C ϕϕ þ σ M C ϕα dt ∂ϕ ∂α dα ∂F ∂F ¼ C αϕ þ σ M C αα ; dt ∂ϕ ∂α (7) where σ M is the mean stress. Each term in the kinetic equation (7) represents a thermodynamic force associated with energy gradient due to damage (∂F/∂α) and porosity (∂F/∂ϕ) changes. The kinetic coefficients, Cαα, Cϕϕ , Cαϕ , and Cϕα have different expressions and values for different deformation regimes, but general structure of the kinetic equation (7) remains the same. The onset of damage accumulation (∂α/∂t = 0) defines the transition between deformation Regimes I and II (yield cap in Figure 1). Since the free energy of the medium is defined in terms of strains, the yield cap is pffiffiffiffi expressed using volumetric strain, εv = I1, second strain invariant, I2, and the strain invariant ratio, ξ ¼ I1 = I2. This ratio is equivalent to the ratio between differential and normal stress components in the stress-space formulation. Substituting energy function (1) into (7) leads to the following equation for the yield cap: Dðζ ÞðI1 Þnþ1 þ I2 ðξ ξ 0 Þ ¼ 0; (8) where the ratio between two kinetics coefficients, D(ζ ) = Cαt/Cαα, is assumed to be a function of the total porosity (water content), ζ , constant power index n > 1, and critical strain invariant ratio, ξ 0. The shape and size of the yield cap are defined by power index n and coupling coefficient D(ζ ). With these assumptions the shape of the yield cap remains the same (self-similar) during its growth provoked by porosity reduction similar to the DiMaggio and Sandler [1971] model that was shown to be realistic [Baud et al., 2006; Tembe et al., 2007; Wong and Baud, 2012]. For the low porosity, crystalline rocks the coupling coefficient vanishes (D(ζ ) → 0) and equation for the yield cap (8) reduces to ξ ¼ ξ0: (9) This condition defines the straight-dotted line in Figure 1 corresponding to the compaction-dilation transition. Lyakhovsky et al. [1997] proposed this equation for the onset of damage in their model, which ignores porosity effects for brittle crystalline rocks. They connected the critical strain invariant ratio, or modified internal friction, ξ 0, with internal angle of the Byerlee [1978] friction law for rocks. In this study we utilize the condition (9) for the transition between deformation Regimes II and III. Below we specify the kinetic coefficients, Cαα, Cϕϕ , Cαϕ , and Cϕα in (7) for each deformation regime. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 5 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 2.2.1. Regime I Deformation Regime I occurs for the loading conditions beneath the yield cap (8) or for the strain invariant values below a critical value: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 1 þ 4Dðζ Þ ðI1 Þn1 ξ 0 ξ ≤ ξ cr ðI1 ; ζ Þ ¼ ; (10) 2Dðζ Þ ðI1 Þn1 In this regime, material strengthening and compaction take place. The rate of material recovery depends exponentially on the damage and porosity: dα α ϕ exp Dðζ Þ ðI1 Þnþ1 þ I2 ðξ ξ 0 Þ : ¼ H1 exp (11) dt H2 H3 Motivated by the observed logarithmic increase of the static coefficient of friction with time, Lyakhovsky et al. [1997] used an exponential function for the kinetics of material strengthening. The analytical and numerical results of Lyakhovsky et al. [2005] provide quantitative connections between the parameters (H1, H2) and the parameters of rate-and-state-dependent friction. Laboratory experiments [Olsen et al., 1998; Tenthorey et al., 2003] reveal a strong connection between material strengthening, grain crushing, and compaction in high-porosity rocks. To account for this porosity-related enhancement of material recovery, Lyakhovsky and Hamiel [2007] introduced additional term with coefficient H3. The material compaction (porosity reduction) is associated with two processes. The pressure-driven compaction or porosity reduction to certain pressure-dependent equilibrium porosity known as the Athy law [Athy, 1930]: Pe ϕ eq ¼ ϕ f þ A1 exp ; (12) A2 where ϕ f is the final minimum porosity and A1, A2 are material properties defined from the experimentally measured and borehole-observed compaction curves. This process is dominant in Regime I but may continue in other regimes in the case of relatively low confining pressures or very fast loading. Another type of compaction is associated with material recovery. Overall porosity reduction is a superposition of both processes: dϕ α ϕ γ exp Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ: ¼ C ϕ Pe ϕ ϕ eq H1 exp (13) dt H2 H3 K The positive constant Cϕ defines the rate of pressure-driven porosity reduction for (ϕ > ϕ eq) and is equal to zero for lower porosity values. While parameters defining pressure-driven compaction (first term in (13)) could be constrained from hydrostatic load experiments [e.g., Shalev et al., 2014], the constrain for H1, H2,and H3 could hardly be achieved, since direct observations of porosity reduction due to material recovery are not available. 2.2.2. Regime II Beyond the yield cap and beneath the C-D transition (ξ cr(I1, ζ ) ≤ ξ ≤ ξ 0) deformation, Regime II takes place, and damage increases as dα ¼ C d Dðζ Þ ðI1 Þnþ1 þ I2 ðξ ξ 0 Þ ; dt (14) where Cd is the kinetic coefficient of damage accumulation. Lyakhovsky et al. [1997] and Hamiel et al. [2004b] used results of laboratory experiments to constrain the coefficient of damage rate, Cd. They found that constant value of the order of 0.5 to 10 s1 provides a good fit to observed laboratory data for damage increase under moderate and high confining pressures. Lyakhovsky et al. [2005] noted that the rate of damage accumulation is enhanced under low pressures of a few tens of MPa. We keep this in mind when we compare model predictions and laboratory observations in the next section. The porosity reduction (compaction) is expressed by the same functional relation of Regime I (equation (13)), but with a different kinetic coefficient and factor χ > 1 defining the enhanced compaction rate in Regime II. dϕ γ ¼ C ϕ Pe ϕ ϕ eq χ C d Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ: dt K LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. (15) 6 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 Figure 2. Rate of porosity and damage changes. Porosity change is mostly affected by the compaction-dilation (C-D) transition, and damage rate is mostly affected by the yield cap. Note that the yield cap grows with compaction during loading and dynamically changes damage rate. The damage-related compaction is significant under high confining conditions. High-pressure experiments enable constraining χ value. Following the onset of positive damage evolution, the cataclastic flow is modeled as gradual accumulation of irreversible strain components, εirij . Comparisons between theoretical predictions, observed deformation, and acoustic emission from laboratory experiments in granites and sandstones led Hamiel et al. [2004b] to suggest that the accumulation rate of the damage-related irreversible deformation is proportional to the rate of damage increase: dεirij dt ¼ Cv dα τ ij ; dt (16) where Cv is a material constant and τ ij = σ ij σ kkδij/3 is the deviatoric stress tensor. 2.2.3. Regime III In the third regime, beyond the C-D transition (ξ > ξ 0), damage increase is also controlled by equation (14) as in Regime II with the same kinetic coefficient Cd. This process lasts until the final macroscopic failure, stress drop, and a fast conversion of the elastic strain into irreversible plastic strain and porosity. The damage-related porosity change is expressed as ((A20) combined with pressure-driven compaction (A10)) dϕ γ γ ¼ C ϕ Pe ϕ ϕ eq þ C d Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ þ C d 1 ½I2 ðξ ξ 0 Þ2 dt K Pe (17) The kinetics of the compaction-dilation and damage accumulation in Regimes II and III is demonstrated in Figure 2. The compaction rate in Regime II strongly increases with volumetric strain, and it is less sensitive to the differential component (Figure 2a). The rate of damage is roughly proportional to the distance from the yield cap and strongly increases with the differential strain in Regime III (Figure 2b). The distance from the C-D transition, where the critical damage is achieved and macroscopic failure occurs, is a result of competition between rate of the damage accumulation and rate of loading. In the case of fast damage accumulation and slow loading, macroscopic failure occurs just at the C-D transition similar to classical plasticity theory utilizing Coulomb-Mohr failure criterion. Under higher loading rates, failure occurs at larger distances from the yield cap (higher stresses). Such connection between loading rate and measured rock strength was widely observed in rock mechanics experiments [e.g., Paterson and Wong, 2005]. The presented mapping ignores pressure-driven compaction term, since it strongly depends on the initial porosity of the rock sample. Rate of the material recovery and corresponding compaction exponentially depends on both porosity and damage. Therefore, it also cannot be presented in such mapping. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 7 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 During the macroscopic failure at the critical damage value, dynamic Property, Notation, Units Value stress drop occurs leading to fast 3 2,600 Density, ρ (kg/m ) conversion of the elastic into plastic 9 Young modulus, E (Pa) 15 × 10 strain components. We adopt here a Poisson ratio, ν 0.2 9 quasi-static stress drop procedure Biot modulus, M (Pa) 5 × 10 Biot modulus, αB 1 suggested by Lyakhovsky and Ben-Zion Modified internal friction, ξ 0 0.9 [2008], who also discussed the 1 9 9a 0.5 × 10 8 × 10 Damage-related viscosity, Cv (s ) definition of the critical damage value a 1 Damage accumulation, Cd (s ) 5–40 1 17 and its physical meaning. Their Damage decrease, H1 (s ) 2 × 10 procedure of the shear stress Damage decrease, H2 0.02 Damage decrease, H3 0.05 relaxation is modified here to account 1 1 9 1 × 10 Pressure-driven compaction, Cφ, (Pa s ) also for conversion of the volumetric a Enhanced compaction (χ) 1.2–10 strain into plastic components Athy law (porosity change) during plastic flow. ϕf 0.19 (19%) According to the Mises principle, the A1 0.01 (1%) 7 4 × 10 A2, (Pa) plastic flow rules should provide a 2 15 Damage-free permeability, k (m ) 10 maximum dissipation rate of the ¯o Permeability-damage coefficient, b 5 mechanical energy. This may be a These values vary between model runs. achieved by adopting the associated plasticity law of Drucker [1949]. The connection between the Drucker’s associated plasticity law and the Mises principle was widely discussed in context of thermomechanical principles of continua [e.g., Prager, 1959; Mosolov and Myasnikov, 1965, 1981; Ziegler, 1983; Collins and Houlsby, 1997; Collins, 2005]. The associated plasticity law of Drucker [1949] formulated in terms of the incremental stress change implies [e.g., Washizu, 1982] that it is proportional to the gradient of the yield surface formulated in stress-space. We follow the same approach using strain-space formulation with yield surface (8). Table 1. Material Properties and Kinetic Coefficients 2.3. Numerical Implementation The system of equations is solved with Hydro-PED software, which combines two numerical methods: Finite Element Method for the fluid flow (pressure diffusion equation (5)) and Explicit Finite Difference Lagrangian Method for the mechanical equilibrium (equation (1)). The original Hydro-PED software [Shalev et al., 2013; Shalev and Lyakhovsky, 2013] was modified here to account for the coupled damage-porosity evolution. The kinetics equations for every spatial element are solved using explicit time integration. A cylindered domain, with a diameter of 18.4 mm, and length of 40.3 mm, represents the geometry of the laboratory samples used by Zhu and Wong [1997] and Wong et al. [1997]. The numerical domain containing 41,889 tetrahedral elements isused for all simulations. The top and bottom element layers are assigned to not have any damage in order to avoid boundary effect. This makes the effective length of the part of the cylinder to be 38.1 mm where the damage is accumulated and porosity changes. Applying these conditions, we eliminate “conic cracks” propagating from the top and bottom edges of the cylinder. 3. Results In this section first we demonstrate the role of the different processes during loading under intermediate confining pressure. Then we run the model with 20, 50, 170, and 260 MPa confining pressure and 10 MPa fluid pressure and compare calculated stress-strain and porosity-strain diagrams with the experimental results of Wong et al. [1997]. All the model parameters are summarized in the Table 1. Only few of them (χ, Cv, and Cd) were slightly changed between these runs to achieve better fit to laboratory results. In all simulations the domain is initially fractured with small randomly distributed damage, α = 0.2 ± 0.1. We assume that before the onset of axial loading the sample stays enough time under hydrostatic load to approach the equilibrium porosity (equation (12)). The constant strain rate during loading is 5 × 105 s1. The sample porosity is calculated by integrating the total porosity over the entire numerical elements and then normalizing it to the sample volume. Since the boundary condition at the top edge of the cylinder is a constant displacement rate (constant velocity), the axial stress is calculated by integrating vertical stress component over the top boundary and normalizing it to its area. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 8 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 3.1. An Illustrative Deformation Path Initially, the rock is loaded up to target hydrostatic confining pressure that results with only volumetric strain. The sample remainsenough time under these conditions to achieve equilibrium compaction (12). This pressure-driven compaction is not simulated here; the equilibrium porosity is used as the initial porosity of the sample. Thereafter, loading in the axial direction proportionally increase both differential and volumetric strain in the quasi-elastic Regime I (red line, segment 1 in Figure 3a). The slopes of the stress-strain and porosity-strain curves (red lines, segment 1 in Figures 3b and 3c) reflect the elastic moduli of the rocks. During ongoing loading, once strain state reaches the yield cap, deformation regime changes from I to II, damage starts to buildup and promotes compaction (transition from red toblue lines in Figure 3a). The stress-strain (Figure 3b) and porosity-strain (Figure 3c) deviates from their straight directions. The slope of the stress-strain curve decreases because of weakening of the rock (blue segment 2 in Figure 3b). The compaction (blue segment 2 in Figure 3c) is significantly enhanced due to large χ value. With smaller χ values, the compaction curve will practically coincide with the extension of the porosity curve of Regime I. As loading continues, the strain state reaches the C-D transition and deformation regime changes from II to III (green segment 3 in Figure 3a). The rock sample stops to compact and starts to dilate (green segment 3 in Figure 3c). The slope of the stress-strain curve becomes negative (green segment 3 in Figure 3b); this slope is a product of superposition of several processes: material weakening, dilation, and damage-related irreversible strain accumulation. Variation of the damage rate coefficient (Cd) and the damage-related viscosity (Cv) allows fitting the Figure 3. An illustrative deformation path for loading in intersimulated stress-strain curve to the experimental mediate effective pressure: differential strain - volumetric strain observations by changing the shape and size (a), stress - strain (b), and porosity - strain (c). Red line (1) represents of segments 2 and 3. Small Cd values lead to elastic loading. Blue line (2) represents loading in the cataclastic high peak stress, while increased Cv values regime with damage increase that decrease the slope of the flatten the stress-strain curve and extend it to stress-strain curve. Green line (3) represents loading in the brittle regime. Fast loading rate and low damage rate increase the length larger strains. Further loading along with the of the green line. Orange line (4) represents brittle failure with damage accumulation leads one or several dynamic macro-structure stress drop. numerical elements to macroscopic failure conditions; the rock sample experiences simultaneous stress drops in failed elements (yellow segment 4 in Figure 3). Following the stress drop, the postfailure material recovery starts. Because of the strengthening created by the effective porosity-enhanced recovery of the damaged element, the stress may be transferred to surrounding “weaker” elements and distributed failure in the sample volume. This will activate and lock multiple slipping surfaces in a quasi-cyclic LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 9 Axial Stress (MPa) Journal of Geophysical Research: Solid Earth 120 20 100 19.8 80 19.6 B A 60 40 19.2 20 19 0 0 1 2 3 4 A 5 D 6 18.8 7 B C B A 19.4 D C 10.1002/2014JB011805 0 1 2 3 4 C 5 6 7 D 4 2 -1 0 ) X (cm ) X (cm ) X (cm -1 0 0 1 -1 0 Y (cm) 1 -1 0 Y (cm) 20 19 0 1 1 [%] [%] 19 1 -1 0 Y (cm) 1 19 0 1 -1 0 Y (cm) 1 [%] [%] 20.2 1 -1 ) X (cm -1 Z (cm) 3 20.4 19 20.7 Figure 4. Numerical results of triaxial loading of Berea sandstone with confining pressure of 20 MPa (effective pressure 10 MPa) along with laboratory results of Zhu and Wong [1997]. Four snapshots of porosity (A, B) mark the development of dilatant shear bands into a dilation band (C, D). manner. Thereafter, macroscopic sample behavior is cataclastic-like flow together with gradual compaction (black segments in Figures 3b and 3c). The resulting stress-strain curve reflects different processes controlling rock deformation. Among them are (a) the onset of damage at the yield cap, (b) the onset of dilation at the C-D transition, (c) the peak stress controlled by the competition between loading rate, damage accumulation rate, and damage related viscosity, and (d) the softening (minimum stress after peak stress) controlled by the competition between material recovery and loading rates. 3.1.1. Dilatant Shear and Dilation Bands The numerical results (Figure 4) of triaxial loading of Berea sandstone under confining pressure of 20 MPa (effective pressure 10 MPa), reasonably reproduce the laboratory results of Zhu and Wong [1997]. The stress-strain curve follows the experimental data except for some small deviation at initial parts of the loading (initial adjustment) and for some perturbations when the loading was interrupted at the laboratory for the permeability measurements. The simulated porosity curve follows the experimental data until the latest postfailure stage, where the model reproduces the general tendency of the porosity increase. Because the initial effective pressure is low, the pressure-driven compaction is essential during the initial stage of loading and is comparable with the elastic compaction. This is compatible with previous results reported by Shalev et al. [2014], who pointed that pressure-driven compaction is efficient up to ~100 MPa effective pressure. The onset of damage is at ~1% strain when the elastic loading curve is deflected, and the rock undergoes dilation (increase in porosity). The transition from deformation Regimes I to II and then to III is very fast without any meaningful effects associated with damage-enhanced compaction. In the stress-strain and porosity curves (Figure 4), Regime II could hardly be recognized. Porosity increase and shear-enhanced LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 10 Axial Stress (MPa) Journal of Geophysical Research: Solid Earth 300 19 250 18 200 17 16 15 100 0 Axial Stress (MPa) a. 170 MPa 150 14 50 13 0 1 2 3 4 5 6 7 12 300 19 250 18 0 1 2 3 4 5 6 7 b. 260 MPa 17 200 16 150 15 100 14 50 0 Axial Stress (MPa) 10.1002/2014JB011805 13 0 1 2 3 4 5 6 7 12 200 19 150 18.5 100 18 50 17.5 0 0 1 2 3 4 5 6 7 c. 50 MPa 17 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 Figure 5. Numerical and laboratory results of triaxial loading experiments with different confining pressures: 170 MPa, 260 MPa, and 50 MPa. In all cases pore pressure is 10 MPa. dilation overcome the pressure-driven and elastic compaction, resulting with subsequent continuous dilation. Relatively low-peak stress (~100 MPa) is reproduced using fast damage accumulation with Cd = 40 s1, and low damage-related viscosity, Cv = 0.5 109 s1. Such elevated Cd value is expected for low confining conditions [e.g., Lyakhovsky et al., 2005]. Strain localization starts with the onset of damage and is manifested by relatively high angle conjugate shear bands (~45° is the angle between the band itself and the maximum principle stress along the axial direction). Brittle failure along these bands consists of many small failure events which are accompanied with a gradual stress decrease instead significant major stress drop that is common in crystalline rocks (Figure 4). After the peak stress, shear brittle failure tends to localize to a fault plane and dilate, forming oblique dilatant shear bands easily recognized in the snapshot (A) of porosity. Porosity of these bands further increases with ongoing deformation (snapshot B). Efficient material recovery keeps relatively high strength of the band transferring stresses around it that leads to widening of the damage zone. This is reflected by the absence of the LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 11 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 significant macroscopic stress drop typical for low-porosity crystalline rocks. It is also exhibited by relatively high stress and even hardening during the later stages of loading (Figure 4, stages B–D). The deformation bands flatten and merge (snapshot 0.02 C) and form relatively wide horizontal dilatant zone (snapshot D), which is required for observed porosity change. b 3.1.2. Cataclastic Flow 0.01 At high confining pressures, the porous rock a undergoes cataclastic flow with no significant localization. This feature experimentally observed by Wong et al. [1997] for Berea sandstone under 170 and 260 MPa was successfully reproduced by the 0 0.02 0.04 0 model (Figures 5a and 5b). With the onset of damage Elastic Volumetric Strain and transition from Regimes I to II, the slope of the stress-strain diagram abruptly changes due to Figure 6. The growth of the yield cap during the simulation with 260 MPa; (a) initial yield cap and (b) final yield cap. material softening, damage-related accumulation of the irreversible strain, and damage-enhanced compaction (Figures 5a and 5b). Independent observations of the axial strain and porosity reduction allow to distinguish between viscous relaxation and compaction and constrain the model parameters Cv and χ. Cv = 8 × 109 s1 was adopted in both cases, while slightly different damage-enhanced compaction with Differential Strain 0.03 Axial Stress (MPa) 300 18.2 250 18 A 200 17.8 B C 150 C 17.6 D 100 A 17.4 50 0 B D 17.2 0 1 2 3 4 A 5 6 17 7 B 0 1 2 3 4 C 5 6 7 D 4 2 ) X (cm X (cm ) X (cm 0 0 1 -1 0 Y (cm) 1 [%] 18.5 0 1 -1 0 Y (cm) 18 0 1 1 [%] 19.2 -1 0 Y (cm) 1 [%] 19.2 1 -1 ) -1 X (cm -1 ) -1 Z (cm) 3 17 0 1 -1 0 Y (cm) 1 [%] 19.2 13 19.2 Figure 7. Numerical results of loading of Berea sandstone with confining pressure of 100 MPa. Before macroscopic brittle failure the rock deforms homogeneously (A). After brittle failure, localization starts as high-angle shear bands (B, C). With continuous deformation a horizontal compaction band is formed (D). LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 12 Journal of Geophysical Research: Solid Earth 1 0.8 0.6 0.4 0.2 20 15 10 0 1 0.9 200 400 600 800 1000 0.2 0.195 0.19 10.1002/2014JB011805 χ = 1.2 for 170 MPa and χ = 1.3 for 260 MPa provide better fit to porosity. Damage accumulation with Cd = 5 s1 is slow relatively to the loading rate. Only minor amount of distributed damage is accumulated during the sample testing. The loading path in both cases never approaches the C-D transition, and therefore, there is no dilation (path 2 in Figure 1) and the porosity reduces continuously (Figure 5). The functional relation between the coupling coefficient and total porosity, D(ζ), is constrained at the transition to the cataclastic flow to be D(ζ ) = 20 ζ + 4.3. All other points along the loading path in Regime II are beyond the yield cap, providing upper limit for the coupling coefficient. The yield cap grows from the initial size at the transition to the cataclastic flow to its final size at the end of the experiment (~7% axial strain) is shown in Figure 6. 0.8 Laboratory results for confining pressure of 50 MPa (Figure 5c) is similar to the high pressure experiment (Figures 5a 0.7 and 5b). The only minor decrease in the 0.18 differential stress after the peak value is 0.6 reported without any stage which could 0.175 be associated with stress drop. Similarly to high-pressure experiments (Figures 5a 0.5 0.17 and 5b) there is no dilation, and the 725 727 729 731 733 735 porosity reduces continuously (Figure 5c). Such behavior of the Berea sandstone Figure 8. A single-element evolution of damage and porosity. This element under relatively low pressures is not experienced three brittle failure events that resulted with large stress characteristic. Most of the studies [e.g., drops and damage and porosity changes. Immediately after an event, Bernabe and Brace, 1990; Wong and Baud, because of the shear strain changes, there are damage and porosity 2012] observed cataclastic flow only increase (dilation). Subsequently, the element starts to compact once the shear stress is removed and compression still exists. This reduces both above ~150 MPa. Nevertheless, the the porosity and damage. experimental results are fitted here by the model with low rate of the damage accumulation (Cd = 6 s1), high rate of the damage-related irreversible strain accumulation (Cv = 3 × 109 s1) and extremely large damage-enhanced compaction (χ = 10). With these values the loading path remains in Regime II without transition to the brittle behavior. 3.1.3. Shear-Enhanced Compaction and Compaction Bands The modelling of the sample loading in the intermediate effective pressures between brittle failure and ductile behavior (cataclastic flow) is presented in this section using 90 MPa effective pressure and representative model parameters, Cd = 8 s1, Cv = 0.5 × 109 s1, χ = 1. There was no experiment with this pressure in Wong et al. [1997]. The goal of this case is to illustrate the theoretical occurrence of shear enhanced compaction bands that did not occur in the above experiments. Therefore, we chose different rock properties that magnify band formation. In these conditions, the transition from localized shear dilation to homogeneous compaction includes localized compaction [Wong and Baud, 2012]. In this case, there are similarities to both low and high confining pressure end-members, which are best illustrated in the porosity-strain curve (Figure 7). After the initial elastic porosity reduction there is a dilation phase up to point A. The accumulated damage and 0.185 LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 13 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 porosity is distributed over the entire sample (snapshot A in Figure 7). This stage is followed by series of stress drops and more dilation (B). Until this point the behavior is similar to low confining pressure. From this point, compaction dominates and porosity is reduced, similar to the high confining pressure. The spatial distribution of the porosity shows that localization starts with the occurrence of macroscopic failure events (point and snapshot B in Figure 7). At this stage the entire rock is still dilating; however, discrete high angle compaction bands are noticeable. This was defined as shear-enhanced compaction [Baud et al., 2004]. At later stages (points C, D, and corresponding snapshots in Figure 7), the bands become horizontal and may be interpreted as compaction bands. A single element in the cylinder might experience several brittle failure events with significant stress drops followed by fast material recovery and compaction (Figure 8). Once the damage in the element approaches critical value, significant shear dilation increases the porosity. The element fails, releases shear stress (stress drop), and remains with mainly volumetric compression. As a result, the element goes through significant compaction that causes a decrease in both porosity and damage (material recovery). The stress drop in the element occurs under relatively high mean stress above 150 MPa. Large shearing under these pressures changes the internal structure of the material. The grains are crushed, and the material may be compacted to significantly lower porosity values. This more effective pressure-driven compaction is accounted by reducing the final minimal porosity (ϕ f in equation (9)) from 19% to 10%. After the element first fails its porosity gradually decreases to lower values than the prefailure porosity (Figure 8). This assumption is supported by various field and laboratory observations reporting significant reduction in the grain size distribution and porosity values within deformation bands [e.g., Wong and Baud, 2012, and references therein]. Although the end result of this process is shear-enhanced compaction, during its formation, the band goes through significant shear-dilation before it compacts. It is suggested here that the mechanism for the localization of the compaction bands is associated with this initial shear dilation. 4. Discussion Deformation of sandstones, when subjected to compressive loading includes damage accumulation, brittle failure, cataclastic flow, dilation, and compaction. We developed new thermodynamically based damage-porosity rheology model. The model equations account for material weakening-strengthening coupled with compaction-dilation and fluid flow. The model addresses different deformation regimes including elastic deformation and pressure-driven compaction, brittle failure, and cataclastic flow. The modeled yield cap grows with material compaction keeping its shape. Such growth of the yield cap suggested by DiMaggio and Sandler [1971] was shown to be realistic [Baud et al., 2006; Wong and Baud, 2012]. The mathematical formulation is implemented into a three dimensional (3-D) numerical model. The model allows distinguishing the impact of each deformation process which is difficult in laboratory studies, because all these processes occur simultaneously. For example, in the laboratory the total porosity change is measured, and in the model we can separate it to elastic compaction, inelastic pressure-driven compaction, and damage-induced compaction during cataclastic flow. We apply this model to a 3-D cylindered domain and compare the results with triaxial laboratory experiments presented by Zhu and Wong [1997] and Wong et al. [1997]. Matching numerical results with independently measured stress-strain and porosity evolution under wide range of confining pressures provide important information for constraining model parameters. Model parameters controlling brittle deformation and failure are constrained using the low-pressure experiments, whereas high pressure experiments are used to constrain cataclastic flow parameters and the yield cap [Tembe et al., 2007]. The intermediate-pressure experiments are interestingly used to compare the competition between processes. We acknowledge that other combinations of model parameters could give reasonable results. However, we claim that the role of each process will remain proportionally the same. Further modeling compared with laboratory results could provide more constrains on the parameters. Better constrain could be achieved using detailed experimental data including porosity, axial and transversal stress and strain, and acoustic emission measured in experiments with different effective pressures and strain rates. It has been widely accepted that shear loading to brittle failure with dilation, localizes to form a narrow fault plane [e.g., Lockner et al., 1992]. On the other hand, compaction does not require shear stress and traditionally has been described as a diffused homogeneous phenomenon [e.g., Terzaghi, 1941]. The localized compaction bands were found to occur in intermediate effective pressures and suggested to LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 14 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 localize in different mechanisms [Rudnicki, 2007; Das et al., 2013]. Results of this study suggest that the mechanism of compaction localization is the same as brittle shear failure. During the brittle failure, the element dilates and loses shear stress. This process tends to be localized. Upon the release of shear stress, the damaged grains became more compacted (material strengthening) than the prefailure grain packing (Figure 8). This mechanism creates discrete conjugate high-angle shear-enhanced compaction bands (Figure 7b). With increasing loading, the bands connect to form one horizontal compaction band. The three steps in the formation of compaction bands include (1) dilation and subsequent compaction of elements on a localized plane, (2) formation of high-angle shear-enhanced compaction band, and (3) formation of horizontal pure compaction band. Stanchits and Dresen [2003] and Fortin et al. [2006] separated tensile and shear cracks based on acoustic emission and noted that the tensile cracks occur at initial stages of failure. The ability to reproduce the evolution of deformation bands is the main advantage of our model that includes time-dependent formulation, which is not addressed by the bifurcation analysis [e.g., Holcomb et al., 2007]. 5. Conclusions We developed new mathematical formulation and numerical model that accounts for three deformation regimes: (I) quasi-elastic: material strengthening and compaction under the yield cap, (II) cataclastic: damage increase and compaction beyond the yield cap and beneath the compaction-dilation transition, and (III) brittle: damage increase and dilation beyond the compaction-dilation transition. A rock in the brittle regime will eventually experience macroscopic failure and stress drop. Results of numerical simulations demonstrate that at low effective pressures the rock experiences the elastic and brittle regimes where localized shear and dilatant bands may form. At high effective pressures the rock experiences the elastic and cataclastic regimes and deforms homogeneously without strain localization. All deformation features take place at the intermediate effective pressures. This complex behavior results with localization as in the brittle regime and compaction typical for the cataclastic flow. Finally, the shear-enhanced compaction bands and compaction bands are formed. We recognize three stages in the formation of compaction bands: (1) dilation and subsequent compaction of elements on a localized plane, (2) formation of high-angle shear-enhanced compaction band, and (3) formation of horizontal pure compaction band. Appendix A: Constitutive Relations and Kinetics of Damage and Porosity Following Biot’s theory of poroelasticity [Biot, 1941, 1956], the free energy of a poroelastic medium, F, is a sum of the elastic energy, and the poroelastic coupling term of the saturated medium F¼ pffiffiffiffi 1 λðα; ϕ Þ 2 I1 þ μðα; ϕ ÞI2 γðα; ϕ ÞI1 I2 þ M ½αB I1 ðζ ϕ Þ2 ; 2 2 (A1) where εij is the elastic strain tensor, I1 = εkk and I2 = εijεij are the first and second invariants of the strain tensor, ζ is the fluid volume content, ϕ is the material porosity, α is the damage variable, λ, μ, are the Lame drained moduli, γ is the strain coupling modulus M is the reciprocal of the specific storage coefficient at constant strain, and αB is the Biot coefficient. Differentiation of the poroelastic energy (A1) leads to constitutive relations for the stress tensor, σ ij, and fluid pressure, p: ∂F γ ¼ λ I1 δij þ ð2μ γξ Þεij þ αB M½αB I1 ðζ ϕ Þδij ; (A2) σ ij ¼ ∂εij ξ p¼ ∂F ∂F ¼ ¼ M½αB I1 þ ðζ ϕ Þ; ∂ζ ∂ϕ (A3) pffiffiffi pffiffiffi pffiffiffiffi where ξ ¼ I1 = I2 is the strain invariant ratio changing from ξ ¼ 3 for isotropic compaction to ξ ¼ 3 for isotropic dilation. The effect of rock degradation is simulated by making the poroelastic moduli functions of the damage variable. In the general case, the poroelastic moduli are functions of the material porosity, ϕ, and the damage variable, α [Hamiel et al., 2004a]. However, in the simulations presented in this paper, the material porosity change within a few presents (between 18 and 21%) and the direct effect of the porosity on the elastic moduli is negligible. Therefore, the poroelastic moduli μ and γ are assumed to depend only on LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 15 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 thedamage variable α. Following Agnon and Lyakhovsky [1995] and Hamiel et al. [2004a], we use linear approximations for the elastic moduli μ and γ: μðαÞ ¼ μ0 þ ξ 0 γ1 α γðαÞ ¼ γ1 α; (A4) where μ0, γ1, and ξ 0 are constants for each material. The coefficient ξ 0 is related to the friction angle [Agnon and Lyakhovsky, 1995], and the moduli μ0, γ1 are constrained by the conditions for loss of convexity [Lyakhovsky et al., 1997, equations (14) and (15)]. Nonlinear, power law relations between the elastic moduli and the damage are discussed by Hamiel et al. [2004b]. Material porosity and damage evolve with time as a result of the ongoing deformation. Using the balance equations of the energy and entropy, Gibbs relation, fluid mass conservation equation, the definitions of the stress tensor (A2), and fluid pressure (A3), Hamiel et al. [2004a] derived the local entropy production ϕ related to kinetics of porosity and damage (σ M = σ kk/3 is a mean stress): ∂F dα ∂F dϕ þ σM ≥ 0 ϕ¼ ∂α dt ∂ϕ dt (A5) Adopting Onsager [1931] reciprocal relations between thermodynamic forces and fluxes, we write the phenomenological equations for the kinetics of α and ϕ as a set of two coupled differential equations [deGroot and Mazur, 1962; Malvern, 1969]: dϕ ∂F ∂F ¼ C ϕϕ þ σ M C ϕα dt ∂ϕ ∂α dα ∂F ∂F ¼ C αϕ þ σ M C αα ; dt ∂ϕ ∂α (A6) where Cαα, Cϕϕ , Cαϕ , and Cϕα are kinetic coefficients. Each term in the kinetic equations (A6) represents a thermodynamic force; one relates to damage and another to porosity change. The terms with the kinetic coefficients Cαα and Cϕϕ are, respectively, responsible for the evolution of α and ϕ as functions of damage- and porosity-related forces, while the coefficients Cαϕ and Cϕα are responsible for the coupling between the kinetics of damage and porosity. Substituting (A6) back into (A5), using fluid pressure equation (A3) and introducing the effective pressure Pe = σ M p (mean stress minus fluid pressure), the local entropy production (A5) is ϕ ¼ C αα ∂F ∂α 2 ∂F ≥0 þ C ϕϕ P2e þ C αϕ þ C ϕα Pe ∂α (A7) The first term is related to kinetics of damage, the second term is related to porosity change, and the third term appears in the case of coupled damage-porosity evolution. In general, the dissipated mechanical energy (A7) is converted into heat; however, efficiency of such conversion may be less than 100%. Some thermo-mechanical models introduce interstitial working [e.g., Dunn and Serrin, 1985]. Following original study by Taylor and Quinney [1934], in many thermomechanically coupled models the heat power is calculated by using a constant Taylor-Quinney coefficient, which describes the efficiency of converting mechanical work into heat. In this study we ignore thermal effects assuming isothermal conditions and should assure nonnegativity of the local entropy production defining the sign and values of the kinetic coefficients. Nonnegative dissipation without coupling between damage and porosity evolution is guaranteed by positive sign of the diagonal coefficients Cαα and Cϕϕ . The processes associated with coupled damage-porosity evolution (third term in (A7)) may serve as an energy source or sink. For many practical applications it is assumed that these processes eliminate the energy dissipation, which leads to condition (Cαϕ + Cϕα) = 0. However, this mathematical condition adopted by Hamiel et al. [2004a] in their model formulation is not a necessary condition. For example, under compacting stresses (Pe > 0) and moderate shear load corresponding to the deformation Regime II with ξ > ξ 0 and ∂F/∂α > 0, the sum (Cαϕ + Cϕα) may be any positive value. This assumption means that part of energy released due to compaction is transferred into fracture energy driving LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 16 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 the damage increase. Detail discussion of the functional form of the kinetic coefficients Cαα, Cϕϕ , Cαϕ , and Cϕα for different deformation regimes is presented below. Kinetic relations (A6) are the most important equation of the porodamage model controlling the coupled evolution of the state variables. Each term represents different mechanism of porosity and damage evolution and competition between these terms reproduces the different deformation regimes of porous rocks. For simple case neglecting damage-porosity coupling the kinetic equations (A6) are reduced to dϕ ¼ C ϕϕ Pe ; dt (A8a) dα ¼ C αα γ1 I2 ðξ ξ 0 Þ; dt (A8b) The first equation (A8a) accounts for the pressure-driven compaction or porosity reduction to certain pressure-dependent equilibrium porosity known as the Athy law [Athy, 1930]: Pe ϕ eq ¼ ϕ f þ A1 exp ; A2 (A9) ϕ f final minimum porosity; A1, A2 are material properties defined from the experimentally measured and borehole-observed compaction curves. The rate of pressure-driven compaction is not well resolved in the laboratory experiments and cannot be constrained by any field (borehole observation). Here we suggest that the compaction rate coefficient, Cϕϕ , is proportional to the difference between the porosity and the equilibrium porosity: Cϕϕ = Cϕ (ϕ ϕ eq). Finally, kinetics of the compaction under confining pressure is dϕ ¼ dt ( C ϕ Pe ϕ ϕ eq ; 0; for ϕ > ϕ eq for ϕ ≤ ϕ eq ; (A10) Under constant effective pressure this gives exponential compaction: ϕ ðtÞ ¼ ϕ 0 ϕ eq exp C ϕ Pe t þ ϕ eq ; (A11) The second equation (A8a) defines the damage evolution with different kinetic coefficient for the damage increase and decrease 8 C d I2 ðξ ξ 0 Þ; dα < ¼ α ϕ dt : H1 exp exp I 2 ðξ ξ 0 Þ H2 H3 for ξ ≥ ξ0 for ξ < ξ0 (A12) where Cd is the damage rate coefficient, while H1, H2, and H3 are constants describing the rate of material recovery, which depends exponentially on the damage and porosity values. The evolution of damage and porosity described by equations (A10) and (A12) ignores coupling of the damage and porosity evolution. Such coupling comes into consideration by nonzero coefficients Cαϕ and Cϕα in (A6). Competition between porosity-driven damage increase and pressure-driven healing leads to different modes of failure in porous rocks. Following the Hertzian contact theory, Hamiel et al. [2004a] suggested that Cαϕ is a power law expression of the effective pressure, Cαϕ ~ Pen. They demonstrated that the transition from positive to negative values of the slope of the yield cap is a general feature of the model for n > 1. With this form of the kinetic coefficient, their model defines the yield condition in stress-space being function of the effective pressure. Their formulation actually mixed stress- and strain-space formulation for the kinetics of damage. Here we define all conditions completely in the strain-space: 8 > > < C αϕ LYAKHOVSKY ET AL. C d Dðζ Þ ðI1 Þn =K; ¼ α ϕ > > : H1 exp exp Dðζ Þ ðI1 Þn =K; H2 H3 ©2015. American Geophysical Union. All Rights Reserved. dα >0 dt : dα <0 for dt for (A13) 17 Journal of Geophysical Research: Solid Earth where K is a bulk modulus and D(ζ ) is either constant or function of the total water content (total porosity) constrained using strain-defined yield cap (see next section). The Cαϕ coefficient is defined only for compaction (I1 < 0). With this assumption, the kinetic for damage accumulation (dα/dt > 0) is 3 e i ttl br .5 D=6 D= 7 .5 D=7 n=2 D=8 Differential strain (%) a. 2 dα ¼ C d Dðζ ÞðI1 Þnþ1 þ I2 ðξ ξ 0 Þ ; dt 1 0 0 1 2 3 4 5 Volumetric strain (%) 3 n=2.5 e n=1.5 ittl and for material recovery (dα/dt < 0): dα α ϕ ¼ H1 exp exp dt H2 H3 Dðζ ÞðI1 Þnþ1 þ I2 ðξ ξ 0 Þ : Dðζ Þ ð I1 Þn1 ξ 2 þ ξ ξ 0 ¼ 0; 1 0 0 1 2 3 (A14a) (A14b) Damage kinetic equation defines yield surface (onset of damage accumulation) in strain-space. Substituting dα/dt = 0 into (A14a) and (A14b) leads to n=2 br Differential strain (%) b. 2 10.1002/2014JB011805 4 Volumetric strain (%) (A15) which is a quadratic equation for the strain invariant ratio. The solution is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 1 þ 4Dðζ Þ ðI1 Þn1 ξ 0 ξ cr ðI1 ; ζ Þ ¼ ; (A16) 2Dðζ Þ ðI1 Þn1 Figure A1. (a) Self-similar increase of the yield cap with decreasing D value and (b) change of the yield cap shape for different power index n. The plus sign before the radical is chosen to provide the solution ξ cr = ξ 0 for small volumetric strains (I1 → 0) and vanishing coupling term (D → 0). The function D(ζ ) controls the shape of the yield cap and its change from the straight line (ξ cr = ξ 0) for low-porosity crystalline rocks with D = 0 to yield cap with maximal compaction (critical effective pressure in stress-space formulation) corresponding to onset of grain crushing under hydrostatic load. In this case ξ cr ¼ pffiffiffi 3, grain crushing occurs at the volumetric compaction strain equal to pffiffiffi 1 =n 1 3 þ ξ0 ðcrÞ . (A17) I1 ¼ 3Dðζ Þ Equations (A16) and (A17) define the shape and size of the yield cap shown in Figure A1. Colored lines in Figure A1a shows self-similar increase of the yield cap with decreasing of the D value from D = 6.5 to 8.0. The constant power index n = 2 is used in this case. Under low confining conditions, all the yield caps overlap with the C-D transition (dash line for ξ 0 = 0.9). Under high hydrostatic load the yielding occurs at volumetric strain ~3.5% for the high D = 8.0 (red line) and increases to ~4.3% for D = 6.5. These strain values correspond to critical confining pressures for grain crushing between 350 and 430 MPa for porous rocks with bulk modulus of the order of 10 GPa. Figure A1b shows how shape of the yield cap changes for different power index n. For each n the D value were calculated using (A17) using 4% volumetric strain (400 MPa pressure) and same ξ 0 = 0.9. With increased n values yield cap follows the brittle conditions at relatively high loading and the negative slope of the yield cap is significantly steeper toward grain crushing under hydrostatic conditions. The term responsible for the damage-related porosity change in the kinetic relations (A6) with coefficient Cϕα is different for different deformation regimes. Under relatively low shear beneath the yield cap (Regime I) it is difficult to distinguish between pressure-driven compaction (equations (A10) and (A11)) and damage-related compaction. For the simplicity we adopt the condition Cϕα = Cαϕ , and the second term in (A6) becomes dϕ α ϕ γ exp Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ: ¼ H1 exp (A18) dt H2 H3 K LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 18 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 With onset of damage accumulation at the yield cap (Regime II), the damage-related compaction becomes more significant. In order to account for the enhanced compaction, we use Cϕα = χ Cαϕ , 1, χ > 1 and the equation for compaction (A18) becomes dϕ γ ¼ χ C d Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ: dt K (A19) Larger values of the factor χ > 1 mean that compaction in Regime II might be significantly higher than it was predicted by Hamiel et al. [2004a] model, which adopted the condition Cϕα = Cαϕ for all deformation regimes. The sign of the term responsible for the damage-related porosity evolution is changed with transition to the brittle Regime III beyond the C-D transition. Instead of compaction, it provokes the dilation. The maximal dilation rate may be obtained if all the mechanical energy released due to the damage increase (material weakening) is transferred into the pore space surface energy. The balance between damage-related energy source and porosity-related sink leads to the more complex equation for the dilation: dϕ γ γ ¼ C d Dðζ Þ ðI1 Þn 1 I2 ðξ ξ 0 Þ þ C d 1 ½I2 ðξ ξ 0 Þ2 : dt K Pe (A20) The first term in (A20) is almost the same as in (A18), but the additional second term appears due to the energy balance considerations. The equations A18, A19, A20 together with pressure-driven compaction (A10) define evolution of the material porosity along every loading path crossing different deformation regimes. Acknowledgments This work was funded by the Israel Science Foundation (ISF 958/13). Partial support by the U.S. Department of Energy, the Office of Basic Energy Sciences (BES), under award DE-FG-0207ER15916 (UMD) is gratefully acknowledged. In this study, all data are available on request at the following e-mail address: [email protected]. We appreciate the comments from two anonymous reviewers that help to improve the presentation. LYAKHOVSKY ET AL. References Agnon, A., and V. Lyakhovsky (1995), Damage distribution and localization during dyke intrusion, in Physics and Chemistry of Dykes, edited by G. Baer and A. Heimann, pp. 65–78, Balkema, Rotterdam, Netherlands. Antonellini, M., and A. Aydin (1994), Effect of faulting on fluid flow in porous sandstones: Petrophysical properties, Am. Assoc. Petrol. Geol. Bull., 78, 355–377. Antonellini, M. A., A. Aydin, and D. D. Pollard (1994), Microstructure of deformation bands in porous sandstones at Arches National Park, Utah, J. Struct. Geol., 16, 941–959. Athy, L. F. (1930), Density, porosity, and compaction of sedimentary rocks, Am. Assoc. Pet. Geol. Bull., 14, 1–24. Aydin, A., and R. Ahmadov (2009), Bed-parallel compaction bands in aeolian sandstone: Their identification, characterization and implications, Tectonophysics, 479, 277–284. Aydin, A., and A. M. Johnson (1983), Analysis of faulting of porous sandstone, J. Struct. Geol., 5, 19–31. Baud, P., E. Klein, and T.-F. Wong (2004), Compaction localization in porous sandstones: Spatial evolution of damage and acoustic emission activity, J. Struct. Geol., 26, 603–624. Baud, P., V. Vajdova, and T.-F. Wong (2006), Shear-enhanced compaction and strain localization: Inelastic deformation and constitutive modeling of four porous sandstones, J. Geophys. Res., 111, B12401, doi:10.1029/2005JB004101. Baud, P., P. Meredith, and E. Townend (2012), Permeability evolution during triaxial compaction of an anisotropic porous sandstone, J. Geophys. Res., 117, B05203, doi:10.1029/2012JB009176. Bernabe, Y., and W. F. Brace (1990), Deformation and fracture of Berea sandstone, in The Brittle-Ductile Transition in Rocks, Geophys. Monograph., vol. 56, pp. 91–101, AGU, Washington, D. C. Besuelle, P. (2001), Compacting and dilating shear bands in porous rock: Theoretical and experimental conditions, J. Geophys. Res., 106, 13,435–13,442, doi:10.1029/2001JB900011. Biot, M. A. (1941), General theory of three-dimensional consolidation, J. Appl. Phys., 12, 155–164. Biot, M. A. (1956), General solutions of the equations of elasticity and consolidation for a porous material, J. Appl. Mech., 78, 91–96. Biot, M. A. (1973), Nonlinear and semilinear rheology of porous solids, J. Geophys. Res., 78, 4924–4937, doi:10.1029/JB078i023p04924. Byerlee, J. D. (1978), Friction of rocks, Pure Appl. Geophys., 116, 615–626. Cai, Z., and D. Bercovici (2014), Two-phase visco-elastic damage theory, with applications to subsurface fluid injection, Geophys. J. Int., 199(3), 1481–1496. Carroll, M. M. (1991), A critical state plasticity theory for porous reservoir rock, Am. Soc. Mech. Eng., 117, Book G00617-1991. Casey, J., and P. M. Naghdi (1983), On the nonequivalence of the stress space and strain space formulations of plasticity theory, J. Appl. Mech., 50, 350–354. Charalampidou, E. M., S. A. Hall, S. Stanchits, G. Viggiani, and H. Lewis (2014), Shear-enhanced compaction band identification at the laboratory scale using acoustic and full-field methods, Int. J. Rock Mech. Min. Sci., 67, 240–252. Chemenda, A. I. (2011), Origin of compaction bands: Anti-cracking or constitutive instability?, Tectonophysics, 499, 156–164. Chemenda, A. I., C. Wibberley, and E. Saillet (2012), Evolution of compactive shear deformation bands: Numerical models and geological data, Tectonophysics, 526–529, 56–66. Collins, I. F. (2005), The concept of stored plastic work or frozen elastic energy in soil mechanics, Geotechnique, 55, 373–382. Collins, I. F., and G. T. Houlsby (1997), Application of thermo-mechanical principles to the modelling of geotechnical materials, Proc. R. Soc. London, Ser. A, 453(1964), 1975–2001. Coussy, O. (1995), Mechanics of Porous Continua, Wiley, Chichester, U. K. Das, A., G. D. Nguyen, and I. Einav (2013), The propagation of compaction bands in porous rocks based on breakage mechanics, J. Geophys. Res. Solid Earth, 118, 2049–2066, doi:10.1002/jgrb.50193. deGroot, S. R., and P. Mazur (1962), Non-Equilibrium Thermodynamics, North-Holland Co., Amsterdam. DiMaggio, F. L., and I. S. Sandler (1971), Material model for granular soils, J. Eng. Mech. Div. ASCE, 97, 935–950. ©2015. American Geophysical Union. All Rights Reserved. 19 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 Drucker, D. C. (1949), Some implications of work-hardening and ideal plasticity, Quart. Appl. Math., 7, 411–418. Du Bernard, X., P. Eichhubl, and A. Aydin (2002), Dilation bands: A new form of localized failure in granular media, Geophys. Res. Lett., 29(24), 2176, doi:10.1029/2002GL015966. Dunn, J. E., and J. Serrin (1985), On the thermodynamics on interstitial working, Arch. Ration. Mech. Anal., 88, 95–133. Eichhubl, P., J. N. Hooker, and S. E. Laubach (2010), Pure and shear-enhanced compaction bands in Aztec Sandstone, J. Struct. Geol., 32, 1873–1886. Einav, I. (2004), Thermo-mechanical relations between stress-space and strain-space models, Geotechnique, 54(5), 315–318. Einav, I. (2005a), Energy and variational principles for piles in dissipative soil, Geotechnique, 55, 515–525. Einav, I. (2005b), A second look at strain space plasticity and latest applications, 18th Australasian Conference on the Mechanics of Structures and Materials, Perth, (1), 225–231. Einav, I. (2007a), Breakage mechanics—Part I: Theory, J. Mech. Phys. Solid., 55, 1274–1297. Einav, I. (2007b), Breakage mechanics—Part II: Modeling granular materials, J. Mech. Phys. Solid., 55, 1298–1320. Fortin, J., S. Stanchits, G. Dresen, and Y. Guéguen (2006), Acoustic emission and velocities associated with the formation of compaction bands in sandstone, J. Geophys. Res., 111, B10203, doi:10.1029/2005JB003854. Fossen, H., R. A. Schultz, Z. K. Zhipton, and K. Mair (2007), Deformation bands in sandstone: A review, J. Geol. Soc. London, 164, 755–769. Fossen, H., R. A. Schultz, and A. Torabi (2011), Conditions and implications for compaction band formation in the Navajo Sandstone, Utah, J. Struct. Geol., 33, 1477–1490. Foster, C. D., R. A. Regueiro, A. F. Fossum, and R. I. Borja (2005), Implicit numerical integration of a three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials, Comput. Meth. Appl. Mech. Eng., 194(50), 5109–5138. Fredrich, J. T., B. Evans, and T.-f. Wong (1989), Micromechanics of the brittle to plastic transition in Carrara Marble, J. Geophys. Res., 94, 4129–4145, doi:10.1029/JB094iB04p04129. Grueschow, E., and J. W. Rudnicki (2005), Elliptic yield cap constitutive modeling for high porosity sandstone, Int. J. Solids Struct., 42(16), 4574–4587. Hamiel, Y., V. Lyakhovsky, and A. Agnon (2004a), Coupled evolution of damage and porosity in poroelastic media: Theory and applications to deformation of porous rocks, Geophys. J. Int., 156, 701–713. Hamiel, Y., Y. Liu, V. Lyakhovsky, Y. Ben-Zion, and D. Lockner (2004b), A viscoelastic damage model with applications to stable and unstable fracturing, Geophys. J. Int., 159, 1155–1165. Hamiel, Y., V. Lyakhovsky, and A. Agnon (2005), Poroelastic damage rheology: Dilation, compaction, and failure of rocks. Geochem. Geophys. Geosyst., 6, Q01008, doi:10.1029/2004GC000813. Han, D. J., and W. F. Chen (1986), Strain-space plasticity formulation for hardening-softening materials with elastoplastic coupling, Int. J. Solids Struct., 22, 935–950. Holcomb, D., J. W. Rudnicki, K. A. Issen, and K. Sternlof (2007), Compaction localization in the Earth and the laboratory: State of the research and research directions, Acta Geotech., 2(1), 1–15. Issen, K. A., and J. W. Rudnicki (2000), Conditions for compaction bands in porous rock, J. Geophys. Res., 105, 21,529–21,536, doi:10.1029/ 2000JB900185. Kachanov, L. M. (1958), On the time to rupture under creep conditions [in Russian], Izv. Acad. Nauk SSSR OTN 8, 26–31. Karrech, A., C. E. Schrank, R. Freij-Ayoub, and K. Regenauer-Lieb (2014), A multi-scaling approach to predict hydraulic damage of poromaterials, Int. J. Mech. Sci., 78, 1–7. Lehane, B. M., and B. Simpson (2000), Modelling glacial till under triaxial conditions using a BRICK soil model, Can. Geotech. J., 37, 1078–1088. Lockner, D. A., J. D. Byerlee, V. Kuksenko, A. Ponomarev, and A. Sidorin (1992), Observations of quasi-static fault growth from acoustic emissions, in Fault Mechanics and Transport Properties of Rocks, edited by B. Evans and T.-F. Wong, pp. 3–32, Academic Press, San Diego, Calif. Lyakhovsky, V., and Y. Ben-Zion (2008), Scaling relations of earthquakes and aseismic deformation in a damage rheology model, Geophys. J. Int., 172, 651–662. Lyakhovsky, V., and Y. Hamiel (2007), Damage evolution and fluid flow in poroelastic rock, Izv. Phys. Solid Earth, 43, 13–23. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon (1997), Distributed damage, faulting, and friction, J. Geophys. Res., 102, 27,635–27,649, doi:10.1029/ 97JB01896. Lyakhovsky, V., Y. Ben-Zion, and A. Agnon (2005), A viscoelastic damage rheology and rate- and state-dependent friction, Geophys. J. Int., 161, 179–190. Malvern, L. E. (1969), Introduction to the Mechanics of a Continuum Medium, Prentice-Hall, Inc., Englewood Cliffs, N. J. Mosolov, P. P., and V. P. Myasnikov (1965), Variational methods in flow theory of visco-plastic media [in Russian], Appl. Math. Mech., 29, 468–492. Mosolov, P. P., and V. P. Myasnikov (1981), Mechanics of Rigid Plastic Media [in Russian], Nauka, Moscow. Naghdi, P. M., and J. A. Trapp (1975), The significance of formulating plasticity theory with reference to loading surfaces in strain space, Int. J. Eng. Sci., 13, 785–797. Olsen, M. P., C. H. Scholz, and A. Leger (1998), Healing and sealing of a simulated fault gouge under hydrothermal conditions: Implications for fault healing, J. Geophys. Res., 103, 7421–7430, doi:10.1029/97JB03402. Olsson, W. A. (1999), Theoretical and experimental investigation of compaction bands in porous rock, J. Geophys. Res., 104, 7219–7228, doi:10.1029/1998JB900120. Onsager, L. (1931), Reciprocal relations in irreversible processes, Phys. Rev., 37, 405–416. Paterson, M. S., and T.-F. Wong (2005), Experimental Rock Deformation—The Brittle Field, 2nd ed., 347 pp., Springer, Berlin. Prager, W. (1959), An Introduction to Plasticity, Addison-Wesley, Amsterdam. Puzrin, A. M., and G. T. Houlsby (2001), Fundamentals of kinematic hardening hyperplasticity, Int. J. Solids Struct., 38, 3771–3794. Ricard, Y., and D. Bercovici (2009), A continuum theory of grain size evolution and damage, J. Geophys. Res., 114, B01204, doi:10.1029/ 2007JB005491. Rudnicki, J. W. (2007), Models for compaction band propagation, in Rock Physics and Geomechanics in the Study of Reservoirs and Repositories, edited by C. David and M. LeRavalec, Geol. Soc. London Spec. Publ., 284, 107–126. Rudnicki, J. W., and J. R. Rice (1975), Conditions for the localization of deformation in pressure-sensitive dilatant materials, J. Mech. Phys. Solid., 23, 371–394. Shalev, E., and V. Lyakhovsky (2013), The processes controlling damage zone propagation induced by wellbore fluid injection, Geophys. J. Int., 193, 209–219. Shalev, E., M. Calò, and V. Lyakhovsky (2013), Formation of damage zone and seismic velocity variations during hydraulic stimulation: Numerical modeling and field observations, Geophys. J. Int., 195, 1023–1033. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 20 Journal of Geophysical Research: Solid Earth 10.1002/2014JB011805 Shalev, E., V. Lyakhovsky, A. Ougier-Simonin, Y. Hamiel, and W. Zhu (2014), Inelastic compaction, dilation and hysteresis of sandstones under hydrostatic conditions, Geophys. J. Int., 197, 920–925, doi:10.1093/gji/ggu052. Stanchits, S., and G. Dresen (2003), Separation of tensile and shear cracks based on acoustic emission analysis of rock fracture, paper presented at the symposium Non-Destructive Testing in Civil Engineering 2003 (NDT-CE), Deutche Gesellschaft fur Zerstorungsfreie Prufugn e.V. (DGZfP), Berlin. Sternlof, K. R., J. W. Rudnicki, and D. D. Pollard (2005), Anticrack-inclusion model for compaction bands in sandstone, J. Geophys. Res. 110, B11403, doi:10.1029/2005JB003764. Sternlof, K. R., M. Karimi-Fard, D. D. Pollard, and L. J. Durlofsky (2006), Flow and transport effects of compaction bands in sandstone at scales relevant to aquifer and reservoir management, Water Resour. Res., 42, W07425, doi:10.1029/2005WR004664. Sun, W., J. E. Andrade, J. W. Rudnicki, and P. Eichhubl (2011), Connecting microstructural attributes and permeability from 3D tomographic images of in situ shear-enhanced compaction bands using multiscale computations, Geophys. Res. Lett., 38, L10302, doi:10.1029/ 2011GL047683. Sun, W., Q. Chen, and J. T. Ostien (2014), Modeling the hydro-mechanical responses of strip and circular punch loadings on water-saturated collapsible geomaterials, Acta Geotech., 9(5), 903–934. Taylor, G. I., and H. Quinney (1934), The latent energy remaining in a metal after cold working, Proc. R. Soc. London, Ser. A, 143, 307–326. Tembe, S., V. Vajdova, P. Baud, W. Zhu, and T.-f. Wong (2007), A new methodology to delineate the compactive yield cap of two porous sandstones under undrained condition, Mech. Mater., 39, 513–523. Tenthorey, E., S. F. Cox, and H. F. Todd (2003), Evolution of strength recovery and permeability during fluid-rock reaction in experimental fault zones, Earth Planet. Sci. Lett., 206, 161–172. Terzaghi, K. (1941), Theoretical Soil Mechanics, Wiley, New York. Wang, H. F. (2000), Theory of Linear Poroelasticity With Applications to Geomechanics and Hydrogeology, Princeton Univ. Press, Princeton, N. J. Washizu, K. (1982), Variational Methods in Elasticity and Plasticity, 3rd ed., 630 pp., Pergamon Press Ltd, Oxford. Wong, T.-F., and P. Baud (2012), The brittle transition in rocks: A review, J. Struct. Geol., 44, 25–53. Wong, T.-F., C. David, and W. Zhu (1997), The transition from brittle faulting to cataclastic flow in porous sandstones: Mechanical deformation, J. Geophys. Res., 102, 3009–3025, doi:10.1029/96JB03281. Yoder, P. J., and W. D. Iwan (1981), On the formulation of strain space plasticity with multiple loading surfaces, J. Appl. Mech., 48, 773–778. Zhu, W., and J. Walsh (2006), A new model for analyzing the effect of fractures on triaxial deformation, Int. J. Rock Mech. Min. Sci., 43, 1241–1255. Zhu, W., and T.-F. Wong (1997), The transition from brittle faulting to cataclastic flow: Permeability evolution, J. Geophys. Res., 102(B2), 3027–3041, doi:10.1029/96JB03282. Zhu, W., and T.-F. Wong (1999), Network modeling of the evolution of permeability and dilatancy in compact rock, J. Geophys. Res., 104(2), 2963–2971, doi:10.1029/1998JB900062. Ziegler, H. (1983), An Introduction to Thermomechanics, 2nd ed., 358 pp., North-Holland Co., Amsterdam. LYAKHOVSKY ET AL. ©2015. American Geophysical Union. All Rights Reserved. 21
© Copyright 2024