# Micromechanics-Based Permeability Evolution in Brittle Materials at High Strain Rates

Thibaut Perol, Harsha S. Bhat

► To cite this version:

Thibaut Perol, Harsha S. Bhat. Micromechanics-Based Permeability Evolution in Brittle Materials at High Strain Rates. Pure and Applied Geophysics, 2016, 173, pp.2857-2868. 10.1007/s00024-016-1354-4. insu-03581305

HAL Id: insu-03581305

https://insu.hal.science/insu-03581305

Submitted on 19 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Micromechanics-Based Permeability Evolution in Brittle Materials at High Strain Rates

THIBAUT PEROL<fnref n="1" /> and HARSHA S. BHAT<fnref n="2" />,<fnref n="3" />

*Abstract*—We develop a micromechanics-based permeability evolution model for brittle materials at high strain rates (≥ 100 s⁻<fnref n="1" />). Extending for undrained deformation the mechanical constitutive description of brittle solids, whose constitutive response is governed by micro-cracks, we now relate the damage-induced strains to micro-crack aperture. We then use an existing permeability model to evaluate the permeability evolution. This model predicts both the percolative and connected regime of permeability evolution of Westerly Granite during triaxial loading at high strain rate. This model can simulate pore pressure history during earthquake coseismic dynamic ruptures under undrained conditions.

**Key words:** Micromechanics, permeability, poroelasticity.

## 1. Introduction

Many geomechanical phenomena, both natural and anthropogenic, involve a strong coupling between mechanical, hydraulic, thermal and chemical processes in the host medium and occur at various lengths and time scales with respect to the said phenomenon. For example, dynamic earthquake ruptures hosted on faults result in high strain rate fracturing of the surrounding medium leading to enhanced permeability and rapid fluid discharge resulting in reduction of pore fluid pressure (Sibson 1994). Simultaneously, rapid heating of the sliding interface will lead to temperature rise in the surrounding medium and due to mismatch in the thermal expansion of the host rock and fluids, the so-called

thermal pressurization mechanism, will lead to increase in pore pressure acting on the fault surface. This will in turn rapidly reduce the frictional resistance, via Terzaghi effect, of the fault (Rice 2006). Clearly one expects a competition between pore pressure loss due to permeability enhancement and gain due to thermal pressurization. In addition to this (Miller et al. 2004) also raised the possibility that aftershocks in large earthquakes, in fluid-rich geologic environment, may be driven by high-pressure fluids propagating through fracture damaged zones created by main shocks. They showed that the 1997 Umbria-Marche seismic sequence propagated at speeds of the order of 1 km/day, which is faster than diffusion timescale. It was shown by Rice (1992) that soliton-like solutions may exist for fluid flow problems that show strong nonlinear dependence on hydraulic and mechanical properties. Murphy et al. (2004) showed that flow in fractures with pressure-dependent aperture results in a pressure field expanding in a shock-wave manner rather than in a smooth diffusive way.

Zoback and Gorelick (2012) and Cox (2016) argue that due to the critically stressed nature of the crust, fluid injection in deep wells can trigger earthquakes when the injection increases pore pressure in the vicinity of pre-existing potentially active faults. This effect was first documented by Healy et al. (1968) in the 1960s in Denver, Colorado, when injection into a 3-km-deep well at the nearby Rocky Mountain Arsenal triggered earthquakes and most recently by Keranen et al. (2013) potentially linking wastewater injection and the 2011 M<sub>w</sub> 5.7 earthquake sequence in Oklahoma, USA. Fluid injection could be associated with waste water injection in shale gas extraction, large-scale geological sequestration of CO<sub>2</sub> or in enhanced geothermal systems. Not just injection but also rapid extraction of groundwater

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA, USA. E-mail: tperol@seas.harvard.edu</footnote>

<footnote n="2">² Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, UMR 7154 CNRS, 75005 Paris, France.</footnote>

<footnote n="3">³ *Present Address:* Laboratoire de Géologie, CNRS UMR 8538, École Normale Supérieure, 75005 Paris, France.</footnote>
from a basin aquifer can lead to potentially damaging earthquakes as shown by González et al. (2012).

Permeability evolution in almost all of the above situations is expected to be non-negligible, if not dramatic. Mitchell and Faulkner (2008) show at least two orders magnitude of increase in permeability in crystalline rocks under triaxial loading conditions at low strain rates. They also noticed that initially permeability and pore volume decreased with loading before permeability started to increase. Evans et al. (2005) and Evans (2005) showed permeability increase of two orders of magnitude associated with reservoir stimulation exercise conducted at Soultz-sous-Forêt. Both the observations above were associated with low loading rates. Morton et al. (2012) conducted split Hopkinson bar experiments to induce high strain rate damage, associated with earthquake ruptures, in crystalline rocks and found permeability jump up to seven orders of magnitude.

All of the above subset of observations, for a variety of geomechanical conditions, highlight that the nonlinear hydro-mechanical coupling cannot be ignored when modeling such phenomenon and call for a fully coupled hydro-mechanical model that takes into account changes in mechanical and hydraulic properties of fractured media and the interplay between them. Typically this coupling is done by relating permeability to confining stress or to Terzaghi effective stress. Other variants exist like Cappa and Rutqvist (2011a, b) who developed a permeability model in which permeability changes with volumetric plastic strain and porosity in a Mohr–Coulomb elastic–plastic medium. As far as the authors are aware no micromechanical model exists that attempts to combine the strain rate-sensitive nature of brittle materials at high strain rates and the corresponding evolution of both mechanical and hydraulic properties. In most brittle materials, such as rocks and ceramics, failure occurs by growth of micro-cracks, that either nucleate or propagate from pre-existing ones, which coalesce leading to macroscopic failure. Bhat et al. (2012) showed that rate sensitivity of micro-crack growth at dynamic loading rates translates into the macroscopic strain rate sensitivity of compressive failure strength. Tensile cracks induced by grain-boundary sliding and opening, commonly characterized as damage, are the

principal source of inelastic deformation in these materials. As these cracks grow the hydraulic connectivity increases and make the rock more permeable. Understanding the hydromechanical coupling taking place during mechanism leading to failure, e.g., hydraulic fracturing, has a significant impact when trying to forecast the dynamic evolution of oil reservoirs and underground waste storage. Quasi-static models of permeability in fractured media are physically insufficient to describe the evolution of conductivity of the network during these processes.

