Microphysics of Inelastic Deformation in Reservoir Sandstones from the Seismogenic Center of the Groningen Gas Field

Physics-based assessment of the effects of hydrocarbon production from sandstone reservoirs on induced subsidence and seismicity hinges on understanding the processes governing compaction of the reservoir. Compaction strains are typically small (ε < 1%) and may be elastic (recoverable), or partly inelastic (permanent), as implied by recent experiments. To describe the inelastic contribution in the seismogenic Groningen gas field, a Cam–clay-type plasticity model was recently developed, based on the triaxial test data obtained for sandstones from the Groningen reservoir (strain rate ~ 10−5 s−1). To underpin the applicability of this model at production-driven strain rates (10−12 s−1), we develop a simplified microphysical model, based on the deformation mechanisms observed in triaxial experiments at in situ conditions and compaction strains (ε < 1%). These mechanisms include consolidation of and slip on µm-thick clay films within sandstone grain contacts, plus intragranular cracking. The mechanical behavior implied by this model agrees favourably with the experimental data and Cam–clay description of the sandstone behavior. At reservoir-relevant strains, the observed behavior is largely accounted for by consolidation of and slip on the intergranular clay films. A simple analysis shows that such clay film deformation is virtually time insensitive at current stresses in the Groningen reservoir, so that reservoir compaction by these mechanisms is also expected to be time insensitive. The Cam–clay model is accordingly anticipated to describe the main trends in compaction behavior at the decade time scales relevant to the field, although compaction strains and lateral stresses may be slightly underestimated due to other, smaller creep effects seen in experiments.


