Rock Mech Rock Eng DOI 10.1007/s00603-016-0996-y

ORIGINAL PAPER

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability Versus Diagenesis Sotiris Alevizos1 • Thomas Poulet1,2 • Mustafa Sari1 • Martin Lesueur1,2 Klaus Regenauer-Lieb1 • Manolis Veveakis1,2

•

Received: 5 November 2015 / Accepted: 25 April 2016  Springer-Verlag Wien 2016

Abstract Understanding the formation, geometry and fluid connectivity of nominally impermeable unconventional shale gas and oil reservoirs is crucial for safe unlocking of these vast energy resources. We present a recent discovery of volumetric instabilities of ductile materials that may explain why impermeable formations become permeable. Here, we present the fundamental mechanisms, the critical parameters and the applicability of the novel theory to unconventional reservoirs. We show that for a reservoir under compaction, there exist certain ambient and permeability conditions at which diagenetic (fluid-release) reactions may provoke channelling localisation instabilities. These channels are periodically interspersed in the matrix and represent areas where the excess fluid from the reaction is segregated at high velocity. We find that channelling instabilities are favoured from pore collapse features for extremely low-permeability formations and fluid-release diagenetic reactions, therefore providing a natural, periodic network of efficient fluid pathways in an otherwise impermeable matrix (i.e. fractures). Such an outcome is of extreme importance the for exploration and extraction phases of unconventional reservoirs. Keywords Shale gas  Diagenetic reactions  Reservoir compaction  Channelling localisation instability

& Manolis Veveakis [email protected] 1

School of Petroleum Engineering, UNSW Australia, Sydney, Australia

2

CSIRO Mineral Resources, North Ryde, NSW, Australia

List of symbols b Arrhenius number for micromechanical processes ‘H McKenzie’s compaction length (m) k Hydro-mechanical coefficient g Chemo-mechanical coefficient a Thermal conductivity (W m1 K1 ) Dh Enthalpy of the reaction (J mol1 ) ij Strain tensor _ref Reference strain rate (s1 ) vp _d Deviatoric visco-plastic strain rate (s1 ) vp _v Volumetric visco-plastic strain rate (s1 ) ks Thermal expansion coefficient of solid (K1 ) lf Fluid viscosity (Pa s) q, [qs , qf ] Density [solid, fluid] (kg m3 ) rij , [r0ij ] Stress tensor [effective stress] (Pa) / Porosity e Cijkl Elasticity tensor (Pa) C Specific heat capacity (m2 K1 s2 ) E Activation energy (J mol1 ) Kc Ratio of forward over reverse preexponential constants M, [Mi ] Molar mass [of species i] (kg mol1 ) Q, [QF , QR ] Activation enthalpy [forward, reverse] (J mol1 ) R Universal gas constant (J K1 mol1 ) T Temperature (K) Pre-exponential factor [forward, reverse] k, [kþ , k ] (s1 ) kp Permeability (m2 ) p, [pc ] Volumetric mean stress [value at yield] (Pa)

123

S. Alevizos et al.

pf q, [qc ] r, [r þ , r  , Ri ] v, [vs , vf ]

Pore fluid pressure (Pa) Equivalent stress [value at yield] (Pa) Molar reaction rates [forward, reverse, species i] (mol m3 s1 ) Velocity [solid, fluid phase] (m s1 )