We develop a coupled hydro-mechanical description of a brittle material by modeling the interplay between constitutive response (stress–strain behavior) and permeability evolution both due to external loads and internal configurational changes (evolution of the fracture network and pore fluid pressure). For this we use a micromechanics-based model of brittle failure developed and validated against experimental results for both quasi-static (Ashby and Sammis 1990) and dynamic strain rate regimes (Deshpande and Evans 2008; Bhat et al. 2012). From the mechanical constitutive relationships we derive the expressions for the aperture of the micro-cracks and use the permeability relation given by Dienes (1982) and Guéguen and Dienes (1989).

## 2. Micromechanics-Based Constitutive Model for Undrained Deformation

In this section, we extend the micromechanical model from Bhat et al. (2012) for brittle materials whose constitutive response is governed by micro-cracks to incorporate the effect of pore fluid pressure into the constitutive relationship following linear poroelasticity theory. This wing crack model, initiated by Ashby and Sammis (1990), has proved to capture the mechanical response of rocks during both creep experiments at low strain rates ($\le 10^{-4} s^{-1}$) Brantut et al. (2012); Mallet et al. (2015) and Hopkinson bar experiments at high strain rates ($\ge 1 s^{-1}$) (Deshpande and Evans 2008; Bhat et al. 2012). In this manuscript, we use the wing crack model to predict permeability evolution during triaxial loading at high strain rates ($\ge 100 s^{-1}$) that
occurs under undrained conditions. Unfortunately, no continuous measurements of permeability at high strain rates are available yet.

## 2.1. Wing Crack Model Ashby and Sammis (1990)

We consider an isotropic elastic solid that contains an array of penny-shaped cracks with radius $a$ (approximately the grain size). In addition, we assume that all the cracks are aligned at the same angle $\Psi$ to the largest (most negative) remote compressive stress $\sigma_1$ (see Fig. 1); $\sigma_3$ is the least compressive stress and the stresses are considered positive in tension. Bhat et al. (2011) showed that the inclusion of size distribution of cracks and allowance for multiple crack orientations do not significantly affect the failure strength predicted by the model used here. Therefore, we retain the single-crack size and single-orientation assumption. By aligning all the cracks at an angle $\Psi$, we only consider the cracks that are optimally oriented from a Coulomb friction perspective for a given stress state. For the loading depicted in Fig. 1, $2\Psi = \tan^{-1}(1/f)$ with $f$ the friction coefficient. We also assume that the population of cracks that exist prior to the loading has the volume density, $N_v$, and remains fixed during loading (i.e., no nucleation of new cracks). The size and the density of these initial flaws are characterized by a dimensionless scalar defined as,

$$D_0 = (4\pi/3)N_v(\alpha a)^3, \quad (1)$$

where $\alpha a$ is the projection of the crack radius in a vertical plane parallel to the direction of $\sigma_1$, $\alpha = \cos \Psi$. Once the loading is sufficiently large to induce inelastic deformation, all additional damage in the solid is in the form of tensile wing cracks that nucleate at the tips of the starter penny-shaped flaws. Wing cracks grow parallel to $\sigma_1$ axis and open in the direction of the least compressive stress $\sigma_3$ (cf. Fig. 1). A scalar damage parameter $D$ accounts for the current size of cracks and the constant number of flaws per unit volume,

$$D = (4\pi/3)N_v(\alpha a + l)^3. \quad (2)$$

This measure of damage simply quantifies the fraction of the volume occupied by the micro-cracks that

have a significant mechanical influence. As $D$ approaches 1, neighboring wing cracks get closer to coalescence that will potentially lead to the macroscopic fracture of the solid.

## 2.2. Gibbs Free Energy of Damaged Solids

To construct a realistic constitutive description of the micro-cracked solid we consider the energetics of such a solid following the approach of Rice (1975). For a given stress state $\sigma_{ij}$ and damage state the Gibbs free energy density of the solid, $W(\sigma_{ij}, D)$, is the sum of the elastic strain energy density for this solid without cracks inside and deforming purely elastically and the inelastic contribution of $N_v$ cracks per unit volume. We ignore Mode II and Mode III contributions at the incipience of the wing cracks as the wing cracks grow quickly, under dynamic loading rates, to become purely tensile. The Gibbs free energy density of the micro-cracked solid is Bhat et al. (2012)

$$W(\sigma_{ij}, D) = \frac{1}{4\mu} \left[ 2\tau^2 + \frac{3(1 - 2\nu)}{1 + \nu} \sigma^2 \right] + N_v \frac{1 - \nu^2}{E} \int_s (K_I^2(\sigma_{ij}, D) - 2\gamma_{\text{surf}}) ds, \quad (3)$$

where $\mu$ is the shear modulus, $\nu$ is the Poisson ratio, $E$ is the Young's modulus, $s$ is the wing crack surface, $\gamma_{\text{surf}}$ is the surface energy and

$$\sigma = \sigma_{kk}/3 \quad \text{and} \quad \tau = (S_{ij}S_{ij}/2)^{1/2}, \quad (4)$$

where $S_{ij} = \sigma_{ij} - \sigma \delta_{ij}$.

We delineate two macroscopic deformation regimes. Regime I concerns the elastic deformation of the solid at a fixed state of damage. In this regime the stresses are not sufficient to induce sliding or opening of the existing cracks and the solid is assumed to behave like an isotropic linear elastic solid. This ignores macroscopic nonlinear elastic behavior of rocks from closing and frictional sliding of pre-existing penny-shaped cracks shown by experimental data Walsh (1965), Kachanov (1982b), David et al. (2012). In Regime II the remote compressive loading overcomes the frictional resistance acting on the penny-shaped cracks. The shear
<figure><img src="image_4.png" /><figcaption>Figure 1</figcaption></figure>

Geometry of the micromechanical damage mechanics model. a Grain boundary flaws are represented by a volume density of $N_v$ wing cracks opened by sliding over penny-shaped crack of radius a. b Upper view of a initial flaw and its associated wing crack, cross section of a wing crack like in a and effective crack used for hydraulic properties calculations.