Introduction
A physics-based assessment of the short-and long-term effects of hydrocarbon production from sandstone reservoirs on induced subsidence and seismicity requires an understanding of the processes operating in the reservoir system (Mallman and Zoback 2007;Van Wees et al. 2018). It is generally understood that the increase in effective stress accompanying pore pressure (P p ) reduction leads to elastic and possibly inelastic compaction of the reservoir, resulting in small vertical strains (ε), totaling 0.1-1.0%, derived from surface subsidence (Morton et al. 2001;Mallman and Zoback 2007) or from in situ compaction measurements (Cannon and Kole 2017). Even though these strains are small, they may induce sufficiently high shear tractions on pre-existing faults to trigger seismogenic rupture (Mulders 2003;Bourne et al. 2014;Buijze et al. 2017). Constraining the magnitude of the in situ reservoir strain and particularly of the inelastic contribution is crucial for realistic geomechanical modeling of induced subsidence and seismicity Buijze et al. 2017;Van Wees et al. 2018;Candela et al. 2019;Smith et al. 2019). However, any inelastic sandstone deformation at the small strains (ε < 1%) relevant to induced reservoir compaction is frequently overlooked, since previous research (Wong and Baud 2012-and references therein) mainly focused on the deformation behavior seen 1 3 at (much) higher strains (ε of 1-10%). Moreover, extrapolation of lab-derived constraints on the reservoir compaction behavior and inelastic-elastic strain partitioning (Schutjens et al. 1995;Hol et al. 2018;Pijnenburg et al. 2018Pijnenburg et al. , 2019a, to field time scales and conditions, requires underpinning by a much better understanding of the governing processes (Spiers et al. 2017).
Experimental studies investigating the compaction behavior of sandstone widely report a characteristic, trimodal sequence in stress-strain behavior, usually expressed using plots of mean effective stress (here denoted P = t 1 + t 2 + t 3 ∕3−P p ) versus total porosity reduction (Δφ t ) (Wong and Baud 2012-and references therein), where σ 1 t , σ 2 t and σ 3 t represent the maximum, intermediate and minimum principal total stresses, respectively. Note that the notation used here is chosen to avoid cumbersome subscript stacking in later derivations. This trimodal sequence ( Fig. 1) typically includes an initial stage of nonlinear, concave-up P-Δφ t behavior (stage 1), which at high confining pressures (≥ 20 MPa) is followed by near-linear behavior (stage 2) and finally nonlinear, concave-down behavior, characterized by compaction (stage 3c). Alternatively, at low confining pressures (< 20 MPa) dilation (stage 3d) may set in after stage 1.
In previous work, the onset of stage 3d is often delineated by a yield envelope that is virtually insensitive to porosity , describing an approximately linear, Coulomb-type dependence of the differential stress (Q = σ 1 − σ 3 ) on P (Paterson and Wong 2005). In contrast, at the onset of stage 3c, the corresponding differential stresses are inversely dependent on P, strongly sensitive to porosity, and delineated by a yield envelop which may be either elliptical Baud et al. 2004;Skurtveit et al. 2013) or linear in a Q-P plot (Fortin et al. 2006;Guéguen and Fortin 2013). Stages 3c and 3d are both associated with brittle deformation behavior, involving localized shear fracturing (stage 3d), or pervasive intra-or transgranular cracking (stage 3c) (Wu et al. 2000;Baud et al. 2004;Fortin et al. 2005;DiGiovanni et al. 2007).
Following the main focus of the previous experimental studies, many microphysical models developed for inelastic deformation of sandstone are typically based on the brittle phenomena observed during the high strain (ε = 1-15%) stages 3d and 3c. Such models invoke combined frictional slip plus intergranular or intragranular cracking to describe dilatant stage 3d behavior (Ashby and Sammis 1990;Guéguen and Fortin 2013;Baud et al. 2014) or intragranular cracking in the case of stage 3c (Sammis and Ashby 1986;Zhang et al. 1990;Wong et al. 1997;Einav 2007a;Guéguen and Fortin 2013). Here, cracks are assumed to emanate from stress concentrators within the microstructure, such as grain boundaries, pores and other pre-existing flaws (e.g., Sammis and Ashby 1986), or from flaws at the periphery of assumed Hertzian grain contacts (Zhang et al. 1990;Wong et al. 1997;Brzesowsky et al. 2014). Although these models significantly enhance mechanistic understanding of the high strain stages 3d and 3c, they may not apply to the lower strains more relevant to producing reservoirs (ε ~ 0.1-1.0%). In stages 1 and 2, shear fracturing is not observed, while intragranular cracking is either absent, or far less common (Menéndez et al. 1996;Wu et al. 2000;Pijnenburg et al. 2019a).
In recent years, seismic activity in the vast (30 by 30 km) and densely populated Groningen gas field in the Netherlands has led to an urgent need to understand better the compaction behavior of the Slochteren reservoir sandstone at production-relevant strains (ε ≤ 1%) (Van Spiers et al. 2017;Van Wees et al. 2018). The latest experiments performed on the Slochteren sandstone (Pijnenburg et al. 2019a) demonstrate similar stages 1, 2, 3c and 3d behavior as typically reported for sandstone (cf. Fig. 1). The small-strain behavior (ε ≤ 1%) seen during stage 2 was found to be the most relevant to in situ compaction. At these small stage 2 strains, one-third, or up to one half of the total strain was found to be inelastic (Hol et al. 2018;Pijnenburg et al. 2018Pijnenburg et al. , 2019a. The stage 2, inelastic deformation behavior was shown to be well-described by Fig. 1 Schematic diagram showing the different stages in mean effective stress (P) versus total porosity reduction (Δφ t ) behavior typically seen in mechanical tests on porous sandstones. Note the effects of porosity and confining pressure on stage 3 behavior (dilation-suffix "d" versus compaction-"c"). The onset of stages 3d and 3c is indicated with "d" and "c", respectively. The deformation mechanisms inferred in the literature are shown per stage, as are typical values of Δφ t measured at the end of each stage a constitutive, Cam-clay-type plasticity model, implying isotropic strain-hardening accompanying inelastic porosity reduction, i.e. inelastic densification (Pijnenburg et al. 2019a). Application of this constitutive model to depletion of the Groningen field yielded similar stress evolution estimates to the values measured over the past 30 years in the reservoir, suggesting a controlling role of inelastic deformation on the in situ stress evolution (Pijnenburg et al. 2019a). However, a mechanistic basis for this empirical model and the underlying experimental data is still lacking, which hampers confident extrapolation of the model results to the decade time scales relevant to field production.
A recent microstructural study by Pijnenburg et al. (2019b) showed that stage 2, inelastic deformation of Slochteren sandstone is accompanied by consolidation of and slip on µm-thick, (intergranular) clay films present on the surfaces of and within the contacts between the sandstone grains, plus minor intragranular cracking. A rough estimate of the inelastic strain magnitude associated with clay film deformation (ε ≈ 0.3%), at the experimentally explored stress changes (σ 1 from 41 to 91 MPa, σ 3 = 41 MPa, P p = 1 MPa), showed similar results to the values measured during testing (0.4%). However, a more thorough microphysical description of the observed behavior is still lacking. Moreover, while the criteria for the onset of pervasive intragranular cracking have been established (Zhang et al. 1990;Wong et al. 1997;Guéguen and Fortin 2013), the role played by intragranular cracking in determining the σ-ε behavior seen in experiments has only been sparsely investigated. Where the effect of intragranular cracking on the overall σ-ε behavior has been explored (Einav 2007a;Marketos and Bolton 2009;Tengattini et al. 2014), the main focus again dominantly lies on describing the high strain behavior not relevant to the Groningen gas reservoir (and other, similar reservoirs), while any influence on this of the intergranular clays present in the Slochteren sandstone has not yet been analyzed.
In this study, we develop and test a series of microphysical models addressing the inelastic deformation processes observed in the conventional triaxial compression experiments on Slochteren sandstone reported by Pijnenburg et al. (2019a). First, we evaluate more rigorously whether intergranular clay consolidation and slip can quantitatively account for the inelastic behavior governing the stage 1 and stage 2 behavior seen in these experiments. The possible effects of the presence of the intergranular clays on the higher strain stage 3d (i.e., dilatant, intergranular slip) and stage 3c (intragranular cracking) is also considered. We investigate model behavior for Slochteren sandstones with porosities of 13.4, 21.5 and 26.4%, as tested in the experiments.
This series of simplified models represents a first exploratory attempt to test whether the microphysical mechanisms described by Pijnenburg et al. (2019b) can (semi-)quantitatively account for the inelastic behavior of the Slochteren sandstone seen at the sample scale in our triaxial deformation experiments. Our hope is that, if so, the underlying physical processes will be explored in further detail in the future, using more rigorous approaches for linking grainscale to aggregate-scale behavior, i.e., by specialists in homogenization treatments (Bardet and Vardoulakis 2001;Fortin et al. 2003), granular mechanics modeling (Einav 2007a, b;Tengattini et al. 2014) and the discrete element method (Kawamoto et al. 2016;Van den Ende et al. 2018). Alternatively, if the processes explored here cannot account for the experimentally observed deformation behavior, then they should be rejected in playing a dominant role in the compaction behavior seen in experiments, and likely in the Groningen and perhaps other, similar reservoir rocks. In line with this, the present paper attempts to assess the physical mechanisms underlying the empirical elastic plus inelastic (Cam-clay type) compaction model proposed for the Groningen reservoir by (Pijnenburg et al. 2019a) as well as the applicability of this model to the reservoir conditions and time scales that apply in Groningen and possibly elsewhere.

Decoupling of Elastic and Inelastic Deformation
We adopt the convention that compressive stresses and strains are positive. In line with our definition of mean effective stress P above, all stresses (denoted σ) supported by the sandstone framework are henceforth understood to be Terzaghi effective values, i.e., total stress (σ t ) minus the pore pressure, to avoid unworkably cumbersome stacking of subscripts in derivations to come. A list of the symbols used is given in Table 1. Assumed values for several of these quantities are listed in Table 2. In general, the principal strain increments (dɛ i t , i = 1, 2 or 3) resulting from changes in the principal effective stress components (dσ i ) and/or time-dependent deformation are given as the sum of the elastic (dɛ i el ) and inelastic (dɛ i inel ) components, as: where d inel i includes both time-independent and timedependent inelastic strain contributions. In this study, we consider dɛ i inel to be time-independent and hence determined by changes in stress only. Since we consider only small, coaxial strains relevant for compaction of producing reservoirs, as studied in conventional triaxial tests (i.e., of order 0.1-1%) and measured with respect to well-defined zerostrain reference states, Eq. 1 also applies in the finite strain form t i = el i + inel i . In using both forms of Equation 1, we assume that the elastic and total inelastic contributions are fully decoupled, so that each can be evaluated independently of the other.

Quantifying the Elastic Contribution
The elastic strain contribution is quantified here using poroelasticity theory (Wang 2000). We assume that the compressibility of the grains constituting the Slochteren sandstone is negligible as compared to that of the grain framework, i.e., the Biot coefficient is assumed close or equal to 1, as estimated by Hettema et al. (2000) and  Pijnenburg et al. (2019a). We further assume that the sandstone can be treated as elastically isotropic. Poroelasticity then prescribes that the principal elastic strain increments (d el i ) are related to incremental changes in the principal effective stresses (dσ i ) and mean effective stress (dP = [dσ 1 + dσ 2 + dσ 3 ]/3) via: Here, G is the shear modulus and K is the bulk modulus of the sandstone. In the previous work (Pijnenburg et al. 2019a), it was shown that the elastic behavior of the Slochteren sandstone is nonlinear and that the values of K and of the Young's modulus (E) obtained for samples with initial porosity φ 0 = 13.4, 21.5 and 26.4%, increase with decreasing φ 0 and with increasing P, up to P ≈ 50 MPa, beyond which they are roughly constant. At low P < 30 MPa, the elastic behavior was shown to be slightly anisotropic, but this anisotropy was not apparent at higher P. For present purposes, we obtained an average value of K at P in the range of 5-50 MPa, using the following, linear relation between K and P: Here, C K 1 and C K 2 are constants obtained for each sample set tested (i.e., each initial porosity φ 0 ), through linear regression using the K-P data presented in Figure 5 of Pijnenburg et al. (2019a). The values obtained are listed in Table 2. The E values measured at each P and φ 0 were combined with Eq. 3a to obtain the corresponding shear moduli, using G = 3KE/(9K − E) (Wang 2000). This showed a similar P-and φ 0 dependence of G, as shown by K. For 5 ≤ P ≤ 50 MPa, the shear moduli are described through: where C G 1 and C G 2 are porosity-sensitive constants obtained through linear regression of the G-P data derived as described above. For values obtained, see Table 2. At P ≥ 50 MPa, the experimentally observed K and G values were found to be roughly constant (within 0.5 GPa), taking values of K ≈ 11.7, 9.4 and 6.1 GPa and G ≈ 7.5, 6.0 and 4.5 GPa, for initial porosities of 13.4, 21.5 and 26.4%, respectively.

Key Microstructural Observations
We start model development by summarizing key microstructural features observed in Slochteren sandstone samples from the Groningen field. These features are illustrated in Fig. 2, using scanning electron microscope graphs of samples used in the previous experimental studies on the Slochteren sandstone (Pijnenburg et al. 2019a, b). Porosities of undeformed Slochteren sandstone range between 10 and 27% (NAM 2013). The sandstone consists predominantly of quartz (72-90 vol%), with lesser amounts of feldspar (8-25 vol%) and clay (0.5-5.5 vol%) (Waldmann et al. 2014;Waldmann and Gaupp 2016). Quartz and feldspar grains are typically 200 ± 50 µm in size (Hol et al. 2015;Pijnenburg et al. 2018) and are flattened at grain-to-grain contacts. The clay is predominantly present in the pore space ( Fig. 2a) and is further found as thin films coating the pore walls and a high proportion of the grain contacts, where they constitute part of the load-bearing framework (Gaupp et al. 1993;Waldmann 2011). The clay present within grain contacts consists predominantly of illite and is inferred to have been largely present already prior to burial (Gaupp and (3a)  φ 0 denotes the initial porosity; z denotes the characteristic initial clay film thickness at each of these porosities, as was estimated from corresponding micrographs; C K 1 , C K 2 , C G 1 and C G 2 are empirically obtained constants describing the P-sensitivity of the bulk modulus (superscript K) or of the shear modulus (superscript G), shown in experiments (Pijnenburg et al. 2019a); r is the grain contact radius (see Eq. 4); R and R c refer to the radii of the quartz grain measured from grain center to pore wall and the grain contact, respectively. d denotes the assumed grain contact flaw spacing in 2-D, chosen to be equivalent to the characteristic interspacing of grain contact surface irregularities seen in micrographs; µ denotes the assumed friction coefficient of the intergranular illite; B is a constant describing the consolidation behavior of illite; E q is the Young's modulus of quartz grains; Y is the fracture energy of quartz and φ cr is the intragranular porosity fraction assumed to prevail after multi-edge cracking  Pijnenburg et al. (2019a, b) shows that these graincontact clay films are thicker in the more porous sandstones. Roughly estimating these thicknesses from the micrographs  Fig. 1). Both sample types were tested at a P p of 0.1 or 1 MPa, a strain rate of 10 −5 s −1 and/or a rate of P increase of 0.1 MPa/s, up to the maximum stress conditions indicated. a Clays are present in the pore space, on the pore walls and within graincontacts; b, c Clay films, including those within the grain contacts (inset) appear thicker at higher sample porosity, than at low porosity; d, e During stage 1, inelastic deformation is mostly accommodated by compression of and local shear on clay films. Similar, although more prominent clay deformation is seen in stage 2; f After stage 3d, sample-scale, shear fractures are often seen, which at low P c and high φ 0 (e.g., here) are mostly intergranular; g after stage 3c, intragranular cracks are ubiquitous, emanating from grain-to-grain contacts. h Cracked grains are either split by a single, meridional crack (yellow box), or else more pervasively crushed through multiple, roughly parallel, or convergent cracks (red arrow), referred to here as multi-edge cracks. Cracks typically open towards the pore space, while in more densely-packed (parts of) grains (white arrow), cracks are absent or closed; i within grain contacts, cracks are typically seen to emanate from µm-scale asperities; j close-up of h. The boundaries of the split grain indent into the clay film coating the contact (color figure online) recovered from the undeformed counterparts of the samples mechanically in Pijnenburg et al. (2019a), having porosities of 26.4% (Fig. 2b), 21.4-21.5% and 13.4% (Fig. 2c), reveals that these typically range between 6 ± 2 µm, 4 ± 2 µm and 2.5 ± 1.5 µm, respectively. Further quantification of (any porosity sensitivity of) the clay film thickness would be desirable, but is not attempted here, since any statistically meaningful analysis would require an elaborate, 3D, and high resolution (< 1 µm) microstructural study of a large number of grain contacts (e.g., Desbois et al. 2016), which is not feasible in the present study.
After experimental deformation of the three porosities of samples, a variety of permanent, inelastic microstructural changes were observed, the nature of which were shown to depend on the stage of P-Δφ t behavior explored in the experiment (refer Fig. 1). During stages 1 and 2, inelastic deformation was found to develop mostly through compression (i.e., consolidation) of, and local slip on the intergranular clay films, as is illustrated for stage 1 in Fig. 2d, e. Intragranular cracking played no role during stage 1, and only a small role during stage 2 (Pijnenburg et al. 2019a). Samples deformed into stage 3d often revealed one or more sample-scale fractures, oriented at an angle between 10 and 45 degrees to the axial compression direction (Fig. 2f). These fractures are dilated at the time of inspection, and are predominantly intergranular at low confining pressure (< 20 MPa) and high φ 0 (> 20%). At higher P and lower φ 0 , they are typically surrounded by a damage zone up to 1 mm wide, in which intragranular cracks are abundant. In samples deformed into stage 3c, such dilatant fractures are absent, while in these samples, pervasive intragranular cracking occurs within 1-10 mm wide bands, oriented at 45°-90° to the main compression direction. Where observed, intragranular cracks emanate from grain-to-grain contacts (Fig. 2g, h), frequently from, or opposing micrometer-sized contact asperities, which are typically interspaced by 3 ± 2 µm ( Fig. 2i, white arrows). Cracked grains are either split by a single, meridional crack ( Fig. 2h, yellow box), or else fragmented by multiple, roughly parallel, or convergent cracks, interspaced by several micrometers (Fig. 2g, h, red arrow). The latter type will be referred to as multi-edge cracks. Cracks are typically opened by 0.1-10 µm. Where markedly opened, the boundaries of the split grain, or grain fragments can be seen to indent into the clay film coating the contact (Fig. 2j). Crack opening is seen more readily in the direction of surrounding pores ( Fig. 2 h, red arrow), whereas fewer cracks and reduced crack opening are seen in more densely packed parts of the microstructure (Fig. 2h, white arrow).

Idealized Microstructure and Microstructural Unit Cell
We assume an idealized and highly simplified model microstructure for Slochteren sandstone undergoing inelastic deformation, as shown in Fig. 3. In the σ 1 -σ 3 plane, we assume a simple, 2D hexagonal pack of quasi-spherical quartz grains ( Fig. 3a) with flat, or truncated, circular grain contacts, coated by clay films (Fig. 3b, c-colored red). In the third dimension, or σ 2 direction, this 2D pack is assumed to be stacked grain-on-grain (Fig. 3a). The grain radius to the pore walls (R) is taken to be 100 µm ( Table 2). The grain contact radius (r) is approximated for a given initial porosity (φ 0 ) through (Renard et al. 1999;Gundersen et al. 2002): where q is a constant assumed to be equal to 0.8 (cf. Spiers et al. 2004), and k is the grain coordination number, here equal to 8. The distance between grain centers and flattened contacts is then given by R c = √(R 2 − r 2 ) (see Fig. 3c). We ignore clays not present within grain contacts, as these fall outside the load-bearing framework. Grain contact clay films are taken to consist of illite (see Sect. 4.1), and are discshaped, having the same top and bottom surface area as that of the contact (πr 2 ) and a thickness z of 3, 4 and 6 µm for porosities of 13.4, 21.5 and 26.4%, respectively (Table 2), following the clay film thicknesses estimated earlier for these porosities (Sect. 4.1). The unit cell of this regular pack is one grain deep (i.e., ~ 2R c ) and shown in Fig. 3. It incorporate four sets of two parallel grain contacts, referred to as grain contacts A, B, C and D (see Fig. 3b).
In the following, we first give general relations for the stresses acting on the unit cell and grain contact scales. We then develop a series of model components describing stress criteria and accompanying stress-strain relations for the above processes assumed to dominate inelastic deformation during stages 1 and 2 (Sect. 5), stage 3d (Sect. 6) and stage 3c (Sect. 7). Finally (Sect. 8), we integrate the different model components to evaluate the total inelastic deformation behavior that results.

Stresses at the Unit Cell and Grain Contact Scales
The maximum (σ 1 ) and minimum (σ 3 ) principal stresses act at an angle θ to the normal of the unit cell top and side, Fig. 3 Assumed model microstructure for inelastic deformation of the clay-bearing Slochteren sandstone undergoing triaxial compression. The microstructure constitutes a simple hexagonal pack of quasi-spherical quartz grains with flattened contacts, each covered with illite clay films (indicated red). a The simple hexagonal grain pack. The basal plane is oriented perpendicular to σ 2 . b Basal planeparallel view on the representative, unit cell volume, taken to be one grain (2R c ) deep. The unit cell contains grain contact types A, B, C and D; c The maximum vertical stress (σ 1 ) acts at an angle θ to the top of the unit cell, implying (external) normal and shear stresses on the top (σ top , τ top ) and side (σ side , τ side ) of the unit cell. These external stresses imply normal (̃) and shear stresses (̃) on grain contacts A, B and C which, if sufficiently large, may cause consolidation of-and/ or slip on the intergranular clay (thickness z) and/or cracking of the quartz grains. Grain contact D supports only ̃ , no ̃ . R and R c refer to the radii of the quartz grain measured from grain center to pore wall and the grain contact, respectively, while r is the grain contact radius (color figure online) respectively, while the intermediate principal stress (σ 2 ) acts perpendicular to σ 1 and σ 3 , hence perpendicular to grain contact D (Fig. 3a, b). By varying θ between 0° and 30°, grain contacts A, B and C are rotated such that all possible grain contact orientations (0°-90°) with respect to σ 1 and σ 3 are described. In a representative volume of a real, reasonably isotropic sandstone, all of these grain contact orientations, hence all angles θ, will be present in different parts of the microstructure. To account for this, we explore the deformation behavior of the unit cell for the two end-member angles θ of 0° and 30° (Fig. 3c). In doing so, we effectively explore the two extremes in mechanical behavior expected, where the behavior at intermediate angles is anticipated to be similarly intermediate. In addition, to explore the onset of dilatant intergranular slip and hence dilatant shear failure in a representative volume of isotropic sandstone (stage 3d), we examine the response of the unit cell for the full range of 0° ≤ θ ≤ 30°, since the value of (σ 1 -σ 3 ) required to cause shear slip on clay-filled grain contacts with a specific friction coefficient, at a given value of σ 3 , is strongly orientationdependent (cf. Jaeger et al. 2007). Assuming a homogeneous stress distribution throughout the sandstone microstructure, then at a given angle θ, the normal (σ top ) and shear stresses (τ top ) acting on the top and bottom planes of the unit cell ( Fig. 3c), as well as the normal (σ side ) and shear stresses (τ side ) acting on the lateral planes of the unit cell, are related to the principal stresses at the sandstone sample scale (Jaeger et al. 2007) through: These stresses, and the intermediate principal stress σ 2 are transmitted onto grain contacts A, B, C and D, producing shear (̃) and/or effective normal (̃) stresses on each contact. Ignoring the thickness of the intergranular clay films (i.e., taking z << R, as observed), force balance considerations imply that the stresses acting on the unit cell and grain contacts (Fig. 3c) are related by: Here, the subscripts to ̃ and ̃ denote the specific grain contact stress components (see Fig. 3c). We assume that these force balance relations are satisfied within each cell-sized grain cluster during deformation. Note that at a given external deviatoric stress state acting on the unit cell, Eqs. 6a-6e contain more unknown contact stresses than we have equations to prescribe them (7 unknowns, 5 equations). To solve for each grain contact stress, at a given external stress state, additional constraints specifying the stresses that can be supported by the clay films present within grain contacts are needed. Thus, the magnitudes of the contact stresses, for any given cell orientation θ, depend on the external loading conditions (σ i ), and on the capacity of the intergranular illite films to support shear and normal load. They will be evaluated in the following.

Consolidation Behavior of Illite Within Grain Contacts
The consolidation (i.e., compaction) behavior of illite powder has been investigated by Brown et al. (2017). These authors performed long-term (months), drained, 1D consolidation tests on discs (1-2 mL in volume) of brine-saturated illite powder (purity 65-85%) at room temperature, exploring effective axial stresses in the range of 5-180 MPa. They found that the void ratio hence the porosity of the illite decreased approximately logarithmically with increasing effective normal stress. Similar behavior has been reported for many other clays (Allman and Atkinson 1992;Smith et al. 1992;Marcial et al. 2002) and shales Ferrari et al. 2016), although generally with uncertainties in degree of drainage versus test duration. Assuming constant pore fluid density, the total pore volume reduction is in 1D compression equal to the fluid volume expelled, and thus gives axial strain when normalized with respect to a given reference disc thickness. We estimate the inelastic contribution to these data by subtracting the axial elastic strain expected for the corresponding changes in effective stress (i.e., Δ E ). We assume a constrained (oedometric) modulus E c of 10 GPa, as measured at an effective axial stress of 100 MPa, in wet Opalinus shales with porosities in the appropriate range (10-20%; Ferrari et al. 2016), and a Biot coefficient of 1. After subtracting the elastic strain contribution (< 10% of total), the residual inelastic strain due to clay consolidation is found to increase again roughly logarithmically with increasing effective normal stress (see: Supporting Information, Figure S1). Applying these data in our model, finite inelastic normal strain (̃c lay ) occurring within a grain contact illite film due to 1D (film normal) consolidation accompanying an increase in grain contact normal stress, from a reference value ̃0 to a value ̃ , is described by: which is, of course, the integral form of the incremental relation d̃c lay = B d̃ . Here, B is a constant estimated from the data treatment illustrated in the Supporting Information ( Figure S1) to be equal to about 0.05. The use of Eq. 7 assumes that inelastic compression of illite films in quartz grain contacts can be treated as 1D, drained consolidation behavior. Given the small radius of grain contacts (10-50 µm in the Slochteren sandstone), hence very short fluid expulsion/drainage path, and the frictional constraints imposed on the clay film surfaces by the enclosing quartz grain contacts, we consider this assumption to be reasonable.

Clay Film Response to Hydrostatic Loading and Effects at Unit Cell Scale
Under hydrostatic effective stress conditions, σ 1 = σ 2 = σ 3 = σ top = σ side = the effective confining pressure P c , while all shear stresses at the unit cell boundaries and grain contacts are zero. In this case, Eqs. 6a and 6c imply that the normal stresses acting on grain contacts A, B and C in Fig. 3c are given: At the same time, Eq. 6e implies that: which highlights contact loading anisotropy in the model geometry chosen. Despite this anisotropy, Eq. 7 plus 8a or 8b show that any increase in isotropic confining pressure to a value P c , from a reference confining pressure P c 0 results in the same inelastic strain response Δ̃c lay A,B,C,D at all contacts, i.e., at A through D. At the scale of the unit cell (Fig. 3c) and above, the principal inelastic strains developing due to clay consolidation during hydrostatic compression can now be written as: where z is the grain contact clay film thickness. Combining these relations with Eqs. 7, 8a and 8b, we obtain the following expressions for the inelastic strains produced by consolidation during an increase in isotropic effective confining pressure (Fig. 4a) from a value P c 0 to a value P c :

Criterion for the Onset of Serially Coupled, Intergranular Slip and Clay Film Consolidation
During deviatoric compression at θ = 0°, the model geometry dictates that any consolidation of the clays present at contact type B must be accompanied by slip and consolidation of the clay films present at A and C (Fig. 4b). During deviatoric loading at θ = 30°, consolidation of the clays at A and B requires slip along both of these contacts (Fig. 4c). Hence, during deviatoric loading at both orientations, consolidation and slip are serially coupled deformation processes, meaning that one process cannot occur without the other. Ignoring cohesion in the clay films, and noting the symmetry of the unit cell at θ = 0°, slip on grain contacts A and C during deviatoric compression at θ = 0° (cf. Fig. 4b) occurs when the criterion for frictional slip is obeyed, i.e., when: Here, µ is the friction coefficient of illite. Similarly, during deviatoric loading at θ = 30° (cf. Fig. 4c), A and B slip when: Wet shear tests performed on layers of relatively pure (95 vol%) illite clay powder (simulated frictional wear material or "gouge" found in upper crustal faults) at room temperature, an effective normal stress of 40 MPa, and a shear velocity of 1 µm/s have shown that the value of µ is strain dependent, increasing roughly linearly from 0.22 to 0.30 as the shear strain increases from 0.4 to 8.0, likely reflecting hardening due to densification and associated planar fabric development (Tembe et al. 2010). For present purposes, we use an intermediate µ value of 0.26, reported by Tembe et al. (2010) at a shear strain of 2.3. We later discuss to what extent the model behavior changes if a variable, shear strain-or consolidationstrain-dependent µ value is assumed.
At deviatoric stresses where Equations 10a and b are not satisfied, i.e., where ̃x <̃x, no clay film consolidation and no inelastic deformation of the unit cell can occur, at least until intragranular cracking is activated. Since the elastic response was decoupled from the present inelastic model, the distribution of the stresses acting within the σ 1 -σ 3 plane on the unit cell boundaries and on the grain contacts is not rigorously constrained under these "inelastically rigid" conditions, and is here assumed to be uniform, following Guéguen and Fortin (2013). This assumption implies an isotropic stress state within the unit cell, such that the normal and shear forces supported on grain contact B are one-third of the normal and shear forces supported on the top and bottom of the unit cell, and hence given by 4R c 2 σ top /√3 and 4R c 2 τ top /√3, respectively (see Fig. 3b). The corresponding stresses are obtained by dividing by the grain contact area (πr 2 ). Combining these relations with Eqs. 5a and 5b, we obtain descriptions of ̃B and ̃B as these develop during deviatoric compression prior to slip, in terms of 1 , 3 and θ, yielding: When ̃x <̃x, the normal and shear stresses acting on A and C can then similarly be described in terms of the applied 1 , 3 and θ (in °), by: Note that for the hydrostatic case when σ 1 = σ 3 = P c , then ̃A , ̃B and ̃C are all zero, and the values of ̃A , ̃B and ̃C given by Eqs. 11a, 11c and 11e are equivalent to those given by Eq. 8a.
At θ = 0°, intergranular slip on A and C, and accompanying serial consolidation at A, B and C (cf. Fig. 4b) occurs when ̃A | |no slip =̃C | |no slip satisfy the criterion given in Eq. 10a. Hence, for this orientation, the criterion of serially coupled intergranular slip plus consolidation is obtained by inserting Eqs. 11c-11f into Eq. 10a, which gives: Similarly, inserting Eqs. 11a-11d into Eq. 10b gives the criterion for slip on and serial consolidation at A and B when θ = 30° (Fig. 4c): For values of µ in the range of 0 to 1/√3 (≈ 0.58), including the presently assumed value of 0.26, Eqs. 12a and 12b show that , indicating that for these values of µ, slip plus consolidation is easier at θ = 0°, than it is at . θ = 30°, provided that consolidation is the easier serial step. The inelastic strains developing by these processes at θ = 0° will dominate over the corresponding strain contributions at θ = 30°, as these develop earlier in the former orientation. Therefore, we will focus on describing the stress-strain behavior associated with serially coupled consolidation and slip at θ = 0° (cf. Fig. 4b) only.

Deviatoric Stress vs. Strain Behavior Due to Serially Coupled Consolidation and Slip at θ = 0°W
hen the vertical stress 1 applied during deviatoric compression at θ = 0° exceeds (see Eq. 12a), the unit cell deforms by slip and consolidation of A and C, and by consolidation of B (Fig. 4b). The normal stress acting on B is then obtained by combining Eqs. 5a, 5b, 6a, 6c and 10a, which gives: Note that during deviatoric compression at constant σ 3 , ̃B is larger during slip and consolidation of A and C (Eq. 13), than it is when A and C do not slip and consolidate (Eq. 11a).
With reference to Fig. 3c, the kinematic constraints imposed by the unit cell geometry at θ = 0° require that the inelastic grain contact normal strains (̃) developing due to consolidation of the clay films at A, B and C are related to the vertical principal strain (ε 1 ) and the horizontal principal strain (ε 3 ) of the unit cell, via: In turn, Eq. 7 can be employed to evaluate the magnitudes of ̃c lay A , ̃c lay B and ̃c lay C developing during deviatoric compression, in terms of the corresponding normal contact stress ( ̃A , ̃B or ̃C ). This is best done by normalizing the contact stresses with respect to a reference contact stress ( ̃h ) acting on grain contacts A, B and C at hydrostatic conditions (σ i = P c ) before the onset of deviatoric compression. Hence, ̃h is equal to ̃A =̃B =̃C at the σ 1 = σ 2 = σ 3 = P c conditions applied before deviatoric loading, and is given ̃h = 4R 2 c √ 3 r 2 3 (cf. Eqs. 8a, 8b). Replacing the initial reference stress ̃0 used in Eq. 7 by ̃h , and by combining this equation with Eqs. 14a, 14b, and subsequently with Eqs. 5a, 5b, 6a, 6c and 10a, we obtain relations describing the changes in maximum and minimum principal inelastic strains (Δε 1 and Δε 3 ), occurring in the unit cell in response to increasing σ 1 up to and beyond slip 1 . The result obtained for serially coupled slip plus consolidation at θ = 0° (Fig. 4b), is given: and: where Δ 1 and Δ 3 are measured relative to the hydrostatic reference state σ i = P c .

Stage 3d: Intergranular Slip-Assisted Dilation
Having explored the unit cell behavior during intergranular slip, we now constrain the conditions for intergranular slipassisted dilation, which was assumed to mark the onset of stage 3d (Fig. 4d, e). For deviatoric loading at θ = 0°, geometry implies that the unit cell microstructure cannot dilate. By contrast, deviatoric compression at θ = 30° may lead to dilation, with slip on A and B causing opening of grain contact C (Fig. 4d). This occurs when the normal stress acting on C (̃C) is reduced to zero, or perhaps to a small, negative (tensile) value corresponding to the tensile strength of the intergranular clay films. For deviatoric stresses driving slip and consolidation of A and B at θ = 30° (i.e., those satisfying Eq. 12b), the normal stress acting on C is obtained by combining Eqs. 5a-5c, 6a-6d and 10b. This gives: Ignoring any tensile strength of the intergranular illite, and setting ̃C = 0 in accordance with the onset of dilation caused by slip on A and B, plus opening of C now yields the following dilation criterion for θ = 30°: For our chosen value of µ = 0.26, this equation reduces to: The stress conditions at which dilation is the easiest at given µ can be more generally considered for θ in the range between 0° and 30°, by considering opening of C due to slip on A only (Fig. 4e). For the contact stresses corresponding to this scenario, i.e., for ̃A =̃A,̃C = 0 and ̃C = 0 , Eqs. 6a-6d deliver the criterion for the onset of dilation due to slip on A only, plus opening of C, for 0 < θ ≤ 30°: . For µ = 0.26, differentiation of Eq. 18 with respect to θ, to obtain the minimum value of dil 1 , demonstrates that dilation is the easiest at θ = 23°, with dil ≈ 5.7 3 . Hence, this orientation is most favorable for slip-coupled dilation in our model, such that Eq. 18 and θ = 23° (cf. Fig. 4e) will be used in subsequent model integration (Sect. 8) to describe the onset of dilation.

Stage 3c: Intragranular Splitting
We now derive descriptions of the stress conditions and the instantaneous strain increment associated with the onset of intragranular cracking (Fig. 4f-k), which was assumed to mark the end of stage 2 and the onset of stage 3c.

Assumed Grain Contact Properties and Behavior at the Onset of Stage 3c
As a first approximation, we base this analysis on the following assumptions and simplifications: 1. Pre-existing stress concentrators and flaws are assumed present at the microscale within the flat grain contacts featured in our microstructural model. In 2D (i.e., in the σ 1 and σ 3 plane) these flaws have a characteristic mean spacing (d) of 3 µm (Table 2), corresponding to the typical wavelength of grain contact asperities and irregularities from which intragranular cracks were observed to emanate in our experiments (Fig. 2i). 2. Because the amplitude of grain contact asperities seen in the Slochteren sandstones is small (< 1% of the grain diameter), the above grain contact roughness and flaw distribution is assumed to have negligible influence on the state of stress in the bulk of the grain, such that . the average grain contact stresses are sufficiently welldescribed in terms of total traction acting over the presently assumed, flat grain contact geometry. 3. The criterion for propagation of grain contact flaws as cracks traversing the full grain diameter is assumed to be controlled by the normal stress (̃) supported on grain contacts, in line with the particle splitting model of Kendall (1978), as applied to grain failure/cracking by Sammis and Ben-Zion (2008). Any effect of the grain contact shear stress (̃) is to be negligible, as is any effect of lateral (confining) stresses acting normal to the crack propagation direction. The latter assumption is based on the notion that crack opening will preferentially occur towards unloaded pore walls (cf. Fig. 2h-red arrow). 4. Because the highest normal load will be supported on grain contacts oriented perpendicular to σ 1 , intragranular cracking will generally be the easiest at grain contact B when θ = 0°. Therefore, intragranular cracking is evaluated for this orientation only. 5. We ignore the anisotropy implied by the stacked, 2D model geometry (see Sect. 5.2). Hence, during hydrostatic compression, ̃D is taken to be equal to ̃A , ̃B and ̃C , thus being described by Eq. 8a.
We note that excluding an effect of a lateral confining stress (Assumption 3) effectively means that we only consider intragranular cracking in grains that are not loaded by other grains in at least one direction, which may lie out of the σ 1 and σ 3 plane. In this sense, the crack opening represented in Fig. 4f-k is schematic. Micrographs of our experimentally deformed samples indeed show that intragranular cracks open more readily towards open pore spaces (Fig. 2h-red arrow), while remaining closed in more densely packed parts of the microstructure (Fig. 2h-white arrow). Still, Assumption 3 is likely only reasonable at earlier stages of intagranular cracking (early stage 3c). Kendall (1978) developed a criterion for meridional, mode I crack propagation of an incipient crack present at the narrow, flat end of a prismatic-shaped grain. The normal grain contact force (F*) required for such crack propagation is given:

Grain Splitting Criterion at the Grain Scale
where E q is the Young's modulus of the grain, Y is the fracture energy, D and b are the width and depth of the crosssection of the thick part of the grain, respectively and w is the width of the rectangular contact region. We take E q to be equal to 95 GPa for quartz grains (cf. Wong and Wu 1995), while Y is taken to be 5 J/m 2 (Table 2), which lies midway in the experimentally obtained range of 0.1-11.5 J/m 2 reported for quartz by Atkinson and Meredith (1987). Equation 19 is applied to the truncated spherical grains represented in our model microstructure, taking D ≈ 2R, b ≈ 2R c and w ≈ 2r (refer Fig. 3c). We thus obtain the following expression for the effective normal contact stress (̃ * ) required for mode I, meridional crack propagation from a flaw near the grain contact center (Fig. 4f-h):

Grain Splitting Criterion in Terms of Stresses at the Unit Cell Scale
During deviatoric compression at θ = 0°, prior to slip on A and C 1 < required for intragranular crack propagation at grain contact B (i.e., when ̃B =̃ * ; Fig. 4g) is obtained by combining Eqs. 11a and 20, which gives: Note that for hydrostatic compression, the confining pressure (P * c ) required for splitting of grain contacts A, B, C and/ or D (Fig. 4f) is equivalent to * 1 . At deviatoric stresses driving slip on A and C and consolidation on A, B and C , the vertical stress ( * 1 | | |slipAC ) required for intragranular cracking at grain contact B (Fig. 4h) is obtained by combining Eqs. 13 and 20. This gives:

Model Development
When either one of the above criteria for grain splitting is reached, transgranular cracks are produced. In the deviatoric case, these cracks propagate diametrically from B to B (Fig. 4g, h), in reality opening preferentially towards open pore space (cf. Fig. 2h-red arrow). In the hydrostatic case, all grain contact types (i.e., A, B, C and D) are equally stressed, while transgranular cracks will form normal to any one pair of grain contacts of the unit cell, i.e., normal to A, B, C or D, assuming that one of these fails slightly earlier than the others (Fig. 4f). Here, we explore the magnitude of crack opening (2c) for the deviatoric case, the strain increment that may develop because of it, and the effect it has on the subsequent loading behavior and intragranular crack evolution. The hydrostatic scenario can be evaluated in a similar way, as will be briefly shown at the end of this section.

Fig. 5
Diagrams illustrating the geometry and geometrical parameters considered in our analysis of the inelastic strain increment associated with intragranular splitting. a During deviatoric compression, splitting at grain contacts B leads to transgranular cracking normal to these contacts. Crack opening is schematically represented here as occurring normal to σ 3 , but will in reality occur towards open pores (see Fig. 2h-red arrow). Note the wedge-shaped opening defined by each intragranular crack; b On the one hand, such opening causes expansion (u exp ) at the quartz-clay interface. On the other hand, it leads to contractions of the clay films present at: c the opened end of the crack (u o cont ) ; and at: d the closed end of the crack (u c cont ) . The factors used in obtaining these contractions are indicated. e Top view on contact B at either end of the crack. By approximating the circular grain contact to be hexagonally shaped, the crack-parallel, grain contact width (L) is simplified to a bi-linear function of the distance towards the crack (x) (color figure online)

3
The assumed geometry of crack opening, and the geometrical factors considered in our analysis are defined in Fig. 5. On the one hand, crack opening by an amount 2c (c = half opening) will cause vertical expansion (u exp ) normal to the failing contact by opening space between the contact and the overlying clay film (Fig. 5b). This expansion is given by: u exp = − sin(β)r = − cr/(2R c ) (Fig. 5b). On the other hand, vertical contraction occurs, at the same site, due to indentation of the split grain halves into the overlying clay films, as seen in micrographs of lab-deformed, Slochteren sandstone (Fig. 2j). At the opened end of the crack (Fig. 5c), this contraction is given by: u o cont = sin(β)(r − m). At the opposite (closed) end of the crack (Fig. 5d), indentation of the tilted contact margin into the contact clay film yields a contraction: u c cont = sin(β)n. The net contraction (u net ) parallel to a single meridional crack, due to its opening, is then the sum of these contraction and expansion effects, given: (n − m) . In the limit of small c, the resultant vertical strain increment of the unit cell due to intragranular splitting (d cr,s 1 ) can hence be written as: Because intragranular splitting is considered for laterally unconfined grains (see Sect. 7.1), opening will not induce any lateral strains, so that d cr,s 3 = 0. The values of c, m and n in Eq. 22 are found by balancing the forces supported by the grain contacts at both ends of the crack prior to opening, versus those supported after opening, i.e., by the portions of the clay contact films that are indented by the split grain halves. The locally increased clay film consolidation strains implied by the indentations shown at both the opened and closed cracks ends in Fig. 5c, d are additive to the reference consolidation strain accumulated in the clay films at the B contacts up to the point of intragranular splitting. This reference strain is given: ̃ * = B ln(̃ * ∕̃h), while the corresponding reference clay film thickness is: z * = z(1−̃ * ) . Assuming crack opening is symmetric on both sides of each failing grain contact, only half of the clay film thickness is indented on each side (i.e., z*/2-see Fig. 5c, d). At the opening end of the crack, the contraction of the half clay film at a distance x from the contact margin is: Δl o (x) = (x − m) sin(β) (see Fig. 5c). This implies that at this end, the increase in normal strain (̂o(x)) within the half film at location x is given: Similarly, at the closed end of the crack, the increase in normal strain in the half clay film at distance x from the contact margin is: Combining Eqs. 23a and 23b with Eq. 7, and inserting Eq. 20) as the reference normal contact stress ̃0 , we obtain the grain contact normal stresses supported by the clay films present at the opened (̂o(x)) and closed ends (̂c(x)) of the crack, as a function of x: In the clay film present at the opened end of the crack, the normal force acting on a constant element dx, at position x is given: dF o (x) = S(x)̂o(x)dx . Here, S(x) is the crack-parallel width of the circular grain contact, at x (see Fig. 5e). At the closed end, the normal force contribution at x is dF c (x) = S(x)̂c(x)dx . The normal force (̃ * r 2 ) supported on grain contact B upon splitting must be equal to the force supported after opening, given the assumed microstructure. At the opening end of the crack (Fig. 5c), this force balance means: while at the closed end (Fig. 5d), we obtain: For easy integration, we approximate the circular grain contact area as being hexagonally shaped (see Fig. 5e), such that for x ≤ r/2, then S(x) = 4x, while for r/2 < x ≤ r, S(x) = 2r. The hexagonal area is similar to the circular solution within 5%. The description of S(x) in Eqs. 25a and 25b thus depends on whether m and n are larger or smaller than r/2, e.g., if 0 < m < r/2, then S(x) is equal to 4x within the bounds of m ≤ x < r/2, and to 2r when r/2 ≤ x < r. If m ≥ r/2, then S(x) = 2r. For the opened end of the crack, we will assume, for now, that 0 < m ≤ r/2 (see Fig. 5e). Later calculations demonstrate that this is indeed the case. The force balance for the grain contact at the open end of the crack is then obtained by inserting Eq. 24a along with the two simplified descriptions given above for S(x), into Eq. 25a, which gives: where the first integral term represents the force contribution within the bounds of m ≤ x < r/2, and the second integral term represents the force contribution along r/2 ≤ x < r. After integration, we obtain: Similarly, assuming that at the closed end of the crack, r/2 ≤ n < r, and inserting Eq. 24b and the above two simplified descriptions of S(x) into Eq. 25b leads to: which on integration gives: At the opened end of the crack (Fig. 5c), the work done to achieve clay indentation over a clay film strip of width dx (length S(x)), at location x is given by the product of the local force acting on this strip at x (i.e., S(x)̂o(x)dx ), and the c o r r e s p o n d i n g i n d e n t a t i o n d i s p l a c e m e n t (Δl o (x) =̂o(x)(z * ∕2)), h e n c e b y t h e p r o d u c t : Similarly, at the closed end of the crack (Fig. 5d), the work done by clay indentation at location x is given: Recalling that two intragranular cracks are formed in each unit cell (Fig. 5a), each indenting the clay on both sides, and at both ends, the corresponding internal work per unit volume dW int is given: Here 8√3R c 3 is the volume of the unit cell, while the first and second integral terms reflect the internal work contributions due to indentation of the clay films at the opening and closed ends of the two cracks, respectively. Inserting Eqs. 23a, 23b, 24a and 24b and the relations given above for S(x) , in Eq. 29a, and integrating over x gives: Equations 26b and 27b can be numerically solved to yield the values of m and n, corresponding to an arbitrary half crack opening c. An additional constraint is provided by the fact that c, and the corresponding values of m and n, must be such that the external work increment per unit volume (dW ext ) done through deforming the unit cell by intragranular splitting is approximately equal to the internal work increment per unit volume (dW int ) done (dissipated) by indenting all B-type clay films within the cell upon crack opening, i.e., assuming that the energy required to create new crack surface is provided by elastic energy released as described by the Kendall criterion (Kendall 1978). For the presently considered, deviatoric case with zero d cr,s 3 , dW ext is given: (28) dW ext = d cr,s 1 *

