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

We develop a micromechanics-based permeability evolution model for brittle materials at high strain rates (≥100\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ge 100$$\end{document} s-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}). 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.


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 expansivity 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 pressuredependent 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 w 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 2 or in enhanced geothermal systems. Not just injection but also rapid extraction of groundwater 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 Soultzsous-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).

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 microcracks 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 ( 10 À4 s À1 ) Brantut et al. (2012); Mallet et al. (2015) and Hopkinson bar experiments at high strain rates ( ! 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 ( ! 100 s À1 ) that occurs under undrained conditions. Unfortunately, no continuous measurements of permeability at high strain rates are available yet. 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 W to the largest (most negative) remote compressive stress r 1 (see Fig. 1); r 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 W, 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, 2W ¼ 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,

Wing Crack Model
where aa is the projection of the crack radius in a vertical plane parallel to the direction of r 1 , a ¼ cos W. 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 r 1 axis and open in the direction of the least compressive stress r 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, 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.

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 r ij and damage state the Gibbs free energy density of the solid, Wðr 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) where l is the shear modulus, m is the Poisson ratio, E is the Young's modulus, s is the wing crack surface, c surf is the surface energy and where S ij ¼ r ij À rd 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 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 r 1 and opens perpendicular the least compressive direction, r 3 . As in Ashby and Sammis (1990) we account for the interaction between neighboring wing cracks. This interaction induces a tensile stress r ð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 r and s. 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) with and 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.

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 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 where n 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 From Eqs. (3) and (5) we derive the stress-strain relationship in Regime II (details in 8) where A 1 ðDÞ ¼ AðDÞ B 1 ðDÞ ¼ BðDÞ and the conjugate strains are defined as with e ij ¼ ij À d ij .

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 where dl=dt v 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 ( ! 1 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 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.

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 m 2 s À1 (consistent with our simulation results), the diffusive length scale for applied strain rates of 1000 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 where K u is the undrained bulk modulus defined as with U 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 _ m ¼ 0 for undrained deformation we obtain the rate of change in pore fluid pressure, 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.

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, U 0 , via w ¼ U 0 =ðN v pa 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 vl where v is the average wing crack opening This wing crack opening is averaged over the length l and over the angular range, 2h m , on the shear crack over which tensile crack growth is active; see Fig. 1. The exact solution for K I ðhÞ 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), 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, dh m =dt ! 0. We now hypothesize that the volumetric inelastic strain in ¼ in kk 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 where a cos h is the chord length. Hence the average opening of the tensile wing crack in this regime is 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 U 0 as w ¼ U 0 =ðN v pa 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.

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 pennyshaped 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; v; h m and w) to the parameters required to calculate the permeability tensor (c the radius of the equivalent circular crack and x 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 WÞ= sin g and tan g ¼ ðl sin WÞ=ða þ l cos WÞ. This gives 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 pc 2 x whereas the volume inside a shear and wing cracks written in terms of the maximal angle over which tensile cracks grow is, 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, Using the apertures worked out in the previous section, the equivalent aperture, x, in Regime I is and, in Regime II, Note that in the elastic regime, inelastic deformation is absent, D ¼ D 0 and h m ¼ 0. Following Dienes (1982) and Guéguen and Dienes (1989), the isotropic permeability of a cracked rock is where A ¼ x=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, 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.

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 À1 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 2 .
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 (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 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) Numerically predicted permeability evolution (a) and crack aspect ratio A (b) with damage in granite during triaxial loading at P diff ¼ 20 MPa. The gray and black curves are for an applied strain rate of 300 and 500 s À1 , respectively.

Figure 5
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 conf .
hence permeability slightly decreases with damage (see Fig. 4). A connected regime appears when neighboring cracks are sufficiently close. Interactioninduced 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 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 interactioninduced 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).

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 ABA-QUS-Explicit does not allow simulations for applied strain rates as low as 10 À5 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.

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). 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 Note that the constant contribution from the surface energy is ignored. Using we obtain Recasting the free energy expression in terms of strain invariants and c, we obtain Solving for r and s from above one obtains, where and Since the Free Energy (Gibbs or Helmholtz) is quadratic in form, we assume the Helmholtz Free energy is of the form This then implies (by comparison with Eq. 41), T. Perol, H. S. Bhat Pure Appl. Geophys. and thus, L 3 ¼b 1 : Therefore, rewritten as, We can now evaluate add the effect of pore fluid pressure following linear poroelasticity theory and retrieve the stress-strain relationship in Regime II written in Eq. (14).