stress applied to crack's faces in excess of the frictional resistance creates a wedging force that opens the mouth of each wing crack. Along the activated part of the rim of the penny-shaped cracks, $K_I > 0$, a tensile crack nucleates and grows in the direction of the most compressive principal stress $\sigma_1$ and opens perpendicular the least compressive direction, $\sigma_3$. As in Ashby and Sammis (1990) we account for the interaction between neighboring wing cracks. This interaction induces a tensile stress $\sigma^{(i)}$ that increases as the neighboring wing cracks approach each other.

The stress intensity factor of the wing crack in Regime II, $K_I$, thus (1) increases with wedging force, (2) decreases with crack length, (3) decreases with confining pressure and (4) increases with interaction between neighboring cracks. As in Bhat et al. (2012) we rewrite the stress intensity factor in terms of the scalar damage variable D and thus effectively homogenize the problem. One may look at this as the effective or homogenized stress intensity factor of a unit volume containing $N_v$ cracks. Furthermore, to

generalize the problem for arbitrary stress state we rewrite the normal and shear stresses acting on the crack faces with their corresponding invariant measures, i.e., by $\sigma$ and $\tau$. This is akin to the generalization of Mohr-Coulomb plasticity introduced by Drucker and Prager (1952) and Rudnicki and Rice (1975). The stress intensity factor in Regime II is Bhat et al. (2012)

$$K_I(\sigma_{ij}, D) = \sqrt{\pi a} [A(D)\sigma + B(D)\tau], \quad (5)$$

with

$$A(D) = f c_1(D) + c_3 (f c_2 + 1), \quad (6)$$

$$B(D) = c_1(D) + c_2(D)c_3(D), \quad (7)$$

and

$$c_1(D) = \sqrt{1 - \alpha^2} \pi^{-1} \alpha^{-3/2} \left[ (D/D_0)^{1/3} - 1 + \beta/\alpha \right]^{-3/2}, \quad (8)$$

$$c_2(D) = \left( \frac{\sqrt{1 - \alpha^2}}{\alpha} \right) \left( \frac{D_0^{2/3}}{1 - D^{2/3}} \right), \quad (9)$$
$$c_3(D) = \frac{2\sqrt{\alpha}}{\pi} \left[ (D/D_0)^{1/3} - 1 \right]^{1/2}, \quad (10)$$

with $\beta = 0.1$. Note that Eq. (6) corrects a small typographical error in relation (15) of Bhat et al. (2012).

## 2.3. Constitutive Relationship of Damaged Solids in Undrained Conditions

Finally the strain-stress relationship under isothermal conditions for a given damage and stress state in the material is derived using

$$\epsilon_{ij} = \partial W(\sigma_{ij}, D) / \partial \sigma_{ij}. \quad (11)$$

After inverting for the stress-strain relationship we incorporate the effect of pore fluid pressure in microcracks following linear poroelasticity Biot (1941), Nur and Byerlee (1971), Rice and Cleary (1976), Wang (2000) and ignore the influence of deviatoric stress on pore pressure suggested by Skempton (1954). Therefore the constitutive relationships in Regime I is

$$\sigma_{ij} = 2\mu \left( \epsilon_{ij} + \frac{\nu}{(1-2\nu)} \epsilon \delta_{ij} \right) - \xi p \delta_{ij}, \quad (12)$$

where $\xi$ is the Biot coefficient that may expressed in terms of the drained bulk modulus of the porous material (i.e., at constant pore fluid pressure $p$), $K_d$, and the bulk modulus of the solid, $K_s$, as

$$\xi = 1 - \frac{K_d}{K_s}. \quad (13)$$

From Eqs. (3) and (5) we derive the stress-strain relationship in Regime II (details in 8)

$$\sigma_{ij} = \frac{\mu}{\Gamma} \left[ \begin{aligned} & \left( \frac{3(1-2\nu)}{(1+\nu)} + A_1^2 - \frac{A_1 B_1 \epsilon}{\gamma} \right) \epsilon_{ij} \\ & + \left( \frac{3\nu}{(1+\nu)} + \frac{B_1^2}{2} - \frac{A_1^2}{3} + \frac{A_1 B_1 \epsilon}{3\gamma} \right) \epsilon \delta_{ij} \\ & - \left( \frac{A_1 B_1}{2} \right) \gamma \delta_{ij} \end{aligned} \right] - \xi p \delta_{ij}, \quad (14)$$

where

$$A_1(D) = A(D) \sqrt{\frac{\pi D_0 (1 - \nu)}{\alpha^3}}, \quad (15)$$

$$B_1(D) = B(D) \sqrt{\frac{\pi D_0 (1 - \nu)}{\alpha^3}}, \quad (16)$$

and the conjugate strains are defined as

$$\epsilon = \epsilon_{kk}, \quad (17)$$

$$\gamma = \sqrt{2e_{ij}e_{ij}}, \quad (18)$$

with $e_{ij} = \epsilon_{ij} - \epsilon\delta_{ij}$.

## 2.4. Crack Growth Law

To complete the above energy-based description of the coupled hydro-mechanical properties of the brittle material we need to prescribe the damage evolution law and the evolution of the pore fluid pressure. Differentiating $D$ with respect to time we obtain

$$\frac{dD}{dt} = (3D^{2/3}D_0^{1/3}/\alpha a)\frac{dI}{dt}, \quad (19)$$

where $dI/dt \equiv \nu$ is the wing crack tip speed. We use the experimentally motivated crack growth law developed by Bhat et al. (2012), which is valid over a wide range of high strain rates ($\ge 1 \text{ s}^{-1}$) and has the capability of suddenly accelerating or decelerating/stopping a propagating crack. The crack-tip speed is obtained by solving the nonlinear equation

$$\frac{K_I(1 - \nu/c_R)}{\sqrt{1 - \nu/c_p}} = K_{IC}^{SS} \left[ \frac{1 + (\nu/v_m)^5}{\sqrt{1 - \nu/c_p}} \right], \quad (20)$$

where $c_R$ and $c_p$ are the Rayleigh wave speed and P-wave speed of the material, respectively, and $v_m$ the crack branching speed. Once the instantaneous crack-tip speed is determined the state of damage is updated by integrating the above equation for $dD/dt$ and consequently the mechanical and hydraulic states are updated as well.

## 2.5. Pore Fluid Pressure Evolution in Undrained Conditions