.
For the hydrostatic case, transgranular cracks may form normal to any pair of the contact types A, B, C or D, all of which are equally likely to fail (Fig. 4f). The external and internal work increments per unit volume associated with crack opening plus indentation at A, C and D will be the same to the corresponding work increments associated with failure of B described above, while the orientation of the strain increments developing in the unit cell due to cracking will depend on which of the contacts fail. Assuming that at the sample scale, all contacts fail at roughly equal proportions, the average, principal, inelastic strain increments d cr,s i (i = 1,2,3), developing during hydrostatic compression are approximately uniform, and each roughly equal to a third of the vertical strain increment given for the deviatoric case by Eq. 22.

Interim Model Application to Estimate Strain Magnitude Associated with Grain Splitting
We now briefly apply the above equations to the case of deviatoric loading at an effective confining pressure of 40 MPa, to explore whether our analysis implies realistic magnitudes of crack opening (2c), and what increments in vertical strain (d cr,s 1 ) should be associated with such opening. We first employ Eqs. 26b, 27b, 28 and 29b, to numerically yield the values of m, n, dW ext and dW int corresponding to a range of c values between 10 nm and 10 µm. We use the input parameters corresponding to initial porosities of 13.4, 21.5 and 26.4%, listed in Table 2. For these initial porosities, the values obtained for half opening c, at which the external and internal work increments per unit volume are balanced (dW ext = dW int ) yield 0.1, 0.3 and 0.5 µm, respectively. These values of c imply full crack openings (2c) between 0.2 and 1.0 µm, which are similar, although perhaps on the low side, of the crack openings typically seen in micrographs (0.1-10 µm; Fig. 2g-j). The corresponding, splitting induced, vertical strain increments (d cr,s 1 ) ranged between 0.01 and 0.03%, being again higher with increasing porosity. We further test whether the values calculated for the indentation limits m and n indeed fall within the ranges of 0 < m < r/2, and r/2 ≤ n < r, as previously assumed in deriving Eqs. 26a, 26b, 27a and 27b. For the full range in stress changes considered here, i.e., for deviatoric compression at P c = 5-80 MPa and for hydrostatic compression up to P * c , the obtained magnitude of m always fell in between 0 and r/2, while the resultant values of n ranged between r/2 and r. This validates the previously assumed bounds on m and n. Based on these results, we infer that our analysis of the inelastic strain associated with grain splitting, as outlined in Sect. 7.4.1, leads to representative magnitudes of crack opening, and hence of the strain increments developing because of it.