1 Introduction The ever increasing access to energy and mineral resources provides the foundation of our current market economy. All around the world, we are experiencing a shift from a coal-based to a gas-based primary energy supply. The International Energy Association (IEA) has therefore nominated this century as the golden age of gas (http:// www.worldenergyoutlook.org/goldenageofgas/). However, the extraction of conventional gas resources is limited compared to unconventional resources like shale gas in deep basins around the world. The expectation is that the classical geomechanics approach that has been successfully used to unlock the US shale gas in the past will be successful in other setting like Australia, for example. Unfortunately, in Australia there have been setbacks and a number of exploration wells have proven to be difficult. A prominent example is the Whicher Well in Western Australia where renewed activity has been started although the complexities and difficulties in this field essentially resulted in it being stranded since it was discovered in 1968. In the Cooper Basin, a successful unconventional shale gas plant is currently operating. The operators note the differences in shale composition (‘‘clay rich shallow overfilled lacustrine basin setting’’) and permeability network (‘‘different poroperm network dominated by diagenetically altered clays’’). However, the basis for the permeability network in such a high temperature, high pressure setting remains obscure. Also commonly accepted is the difference in stress field between Australia (compressive) and the USA (nominally neutral-extensional), putting the Australian plays squarely outside the conventional approach based on the US experience. As a consequence, conventional methods constitute a poor approach for tackling the fundamental problem of unconventional energy and mineral resources. The traditional geomechanics workflow is not sufficient to model processes occurring outside of the temperature, pressure and timescale window accessible to the laboratory (Kohlstedt and Holtzman 2009). The limited success of traditional stimulation techniques (hydraulic fracturing) applied to an unconventional geothermal reservoir has led to the drastic decline of the geothermal industry in Australia (Regenauer-Lieb et al. 2015). Surprisingly enough, even in these hostile settings where new fracture networks cannot be easily created, the empirical approach of the

123

industry revealed that the pre-existing fractures of the rock formation are accessible and can act as high-pressure pathways that allow for the exploitation of unconventional reservoirs. The reason for the presence of these fractures, as well as their performance and response to human intervention, is however poorly understood. Recently, two fundamental fluid-release instabilities have been discovered in geomaterials that are deforming in a ductile manner. The first instability is a volumetric collapse with a periodic feature of high-pressure fluid channels, so-called cnoidal waves (Veveakis and RegenauerLieb 2015; Veveakis et al. 2014b; Regenauer-Lieb et al. 2016). The second instability is an episodic fluid-assisted seismic slip event that is caused by shear deformation (Alevizos et al. 2014; Veveakis et al. 2014a; Poulet et al. 2014a, b). Both features can form in an otherwise impermeable ductile setting and create efficient, high-pressure fluid pathways. They form as material instabilities not dissimilar to the shear and compaction bands encountered in classical (brittle) solid mechanics (Rudnicki and Rice 1975; Issen and Rudnicki 2000). This work extends these previous studies to model the unconventional behaviour of geomaterials prone to fluid-release reactions under high pressure and temperature. Such fluid-release reactions can be encountered in shales in their diagenetic temperature window (100–200 C) (Fowler and Yang 2003). This temperature range corresponds to the depth that shale reservoirs are commonly located in Australia (3–5 km). In this temperature regime and particularly under compressive stress conditions, the dissolution of hydrated clay minerals (smectite) to anhydrous forms (illite) releases pore water that can lead to excess pore pressures (Fig. 1). The main physical change during this process is compaction, which drives pore water out of the material and therefore controls the permeability of the formation (Bachrach 2010) . Under normal circumstances in the absence of tectonic driving stress, the formation will tend towards thermodynamic equilibrium after compaction and the preferred mineral reaction is that of the reverse reaction, namely precipitation, also known as cementation. This reaction is exothermic hence does not need a driving stress. The normal situation therefore is a cementation of the reservoir as illustrated in Fig. 1. The wellknown cementation reaction renders nominally permeable formations into impermeable ones as shown in the example of sandstone in Fig. 1. However, this is only one branch of the dissolution/precipitation scheme which assumes equilibrium. The other, more interesting branch is for far-fromequilibrium conditions where a tectonic driving stress provides a steady energy supply for the endothermic dissolution reaction which releases fluids. The subject of this paper is the interplay between this tectonic environment

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability… Fig. 1 Effect of diagenesis in shales (left) and sand (middle column). The initial phases of deposition and packing is similar to both materials. Diagenesis occurs at the third step, with shales releasing water from their interlayer structure, thus maintaining locally increased porosity. Sand diagenesis on the other hand cements the matrix sealing its porosity (right column). Figure reproduced from Asveth (2003)

and a dissolution–precipitation reaction, and its role for the genesis of gas/oil-producing unconventional shale gas/oil reservoirs. In that particular context, diagenesis is of crucial importance for the oil and gas industry, and in nominally impermeable materials like shales it may offer new fluid pathways that are traditionally unexplored. Building on previous studies on the kinetics of shale diagenesis (Fowler and Yang 2003), we show here how a realistic model of diagenesis can be included in a compactive model of the shale formation, leading to channelling fluid instabilities at critical conditions where diagenesis under tectonic loading does not equilibrate in the cementation direction but rather the fluid-release direction.

2 Upscaled Model for Dissolution/Precipitation of Shales

Fowler and Yang (2003) highlighting the role of potasium feldspar KFs and its dissolution product to provide potassium ions combined with the aqueous silica compound formed from the smectite dissolution in pore water. A simple reaction for this scheme is to write (Fowler and Yang 2003): r

ð1Þ

2SðsolidÞ þ KFsðsolidÞ ! 2I þ 4QzðfluidÞ þ 9H2 O

where QzðfluidÞ represents the aqueous quartz (i.e. in solution) and H2 O the excess pore water. However, Eq. (1) is not a single reaction, but a compound one consisting of the processes of smectite dissolution, illite precipitation, Kfeldspar dissolution and quartz precipitation. They can be written as the following sequence of reactions (Fowler and Yang 2003): R1

ð2Þ

SðsolidÞ ! XðfluidÞ þ nH2 O R2

By way of introduction into our upscaling model, we recapitulate first some of the concepts of shale diagenesis presented in the fundamental work of Fowler and Yang (2003). In their analysis of dissolution/precipitation of shales, Fowler and Yang (2003) mention that diagenesis occurs in sedimentary basins when hydrated clay minerals such as smectite are converted to illite by a dehydration reaction. The process is facilitated by high temperatures and pressures, and the resultant water which is released invades the pore space. One consequence of this is that pore water pressures can be increased, thus leading to the phenomenon of over pressuring, which is of concern in oil drilling operations. The far-from-equilibrium dissolution (dehydration) reaction of smectite to illite has been thoroughly reviewed

1 KFsðsolidÞ ! Kþ fluid þ Al(OH)4ðfluidÞ þ sSiO2ðfluidÞ

ð3Þ

R3

1 Kþ ðfluidÞ þ Al(OH)4ðfluidÞ þ f Xfluid ! f Isolid þ SiO2ðfluidÞ

ð4Þ Rþ 4

SiO2ðfluidÞ  Qz R4

ð5Þ

The first reaction is the smectite dissolution, where S is smectite, X is an aqueous silica combination, like Si4 O10 (OH)2 ; the second reaction is the dissolution of potassium feldspar (KFs) to potassium and aluminium hydroxyl ions and aqueous silica SiO2 ; the third step represents the illite precipitation through combining the silica compound with potassium and aluminium hydroxyl ions; and the final step is the precipitation/dissolution of quartz

123

S. Alevizos et al.

(Qz) from/to the aqueous silica phase. The coefficients n, s, f are stoichiometric coefficients and Ri are the reaction rates of the various reactions. In order to describe the evolution of the species concentrations, we follow Fowler and Yang (2003) and denote by S, X, I, K, A, L, Q, W and F the concentrations of smectite (S), aqueous silica compound (X), illite (I), potassium ions (K), aluminium hydroxyl ions (A), aqueous silica (L), quartz (Q), water (W) and feldspar (F). The rates of concentration variation of the species are:

d½W ¼ n½SR2 dt

ð20Þ

For the macroscopic modelling purposes of the present work such detailed information is not required, and an equivalent (upscaled) single-step first-order dissolution/precipitation reaction needs to be obtained, being of the type rþ

ð21Þ

ABðsolidÞ

AðsolidÞ þ BðfluidÞ :  r

The forward reaction corresponds to the dissolution (dehydration) part, whereas the reverse corresponds to the precipitation (cementation) part of the smectite to illite transition. In this reaction, AB is incorporating all the solid reactants (like S and KFs), A is including all the solid products (like I) and B all the fluid products (like H2 O and Qz). In this case, the equivalent reaction rates would be

d½S ¼ ½SR1 dt

ð6Þ

d½X ¼ ½SR1 dt

ð7Þ

d½I ¼ f ½I½LR3 dt

ð8Þ

d½K ¼ ½FR2  ½K½A½Xf R3 dt

ð9Þ

d½A ¼ ½FR2  ½K½A½Xf R3 dt

ð10Þ

r  ¼CA CB k ðTÞ ¼ 2nf ½I½S½F2

d½L þ ¼ ½FR2 þ ½K½A½Xf R3 þ ½LðR 4  R4 Þ dt

ð11Þ

d½Q  ¼ ½LðRþ 4  R4 Þ dt

ð12Þ

d½F ¼ ½FR2 dt

ð13Þ

d½W ¼ n½SR2 dt

ð14Þ

where CAB ; CA ; CB the upscaled concentrations of AB, A, B. Therefore, the equivalent reaction rate of Eq. (21) can be obtained from the micromechanisms of shale diagenesis. Fowler and Yang (2003) analysed this system of reactions and provided analytical expressions for the Ri ’s. They involve complicated expressions of the physics and structure of the grain interface at the microscale, where the reactions take place (mineral surface area, grain size, volume fractions, etc.), as well as the relative reaction rates at this scale. In this work, we will assume equivalent (upscaled) reaction constants of the form   QF þ þ k ðTÞ ¼ R1 R2 ¼ k0 exp  ð24Þ RT   R3 R3 QR ð25Þ k ðTÞ ¼ 2nf þ 2  ¼ k0 exp  RT ðR4  R4 Þ

Note that in the limit where mass flux amongst the species ðaÞ

is not permitted, d=dt ¼ o=ot þ vk o=oxk represents the material time derivative of the corresponding a-species (a is any of S, X, I, K, A, L, Q, W and F). Invoking the steadystate assumption for the intermediate species (Law 2006), we set d½X d½A d½L ¼ ¼ ¼0 dt dt dt

ð15Þ

in the set of Eq. (16). We then solve for [L], [A] and substitute them to obtain d½S ¼  ½SR1 dt

ð16Þ

d½F ¼  ½FR2 dt

ð17Þ

d½I R2 R3 ¼ f ½I½F þ dt ðR4  R 4Þ

ð18Þ

d½Q ¼ 2½FR2 dt

ð19Þ

123

r þ ¼CAB kþ ðTÞ ¼ ½S½FR1 R2

ð22Þ R32 R3  R 4Þ

ðRþ 4

ð23Þ

where QF ; QR are the activation enthalpies of the forward (dissolution) and reverse (precipitation) reactions. Note that Qi ¼ Ei þ pViact with E the activation energy, p the mean pressure and V act the activation volume for internal volume changes. Similarly, the activation energy Ei incorporates the shear-induced microstructural mechanisms required for the macroscopic shear failure to initiate and can be further expanded as Ei ¼ Ei0 þ qX, where q is the deviatoric invariant of the stress tensor and X the activation volume for shear changes (Rice et al. 2001). It is therefore straightforward to define QF ¼ EF0 þ r0ij act ij :

ð26Þ

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability…

where the quantity r0ij act ij is the scalar product (double dot product) of two tensors, allowing for different activation volumes act for deviatoric- and volumetric-triggered ij reactions. Therefore, although not explicitly using the microscopic considerations of Fowler and Yang (2003), we consider an equivalent formulation and provide a direct link to the microstructure, showcasing the processes required to be quantified for evaluating the present model. Up to this point, we have not discussed the mechanical interactions but have laid down an upscaled version of the model of Fowler and Yang (2003) which provides an effective formulation of the chemical reactions. We emphasise that the upscaled formulation is valid for any set of chemical reactions involving dissolution/precipitation mechanisms.

3 Poromechanics We consider a 1D model formulation based on the method of characteristics developed in the general theory of plasticity (Hill 1950). This method reduces a generalised geometrical problem in two dimensions into an equivalent 1D failure line along a marching coordinate system n. These failure lines are constructed in a stress space that allows a decomposition of purely volumetric and purely deviatoric components acting across and along them, respectively. For more details on the rescaling along the characteristics refer to Veveakis and Regenauer-Lieb (2015). In triaxial or oedometric conditions for example (where in an x  y  z Cartesian coordinate system xx ¼ zz ¼ 0 and yy ¼ 6 0), the direction of n coincides with the purely

Fig. 2 Direction of the marching coordinate system n in a a triaxial test, following the maximum compressive direction of the y-axis perpendicular to the failure lines (dark lines); b in the problem of the expansion of a thick cylinder where the n direction is a logarithmic

compactive/dilatational y-direction (see Fig. 2a). In the problem of the expansion of a thick cylinder the slip lines are logarithmic spirals, and the n direction is also a logarithmic spiral (see Fig. 2b). Volumetric instabilities are therefore expected in the same orientation as the classical approach by Arthur et al. (1977) and Vardoulakis (1980). For a reservoir under compaction, we therefore expect similar response as a material with triaxial loading at the compressive cap of the yield surface (see Fig. 2a) or the logarithmic spirals emerging during the expansion of a thick cylinder (see Fig. 2b). In this contribution, we investigate the influence of the volumetric mechanism acting along the axis normal to the instability (n) and solve for the distribution of the mean effective stress p0 inside the compressed specimen of height 2H, defined as the area that reaches the yield stress. For mathematical simplicity, we assume that all the material properties, including solid viscosity and host permeability, are constant. Based on this approach, we position ourselves in the plastified area of height 2H where the initial elasto-plastic loading has occurred and therefore all the stress quantities are exceeding the yield regime, in a classical setting of overstress visco-plasticity (Perzyna 1966). Stress equilibrium in the n-direction combined with the stress decomposition p ¼ p0 þ pf (pf is the pressure of the fluid and p0 is the mean effective stress, using the convention of positive compressive stresses) provides op0 opf ¼ on on

ð27Þ

When this equation is combined with the mass balance considerations for the fluid–solid mixture, the following equations is obtained (Veveakis et al. 2014b)

spiral running perpendicular to the failure lines. Figures reproduced from Baxevanis et al. (2006) and Papamichos et al. (2001), see also Veveakis and Regenauer-Lieb (2015)

123

S. Alevizos et al.

  k p o2 pf 1 1 þ _ ¼    r V qf qs lf on2

ð28Þ

where kp is the constant permeability, lf the fluid viscosity, qf , qs the densities of the fluid and solid constituents, respectively, and r the rate of mass exchange of the diagenetic reaction given by Eq. (22). The theory has been detailed in Veveakis et al. (2014b) for melt as the saturating fluid, and for convenience it is summarised in Appendix 1. As already discussed in the previous sections, we decompose the volumetric strain rate into an elastic and an inelastic part, and consider isothermal conditions. Under these conditions, the problem is formulated as, _v ¼ _ev þ _vp v

p_ 0 ¼ þ _vp v K

ð29Þ

where K is the elastic bulk modulus. For the visco-plastic part of the strain rate a standard overstress plasticity approach is followed (Perzyna 1966). To this end, a power law rheology is assigned (see also Hickman and Gutierez 2007; Oka et al. 2011 and references therein):  0  p  p0c m _ ð30Þ _vp ¼  ; 0 v rref

bifurcation (Hill 1962,) we will hereby emphasise on the stationary attractor of this elasto-visco-plastic problem (i.e. at the limit of rigid visco-plastic material) and investigate the response of the matrix to volumetric deformation (Veveakis and Regenauer-Lieb 2015; Veveakis et al. 2014b, for further details). To this end, in the forthcoming sections we will perform a linear stability analysis of the problem at hand, before studying in depth the possible bifurcations from the homogeneous deformation state at the steady-state limit of the governing equations.

4 Analysis of the Main Model Parameters Following the previous considerations, we may combine Eqs. (27), (28), (29) and (31) to obtain the final equation (see also Sect. 1):   1 op0 kp o2 pf 1 1 þ ð32Þ _ ¼   þ  r V K ot qf qs lf on2 Equation (32) is brought into a dimensionless form: or0 d2 r0 0 ¼ 2  kr0m þ gebr ; os dz

ð33Þ

H

where _0 ¼ _ref eRT is the rate coefficient, with H the activation enthalpy for any micromechanical processes, here taken constant because of the isothermal assumption invoked (see also Poulet and Veveakis 2016). In the present work, since we restrict the analysis on the isothermal limit, we will assume that the rate coefficient is constant. Furthermore, the loading strain rate at the boundary where h 0 0 im p p equivalent loading p0n is applied, is thus _n ¼ _0 nrref c . If we accept rref ¼ p0n  p0c , then _0 ¼ _n and the constitutive model is written as  0  p  p0c m vp ð31Þ : _v ¼ _n 0 pn  p0c where we assume that the Macauley brackets are active, i.e. the material is in the post-yield regime. Note that the kinematic quantities _n and _v , as well as the static fields p0 and p0n must always have the same sign, in order to ensure positive mechanical work and satisfy the second law of thermodynamics. This means that the present formulation is invariant of the sign convention of the static and kinematic fields and applies to purely volumetric compaction and dilatational problems. A complete solution to the problem considers the strain rate to be decomposed into elastic, plastic and viscoplastic components and makes it necessary to calculate the elastic transients of the stress and porosity wave propagation. Following the classical approach of material

123

by assuming that the mean effective yield stress is constant and considering s ¼ Hkp2 Kl t, z ¼ Hn (H is the length of the f

p0 p0

reservoir in the compression direction), r0 ¼ p0 pc0 and n

 2 H _n lf 2 k¼ H ¼ ; ‘H kp p0n   þ 1 1 k0 lf H 2 QF0 g ¼ qAB e RT ;  qf qs MAB kp jp0n j QF0 ¼ b¼ 

c

ð34Þ

EF0 þ p0c act V ; 0 act jpn jV RT

Note that all the dimensionless groups are normalised so that they are not influenced by the stress sign convention used. The material parameters driving the response of the system are therefore the rate sensitivity exponent m and the three dimensionless groups of Eq. (34), namely k, g and b. Note that k scales with a modified compaction length that 0

p includes the adjusted matrix viscosity ls ¼ _nn , sffiffiffiffiffiffiffiffiffi kp ls : ‘H ¼ lf

ð35Þ

4.1 The Dimensionless Groups The parameter k is a hydro(poro)-mechanical group, representing the ratio of the mechanical matrix deformation

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability…

diffusivity over the internal mass diffusive transfer of the fluid. This parameter has a profound physical meaning as it states that volumetric instabilities are a result of competition of two time-dependent processes, which are the mechanical deformation of the matrix and the internal response of the embedded fluid phase. We may anticipate that stable deformation occurs when k  1, i.e. when the matrix deformation is much slower than the fluid diffusion rate and the specimen has the time to diffuse away any fluid pressure variations induced by the loading conditions. At the other extreme, when k  1, the loading rate is much faster than the fluid diffusion rate and internal mass variations cannot be equilibrated. In that case, coupled poromechanical instabilities are expected. As an illustrative example, we consider here the case of a cylindrical reservoir (2 km height) at an intermediate depth of 5 km in a tectonically deforming triaxial-like compressive environment along the direction of the maximum (absolute) principal stress r1 in the direction of the cylinder’s axis. We assume a tectonic loading of approximately 10 cm year1, resulting to a loading rate of _n  5 1012 s1 in the r1 direction. For simplicity in the rational, we assume a shale material that is being adequately represented by a modified cam 2 clay yield envelope, Mq þpðp  p0 Þ ¼ 0 with M ¼ 0:6, and a pre-consolidation pressure of around p0 = 70 MPa. The case where r1 = 110 MPa and r2 ¼ r3 = 70 MPa is a scenario where the reservoir experiences compressive tectonic loading that deviates it from the isotropic case r1 ¼ r2 ¼ r3 = 75 MPa, for an indicative density of 2500 kg m3 . In this scenario, the loading pressure is pn = 120 MPa. For a fluid density of 1000 kg m3 the applied effective pressure is 50 MPa, providing an effective loading stress p0n = 70 MPa and the experienced overstress at this stress state p0n  p0c = 10 MPa. These values, together with a standard fluid viscosity for the gas at this depth of lf ¼ 105 Pa s, allow to estimate the value of k as k

3 1018 kp

ð36Þ

We therefore deduce that the regime k  1 in which instabilities are expected is promoted only under extremely low-background-permeability formations, putting shale plays where the permeability is estimated well below 1018 m2 (i.e. at the nano-darcy level) as potential candidates for the interplay described in the following sections to operate. Similarly to k, g is a hydro-chemical group, representing the ratio of the chemical diffusivity with the diffusivity of the fluid in the porous matrix. It can be correlated to k through the following considerations:

g ¼ g0 k g0 ¼

  qAB 1 1 þ QF0  k0 _n e RT MAB qf qs

ð37Þ

Finally, the dimensionless group b represents the extend of the internal mass production mechanism under the given external loading. In this case where the mass production mechanism is a diagenetic reaction, b is a chemo-mechanical parameter determining the interplay between the tectonic loading and the mass exchange rate. Note that b includes the effect of the activation volume of the forward reaction, act V . However, unlike the hydro-mechanical parameters that can be relatively well constrained, the chemical parameters may vary quite significantly in the literature. Both the reaction rate constant k0þ (usually between 105 –1015 s1 ) and the activation enthalpy QF0 (usually between 200– 500 kJ mol-1) can vary several orders of magnitude, making the estimation of the value of g difficult. However, we may notice from Eq. (37) that within these ranges of parameters g can be significantly smaller, comparable or larger than k. For the indicative values presented here, g0 is extremely small, ranging as g0 104 to 1015 . It therefore makes sense to consider the reaction term negligible, setting g ¼ 0. However, recent studies (Stefanou and Sulem 2014) have emphasised the importance of this term, the presence of which determines the thickness of the localisation features (in this case compaction bands or fluid channels), providing the necessary internal length to regularise an otherwise ill-posed problem. The role of g therefore needs to be thoroughly investigated, and such an analysis is the subject of the remaining of the paper. 4.2 Linear Stability Analysis We start by studying the effect of the different physics on the quality of a steady-state solution r00 of Eq. (33) which is constant. As it is evident, such a solution exists only when g [ 0. As a next step, we perturb the stress around this solution, by inserting in Eq. (33) the expression r0 ¼ r00 þ r~0 . In this expression, r~0 is a stress perturbation, usually obeying a Lyapunov type law of the form r~0 ¼ Reikzþss

ð38Þ

where R is a constant, k is a wave number and s is the Lyapunov exponent (i is the imaginary unit). By linearising Eq. (33) around the steady state r00 and inserting the perturbation of Eq. (38), we obtain the characteristic polynomial of the equation: 0

0

s ¼ k2 þ mkr0m1  gbebr0

ð39Þ

The Lyapunov exponent s determines the stability of the system, maintaining stability when s\0. Since the sign of

123

S. Alevizos et al.

the Arrhenius number for the micromechanical processes b in not restricted, the exponent s can in principle admit positive values. The system becomes critically unstable when s ¼ 0, which corresponds to positive wave numbers of the form: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 0 ð40Þ k ¼ mkr0m1 þ gbebr0 Since k [ 0 we notice that the wave number can become 0 imaginary for small values of the term gbebr0 leading the problem to mathematical ill-posedness, since the absence of positive real wave numbers (i.e. wavelengths) signifies that the instability cannot be arrested in finite widths. This 0 is prevented when the reaction term gbebr is positive and significant to overcome the effect of the mechanical deformation. In this case, the instability is arrested and localised zones of finite width emerge, in accordance with the results of Stefanou and Sulem (2014). It therefore seems that despite being a small term, g—and therefore the corresponding physical mechanism of the diagenetic reaction—can provide the necessary internal length scale to regularise the problem. The above analysis was performed for an arbitrary constant stress profile across the sample, in order to highlight the regularising effect of the diagenetic reaction. As such, it does not provide information on the location and the pattern of potential localisation instabilities. This information is obtained by studying the stationary wave limit, i.e. the steady state of the governing Eq. (33) in the case where variations in space are allowed. Because of the small magnitude of the mass production term, in the following sections we will start by indeed neglecting it, and study the steady-state response of the system. In a later

(a) Fig. 3 Distribution of the normalised effective stress insidepffi halfpffiof kn  kn pffi ) the specimen, by symmetry, a in the case of m ¼ 1 (r0 ¼ eepffik þe þe k for various values of k. As anticipated, no singularities are observed. b in the case of m ¼ 2, for 2 values of k above kcr . This case presents similar response to the m ¼ 3 one, as anticipated, since the

123

step, the mass exchange rate will be re-introduced in the equation and its role in the system’s response will be thoroughly discussed.

5 First Approach: Weak Mass Exchange Rate We start by assuming that the mass production rate of Eq. (33) is small compared to the mechanical response, g ! 0, to obtain the final equation d2 r0 0  kr m ¼ 0; 2 dz

ð41Þ

5.1 Analytical Solution The solution of Eq. (41) depends on the value of the rate sensitivity coefficient m. For m ¼ 1 the system degenerates into the classical McKenzie equation (McKenzie 1984) without instabilities (Fig. 3a). For all m [ 1 (nonlinear cases), the solution for r0 is nontrivial, as it is elliptic. For odd values, it is the Jacobi and for even values it is the Weierstrass equation (the two are linearly related, as discussed in Veveakis and Regenauer-Lieb (2015)). This behaviour is frequently met in Korteweg–de Vries equation of the shallow water theory, with the solutions presenting periodic singularities that are spectrally stable (Veveakis et al. 2014b). This means that the transient orbits of the system are attracted to the solutions of the steady-state equation (41). For integer values of m\4 this equation admits closedform solutions. Without loss of generality, we focus in this study on the relevant solutions for a representative power law of m ¼ 3 (see also Eq. 1a of Rabinowicz and

(b) singularities increase in number monotonically with k. However, the complexity of the } function allows only for implicit solution with respect to k, and hence, providing a critical value of k for instabilities was not possible. For this reason, the case m ¼ 2 was not further pursued in the present work

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability…

Vigneresse 2004) because this is the most common power law exponent for earth materials. For m ¼ 3 and for con0 stant boundary conditions (r0 ð1Þ ¼ 1 and dr dz jz¼0 ¼ 0), the analytic solution of Eq. (41) is " rffiffiffiffiffiffiffi # ! k Icnð0; ıÞ zþ r0 ¼ C2 sn ð42Þ C2 ; ı 2 C2 where sn and Icn are the Jacobi SN and the inverse Jacobi CN function, and C2 the solution of the transcendental equation r0 ð1Þ ¼ 1. 5.2 Visco-plastic Instability Criterion The profiles of the normalised effective stress depend on the poro-mechanical feedback parameter k. For m ¼ 3 (see Fig. 3 for other cases), Fig. 4 depicts a complex response, providing a multiplicity of singularities for the normalised effective stress as k increases. As expected from the abovediscussed rate competition of diffusion process for small values of k (approximately k\13), the effective stress

(a)

profiles present a smooth solution with a minimum at the origin (z ¼ 0), as shown in Fig. 4a, b. For values k [ kcr  13 the effective stress presents multiple singularities, the number of which increases with k. Note that, for the example considered in Eq. (36), k [ 13 corresponds to formations with permeability range of kp \ 2:5 1019 m2 . Since these singularities are localised in space, they indicate zones of fluid flow barriers. The derivation of Appendix 2, also shown in the next section, shows that these singularities are low-porosity fluid channels, otherwise known as compaction bands. In the case where there is no reaction affecting the volume change of the mixture, we therefore confirm that the mechanics promote pore collapse and localisation of the volumetric deformation in the form of compaction bands. A diagnostic element for geological applications in the field is the spacing h between the compaction bands, as annotated in Fig. 4d. Although not periodically placed across the specimen, the compaction bands (CBs) divide the space into equal layers of distance h. As a consequence, h is then defined as the inverse density of the

(b)

(c)

(e)

(f)

h/H h1/H

h2/H

(d) Fig. 4 Distribution of the normalised effective stress r0 (Eq. 42) inside the specimen for six different values of k. In d the dimensionless distance h/H between the stress singularities is highlighted. In the stationary solution, both signs of the plastic

P-wave (dilational and compactive) can form identical channels of the fluid at approximately right angles to each other, forming a rectangular channel pattern when superposed in the purely volumetric instability case

123

S. Alevizos et al. 10

NC

8

25

N C = 0.26 λ

20

6

h≈ 4 H

2 0

0

500

λ

1000

H = H N C 0.26

=

15 10

kπ ⋅ Δσ 'n εn ⋅μ f

5

1500

0 0.0

Fig. 5 Number of channels NC as a function of k. From the definition of k, we deduce that the spacing between the stress singularities h scales with the characteristic hydraulic length ‘H , providing Eq. (43)

0.2

0.4

0.6

0.8

1.0

Fig. 6 Plot of the stress solution of Eq. (44). The addition of the chemical reaction term regularises the problem by providing fluid channels of finite thickness (highlighted in grey). The result is plotted for k ¼ 100, b ¼ 1 and g ¼ 104

CBs, h ¼ NHC . Through the results of Fig. 5, we obtain pffiffiffi NC ¼ 0:26 k, or rffiffiffiffiffiffiffiffiffiffi H ‘H l ð43Þ h¼  4 kp s : ¼ NC 0:26 lf

6 Strong Mass Exchange Rate: The Effect of Shale Diagenesis The previous considerations lead to a solution that is lacking an internal length regularising the thickness of the channels. This is manifested through the stress singularities of the solution presented in Fig. 4 and the localisation at zones of zero thickness. If the chemical term is included, the steady-state form of the final Eq. (41) reads d2 r0 0 0  kr m þ gebr ¼ 0; 2 dz

ð44Þ

This equation can only be solved numerically, providing an additional length scale that regularises the thickness of the fluid channels according to Eq. (40), as discussed in Section 4.2. The steady-state solution for the stress in this case is shown in Fig. 6. It needs to be noted that the presence of the additional chemical term does not alter the bifurcation criterion k [ kcr , or the scaling laws of the domain outside the channels. In order to showcase this, we perform a series of simulations exploring the range g; b 1. The results of these simulations are summarised in Fig. 7, where we notice that the scaling obtained for the number of channels as a function of k is the same as in Fig. 5 for g ¼ 0. We therefore conclude that the system described involves two internal length scales. The competition of hydro-mechanical diffusion processes—expressed through k—controls the

123

Fig. 7 Number of channels NC as a function of k. The result is plotted for simulations in the range b 1 and g 1. The blue line corresponds to the scaling curve obtained in Fig. 5 for g ¼ 0, pffiffiffi NC ¼ 0:27 k. Therefore, the addition of the mass production term does not affect the scaling of the outer length scale of the domain, the distance between the channels

outer length scale determining the distance between the channels. In addition, the hydro-chemo-mechanical diffusion processes—expressed through g and b—control the inner length scale determining the thickness of the channels and regularising the mathematical problem, as initially suggested recently by Stefanou and Sulem (2014). This result is of particular importance for the modelling community, since it introduces an internal length on the volumetric component of deformation which traditionally cannot admit any. Classical higher-order continua, like the Cosserat or micromorphic, admit this internal length only on the deviatoric component of the constitutive description (Veveakis et al. 2012, 2013; Sulem et al. 2011; Papanicolopulos and Veveakis 2011). For an in-depth discussion

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability…

on the regularisation mechanisms in shear and compaction bands, see Sulem and Stefanou (2016).

7 Porosity Evolution in the Localisation Bands The diagenetic reaction considered in this work (Eq. 21) takes place with a significant volume change, since the decomposition (forward) reaction is producing excess fluid and thus pore space. To model this volume change, we first decompose the total (Eulerian) porosity / to an initial contribution /0 due to the evolution of the interconnected pore volume, and a mechanical (D/mech ) and chemical (D/chem ) contribution. / ¼ /0 þ D/mech þ D/chem ¼

VB ; V

ð45Þ

with VB the volume occupied by fluid B. Following the analysis of Appendix 1, the final expression for the porosity can be obtained by using Eqs. (26), (45) and (58) 0

/ ¼ /0 þ ð1  /0 Þr m Dt  0 þ ð1  /0 Þ 1  r m Dt

0

ð46Þ

gsol eDbr 1 þ gprod eDbr0

where the key dimensionless groups are, gsol ¼ Kc ¼

qAB MB q MA Kc ; gprod ¼ B Kc ; qB MAB qA MB

ð47aÞ

k0þ ðArþDbÞ e k0

ð47bÞ

act with Ar ¼ ðEF  ER Þ=RT, and Db ¼ ðact VF  VR Þ=RT the difference between the activation volumes of the forward and the reverse reactions. We note that gsol and gprod are concentration constants that are of the order of Kc (usually gsol =gprod Oð101 Þ), which is the equilibrium constant of

the diagenetic reaction (ratio of forward to reverse frequency factors). When Kc ! 0 then the rate of the reverse (cementation) reaction is significantly larger than the rate of the forward. In this case, gsol ¼ gprod ¼ 0 and the total porosity is only influenced by the mechanical deformation. In this case, the localisation pattern corresponds to low-porosity compaction bands, dominated by pore collapse. This profile is shown in Fig. 8 for arbitrary values of the parameters. At the other extreme, when Kc ! 1, the rate of the forward (decomposition) reaction is significantly larger than the rate of the reverse. Eq. (46) denotes that the chemical reaction has an opposite effect to the porosity than the mechanical deformation, thus defining two competing mechanisms dominating the porosity distribution. Based on this, we expect the porosity profile to shift from low-porosity compaction bands to open channels under certain values of the equilibrium constant Kc . Figures 9 and 10 verify this claim, showing that with increasing values of Kc the porosity distribution departs from the pore collapse profile of Fig. 8, reaching a nearequilibrium profile for a value of Kc (Fig. 9). Upon further increase of Kc , porosity channels emerge at the same spatial configuration that compaction bands where forming before (Fig. 10). It should be noted that the behaviour of these solutions is affected highly by chemical coefficients. As it is shown in Fig. 11, the value of the critical relative reaction rate Kc that controls the equilibrium between the hydro-mechanical and the chemo-mechanical effects is indeed a function of the relative Arrhenius of the forward and the reverse reaction Db. We note that the thickness of the channels shown in Fig. 11 are predominantly controlled by the relative properties of the chemical reactions (forward and reverse). The above results showcase the importance of diagenetic reactions in a compacting shale (or in a material

0.010

0.012

0.009

0.011

0.008

0.010

0.007

0.009

0.006

0.008

0.005

0.007

0.004

0.006 0.0

0.2

0.4

0.6

0.8

1.0

Fig. 8 Porosity decrease in areas of stress concentration. _V0 Dt ¼ 4 107 , Kc ¼ 0, /0 ¼ 0:01

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 9 Porosity equilibrium in areas of stress concentration. _V0 Dt ¼ 4 107 , Kc ¼ 1014 , Db ¼ 1, /0 ¼ 0:01

123

S. Alevizos et al. 0.35

0.0100

0.30

0.0099

0.25

0.0098

0.20

0.0097

0.15

0.0096

0.10

0.0095

0.05

0.0094

0.00

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

0.0104

Fig. 10 Porosity channels in areas of stress concentration. _V0 Dt ¼ 4 107 , Kc ¼ 104 , Db ¼ 1, /0 ¼ 0:01

prone to diagenesis in general), being able to create fluid channels in an otherwise impermeable matrix even under compressive loading. The spatial distribution of the channels are determined by the hydro-mechanical interplay between mechanical deformation and fluid diffusion (see Sect. 5), and the thickness of the channels together with the magnitude of the fluid release is determined by the chemical interplay between the rates of the chemical reactions (forward vs. reverse).

0.0102

0.0100

0.0098

0.0096

0.35 0.30

8 Discussion and Conclusions Unconventional resources (coalbed methane, shale oil and gas, tight gas, heavy oil tar sands, gas hydrates, mineral deposits. etc.) are considered a promising future for humankind as their global significance surmounts that of the known conventional resources. They are unconventional resources in the sense of being in volatile conditions, deeper and hotter than ever reached before (e.g. deep geothermal, deep in situ leaching, shale gas in Australia) and in an environment that is nominally impermeable and therefore difficult to extract. Traditional modelling techniques of fluid flow through porous rocks are reaching their validity threshold at this range of temperature, pressure and low permeability. Still, the fact that even under these conditions shale gas reservoirs can be operated suggests that there is an overlooked mechanism that governs the response of the fluid flow problem at these conditions. In this contribution, we have extended the theory of fluid flow through porous media to the realm of chemo-plastic, nominally impermeable soft materials. Chemo-plasticity is explicitly studied as two counteracting mechanisms; the visco-plastic mechanical deformation of the medium that tends to reduce the porosity and the diagenetic fluid-release

123

0.25 0.20 0.15 0.10 0.05 0.00

0.0

0.2

0.4

0.6

Fig. 11 Porosity profiles for varying values of the relative Arrhenius number Db. Db ¼ 0:1; 1; 10 from top to bottom, respectively. _V0 Dt ¼ 4 107 , Kc ¼ 1014 , /0 ¼ 0:01

reactions that tend to increase porosity. We show that for a reservoir under compaction, there exist certain ambient and permeability conditions at which periodic stress concentrations can be encountered. In the absence of fluid-release reactions, mechanical deformation and/or cementation reactions tend to translate these stress concentrations into low-porosity compaction bands. However, in the presence of a diagenetic (fluid-release) reaction these stress singularities can provoke high-porosity (channelling) localisation instabilities. These channels are periodically

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability…

interspersed in the matrix and represent areas where the excess fluid from the reaction can be segregated at high velocity. These instabilities are favoured for extremely low-permeability formations and therefore provide efficient fluid pathways in an otherwise impermeable matrix, an outcome that could be very useful for the exploration and extraction phases of unconventional reservoirs. The critical condition for the onset of these periodic stress concentrations is given as a function of the poromechanical parameters of the problem, through the condition k¼

_n lf 2 H  kcr kp p0n

ð48Þ

The emergence of channelling instabilities or compaction bands is determined by the chemical balance between fluid-release and cementation reactions in the medium. It is expected that materials like sandstone that can only experience cementation diagenetic reactions (see Fig. 1) will always admit compaction bands under compression. However, materials like shales, that can experience fluid-release reactions, can admit the opposite manifestation of instability, i.e. channelling patterns. These channelling features are therefore promoted in environments with low permeability and high confining pressures (i.e. at the cap of the yield surface), conditions that are frequently met in shale reservoirs under compaction. This type of instability is classically unexpected and signifies a new response in the fluid–solid interaction that allows for fluid to flow through a nominally impermeable host rock formation. In a typical resource setting where the shale reservoir reached thermodynamic conditions favoured for the generation of the channelling instabilities, i.e. high enough temperature to trigger diagenesis and sufficient pressure and tectonic loading, the first of the reactions will be a fluid-release reaction. This reaction will be maintained as a (pervasive) diffusive front for low k values, which at a critical k value will give rise to channelling instabilities. Upon further strain or reduced loading, or depletion of the reactants, the channels will not be thermodynamically sustainable and the system will favour the reverse reaction of cementation. As explained in the previous section, this will ultimately lead to localised reduced porosity features (i.e. compaction bands) through cementation processes. In conclusion, we have shown the critical role of multiphysical feedbacks between the diagenetic (fluid-release) chemical reactions and the hydro-mechanical response of unconventional reservoirs. Critical channelling instabilities have been shown to emerge under certain chemo-hydromechanical conditions, allowing for effective high-velocity fluid transport through an otherwise impermeable host rock.

Acknowledgments The authors would like to thank the two anonymous reviewers for a thorough and in-depth review from which the manuscript has benefited greatly.

Appendix 1: Mass Balance In the bi-phasic setting considered, we define the partial densities qð1Þ ¼ ð1  /Þqs and qð2Þ ¼ /qf of the solid and fluid phase, respectively, / being the porosity (fluid content), qs and qf the solid skeleton and fluid density, respectively. The mass balance for each of the phases is oqðaÞ v

oqðaÞ ot

ðaÞ

þ on n ¼ 0, where a ¼ 1; 2 are superscripts for the solid and fluid phase, respectively. For incompressible solid and fluid (i.e. qs ; qf ¼ const.), the mass balance equation for each of the phases reduces to: ð1Þ



ð1Þ

ovn o/ o/vn m_ s  þ ¼ ; ot on on qs

ð49Þ

ð2Þ

o/ o/vn m_ f þ ¼ ot on qf

ð50Þ

where m_ f ¼ m_ s ¼ r so that the total mass is reserved (see also Alevizos et al. 2014). By adding them, we obtain the mixture’s mass balance equation   ð2Þ ð1Þ ð1Þ o/ðvn  vn Þ ovn 1 1 ð51Þ þ ¼  m_ f : qf qs on on In this work, we consider the case where porosity generation is large enough to breach the percolation threshold so that the continuum Darcy assumption holds. We therefore accept Darcy’s law for the separation velocity ð2Þ

ð1Þ

ð2Þ

ð1Þ

/ðvn  vn Þ, /ðvn  vn Þ ¼ 

kp opf lf on

ð52Þ

where kp is the constant permeability and lf the fluid ov

ð1Þ

viscosity. The volumetric strain rate is defined as _V ¼ onn . Under these considerations, and by considering constant permeability, we obtain Vardoulakis and Sulem (1995),   k p o2 pf 1 1 ð53Þ  _V þ  m_ f ¼ 0 qf qs lf on2

Appendix Porosity

2:

Chemico-mechanically

Induced

The mechanical porosity evolution can in principle be decomposed into two components, a reversible (elastic) one that is connected to the compressibility and thermal expansion of the mixture’s constituents and an irreversible (visco-plastic) part that is related to the volumetric strain

123

S. Alevizos et al.

rate. In the present framework, based on our emphasis on the initiation of localisation and therefore on the steadystate limit of the equations, we will assume that all mechanically induced porosity variations are due to the irreversible processes that act on the solid phase and thus, vp D/mech  ð1  /0 ÞDvp V ¼ ð1  /0 Þ_V Dt. If we take into consideration the visco-plastic law, Eq. 31, then the mechanically induced porosity variations can be written as:  0  p  p0c m 0 D/mech ¼ ð1  /0 Þ_n 0 Dt ¼ ð1  /0 Þr m Dt 0 pn  pc

where A/ is a coefficient that determines the amount of the interconnected pore volume (porosity) created due to the reaction. We assume that all the fluid generated contributes to the interconnected pore volume, and thus set A/ ¼ 1. We note that when the chemical reaction is endothermic (i.e. the forward decomposition part of (21)), DQ\0 and the system generates porosity. This can correspond to an unusual behaviour where porosity increases locally even at compressive environments, thus leading to what is commonly addressed in the melt segregation community as ‘‘decompaction bands’’ (Rabinowicz and Vigneresse 2004).

ð54Þ 0

p0 p0c p0n p0c

where r ¼ and Dt ¼ Dt _n are the normalised mean effective stress and normalised time increment (cumulative strain), respectively. In addition to the porosity, we define the solid ratio s¼

VA VA ; ¼ Vs ð1  /ÞV

ð55Þ

where V is a representative volume, VA and Vs are the volumes of solid phase A and all solid within V, respectively. The solid ratio is a measure of the extent of reaction (21). We assume the following relations for the partial molar reaction rates of the species involved   qAB xAB ¼  ð1  /Þð1  sÞ k0þ expðQF =RTÞ; ð56aÞ MAB   qA ð1  /Þs kA expðQR =RTÞ; ð56bÞ xA ¼ MA   qB ð56cÞ xB ¼ D/chem kB expðQR =RTÞ; MB where k is the pre-exponential factors (the assumption pffiffiffiffiffiffiffiffiffi k0 ¼ kA kB is used), q the density, M the molar mass and Q the activation enthalpies for the respective entities denoted by the subscripts: A, B, AB, as well as F and R for the forward and reverse reactions. Those reactions rates are linked by the stoichiometry of the considered reaction as xAB ¼ xA ¼ xB :

ð57Þ

From Eqs. (56-57), we may derive the poro-chemical model D/chem ¼ A/

1  /0  D/mech ; A 1 1 þ qqB M MB s

ð58aÞ

A

s¼

xrel ; 1 þ xrel

xrel ¼

  qAB MA DQ Kc exp  ; RT qA MAB

123

ð58bÞ ð58cÞ

References Alevizos S, Poulet T, Veveakis E (2014) Thermo-poro-mechanics of chemically active creeping faults. 1: theory and steady state considerations. J Geophys Res Solid Earth 119(6):4558–4582. doi:10.1002/2013JB010070 Arthur J, Dustan T, Al-Ani Q, Assadi A (1977) Plastic deformation and failure in granular media. Geotechnique 42:395–410 Asveth P (2003) Avo classification of lithology and pore fluids constrained by rock physics depth trends. Lead Edge 22(10):1004–1011. doi:10.1190/1.1623641 Bachrach R (2010) Elastic and resistivity anisotropy of shale during compaction and diagenesis: joint effective medium modeling and field observations. Geophysics 76(6):175–186 Baxevanis T, Papamichos E, Flornes O, Larsen I (2006) Compaction bands and induced permeability reduction in Tuffeau de Maastricht calcarenite. Acta Geotech 1:123–135. doi:10.1007/ s11440-006-0011-y Fowler AC, Yang X-S (2003) Dissolution/precipitation mechanisms for diagenesis in sedimentary basins. J Geophys Res Solid Earth. doi:10.1029/2002JB002269 Hickman R, Gutierez M (2007) Formulation of three-dimensional rate-dependent constitutive model for chalk and porous rocks. Int J Numer Anal Meth Geomech 31:583–605. doi:10.1002/nag.546 Hill R (1950) The mathematical theory of plasticity. Oxford University Press, London Hill R (1962) Acceleration waves in solids. J Mech Phys Solids 10:1–16 Issen K, Rudnicki J (2000) Conditions for compaction bands in porous rock. J Geophys Res 105:21529–21536 Kohlstedt DL, Holtzman BK (2009) Shearing melt out of the earth: an experimentalist’s perspective on the influence of deformation on melt extraction. Ann Rev Earth Planet Sci 37(1):561–593. doi:10.1146/annurev.earth.031208.100104 Law CK (ed) (2006) Combustion physics. Cambridge University Press, Cambridge McKenzie D (1984) The generation and compaction of partially molten rocks. J Petrol 25(3):713–765 Oka F, Kimoto S, Higo Y, Ohta H, Sanagawa T, Kodaka T (2011) An elasto-viscoplastic model for diatomaceous mudstone and numerical simulation of compaction bands. Int J Numer Anal Methods Geomech 35(2):244–263. doi:10.1002/nag.987 Papamichos E, Vardoulakis I, Tronvoll J, Skjærstein A (2001) Volumetric sand production model and experiment. Int J Numer Anal Methods Geomech 25(8):789–808. doi:10.1002/nag.154 Papanicolopulos S, Veveakis E (2011) Sliding and rolling dissipation in cosserat plasticity. Granul Matter 13(3):197–204

A Framework for Fracture Network Formation in Overpressurised Impermeable Shale: Deformability… Perzyna P (1966) Fundamental problems in viscoplasticity. Adv Appl Mech 9:243–377 Poulet T, Veveakis E, Regenauer-Lieb K, Yuen DA (2014a) Thermoporo-mechanics of chemically active creeping faults: 3. The role of serpentinite in episodic tremor and slip sequences, and transition to chaos. J Geophys Res Solid Earth 119(6):4606–4625. doi:10.1002/2014JB011004 Poulet T, Veveakis M (2016) A viscoplastic approach for pore collapse in saturated soft rocks using redback: an open-source parallel simulator for rock mechanics with dissipative feedbacks. Comput Geotech. doi:10.1016/j.compgeo.2015.12.015 Poulet T, Veveakis M, Herwegh M, Buckingham T, Regenauer-Lieb K (2014b) Modeling episodic fluid-release events in the ductile carbonates of the glarus thrust. Geophys Res Lett 41(20):7121–7128. doi:10.1002/2014GL061715 Rabinowicz M, Vigneresse J-L (2004) Melt segregation under compaction and shear channeling: application to granitic magma segregation in a continental crust. J Geophys Res. doi:10.1029/ 2002JB002372 Regenauer-Lieb K, Gaede O, Schrank C (2015) Deep geothermal: the ‘moon landing’ mission in the unconventional energy and minerals space. J Earth Sci 26(1):2–10 Regenauer-Lieb K, Poulet T, Veveakis M (2016) A novel wavemechanics approach for fluid flow in unconventional resources. Lead Edge 35(1):90–97. doi:10.1190/tle35010090.1 Rice JR, Lapusta N, Ranjith K (2001) Rate and state dependent friction and the stability of sliding between elastically deformable solids. J Mech Phys Solids 49(9):1865–1898. doi:10.1016/ S0022-5096(01)00042-4 Rudnicki JW, Rice JR (1975) Conditions for the localization of deformation in pressure sensitive dilatant materials. J Mech Phys Solids 23:371–394

Stefanou I, Sulem J (2014) Chemically induced compaction bands: triggering conditions and band thickness. J Geophys Res Solid Earth 119(2):880–899. doi:10.1002/2013JB010342 Sulem J, Stefanou I, Veveakis E (2011) Stability analysis of undrained adiabatic shearing of a rock layer with cosserat microstructure. Granul Matter 13(3):261–268. doi:10.1007/ s10035-010-0244-1 Sulem J, Stefanou I (2016) Thermal and chemical effects in shear and compaction bands. Geomech Energy Environ. doi:10.1016/j. gete.2015.12.004 Vardoulakis I (1980) Shear band inclination and shear modulus of sand in biaxial tests. Int J Numer Anal Meth Geomech 4(3):119 Vardoulakis IG, Sulem J (eds) (1995) Bifurcation analysis in geomechanics. Blankie Acc. and Professional Veveakis E, Regenauer-Lieb K (2015) Cnoidal waves in solids. J Mech Phys Solids 78:231–248 Veveakis E, Poulet T, Alevizos S (2014a) Thermo-poro-mechanics of chemically active creeping faults: 2. Transient considerations. J Geophys Res Solid Earth 119(6):4583–4605. doi:10.1002/ 2013JB010071 Veveakis E, Regenauer-Lieb K, Weinberg RF (2014b) Ductile compaction of partially molten rocks: the effect of non-linear viscous rheology on instability and segregation. Geophys J Int 200(1):519–523 Veveakis E, Stefanou I, Sulem J (2013) Failure in shear bands for granular materials: thermo-hydro-chemo-mechanical effects. Geotechn Lett 3(2):31–36 Veveakis E, Sulem J, Stefanou I (2012) Modeling of fault gouges with cosserat continuum mechanics: influence of thermal pressurization and chemical decomposition as coseismic weakening mechanisms. J Struct Geol 38:254–264. doi:10.1016/j.jsg.2011. 09.012

123