During deformation at high strain rates, the rate of crack growth is large relative to the rate of fluid flow and we can reasonably expect undrained deformation of the sample during the simulation. Assuming a hydraulic diffusivity of approximately $10^{-4} \text{ m}^2 \text{ s}^{-1}$ (consistent with our simulation results), the diffusive length scale for applied strain rates of $1000 \text{ s}^{-1}$ is 10 cm. This is larger than the crack size and typically the
size of the entire sample. The rate of change in fluid mass, m (mass of fluid per unit bulk volume of porous material), is

$$ \dot{m} = \rho_f \xi \left( \dot{\epsilon} + \frac{\xi}{K_u - K} \dot{p} \right), \quad (21) $$

where $K_u$ is the undrained bulk modulus defined as

$$ K_u = K + \frac{\xi^2 K_f}{\Phi_0 + (\xi - \Phi_0) K_f / K_s}, \quad (22) $$

with $\Phi_0$ the initial porosity and $K_f$ the fluid bulk modulus. Note that in our case the strain is the sum of the elastic and the inelastic (damage-induced) strain. Setting $\dot{m} = 0$ for undrained deformation we obtain the rate of change in pore fluid pressure,

$$ \dot{p} = \dot{\epsilon} (K - K_u) / \xi. \quad (23) $$

From this complete hydro-mechanical description of deformation, we now work out the evolution of micro-crack aperture (opening) during inelastic deformation. Together with micro-crack length evolution we can then work out the net volume available for fluid flow in the fractured medium which will then be used to calculate the permeability of the fractured medium.

### 3. Opening and Sliding of the Micro-cracks

In Regime I the solid is deforming purely elastically. This implies $K_I < 0$. We relate the aperture of the penny-shaped flaws, $w$, to the initial porosity, $\Phi_0$, via $w = \Phi_0 / (N_v \pi a^2)$.

In Regime II tensile wing cracks are wedged open by sliding on the penny-shaped shear cracks. The surface created by two of these wing cracks (associated with each penny-shaped shear crack) is $2\bar{v}l$ where $\bar{v}$ is the average wing crack opening,

$$ \bar{v} = (2\theta_m)^{-1} \int_{-\theta_m}^{\theta_m} l^{-1} \int_0^l v(x) \, dx \, d\theta. \quad (24) $$

This wing crack opening is averaged over the length $l$ and over the angular range, $2\theta_m$, on the shear crack over which tensile crack growth is active; see Fig. 1. The exact solution for $K_I(\theta)$ at the tip of a circular crack loaded in Mode II in 3D is given by Tada et al.

(2000), page 24.21. Following this solution we can derive, as in Kachanov (1982a, b),

$$ \cos \theta_m = \sqrt{3\pi(1 - 2\nu)K_{IC}} / [8(\tau + f\sigma)\sqrt{a}], \quad (25) $$

where $K_{IC}$ is the fracture toughness. To ensure a realistic description of angular range of activation of the tensile wing cracks, that is of ignoring any healing effect of the cracks, we add an additional constraint on this angle, $d\theta_m/dt \ge 0$.

We now hypothesize that the volumetric inelastic strain $\epsilon^{in} = \epsilon_{kk}^{in}$ is a measure of net opening due to all the tensile wing cracks. The volumetric inelastic strain is the volume created by the opening of two growing tensile wing cracks, for each initial shear crack, times the volume density of these mono-sized, mono-oriented flaws. Therefore

$$ \epsilon^{in} = N_v \int_{-\theta_m}^{\theta_m} (2\bar{v}l) a \cos \theta \, d\theta, \quad (26) $$

where $a \cos \theta$ is the chord length. Hence the average opening of the tensile wing crack in this regime is

$$ \bar{v} = \epsilon^{in} / (4N_v l a \sin \theta_m). \quad (27) $$

Note that the volumetric inelastic strain depends both on the invariant shear and normal stresses.

The aperture of the penny-shaped cracks in Regime II is related to the initial porosity $\Phi_0$ as $w = \Phi_0 / (N_v \pi a^2)$. This ignores the surface reduction of circular cracks due to the growth of tensile secondary cracks at their rim. However the volume created by the opening of the wing cracks is considerably larger than the one generated by opening of the circular cracks and hence we retain this approximation. Furthermore, we have assumed that the penny-shaped cracks remain open during sliding and growth of the wing cracks because of the pore fluid pressure.

### 4. Permeability Model

The statistical crack mechanics-based permeability model developed by Dienes (1982) calculates the permeability tensor in the cracked material by summing the contribution of flow in individual penny-shaped cracks. Following Simpson et al. (2001) we
suggest that the fluid flow in a kinked crack can be described by an “equivalent” penny-shaped crack. This allows us to link the parameters of the mechanical model (i.e., l, $\bar{v}$, $\theta_m$ and w) to the parameters required to calculate the permeability tensor (c the radius of the equivalent circular crack and $\omega$ its aperture). We choose to replace the actual penny and wing cracks by an equivalent crack that has the same total volume of fluid inside. From geometrical considerations (cf. Fig. 1), $c = (l \sin \Psi)/\sin \eta$ and $\tan \eta = (l \sin \Psi)/(a + l \cos \Psi)$. This gives

$$c = a\sqrt{1 + \alpha^2 [(D/D_0)^{2/3} - 1]} \quad (28)$$

This latter relation is valid for every regime of deformation. Setting $D = D_0$ we retrieve our initial penny-shaped crack, i.e., $c = a$.

The volume inside the effective circular crack is $\pi c^2 \omega$ whereas the volume inside a shear and wing cracks written in terms of the maximal angle over which tensile cracks grow is,

$$V = wa^2[\pi - 2\theta_m - \sin 2\theta_m] + 4\bar{v}la \sin \theta_m. \quad (29)$$

The first term corresponds to the reduced area of the initial circular crack due to activation of the tensile wing crack on which sliding no longer exists. Therefore the aperture of the effective circular crack is,

$$\omega = \frac{wa^2[\pi - 2\theta_m - \sin(2\theta_m)]}{\pi c^2} + \frac{4\bar{v}la \sin \theta_m}{\pi c^2}. \quad (30)$$

Using the apertures worked out in the previous section, the equivalent aperture, $\omega$, in Regime I is

$$\omega = \frac{\Phi_0}{N_v \pi a^2} \quad (31)$$

and, in Regime II,