Subsequent Intragranular Evolution: Multi-edge Cracking of Split Grains
Micrographs show that cracked grains often host multiple, intragranular cracks, oriented roughly parallel, or convergent unit cell. At a constant loading rate, the subsequent behavior hinges on the capacity of the loading system to adjust to the instantaneous increase in strain, i.e., it depends on the stiffness of the sandstone plus loading system. For the purpose of giving a rough estimate of the strain that would develop upon the formation of multiple intragranular cracks in a single grain in our model for clay-bearing sandstone, we will assume that the stress acting on the unit cell upon initial splitting at θ = 0° remains constant.

Criterion for Edge Cracking
For given contact stress (e.g., at B in the θ = 0° configuration), the normal stresses acting around a newly formed, meridional crack will be elevated where the clay film is indented the most, notably around the opened crack end (Fig. 5c), and near the grain contact periphery at the closed end (Fig. 5d). Near the grain contact periphery at the closed end, the elevated contact stress leads to elastic compression within the bulky volumes constituting the sides of the split grain, and are otherwise likely of little effect. However, at the opened end, the normal stress acting on the most uplifted edges within one flaw spacing width d (3 µm; Table 2) from both sides of the initial meridional crack will be strongly enhanced and will tend to extend the first neighboring contact flaws to produce edge cracks running parallel to the initial crack surface (Lawn 1993;Tada et al. 2000). An estimate of the enhanced edge stress (̃e dge ) acting over an edge width d is obtained by subtracting the normal force acting on one split grain contact half between m < x < r − d, from the nominal force supported by the full grain contact half (i.e., 1/2 ̃B r 2 ). Using the approximate relations for S(x) outlined in Sect. 7.4.1, an expression for the normal stress acting on the most uplifted edge (̃e dge ) alongside the opened end of the meridional crack can hence be written as: which after integration yields: Here, c and m can be obtained at given stress conditions, grain size and porosity by following the approach outlined in Sect. 7.4.1. For the presently considered range in porosity (13.4-26.4%) and confining pressures (5-80 MPa), these enhanced normal stresses are thus estimated to range between 480 and 928 MPa. The critical normal stress ( * edge ), to one another (Fig. 2g, h). The stress-strain behavior, and the associated microstructural evolution following the formation of the initial, single split will depend on the boundary conditions imposed. At constant strain rate, the increase in strain due to initial intragranular splitting ( d cr,s 1 -see Sect. 7.4.1) will lead to a drop in the stress carried by the required for mode I edge crack propagation has been described by Thouless et al. (1987) to be equal to * edge = 0.87 where K Ic is the Griffith critical stress intensity factor for mode I crack propagation. For the E q , Y and d values listed in Table 2, this critical normal stress equals 489 MPa. Hence, for virtually all presently considered porosities and stresses, the enhanced normal stress acting on the edge adjacent the initial meridional crack ̃e dge ≥ * edge , implying that initial meridional splitting and crack opening is spontaneously followed by edge cracking, at assumed constant stress.

