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 ﬂow, 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 ﬂow
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 conﬁning 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 conﬁning 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 ﬂow 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 ﬂow. 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 ﬂow 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 ﬂow [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 signiﬁcant porosity reduction with
concomitant permeability decrease is generally observed in the cataclastic ﬂow 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, conﬁning 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 ﬁeld 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 ﬁeld [Sternlof et al., 2006; Aydin and Ahmadov, 2009; Fossen et al., 2011]. Different types of
deformation bands may either act as ﬂuid conduits or ﬂow barriers [e.g., Zhu and Wong, 1997; Sun et al., 2011;
Baud et al., 2012]. Understanding ﬂuid ﬂow 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 ﬂow 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 ﬂow 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 signiﬁcant 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 deﬁne the brittle failure and by a yield cap criterion that deﬁnes cataclastic
ﬂow [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 signiﬁcant 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 ﬂaws [e.g., Han and Chen, 1986]. As a result, the application of the stress-space formulation encounters
difﬁculties 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
signiﬁcant 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 deﬁned 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 ﬂow (Δεd). The amount of
the growth of the yield cap (f) is proportional to compaction rate.
10.1002/2014JB011805
signiﬁcant 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 difﬁculties 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 ﬁnite
elements are formulated in terms of
displacements (strains), and the
description of the yield cap in the
stress-space is not efﬁcient.
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 deﬁned by two envelopes (Figure 1).
Heavy line schematically deﬁnes 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 deﬁne 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 conﬁning pressures usually treated as cataclastic ﬂow
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 signiﬁcant 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 (conﬁning-pore) pressure (loading path 1 in Figure 1), the
strain state crosses the C-D transition, which overlaps with the yield cap at low conﬁnement, 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 ﬁeld 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 conﬁnement (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 ﬂow reduces
the differential strain (Maxwell relaxation) and also shifts the loading path. Compaction and cataclastic ﬂow
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 ﬂow. The rates of compaction and
cataclastic ﬂow depend on the distance between the strain state and the yield cap. If initially, compaction and
cataclastic ﬂow 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 ﬂuid
ﬂow 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¼
pﬃﬃﬃﬃ 1
λðα; ϕ Þ 2
I1 þ μðα; ϕ ÞI2 γðα; ϕ ÞI1 I2 þ M ½αB I1 ðζ ϕ Þ2 ;
2
2
(1)
where I1 = εkk and I2 = εijεij are the ﬁrst and second invariants of the
elastic strain
tensor, deﬁned 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 ﬂuid 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 speciﬁc storage coefﬁcient at constant strain, and αB is the Biot coefﬁcient.
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 ﬁrst 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 speciﬁc storage coefﬁcient at constant strain. The ﬂuid mass conservation
equation combined with the Darcy law leads to pressure-diffusion equation for saturated isothermal ﬂow
[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 ﬂuid viscosity, ρ is the ﬂuid 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 modiﬁed 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 coefﬁcients, 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) deﬁnes the transition between deformation Regimes I and II
(yield cap in Figure 1). Since the free energy of the medium is deﬁned in terms of strains, the yield cap is
pﬃﬃﬃﬃ
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 coefﬁcients, 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 deﬁned by power index n and coupling coefﬁcient 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 coefﬁcient vanishes (D(ζ ) → 0)
and equation for the yield cap (8) reduces to
ξ ¼ ξ0:
(9)
This condition deﬁnes 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
modiﬁed 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 coefﬁcients, 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:
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
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 coefﬁcient 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 coefﬁcient 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 ﬁnal minimum porosity and A1, A2 are material properties deﬁned 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 conﬁning 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ϕ deﬁnes the rate of pressure-driven porosity reduction for (ϕ > ϕ eq) and is equal to
zero for lower porosity values. While parameters deﬁning pressure-driven compaction (ﬁrst 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 coefﬁcient of damage accumulation. Lyakhovsky et al. [1997] and Hamiel et al. [2004b]
used results of laboratory experiments to constrain the coefﬁcient of damage rate, Cd. They found that
constant value of the order of 0.5 to 10 s1 provides a good ﬁt to observed laboratory data for damage
increase under moderate and high conﬁning 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 coefﬁcient and factor χ > 1 deﬁning 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 signiﬁcant under high conﬁning conditions. High-pressure experiments
enable constraining χ value.
Following the onset of positive damage evolution, the cataclastic ﬂow 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 coefﬁcient Cd. This process lasts until the ﬁnal 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
Modiﬁed internal friction, ξ 0
0.9
[2008], who also discussed the
1
9
9a
0.5 × 10 8 × 10
Damage-related viscosity, Cv (s )
deﬁnition 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 modiﬁed 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 ﬂow.
ϕf
0.19 (19%)
According to the Mises principle, the
A1
0.01 (1%)
7
4 × 10
A2, (Pa)
plastic ﬂow rules should provide a
2
15
Damage-free permeability, k (m )
10
maximum dissipation rate of the
¯o
Permeability-damage coefﬁcient, 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 Coefﬁcients
2.3. Numerical Implementation
The system of equations is solved with Hydro-PED software, which combines two numerical methods: Finite
Element Method for the ﬂuid ﬂow (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 modiﬁed 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 ﬁrst we demonstrate the role of the different processes during loading under intermediate
conﬁning pressure. Then we run the model with 20, 50, 170, and 260 MPa conﬁning pressure and 10 MPa ﬂuid
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 ﬁt 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 conﬁning 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) reﬂect 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 signiﬁcantly 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 coefﬁcient (Cd) and the
damage-related viscosity (Cv) allows ﬁtting 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
ﬂatten
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 conﬁning 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 ﬂow together with gradual compaction
(black segments in Figures 3b and 3c).
The resulting stress-strain curve reﬂects 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 conﬁning 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 efﬁcient up to ~100 MPa effective
pressure. The onset of damage is at ~1% strain when the elastic loading curve is deﬂected, 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 conﬁning 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 conﬁning 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 signiﬁcant 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). Efﬁcient material recovery keeps relatively high strength of the band
transferring stresses around it that leads to widening of the damage zone. This is reﬂected 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
signiﬁcant 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 ﬂatten 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 conﬁning pressures, the porous rock
a
undergoes cataclastic ﬂow with no signiﬁcant
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) ﬁnal 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 conﬁning 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 ﬁt 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 coefﬁcient
and total porosity, D(ζ), is constrained
at the transition to the cataclastic ﬂow 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 coefﬁcient. The
yield cap grows from the initial size at the
transition to the cataclastic ﬂow to its
ﬁnal size at the end of the experiment
(~7% axial strain) is shown in Figure 6.
0.8
Laboratory results for conﬁning 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 ﬂow 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 ﬁtted 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 ﬂow) 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 conﬁning 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 conﬁning pressure. From this
point, compaction dominates and porosity is reduced, similar to the high conﬁning 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 deﬁned 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 signiﬁcant stress drops
followed by fast material recovery and compaction (Figure 8). Once the damage in the element approaches
critical value, signiﬁcant 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 signiﬁcant
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
signiﬁcantly lower porosity values. This more effective pressure-driven compaction is accounted by reducing
the ﬁnal minimal porosity (ϕ f in equation (9)) from 19% to 10%. After the element ﬁrst fails its porosity gradually
decreases to lower values than the prefailure porosity (Figure 8). This assumption is supported by various
ﬁeld and laboratory observations reporting signiﬁcant 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 signiﬁcant
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 ﬂow, 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 ﬂuid ﬂow. The model addresses different deformation regimes including
elastic deformation and pressure-driven compaction, brittle failure, and cataclastic ﬂow. 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 difﬁcult 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 ﬂow. 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
conﬁning 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 ﬂow 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 ﬂow. 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¼
pﬃﬃﬃﬃ 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 ﬁrst and second invariants of the strain tensor, ζ
is the ﬂuid 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 speciﬁc storage coefﬁcient at constant
strain, and αB is the Biot coefﬁcient.
Differentiation of the poroelastic energy (A1) leads to constitutive relations for the stress tensor, σ ij, and ﬂuid
pressure, p:
∂F
γ
¼ λ I1 δij þ ð2μ γξ Þεij þ αB M½αB I1 ðζ ϕ Þδij ;
(A2)
σ ij ¼
∂εij
ξ
p¼
∂F
∂F
¼
¼ M½αB I1 þ ðζ ϕ Þ;
∂ζ
∂ϕ
(A3)
pﬃﬃﬃ
pﬃﬃﬃ
pﬃﬃﬃﬃ
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 coefﬁcient ξ 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, ﬂuid mass conservation equation, the deﬁnitions of the
stress tensor (A2), and ﬂuid 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 ﬂuxes, 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 coefﬁcients. 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 coefﬁcients Cαα and Cϕϕ are, respectively, responsible for the evolution of α and ϕ as functions of
damage- and porosity-related forces, while the coefﬁcients Cαϕ and Cϕα are responsible for the coupling
between the kinetics of damage and porosity. Substituting (A6) back into (A5), using ﬂuid pressure equation
(A3) and introducing the effective pressure Pe = σ M p (mean stress minus ﬂuid pressure), the local entropy
production (A5) is
ϕ ¼ C αα
∂F
∂α
2
∂F
≥0
þ C ϕϕ P2e þ C αϕ þ C ϕα Pe
∂α
(A7)
The ﬁrst 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, efﬁciency 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 coefﬁcient, which describes the efﬁciency 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 deﬁning the sign and values of the kinetic coefﬁcients.
Nonnegative dissipation without coupling between damage and porosity evolution is guaranteed by positive
sign of the diagonal coefﬁcients 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 coefﬁcients 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 ﬁrst 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 ﬁnal minimum porosity; A1, A2 are material properties deﬁned 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 ﬁeld (borehole observation). Here we suggest
that the compaction rate coefﬁcient, Cϕϕ , is proportional to the difference between the porosity and the
equilibrium porosity: Cϕϕ = Cϕ (ϕ ϕ eq). Finally, kinetics of the compaction under conﬁning 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) deﬁnes the damage evolution with different kinetic coefﬁcient 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 coefﬁcient, 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 coefﬁcients 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 coefﬁcient, their model deﬁnes 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 deﬁne 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-deﬁned yield cap
(see next section). The Cαϕ coefﬁcient is deﬁned
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 deﬁnes 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
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
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 ¼
pﬃﬃﬃ
3, grain crushing occurs at the volumetric compaction strain equal to
pﬃﬃﬃ
1 =n 1
3 þ ξ0
ðcrÞ
.
(A17)
I1 ¼ 3Dðζ Þ
Equations (A16) and (A17) deﬁne 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 conﬁning 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 conﬁning 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 signiﬁcantly steeper toward grain crushing under hydrostatic conditions.
The term responsible for the damage-related porosity change in the kinetic relations (A6) with coefﬁcient Cϕα is
different for different deformation regimes. Under relatively low shear beneath the yield cap (Regime I) it is
difﬁcult 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 signiﬁcant. 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 signiﬁcantly 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 ﬁrst 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)
deﬁne 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 Ofﬁce 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 ﬂuid ﬂow 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 identiﬁcation, 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 ﬂuid 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 identiﬁcation at the
laboratory scale using acoustic and full-ﬁeld 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 ﬂuid ﬂow 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 ﬂow 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 signiﬁcance 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 ﬂuid 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 ﬁeld 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 ﬂuid-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 ﬂow 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 ﬂow: 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