$$\omega = \left[ \frac{1}{N_v \pi c^2} \right] \left[ \frac{\Phi_0[\pi - 2\theta_m - \sin 2\theta_m]}{\pi} + \epsilon^{in} \right]. \quad (32)$$

Note that in the elastic regime, inelastic deformation is absent, $D = D_0$ and $\theta_m = 0$.

Following Dienes (1982) and Guéguen and Dienes (1989), the isotropic permeability of a cracked rock is

$$k = \frac{4\pi^2}{15} PA^3 N_v c^5, \quad (33)$$

where $A = \omega/2c$ is the aspect ratio and P is the fraction of cracks that are hydraulically connected. This percolation factor is assumed to be a linear function of D here,

$$P = \frac{D - D_0}{1 - D_0} \quad (34)$$

Using the expression for the apertures, Eqs. (31) and (32), and the expression for c (Eq. (28)) we have an expression for the permeability for a given mechanical and damage state in the material.

## 5. Permeability Evolution During Triaxial Loading at High Strain Rates

We implemented explicit constitutive update in ABAQUS-Explicit by writing a USER MATERIAL (VUMAT) subroutine Bhat et al. (2012). We conducted axisymmetric calculations on cylindrical specimens of height H with a constant velocity V applied on one end resulting in average applied strain rate of $V/H$—a setting similar to a Hopkinson bar experiment except that we were able to impose both confining pressure and the initial pore fluid pressure. Damage and permeability are computed at each grid point. They are homogenized mechanical and hydraulic properties of a unit volume with $N_v$ cracks. For clarity, we display material properties only at one grid point located in the middle of the deformed sample.

We performed numerical simulations for a range of applied strain rates from 100 to 1000 s⁻¹ and various differential pressures $P_{diff}$ (confining pressure minus pore fluid pressure) for Westerly Granite. The parameters used for simulations are displayed in Table 1. In all simulations the pore fluid pressure is set to 50 MPa and we varied the confining pressure. The simulation results shown in Figs. 2, 3 and 4 are done at an initial differential pressure of 20 MPa. In Fig. 5 we show the effect of initial differential pressure on permeability evolution in Westerly Granite. Note that we only show permeabilities above $10^{-23}$ m².

As Guéguen and Schubnel (2003) emphasized, based on experiments, our modeled permeability evolution highlights two regimes: a percolative regime and a connected regime. In the percolative regime permeability increases rapidly by several order of magnitudes
<table><caption>Table 1</caption><caption>Parameter values of Westerly Granite used in the simulations</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Units</th><th>Granite</th></tr></thead><tbody><tr><td>Rock properties</td><td></td><td></td><td></td></tr><tr><td>Young modulus</td><td><i>E</i></td><td>GPa</td><td>70</td></tr><tr><td>Poisson ratio</td><td><i>ν</i></td><td>1</td><td>0.30</td></tr><tr><td>Drained bulk modulus</td><td><i>K</i></td><td>GPa</td><td>27</td></tr><tr><td>Density</td><td><i>ρ</i></td><td>kg m<sup>-3</sup></td><td>2640</td></tr><tr><td>Grain size</td><td><i>a</i></td><td>mm</td><td>0.5</td></tr><tr><td>Initial damage</td><td><i>D</i><sub>0</sub></td><td>1</td><td>0.1</td></tr><tr><td>Friction coefficient</td><td><i>f</i></td><td>1</td><td>0.64</td></tr><tr><td>Fracture toughness</td><td><i>K</i><sub>IC</sub></td><td>MPa m<sup>1/2</sup></td><td>1</td></tr><tr><td>Crack branching speed ratio</td><td><i>ν</i><sub>m</sub>/<i>c</i><sub>R</sub></td><td>1</td><td>0.3</td></tr><tr><td>Initial porosity</td><td><i>Φ</i><sub>0</sub></td><td>%</td><td>0.25</td></tr><tr><td>Sample size</td><td></td><td></td><td></td></tr><tr><td>Sample height</td><td><i>H</i></td><td>cm</td><td>10</td></tr><tr><td>Sample radius</td><td><i>R</i></td><td>cm</td><td>5</td></tr></tbody></table>

All parameters are from Ashby and Sammis (1990) except the initial porosity from Zoback and Byerlee (1975) and the drained bulk modulus from Heard and Page (1982). The bulk modulus of water, $K_f$, is 2 GPA Palciauskas and Domenico (1982)

<figure><img src="image_5.png" /><figcaption>Figure 2</figcaption></figure>

Numerically predicted permeability evolution (a) and damage (b) in Westerly Granite during triaxial loading at $P_{\text{diff}}$ = 20 MPa. The gray and black curves are for an applied strain rate of 300 and 310 s$^{-1}$, respectively.

<figure><img src="image_6.png" /><figcaption>Figure 3</figcaption></figure>

Evolution of the crack aspect ratio A (a) and crack radius c (b) with damage D in Westerly Granite during triaxial loading at $P_{\text{diff}}$ = 20 MPa and at an applied strain rate of 500 s$^{-1}$.

<figure><img src="image_7.png" /><figcaption>Figure 4</figcaption></figure>

Numerically predicted permeability evolution (a) and crack aspect ratio A (b) with damage in granite during triaxial loading at $P_{\text{diff}}$ = 20 MPa. The gray and black curves are for an applied strain rate of 300 and 500 s$^{-1}$, respectively.

<figure><img src="image_8.png" /><figcaption>Figure 5</figcaption></figure>

Numerically predicted permeability evolution in granite during triaxial loading at an applied strain rate of 1000 s$^{-1}$ and for a various confining pressures $P_{\text{conf}}$.