Strain Increment Associated with Multi-edge Cracking
Upon edge cracking, the two grain slices located on both sides of the initial split between r − d < x < r (Fig. 5c) will break off, and be displaced into the gaps opened by the initial, intragranular split (Fig. 5a). If the broken-off slices do not sufficiently fill this gap, they will not support load. The contact normal stress will now be supported on the remaining contact area, where it is largest at the now most-uplifted edge adjacent to the first set of edge cracks (i.e., between r − 2d < x < r − d; Fig. 5c). The normal stress acting on this edge will again exceed * edge , resulting in edge crack propagation from of the next pair of flaws. Thus, at constant stress, an instable sequence of multiple stages of edge cracking is expected, which will progress until the gap opened by the initial split is filled, such that the broken-off slices will start to support load. The strain increments ( d cr,m 1 and d cr,m 3 ) developing up to this point due to multi-edge cracking of grain contacts B are estimated by considering the gap filling process in the 2D plane, where the 2D gap area is given by: cR c + c(r − n) 2 /R c (see Supplementary Figure S2). The 2D porosity fraction φ cr characterizing this 2D gap cross section before the cracked grain starts to support load is taken to be 0.5, i.e., equivalent to a typical porosity fraction of loosely packed, angular sand (0.4-0.5; Mavko et al. 2009). During deviatoric loading, d cr,m 1 is then given by: while d cr,m 3 is again zero. For the values and porosities listed in Table 2, the value of d cr,m 1 obtained for deviatoric compression at P c = 5-80 MPa ranges between 0.08 and 0.40%. During hydrostatic compression, multi-edge cracking occurs at whichever grain contacts A, B, C or D an initial meridional crack was formed. As it was assumed that the unit cells including initial meridional cracks from any one of the contact types A, B, C or D, are equally present throughout the microstructure at the sample scale (cf. Sect. 7.4.1), the average, principal inelastic strain increments ( d cr,m i , i = 1, 2, 3) developing due to multi-edge cracking at the sample scale are then each equivalent to a third of the deviatoric, vertical strain increment d cr,m 1 , described by Eq. 31. Note that the sum of the strain increments due to initial splitting, and those resulting from subsequent multi-edge cracking developing at constant stress, constitutes the total inelastic strain increment developing at the vertical stress ( * 1 ) corresponding to initial intragranular splitting (Eqs. 21a and 21b).

Multimechanism Model Integration
We have developed a series of models for inelastic deformation in Slochteren sandstone, based on the grain-scale processes assumed to operate during each of the main stages 1, 2, 3d and 3c of mean effective stress (P) versus total (Δφ t ) and inelastic (Δφ i ) porosity reduction behavior seen in the experiments of Pijnenburg et al. (2019a). We have obtained equations for the (yield) stresses required to activate these processes, notably for: (1) isotropic consolidation of intergranular clay during hydrostatic compression; (2) serially coupled, intergranular clay consolidation and slip during deviatoric compression (Eq. 12a); (3) intergranular clay slip, leading to dilation (Eq. 18); and (4) intragranular cracking (i.e., meridional splitting, followed by multi-edge cracking at assumed constant stress), whether operating with or without intergranular slip (Eqs. 21a, 21b). The contribution of each of these processes to the inelastic deformation behavior of the sandstone was assumed to be determined by the inelastic response of the unit cell oriented most favorably for this process to occur, i.e., for the value of θ at which the corresponding activation stress is lowest. Most of these processes were found to be the easiest at a unit cell orientation θ of 0°, except for dilatant intergranular slip (Eq. 18), which for µ = 0.26 was shown to be easiest at θ = 23°. For all compactive processes (θ = 0°), approximate stress-strain relations are obtained (Eqs. 9, 15a, 15b, 22 and 31), so that the model hardening behavior due to inelastic densification can be evaluated. Elastic deformation is considered to develop independently (i.e., additively-Eq. 1) and is described through the nonlinear poroelasticity relations given in Eqs. 2, 3a and 3b.
We now integrate the above models by compiling the full set of equations using Matlab software to calculate the implied, integrated stress versus elastic and inelastic strain behavior. The integrated behavior is computed for the initial porosities (φ 0 = 13.4, 21.5 and 26.4%) used in the experimental study (Pijnenburg et al. 2019a), and the corresponding input parameter values listed in Table 2, to allow for a direct comparison between experimental and model results.