(see Fig. 2a). Before onset of damage (axial strain lower than 0.47 % in Fig. 2a) the percolation factor vanishes for intact rocks and, hence, the permeability is null. When cracks growth starts ($D > D_0$, see Fig. 2b), fluid seeps into fresh wing cracks and permeability is no longer zero. A fluid path, although tortuous, allows hydraulic connection between two neighboring microcracks, resulting in permeability increase. As the cracks are growing, the probability of connectivity, P, continues to increase. However, the dominant physical process that dictates permeability evolution is the competition between the volume of fluid present in growing wing cracks and penny-shaped cracks. Reduction of penny-shaped cracks' volume by activation of secondary cracks at their rim together with small wing cracks' aperture introduce a plateau (or slight decrease) in permeability evolution. In this regime damage is increasing but the aspect ratio, of the tensile wing crack, is decreasing (see Fig. 3a) and
hence permeability slightly decreases with damage (see Fig. 4). A connected regime appears when neighboring cracks are sufficiently close. Interaction-induced tensile stress leads to an important augmentation of wing cracks opening resulting in their runaway growth (Fig. 3a, b). Le Ravalec and Guéguen (1994) and Zhu and Wong (1999) also found that, with crack propagation increment, permeability can increase with improved connectivity and decrease when crack aperture compensates for the increase of crack length. For deformation of Westerly Granite at $300 \text{ s}^{-1}$ (Fig. 2a), at 0.51 % of axial strain damage increases dramatically from 0.7 to 1 leading to a jump in permeability by three orders of magnitude. There the rock becomes granular (all the micro-cracks have linked up) and the wing crack model cannot capture the changes in permeability anymore. Presumably, the percolation factor would decrease drastically in this regime because of the tortuosity. However, this is beyond the scope of this paper.

An interesting feature is the effect of strain rate on permeability evolution; increasing strain rate leads to slower permeability increase in the connected regime (see Fig. 2a). This physically means that at low strain rates, because cracks propagation resistance is low, cracks acceleration is more sudden (see Fig. 2b) and hence permeability increases quickly with the runaway growth of cracks. Figure 4a shows that the hydraulic connectivity at low level of damage decreases with increasing strain rate as inertial resistance of the cracks is dominant. The crack aspect ratio at low level of damage decreases with strain rate (Fig. 4b). For larger values of damage interaction-induced tensile stresses build up, which increase with strain rate, which leads to bigger aperture of the wing cracks and hence, larger value of permeability.

Figure 5 shows the permeability evolution of Westerly Granite at various differential pressures. In early stage of deformation the dominant contributions to the Mode I stress intensity factor at the tip of the wing cracks are the wedging force associated with sliding on the penny cracks and the remote confining stress closing the wing cracks. The unbroken ligament between neighboring wing cracks is long and, hence, the interaction-induced contribution is negligible. Increasing the confining pressure delays the onset of crack growth and the initiation of hydraulic

connectivity in the rock. However once damage has developed all permeability curves collapse. The model does not reveal strong dependence on the differential pressure as hypothesized in, e.g., Rice (1992).

## 6. Discussion

The permeability curves presented in this manuscript are simulated for strain rates at which continuous permeability measurements are not available. The explicit time stepping method implemented in ABAQUS-Explicit does not allow simulations for applied strain rates as low as $10^{-5} \text{ s}^{-1}$ at which permeability evolution is typically measured Mitchell and Faulkner (2008). For the sake of simplicity, we ignored any permeability anisotropy in the model that could arise from preferred cracks orientation. We have also simplified the treatment of pore fluid evolution by only modifying the constitutive relationship according to linear poroelasticity. A more careful derivation of the stress-strain relationship would incorporate in the Gibbs free energy density the role of the fluid pressure in the system. The pore fluid pressure would potentially also modify the stress intensity factor expression as formulated by Ashby and Sammis (1990). However, it appears that the pore fluid pressure evolution does not affect significantly the permeability evolution.

## 7. Conclusions and Perspectives

We developed a micromechanics-based model of permeability evolution, extending the micromechanical model developed by Bhat et al. (2012), in brittle materials that captures for the first time both the percolative regime and the connected regime (Guéguen and Schubnel 2003) for simulations at high strain rates.

Using this coupled hydro-mechanical model valid under undrained conditions, we plan to implement the developed stress-strain relationship and the pore fluid pressure evolution into earthquake dynamic rupture codes. With this model, one can simulate earthquake cycles with coseismic dynamic ruptures under undrained conditions and pore fluid diffusion under drained conditions in interseismic periods. Purely accounting for dynamic damage off the fault already
has a significant role on radiated ground motion spectra (Thomas and Bhat 2016). While the dynamic evolution of permeability, and corresponding changes in pore pressure, will not affect thermal pressurization like phenomenon at the rupture tip, we anticipate it to play a role behind the crack tip thus effecting the size of the earthquake pulse. This should play a role in the radiated ground motion.

We also aim to extend this coupled model to the granular regime (when all the micro-cracks have linked up or when the mechanical model suffers from loss of convexity).

### Acknowledgments

T. P.'s research was supported by the GPX master's fellowship and the NSF grant Division for Materials Research 14-20570 to Harvard University. H. S. B. would like to thank Tom Mitchell at U. C. London for inspiring this work and for help in getting it started.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## Appendix: Derivations of the Constitutive Laws for Undrained Deformation

In this section, we derive the stress-strain relationship in Regime II for undrained deformation. Following (Bhat et al. 2012) the Gibbs free energy in Regime II is

$$W = \frac{1}{4\mu} \left[ 2\tau^2 + \frac{3(1-2\nu)}{1+\nu} \sigma^2 \right] + \frac{1}{4\mu} (A_1\sigma + B_1\tau)^2. \quad (35)$$

Note that the constant contribution from the surface energy is ignored. Using

$$\epsilon_{ij} = \frac{\partial W(\sigma_{ij}, D)}{\partial \sigma_{ij}}, \quad (36)$$

we obtain

$$\epsilon_{ij} = \frac{1}{2\mu} \left[ \left( 1 + \frac{A_1 B_1 \sigma}{2\tau} + \frac{B_1^2}{2} \right) \sigma_{ij} - \left( \frac{3\nu}{(1+\nu)} + \frac{A_1 B_1 \sigma}{2\tau} - \frac{A_1^2}{3} + \frac{B_1^2}{2} \right) \sigma \delta_{ij} + \left( \frac{A_1 B_1}{3} \right) \tau \delta_{ij} \right]. \quad (37)$$

Recasting the free energy expression in terms of strain invariants $\epsilon$ and $\gamma$, we obtain

$$\epsilon = \frac{\partial W}{\partial \sigma} = \left[ \frac{3(1-2\nu)}{2\mu(1+\nu)} + \frac{A_1^2}{2\mu} \right] \sigma + \left[ \frac{A_1 B_1}{2\mu} \right] \tau, \quad (38)$$

$$\gamma = \frac{\partial W}{\partial \tau} = \left[ \frac{A_1 B_1}{2\mu} \right] \sigma + \left[ \frac{1}{\mu} + \frac{B_1^2}{2\mu} \right] \tau. \quad (39)$$

Solving for $\sigma$ and $\tau$ from above one obtains,

$$\sigma = \mu(a_1\epsilon + b_1\gamma), \quad (40)$$

$$\tau = \mu(b_1\epsilon + b_2\gamma), \quad (41)$$

where

$$a_1 = \frac{1}{\Gamma} \left( 1 + \frac{B_1^2}{2} \right) \quad (42)$$

$$b_1 = -\frac{1}{\Gamma} \left( \frac{A_1 B_1}{2} \right) \quad (43)$$

$$b_2 = \frac{1}{\Gamma} \left[ \frac{A_1^2}{2} + \frac{3(1-2\nu)}{2(1+\nu)} \right] \quad (44)$$

and

$$\Gamma = \left[ \frac{3(1-2\nu)}{2(1+\nu)} + \frac{3(1-2\nu)B_1^2}{4(1+\nu)} + \frac{A_1^2}{2} \right]. \quad (45)$$

Since the Free Energy (Gibbs or Helmholtz) is quadratic in form, we assume the Helmholtz Free energy is of the form

$$W^H = \mu(L_1\epsilon^2 + L_2\gamma^2 + L_3\epsilon\gamma) \quad (46)$$

This then implies (by comparison with Eq. 41),

$$\sigma = \frac{\partial W^H}{\partial \epsilon} = \mu[2L_1\epsilon + L_3\gamma] \quad (47)$$

$$\tau = \frac{\partial W^H}{\partial \gamma} = \mu[L_3\epsilon + 2L_2\gamma] \quad (48)$$
and thus,

$$L_1 = \frac{a_1}{2} \qquad (49)$$

$$L_2 = \frac{b_2}{2} \qquad (50)$$

$$L_3 = b_1. \qquad (51)$$

Therefore,

$$W^H = \mu \left[ \frac{a_1}{2} \epsilon^2 + \frac{b_2}{2} \gamma^2 + b_1 \epsilon \gamma \right], \qquad (52)$$

rewritten as,

$$W^H = \frac{\mu}{\Gamma} \left\{ \left( \frac{1}{2} + \frac{B_1^2}{4} \right) \epsilon^2 + \left[ \frac{A_1^2}{4} + \frac{3(1-2\nu)}{4(1+\nu)} \right] \gamma^2 - \left( \frac{A_1 B_1}{2} \right) \epsilon \gamma \right\}. \qquad (53)$$

We can now evaluate

$$\sigma_{ij} = \frac{\partial W^H}{\partial \epsilon_{ij}}, \qquad (54)$$

add the effect of pore fluid pressure following linear poroelasticity theory and retrieve the stress-strain relationship in Regime II written in Eq. (14).

REFERENCES

Ashby, M. F., & Sammis, C. G. (1990). The damage mechanics of brittle solids in compression. *Pure and Applied Geophysics*, **133**(3), 489–521.

Bhat, H. S., Sammis, C. G., & Rosakis, A. J. (2011). The micromechanics of westerley granite at large compressive loads. *Pure and Applied Geophysics*, **168**(12), 1–18. doi:10.1007/s00024-011-0271-9.

Bhat, H. S., Rosakis, A. J., & Sammis, C. G. (2012). A micromechanics based constitutive model for brittle failure at high strain rates. *Journal of Applied Mechanics*, **79**(3), 031016.

Biot, M. A. (1941). General theory of three-dimensional consolidation. *Journal of Applied Physics*, **12**(2), 155–164.

Brantut, N., Baud, P., Heap, M. J., & Meredith, P. G. (2012). Micromechanics of brittle creep in rocks. *Journal of Geophysical Research: Solid Earth*, **117**, B08412. doi:10.1029/2012JB009299.

Cappa, F., & Rutqvist, J. (2011a). Impact of CO₂ geological sequestration on the nucleation of earthquakes. *Geophysical Research Letters*, **38**, L17313. doi:10.1029/2011GL048487.

Cappa, F., & Rutqvist, J. (2011b). Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO₂. *International Journal of Greenhouse Gas Control*, **5**(2), 336–346.

Cox, S. F. (2016). Injection-driven swarm seismicity and permeability enhancement: implications for the dynamics of hydrothermal ore systems in high fluid-flux, overpressured faulting regimes-an invited paper. *Economic Geology*, **111**(3), 559–587.

David, E. C., Brantut, N., Schubnel, A., & Zimmerman, R. W. (2012). Sliding crack model for nonlinearity and hysteresis in the uniaxial stressstrain curve of rock. *International Journal of Rock Mechanics and Mining Sciences*, **52**, 9–17. doi:10.1016/j.ijrmms.2012.02.001.

Deshpande, V. S., & Evans, A. G. (2008). Inelastic deformation and energy dissipation in ceramics: a mechanism-based constitutive model. *Journal of the Mechanics and Physics of Solids*, **56**(10), 3077–3100.

Dienes, J. K. (1982). Permeability, percolation and statistical crack mechanics, in *The 23rd US Symposium on Rock Mechanics (USRMS)*, edited by Y. S. Kim, pp. 86–94.

Drucker, D. C., & Prager, W. (1952). Soil mechanics and plastic analysis or limit design. *Quarterly of Applied Mathematics*, **10**(2), 157–165.

Evans, K. F. (2005). Permeability creation and damage due to massive fluid injections into granite at 3.5 km at soultz: 2. critical stress and fracture strength. *Journal of Geophysical Research: Solid. Earth*, **110**, B04204. doi:10.1029/2004JB003169.

Evans, K. F., Genter, A., & Sausse, J. (2005). Permeability creation and damage due to massive fluid injections into granite at 3.5 km at soultz: 1. borehole observations. *Journal of Geophysical Research: Solid. Earth*, **110**, B04203. doi:10.1029/2004JB003168.

González, P. J., Tiampo, K. F., Palano, M., Cannavó, F., & Fernández, J. (2012). The 2011 lorca earthquake slip distribution controlled by groundwater crustal unloading. *Nature Geoscience*, **5**, 821–825. doi:10.1038/NGEO1610.

Guéguen, Y., & Dienes, J. (1989). Transport properties of rocks from statistics and percolation. *Mathematical Geology*, **21**(1), 1–13.