Total and Inelastic Deformation Behavior
The triaxial stress-cycling experiments reported by Pijnenburg et al. (2019a) were performed on Slochteren sandstone samples with initial porosities (φ 0 ) of 13.4, 21.5 and 26.4%. These samples were either deviatorically compressed at a constant effective confining pressure (P c ) of 5 MPa from the outset, or else first hydrostatically compressed from an initial 5 MPa, to 20, 40 or 80 MPa, followed by deviatoric compression at constant P c until stage 3d or stage 3c behavior was seen. The stress conditions used in the model were chosen to match. The modeled mean effective stress (P = [σ 1 + σ 2 + σ 3 ]/3 = [σ 1 + 2σ 3 ]/3) versus inelastic porosity reduction (Δφ i ≈ ε 1 inel +2ε 3 inel ) behavior and the P versus total porosity reduction (Δφ t ≈ ε 1 el + 2ε 3 el + Δφ i ) behavior are shown in Fig. 6 (lines), where they are compared to the corresponding experimental data (circles).
At each of the initial porosities explored, the modeled and experimentally obtained P versus Δφ t and Δφ i curves obtained are typically similar within 0.1-0.2% porosity reduction (Fig. 6). An exception is seen at the lowest P c of 5 MPa tested, where the model significantly overestimates the P-Δφ t behavior by a factor of 1.5-2, particularly for φ 0 ≥ 21.5% (Fig. 6b,  c). The experimental behavior likely reflects the anisotropic elastic behavior shown to occur at these low values of P c (Pijnenburg et al. 2019a), which is at present not included in our model (Sect. 3). Still, during hydrostatic compression up to P c = 20, 40 or 80 MPa, the modeled and experimental P-Δφ t and P-Δφ i are closely similar. During subsequent deviatoric compression, the model initially shows no inelasticity, reflecting the inability of the model to compact by serially coupled clay consolidation and slip, when slip is not yet activated (i.e., when P c < 1 < slip 1 ; Eq. 12a). In our experiments, Δφ i is never fully absent, albeit small at low deviatoric stresses. In our model, the onset of inelastic deformation during deviatoric compression (i.e., when 1 = slip 1 accompanied by an instantaneous increase in Δφ i of about 0.02-0.05% (Fig. 6). This reflects the re-equilibration of the clay films by clay consolidation and slip to the contact stresses prevalent at this point in loading (behavior cf. Eqs. 15a and 15b). During subsequent compression, the model shows more gradual P-Δφ t and P-Δφ i behavior, yielding similar, roughly linear (stage 2) hardening rates to those shown by the experimental data at the same P. Stages 1 and 2 are then followed by stage 3d or stage 3c behavior. The dilatant (suffix d) or compactive (suffix c) nature of these stages, as seen in our experiments (solid arrows in Fig. 6), is in most cases well-predicted by the model (open arrows), except for φ 0 = 13.4%, at P c = 40 MPa and for φ 0 = 26.4%, at P c = 20 MPa. For these two cases, stage 3c is predicted, while stage 3d was observed.
Where described correctly, the predicted mean effective stress marking the onset of stage 3d or stage 3c falls within 10 MPa from the experimentally obtained value. In our experiments, stage 3c is associated with concave-down P-Δφ i behavior, while our model implies an instantaneous increase in Δφ i at constant stress. On the other hand, the inelastic porosity reduction seen in our experiments as stage 3c starts to develop (indicated by c-arrows), is closely similar to the modeled inelastic porosity reduction after initiating stage 3c, i.e., within a relative difference of 5%. Despite the noted, modest discrepancies, the modelled and experimental P-Δφ t and P-Δφ i data show a general agreement, implying that the P-Δφ t and P-Δφ i behavior can be largely accounted for by the mechanisms assumed.

Stress Conditions Required for Stages 3d and 3c Inelastic Deformation
Yield envelopes describing the differential stress (Q = σ 1 − σ 3 ) versus mean effective stress (P = (σ 1 + σ 2 + σ 3 )/3) conditions required to activate each of the assumed, inelastic processes at modelled initial porosities of 13.4, 21.5 and 26.4% are shown in Fig. 7. The yield envelopes obtained for dilatant intergranular slip (i.e., stage 3d) and for serially coupled clay slip plus clay consolidation are insensitive to φ 0 , while the envelopes obtained for intragranular cracking (stage 3c) are markedly φ 0 -sensitive (Fig. 7). Furthermore, the combined outline of the predicted stages 3d and 3c envelopes describes a trilinear, or broadly elliptical shape in Q-P space, transecting the P axis at the origin and at the hydrostatic confining pressure for isotropic grain crushing (P * c ) . This shape and the predicted φ 0 -(in)sensitivities characterizing stages 3d and 3c in our model are in qualitative agreement with the behavior seen in experiments performed on Slochteren sandstones (Pijnenburg et al. 2019a) and on many other sandstones (Wong and Baud 2012-see also our Introduction). Hence, at least for these high strain stages 3d and 3c (ε > 1%), the agreement between predicted and observed stress conditions Fig. 6 Comparison between the mean effective stress (P) versus total-(Δφ t ) and inelastic porosity reduction (Δφ i ) data implied by the model (continuous and dashed lines) and the corresponding experimental data (dots) obtained in our stress-cycling tests reported by Pijnenburg et al. (2019a), for initial porosities (φ 0 ) of: a 13.4%; b 21.5%; and c 26.4%. We explored hydrostatic compression up to the maximum effective confining pressures (max P c ) indicated, followed by deviatoric compression up to the onset of dilation (stage 3d indicated "d"-cf. Fig. 1) or up to and including the onset of enhanced compaction (stage 3c-indicated "c"). The experimental data reflect the maximum P conditions and the corresponding values of Δφ t and Δφ i explored in each stress cycle. The modeled and experimental data show roughly similar, stages 1, 2 and 3c P-Δφ t and Δφ i behavior. Note that at φ 0 ≥ 21.5%, the experimental data show dilation prior to that predicted by the model (see inset figures) ◂ is found to be satisfactory. This implies that the dilatant behavior is well-captured by considering the onset of tension at grain boundaries developing during intergranular slip (see also : Guéguen and Fortin 2013), while stage 3c initiation can be accounted for by considering intragranular cracking, which at Q > 0.6P is made easier by accompanying intergranular slip. We now go on to discuss whether the model captures the inelastic yield envelope expansion implied for stages 1 and 2 of the experiments performed by Pijnenburg et al. (2019a). In that way, we evaluate to what extent the associated hardening behavior due to inelastic densification is accounted for by consolidation of-and slip on the intergranular clay films.

Inelastic Compaction Behavior During Stages 1 and 2: Improved Mechanistic Constraints
To constrain the model hardening behavior (i.e., the yield envelope expansion) due to inelastic porosity reduction, we show contour lines delineating fixed values of inelastic porosity reduction in Fig. 7 (thin green lines; magnitude indicated in %). These data are obtained by interpolating Q-P-Δφ i data rendered in multiple model runs, in which hydrostatic-and subsequent deviatoric compression is simulated at a maximum P c in the range of 5-140 MPa and at a φ 0 value of 21.5%. Similar Δφ i -contour lines were obtained for initial porosities of 13.4 and 26.4%. These are not shown, as the behavior was highly similar to that seen for φ 0 = 21.5%, albeit less-and more compliant, respectively. Overall, the modeled inelastic porosity reduction contours obtained for a φ 0 value of 21.5% at the Q and P conditions enveloped by stage 3d (dilation) and stage 3c (intragranular cracking) are positively and steeply inclined (Fig. 7), indicating that inelastic compaction is primarily sensitive to P, not to Q. The interspacing between the lines decreases with increasing P, reflecting stiffening of the clay films as these consolidate. Recall that throughout hydrostatic compression, the unit cell deforms by inelastic clay consolidation (cf. Fig. 7 Differential stress (Q = σ 1 − σ 3 ) versus mean effective stress (P = [σ 1 + σ 2 + σ 3 ]/3) conditions delineating the yield conditions for each of the (combination of) inelastic processes assumed to govern the indicated stage of P-Δφ t behavior in our model, at initial porosities (φ 0 ) of 13.4%, 21.5% and 26.4%. All processes are described at θ = 0° (refer Fig. 3c), except for dilatant intergranular slip (shown at top-left) which for the assumed intergranular illite friction coefficient µ = 0.26 was found to be the easiest at θ = 23°. Stress conditions delineating fixed values of the inelastic porosity reduction (Δφ i ) com-puted for φ 0 = 21.5% are contoured with thin green lines (Δφ i values in % are indicated). At 0 < Q < 0.6P, no inelastic deformation by intergranular slip and serially coupled clay film consolidation occurs, while fixed Δφ i stress contours have a slope dQ/dP equal to 3, equivalent to that of the stress path implied for deviatoric loading (increasing σ 1 ) at constant σ 2 = σ 3 . Upon intragranular cracking (stage 3c), fixed Δφ i stress contours follow the cracking envelopes towards lower P (color figure online) Eq. 9). However, at deviatoric stress conditions below those required to activate intergranular slip (Q < 0.6P) and intragranular cracking, the unit cell cannot deform inelastically. Therefore, at these deviatoric stresses, corresponding to inelastically rigid behavior, the slope (dQ/dP) of the inelastic porosity reduction contours is equal to 3, i.e., equivalent to the slope dQ/dP implied during deviatoric loading at constant confining pressure. Upon activating intergranular slip plus clay consolidation (Q = 0.6P), all contours shift along the corresponding envelope to lower Q and P, reflecting the instantaneous increase in Δφ i accompanying the onset of serially coupled, intergranular slip and clay film consolidation, as noted earlier (Sect. 8.1.1). At Q > 0.6P, contour lines are steeper (dQ/dP = 5 ± 1), due to inelastic compaction by active clay consolidation plus slip, at these conditions. Where Δφ i contours intersect the intragranular cracking envelopes (Fig. 7), their orientation changes to follow these envelopes towards lower P values, implying negative slopes (dQ/dP) of − 0.5 (if accompanied by intergranular slip), or − 1.5 (no intergranular slip). Hence, the onset of intragranular cracking in our model marks a sharp change from primarily P-sensitive-(due to clay film deformation) to P-and Q-sensitive inelastic porosity reduction behavior (due to intragranular cracking).
In the experiments performed by Pijnenburg et al. (2019a), a similar transition was noted, notably from P-sensitive Δφ i behavior during stage 1 and the initial portions of stage 2, to P-and Q-sensitive behavior, towards the end of stage 2, and into stage 3c. The initial P-sensitive Δφ i behavior of stages 1 and 2 was then speculated to be caused by a yet unidentified, intergranular deformation process. The increased Q-sensitivity towards the end of stage 2 and into stage 3c was tentatively attributed to a gradually increasing role played by intragranular cracking, in accordance with accompanying crack density data and earlier inferences (Curran and Carroll 1979;Wong et al. 1997). The present microphysical model provides a physical basis to underpin these earlier inferences. It confirms that while predominantly P-sensitive hardening behavior can be explained by intergranular clay consolidation and accompanying intergranular slip, Q-sensitive hardening behavior is to be attributed to a role played by intragranular cracking (Shah and Wong 1997).

Model Evaluation and Suggestions for Future Improvement
We have shown that the main experimental trends in Q-P-Δφ i behavior seen in earlier experiments on Slochteren sandstones can be accounted for by the present model. However, discrepancies were also found, including the discontinuous P-Δφ i behavior implied by the model at the onset of serially coupled, intergranular slip plus clay consolidation, and at the onset of intragranular cracking, which was not seen in experiments. These discrepancies mainly arise from the discrete yield conditions implied by the simplifications underlying our approach, i.e., our assumption of a regular 2D array of identical quartz grains and intervening clay films, with uniform dimensions, packing, contact geometry and microscale strength parameters. In natural sandstones, these microscale features and properties will be distributed (McDowell et al. 1996;Brzesowsky et al. 2011;Cook et al. 2015), so that the P-Δφ i behavior is likely to be smoothed. Better quantification of the stress-strain behavior exhibited by Slochteren sandstone requires accounting for such microstructural and microscale strength parameter variations and distributions in the framework of a more rigorous grain-scale to aggregate-scale (representative elementary volume-scale) averaging approach.
To illustrate the effects of micro-scale parameter variability, we note for example that the empirical plasticity model reported by Pijnenburg et al. (2019a) predicts fully balanced inelastic dilation and compaction (i.e., critical state) for φ 0 ≥ 21.5%, at Q = (1.7 ± 0.1)P. This means that in the experiments and at these high porosities, dilation occurs prior to the onset of dilation implied by the present microphysical model at the assumed constant grain contact friction coefficient µ of 0.26 (Q = 1.8P). If the value of µ is taken to vary across the microstructure, for instance due to variations and changes in the consolidation state of and/or shear strain within the intergranular clay (see: Tembe et al. 2010), then the onset of dilation will accordingly be spatially variable. For the smallest µ value of 0.22 within the range described by Tembe et al. (2010) (µ = 0.22-0.30, at shear strains between 0.4 and 8.0-see Sect. 5.3), Eq. 18 implies that dilation occurs at Q = 1.7P and at θ = 24°, hence at similar stress conditions to those implied by the plasticity model at critical state. In the limit of zero intergranular friction (µ = 0), which may locally be the case when grains are unconstrained in the slip direction, Eq. 18 implies dilation when Q = 1.2P and θ = 0°. Note that in this end-member case, the resistance to dilation is controlled by the dilatancy angle implied by the model packing (see Niemeijer and Spiers 2007).
In future, the effects of microstructural heterogeneity within the quartz grain framework, and of varying intergranular clay film thickness and properties on inelastic deformation behavior of sandstones of the type considered here need to be investigated. This can be done by applying a more rigorous 3D averaging procedure to obtain the Q-P-Δφ i behavior of a representative volume of sandstone with distributed microstructural parameters (grain size/shape, grain contact orientation, grain contact area, clay film thickness) and micro-scale strength parameters. As mentioned in the Introduction, routes to achieve this include classical homogenization treatments (e.g., Bardet and Vardoulakis 2001;Fortin et al. 2003), thermodynamically-based granular mechanics modeling methods (e.g. Einav 2007a, b;Tengattini et al. 2014)

Implications: Improved Basis to Assess Experimental Data Applicability at Decade Time Scales
In its present form, the microphysical model contributes to understanding of the processes governing inelastic deformation, particularly at small stresses and strains (stages 1 and 2) relevant for reservoir compaction during gas production. The model can thus be used to underpin the applicability of existing geomechanical models for hydrocarbon productioninduced subsidence and seismicity, at field conditions and time scales. An important finding is that the stage 1 and particularly the stage 2 inelastic behavior most relevant for in situ compaction can be largely accounted for by intergranular clay film deformation, with perhaps a small role played by intragranular cracking towards the end of stage 2 (see Sect. 8.1.3). A simple analysis shows that compaction by intergranular slip plus clay consolidation is virtually time or rate insensitive at effective stresses pertaining to the current state of depletion of the Groningen reservoir (σ 1 = 57 MPa; σ 3 = 27 MPa), so that reservoir compaction by these mechanisms is expected to be similarly time or rate insensitive (Pijnenburg et al. 2019b). In accordance with this, long-term compaction experiments performed on the Slochteren sandstone over weeks (Pijnenburg et al. 2018) to months (Hol et al. 2015) at similar stresses, temperature (100 °C) and chemical conditions (4 M salt brine) to those seen in the Groningen field showed that the bulk of stage 2 inelastic compaction is instantaneous, while decelerating creep deformation contributed a modest 10-20% to the inelastic strain accumulated during active loading. The time-independent plasticity plus poroelasticity model outlined by Pijnenburg et al. (2019a) is accordingly anticipated to capture the main trends of the in situ compaction behavior at the decade time scales relevant to the field, although compaction strains and lateral stresses may be slightly underestimated by 10-20% due to other (decelerating) creep effects seen in long-term (weeks-months) experiments (Hol et al. 2015;Pijnenburg et al. 2018).

Conclusions
In this study, a series of simplified microphysical models has been derived to explain the inelastic deformation behavior shown in triaxial experiments performed under in situ conditions of temperature (100 °C), stress and pore fluid chemistry (~ 4 M brine), on Slochteren sandstone samples from the seismogenic center of the Groningen gas field (Pijnenburg et al. 2019a). In particular, we sought a mechanistic basis for the continuous inelastic (permanent) deformation seen at small strains (ε ≤ 1.0%) relevant to the field, i.e., during stage 2 of the experiments considered (i.e., the near-linear mean effective stress P versus total porosity reduction Δφ t behavior). On the basis of microstructural evidence obtained in previous experiments, inelastic deformation was inferred to be governed by intergranular slip and consolidation of intergranular clay films during stage 1 (concave-up P-Δφ t behavior) and stage 2. The subsequent dilatant (stage 3d) and nonlinear compaction behavior (stage 3c) were inferred to be governed by dilatant, frictional slip and intragranular cracking, respectively. The model microstructure was taken to consist of a regularly packed 2D array of quartz grains (radius 100 µm) with flattened/dissolved contacts and intergranular clay films (thickness 3-6 µm). The integrated model behavior was explored for the same initial sample porosities (φ 0 = 13.4, 21.5 and 26.4%) used in our experiments and for stress conditions covering the bulk of those explored in our tests, i.e., differential stresses (Q) up to 130 MPa and effective confining pressures up to 165 MPa. This allowed for a direct comparison between the model-implied and experimentally obtained behavior. We draw the following main conclusions: 1. Our model captures the main trends in P versus total (Δφ t ) and inelastic porosity reduction (Δφ i ) seen in experiments on Slochteren sandstone, showing a similar progression in stages 1, 2 and 3d or stage 3c P-Δφ t and P-Δφ i behavior. 2. During stages 1 and 2, the model implies a yield envelope that expands with increasing Δφ i from the onset of compression, reflecting hardening due to intergranular clay consolidation. The hardening behavior was found to be primarily P-sensitive when governed by clay consolidation and slip (stages 1 and 2) and markedly Q-sensitive upon subsequent intragranular cracking (i.e., upon stage 3c). This behavior is qualitatively similar to that seen in our experiments. 3. The yield envelope describing stage 3d (dilatant intergranular slip) is φ 0 -insensitive and describes a direct dependence of Q on P (i.e., Q = 1.8P). The yield envelopes obtained for stage 3c (intragranular cracking) are φ 0 -sensitive and show an inverse dependence of Q on P, being shallower when accompanied by intergranular slip (dQ/dP = − 0.5, compared to dQ/dP = − 1.5 absent slip).
The φ 0 -(in)sensitivities and the orientations described in Q-P space are in qualitative agreement with the stage 3d and stage 3c seen in the Slochteren sandstone and in many other sandstones.
4. We infer that the processes assumed in the present, simplified model, i.e., clay consolidation, intergranular slip and intragranular cracking can indeed account for the main trends in inelastic deformation behavior seen in our experiments. The discrepancies between the model predictions and the observations are attributed to the discontinuous inelastic deformation behavior implied by our model, which in reality is likely to be more smoothed out due to distributed values of grain (contact) strength, or due to variations in grain packing. 5. The present model provides mechanistic underpinning for the empirical, constitutive (poroelastic plus plastic) compaction model for the Groningen reservoir (Slochteren) sandstone reported by Pijnenburg et al. (2019a), where at the small strains relevant to reservoir compaction (ε ≤ 1%), the inelastic contribution is largely governed by serially coupled intergranular slip and clay film consolidation. 6. These mechanisms were shown earlier to be virtually rate-insensitive at the decade time scales and in situ effective stress conditions pertaining to the Groningen gas reservoir. Creep effects due to other processes are seen in experiments, although their added contribution to the instantaneous inelastic strain accumulated during loading is modest (10-20% of instantaneous value) (Hol et al. 2015;Pijnenburg et al. 2018). 7. This means that the empirical, constitutive model reported by Pijnenburg et al. (2019a) is expected to capture the main trends of the in situ compaction behavior at the decade time scales relevant to the Groningen field, although compaction strains and lateral stresses may be slightly underestimated due to the above-mentioned, modest creep effects. Hence, this empirical model can be incorporated in geomechanical models that investigate the effects of gas production from the Groningen, where it will improve estimates of the stress evolution in the reservoir and the elastic strain energy budget available in the reservoir for release in induced seismic slip on faults. 8. The previously reported empirical model, underpinned by the present microphysical model, contributes directly to a more realistic, physics-based description of induced seismicity in the Groningen gas field at relevant decade time scales. Still, to fully constrain the long-term compaction behavior, the modest creep effects seen in experiments at reservoir-relevant stresses and strains need to be investigated further in future work. 9. In addition, for further improvement in constitutive modeling efforts, the mechanisms observed in our experiments and inferred from our simple model to be important in the Groningen reservoir sandstone (and probably in other, similar sandstones) should, in future, be incorporated into 3D models based on the more rigorous grain-to aggregate-scale averaging procedures. Suitable candidates here range from classical homogenization and granular mechanics approaches to the discrete element method.
Dr. Ronald P. J. Pijnenburg is a postdoc researcher at the High Pressure and Temperature Laboratory at Utrecht University. He studies the mechanical and mechanistic deformation behavior of sandstone, primarily within the context of anthropogenic use of the subsurface, for energy production or f luid storage.
Prof. Christopher J. Spiers is Professor of Earth Materials and Head of the High Pressure and Temperature Laboratory at Utrecht University. He specializes in research on the mechanical behavior of rocks and faults under the pressure, temperature and chemical conditions that pertain in the Earth's crust and on the controlling microscale processes.