Guéguen, Y., & Schubnel, A. (2003). Elastic wave velocities and permeability of cracked rocks. *Tectonophysics*, **370**, 163–176. doi:10.1016/S0040-1951(03)00184-7.

Healy, J. H., Rubey, W. W., Griggs, D. T., & Raleigh, C. B. (1968). The denver earthquakes. *Science*, **161**(3848), 1301–1310.

Heard, H. C., & Page, L. (1982). Elastic moduli, thermal expansion, and inferred permeability of two granites to 350 °c and 55 megapascals. *Journal of Geophysical Research: Solid Earth*, **87**(B11), 9340–9348. doi:10.1029/JB087iB11p09340.

Kachanov, M. L. (1982a). A microcrack model of rock inelasticity part I: Frictional sliding on microcracks. *Mechanics of Materials*, **1**(1), 19–27.

Kachanov, M. L. (1982b). A microcrack model of rock inelasticity part II: Propagation of microcracks. *Mechanics of Materials*, **1**(1), 19–27.

Keranen, K. M., Savage, H. M., Abers, G. A., & Cochran, E. S. (2013). Potentially induced earthquakes in okahoma, usa: links between wastewater injection and the 2011 mw 5.7 earthquake sequence. *Geology*, **41**(6), 699–702. doi:10.1130/G34045.1.

Le Ravalec, M., & Guéguen, Y. (1994). Permeability models for heated saturated igneous rocks. *Journal of Geophysical Research: Solid Earth*, **99**(B12), 24251–24261. doi:10.1029/94JB02124.

Mallet, C., Fortin, J., Guéguen, Y., & Bouyer, F. (2015). Brittle creep and subcritical crack propagation in glass submitted to triaxial conditions. *Journal of Geophysical Research: Solid Earth*, **120**(2), 879–893. doi:10.1002/2014JB011231.
Miller, S. A., Collettini, C., Chiaraluce, L., Cocco, M., Barchi, M., & Kaus, B. J. P. (2004). Aftershocks driven by a high-pressure co2 source at depth. *Nature*, **427**(6976), 724–727.

Mitchell, T. M., & Faulkner, D. R. (2008). Experimental measurements of permeability evolution during triaxial compression of initially intact crystalline rocks and implications for fluid flow in fault zones. *Journal of Geophysical Research: Solid Earth*, **113**, B114412. doi:10.1029/2008JB005588.

Morton, N., T. M. Mitchell, M. L. Doan, G. Girty, T. K. Rockwell, and J. Renner (2012), Ultra-high permeability induced by seismic shockwaves near the san jacinto fault, 2012 Fall Meeting, AGU, San Francisco, Calif., 3-7 Dec.

Murphy, H., C. Huang, Z. Dash, G. Zyvoloski, & A. White (2004), Semianalytical solutions for fluid flow in rock joints with pressure-dependent openings, *Water Resour. Res.*, **40**(12), W12,506, doi:10.1029/2004WR003005.

Nur, A., & Byerlee, J. D. (1971). An exact effective stress law for elastic deformation of rock with fluids. *Journal of Geophysical Research*, **76**(26), 6414–6419.

Palciauskas, V. V., & Domenico, P. A. (1982). Characterization of drained and undrained response of thermally loaded repository rocks. *Water Resources Research*, **18**(2), 281–290. doi:10.1029/WR018i002p00281.

Rice, J. R. (1975). Continuum mechanics and thermodynamics of plasticity in relation to microscale deformation mechanisms. In A. S. Argon & M. I. T. Press (Eds.), *Constitutive Equations in Plasticity*. Cambridge: Mass.

Rice, J. R. (1992), Fault stress states, pore pressure distributions, and the weakness of the san andreas fault, in *Fault Mechanics and Transport Properties in Rocks*, edited by B. Evans and T. F. Wong, pp. 475–503, Academic Press.

Rice, J. R. (2006). Heating and weakening of faults during earthquake slip. *Journal of Geophysical Research: Solid Earth*, **111**, B05311. doi:10.1029/2005JB004006.

Rice, J. R., & Cleary, M. P. (1976). Some basic stress diffusion solutions for fluid-saturated elastic porous media with

compressible constituents. *Reviews of Geophysics*, 14(2), 227–241.

Rudnicki, J. & Rice, J. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. Journal of the Mechanics and Physics of Solids, 23(6):371–394.

Sibson, R. H. (1994). Crustal stress, faulting and fluid flow. *Geological Society, London, Special Publications*, **78**(1), 69–84. doi:10.1144/GSL.SP.1994.078.01.07.

Simpson, G., Guéguen, Y., & Schneider, F. (2001). Permeability enhancement due to microcrack dilatancy in the damage regime. *Journal of Geophysical Research*, 106(B3), 3999–4016.

Skempton, A. (1954). The pore-pressure coefficients a and b. *Geotechnique*, 4(4), 143–147.

Tada, H., Paris, P. C., & Irwin, G. R. (2000). *The stress analysis of cracks handbook*. New York: ASME Press.

Thomas, M. Y., & H. S. Bhat (2016), *Fault Zone Properties and Processes during Dynamic Rupture*, vol. (in revision), chap. Dynamic off-fault damage and seismic radiation, AGU Monographs.

Walsh, J. B. (1965). The effect of cracks on the uniaxial elastic compression of rocks. *Journal of Geophysical Research*, **70**(2), 399–411. doi:10.1029/JZ070i002p00399.

Wang, H. (2000), *Theory of linear poroelasticity with applications to geomechanics and hydrogeology*, Princeton University Press.

Zhu, W., & Wong, T. (1999). Network modeling of the evolution of permeability and dilatancy in compact rock. *Journal of Geophysical Research: Solid Earth*, **104**(B2), 2963–2971.

Zoback, M. D., & Byerlee, J. D. (1975). The effect of microcrack dilatancy on the permeability of westerly granite. *Journal of Geophysical Research: Solid Earth*, **80**(5), 752–755. doi:10.1029/JB080i005p00752.

Zoback, M. D., & Gorelick, S. M. (2012). Earthquake triggering and large-scale geologic storage of carbon dioxide. *Proceedings of the National Academy of Sciences USA*, **109**(26), 10164–10168. doi:10.1073/pnas.1202473109.

(Received March 14, 2016, revised July 11, 2016, accepted July 12, 2016, Published online July 27, 2016)