Compaction of the Groningen Gas Reservoir Sandstone: Discrete Element Modeling Using Microphysically Based Grain‐Scale Interaction Laws

Reservoir compaction, surface subsidence, and induced seismicity are often associated with prolonged hydrocarbon production. Recent experiments conducted on the Groningen gas field's Slochteren sandstone reservoir rock, at in‐situ conditions, have shown that compaction involves both poroelastic strain and time independent, permanent strain, caused by consolidation and shear of clay films coating the sandstone grains, with grain failure occurring at higher stresses. To model compaction of the reservoir in space and time, numerical approaches, such as the Discrete Element Method (DEM), populated with realistic grain‐scale mechanisms are needed. We developed a new particle‐interaction law (contact model) for classic DEM to explicitly account for the experimentally observed mechanisms of nonlinear elasticity, intergranular clay film deformation, and grain breakage. It was calibrated against both hydrostatic and conventional triaxial compression experiments and validated against an independent set of pore pressure depletion experiments conducted under uniaxial strain conditions, using a range of sample porosities, grain size distributions, and clay contents. The model obtained was used to predict compaction of the Groningen reservoir. These results were compared with field measurements of in‐situ compaction and matched favorably, within field measurement uncertainties. The new model allows systematic investigation of the effects of mineralogy, microstructure, boundary conditions, and loading path on compaction behavior of the reservoir. It also offers a means of generating a data bank suitable for developing generalized constitutive models and for predicting reservoir response to different scenarios of gas extraction, or of fluid injection for stabilization or storage purposes.

Inverse modeling using surface subsidence data, and forward modeling using experimental data, are commonly employed to estimate and predict compaction at the reservoir level (Cannon & Kole, 2017;Dusseault & Rothenburg, 2002;NAM, 2013), often assuming linear poroelastic behavior (Du & Olson, 2001;J. D. Smith et al., 2019;H. F. Wang, 2000). However, the models obtained are of questionable reliability for forecasting compaction, as the partitioning between elastic and any inelastic deformation that may be occurring at the grain-scale is generally unknown. Extrapolation is accordingly unconstrained. Ideally, reservoir-scale models incorporating constitutive relations describing the grain-scale deformation mechanisms that actually operate under in-situ conditions are needed to improve understanding and prediction of reservoir compaction.
In the N.E. Netherlands, gas production from the vast (30 × 30 km) Groningen gas field (Europe's largest) since 1963 has led to ∼34 cm of subsidence at the center of the field and to noticeable seismicity since the early 1990's, culminating in a magnitude 3.6 event in 2012. It is believed that these phenomena are driven by compaction of the Slochteren Sandstone reservoir formation (Mid-Upper Permian), located at ∼3 km depth (Spiers et al., 2017). The reservoir consists of highly porous aeolian and alluvial sandstones (15%-18% average porosity [NAM, 2016]) with a thickness ranging from 300 m in the NNW of the field to 140 m in the SSE (NAM, 2016). If the surface subsidence is representative of compaction at depth, this means that the reservoir has undergone a vertical strain of 0.1%-0.2% in the center of the field, where porosity is highest and seismicity is most marked (NAM, 2013). Experimental studies indicate that much of this strain is inelastic, i.e. permanent (Pijnenburg et al., 2019a). While elastic deformation is generally rather easily predicted from laboratory data on reversible behavior (H. F. Wang, 2000), inelastic deformation of the Groningen reservoir is not, since the grain-scale mechanisms responsible are poorly understood and have only recently been identified (Pijnenburg et al., 2019a).
Much research has been done on understanding the inelastic compaction and failure behavior of sandstones (Wong & Baud, 2012). However, most of this work focuses on large strain deformation (>>1%, Wong & Baud, 2012), whereas production-induced compaction in sandstone reservoirs like Groningen involves much lower strains (<<1%) (Hol et al., 2018;Pijnenburg et al., 2019a). Typical sandstone compaction behavior under first hydrostatic and then deviatoric loading conditions can be divided into three stages (Pijnenburg et al., 2019a). During initial loading (Stage 1), the mean stress versus volumetric strain behavior shows a concave upward trend, often associated with poroelastic compaction and closure of preexisting cracks. Stage 2 is typically associated with (near-) linear behavior, leading to nonlinear Stage 3 behavior, characterized by dilatant shear failure at low confinement and by compactive failure at high confinement. It is Stage 3, where grain failure is a major deformation mechanism especially at high confinement that has been investigated in most previous studies of inelastic deformation. The transition from (near-) linear Stage 2 to nonlinear Stage 3 is conventionally identified as the plastic (compactive) yield point in tests at high confinement (Wong & Baud, 2012). At relatively low confinement, where Stage 3 corresponds to dilatant shear failure (Stage 3), Stage 2 may be absent and dilatant yield is identified as the transition from Stage 1 to Stage 3 (Wong & Baud, 2012). In this paper, alongside these traditional definitions of yield point, as marking the onset of large inelastic strains at the start of Stage 3, we note that significant inelastic deformation can also occur in Stages 1 and 2 of compaction (Pijnenburg et al., 2019a; see also Wong et al., 1992). Therefore, from now on in this paper, we refer to the onset of Stage 3 as the macroscopic yield point, recognizing that inelastic deformation (strain hardening yield) generally occurs continuously from the onset of loading.
For the Slochteren sandstone at the center of the Groningen field, where the mean porosity is >18%, the relevant in-situ strains (0.1%-1.0%) have been shown to correspond to (near-) linear Stage 2 stress-strain behavior, with 30%-50% of the total volumetric strain being irrecoverable (Hol et al., 2018;Pijnenburg et al., 2018). Though this type of (near-) linear Stage two behavior is typically assumed to be purely poroelastic (Baud et al., 2004;Wong & Baud, 2012), in the case of the Groningen reservoir sandstone around 30%-50% of the deformation is inelastic due to the presence of thin intragranular clay films (grain coatings), which accommodate significant inelastic consolidation and shear (Pijnenburg et al., 2019b;Verberne et al., 2020). In Figure 1, clay films around grains are shown for Slochteren sandstone sample. These films 10.1029/2021JB021722 3 of 23 accordingly dissipate 30%-50% of the mechanical work done during compaction, so that less elastic energy will be available for release during fault rupture than if the reservoir were fully elastic (Pijnenburg, 2019). This means that understanding the elastic and inelastic strain partitioning in the reservoir is key for modeling both compaction and induced seismicity in Groningen. However, for a reservoir sandstone of 140-300 m thickness and spanning an area of 30 km by 30 km, the mechanical properties will vary, both vertically and laterally, due to variability in porosity, grain size, clay content, and other microstructural parameters (Hol et al., 2018;Wong & Baud, 2012). Experimental investigation of the deformation behavior of the entire reservoir is therefore not feasible. Instead, (rapid) numerical modeling capability, based on the microphysical processes operating at the grain-scale, is needed to explore the effect of this variability. Only once we can model granular behavior at the lab-sample scale accurately is it possible to upscale laboratory data to describe the average mechanical behavior of the reservoir at the individual borehole or regional scale and extrapolate this behavior in time.
Discrete Element Method (DEM) is a powerful numerical method intended for linking grain-scale properties to macro-mechanical behavior of granular media (Marketos & Bolton, 2009), which may offer an appropriate tool for predicting sandstone behavior. However, while DEM simulations can qualitatively describe the stress-strain behavior of granular aggregates, and can be fitted to specific experimental data sets, popular contact models generally assume (visco-) elastic and frictional interactions at grain contacts. More recently, studies have been reported, taking into account complex contact laws for materials such as concrete (Tran et al., 2011) and clay-coated loose sands (Gilabert et al., 2007;B. Wang et al., 2008), further developments are still needed, especially for highly porous sandstones containing clay minerals. Currently, then, DEM simulations are unable to capture the effects of intergranular clay film deformation, or of grain crushing and pore collapse resulting from microcracking, in a manner that describes the physics of grain-tograin interactions, as observed in cemented sandstones, such as the Slochteren sandstone. To improve the predictive capability of DEM simulations for sandstones under reservoir conditions, it is therefore necessary to incorporate independent and verifiable descriptions of the physical and chemical processes that occur at the grain-contact scale.
In this paper, we take a step toward this aim by incorporating a new contact model describing both nonlinear elastic and inelastic interactions, resulting from compaction of and slip along intergranular clay films (Pijnenburg et al., 2019b). Grain failure is also introduced, based on a grain indentation model leading to the propagation of a penny-shaped microfracture. We calibrate our new contact model against hydrostatic and deviatoric loading data obtained in conventional triaxial compression tests performed on the Slochteren sandstone (by Pijnenburg et al., 2019a). These experiments were performed on sets of samples having porosities of approximately 13.4%, 21.5%, and 26.4%, and having mean grain sizes in the range 150-250 μm (Pijnenburg et al., 2018). The calibrated discrete element model is then used to systematically investigate the compaction behavior of Slochteren sandstone with different porosities, particle size distributions, and clay contents, and under different stress-strain boundary conditions. Specifically, the DEM simulation reproduces the compaction behavior of Slochteren sandstone measured in independent depletion experiments (Hol et al., 2018), performed under in-situ reservoir stress and uniaxial strain conditions. Moreover, our discrete element model is able to predict reservoir compaction in the center of the Groningen gas field, as determined from in-situ measurements within the Zeerijp-3a well, within measurement uncertainties (Cannon & Kole, 2017). For convenience, all symbols in this paper are listed and defined in the Appendix A.

Key Microstructural Characteristics
The Slochteren formation is comprised of fine-to coarse-grained sandstones and locally even conglomerates (Waldmann, 2011). The typical composition of the sandstones is 72%-90% quartz, 8%-25% feldspar 10.1029/2021JB021722 4 of 23 (mainly K-feldspar), 3%-10% lithic fragments, and 0.1%-10% clay minerals (kaolinite and illite) (Pijnenburg et al., 2019a;Waldmann et al., 2014). The grain size ranges from 71 to 978 μm with an average value of 256 μm (Waldmann, 2011). In general, grain-to-grain contact geometries range from point to extended contacts (e.g., Waldmann, 2011). In the Groningen gas field, the average sandstone porosity increases from 12% to 16% at the margins to 18%-22% at the center (NAM, 2016). The higher porosity sandstones in the center of the field contain more clay, both as pore-filling vermicular overgrowths and as thin, grain-coating rims (NAM, 2013;Waldmann, 2011;Waldmann et al., 2014). The lower porosity sandstones are well-cemented with carbonates and quartz (Gaupp et al., 1993). Clay coatings can be up to 15 μm thick (Waldmann, 2011). While kaolinite is the dominant clay mineral in the Slochteren sandstone (up to 10%), it is mainly present as pore-filling clay, resulting from diagenetic processes during burial, such as feldspar dissolution. On the other hand, illite makes up only 1% of the rock mineralogy, but since it was formed prior to burial, it is the main grain-coating clay mineral (Waldmann, 2011) and is widely found in grain contacts in films up to 10 μm thick (Pijnenburg et al., 2019b).

Nonlinear Elastic and Inelastic Behavior
As mentioned in Section 1, nonlinear elastic behavior of sandstones under compression is typically observed during Stage 1 (Khan et al., 1991;Morgenstern & Tamuly Phukan, 1969;Wong et al., 2004;Zimmerman, 1991). Previous studies relating effective hydrostatic or mean stress to stiffness or pore compressibility found exponential (Zimmerman, 1991) or logarithmic (Morgenstern & Tamuly Phukan, 1969) correlations, attributed to pore volume reduction, increasing grain contact areas or the closure of crack-like voids with a high aspect ratio (David et al., 2012;Wong & Baud, 2012). Slochteren sandstone exhibits similar nonlinear elastic behavior during hydrostatic loading (Pijnenburg et al., 2019a), with an increase in confining pressure from 5 to 40 MPa leading to an increase in stiffness of the sandstone by a factor of 1.5-3, depending on initial porosity. For higher confining pressure of 40 MPa (Stage 2), Pijnenburg et al. (2019b) assumed constant stiffness. Although that is an acceptable assumption, in this study we assumed nonlinear elastic behavior at all confining pressures. This is illustrated in Figure 2, which displays the unloading response of nine samples of Slochteren sandstone, with three different porosities (13.4%, 21.5%, and 26.4%), as obtained during incremental, hydrostatic load-cycling tests. Note that the different porosity samples were obtained from three different depths (i.e., horizons) in a single core (Zeerijp-3a), taken from the center of the Groningen gas field, so that for each set of samples with a given porosity, the mineralogy, grain size, and microstructure were virtually the same. Since unloading data are displayed, the response approaches being fully elastic yet displays clear nonlinearity. The corresponding mean effective stress (p) versus elastic volumetric strain ( e v E e) data can be described using least-square quadratic fits. It should also be noted that for material with a 10.1029/2021JB021722 5 of 23 porosity of 26.4%, pressurized to a final effective confining pressure of 80 MPa, sample stiffness decreases once an effective confining pressure of 60 MPa is reached. Therefore, the quadratic fit is applied only to the data obtained up to 60 MPa effective confining pressure, beyond which a decrease in stiffness is observed, presumably due to inelastic effects.
As already discussed, Pijnenburg et al. (2019b) attributed inelastic deformation of Slochteren sandstone under in-reservoir stress conditions (Stage 2) to compaction of and slip along thin clay films (up to 10-15 μm) present between the sand grains. The local inelastic strain showed a positive correlation with clay content. Pijnenburg et al. (2019b) further suggested that clay consolidation could be described using a logarithmic dependence on grain contact normal stress, in line with the log laws generally used to describe clay consolidation under fully drained conditions (P. R. Smith et al., 1992). Kasyap et al. (2021) also formulate the effect of clay consolidation at grain contacts using the Hertzian contact, for which the radius of contact is increasing by consolidation.
Inelastic deformation of Slochteren sandstone at high mean stresses, where grain failure starts to dominate, i.e. in Stage 3 (Pijnenburg et al., 2019b), is more difficult to quantify at the microphysical level. In a sandstone, mechanisms such as elastic mismatch, grain bending, (sub)Hertzian contact indentation, and stress concentrations generated at grain surface and internal flaws (Kemeny & Cook, 1991;Kendall, 1978; can all generate intragranular fractures. Resulting fracture configuration include chipping/spalling, splitting and fragmentation, as observed, for example, by Karatza et al. (2019) using in-situ X-ray computed microtomography of a sand aggregate undergoing uniaxial compaction (see also Hangx et al., 2010;Hol et al., 2018;.

Present Discrete Element Approach and Particle Interaction Laws
DEM is a discontinuum mechanics-based method in which materials are defined as assemblages of rigid or deformable blocks or particles (Cundall, 1971;Cundall & Strack, 1979). Newton's second law of motion is applicable to each particle, with force-displacement laws acting between contacting particles. Unlike other methods, such as Finite Element or Finite Difference Methods, displacement compatibility is not necessary between the particles in DEM. Therefore, discrete element models can undergo large deformations with contacts being continuously made and unmade, with automatic detection, while the calculation progresses.
Since the introduction of "classical" DEM, various reliable discontinuum-based numerical approaches have been developed to model discontinuous materials (Munjiza, 2004;Munjiza et al., 2011), such as combined finite-discrete element (FDEM) (Gao et al., 2018;Munjiza, 2004), molecular dynamics, (Alder & Wainwright, 1959;van den Ende et al., 2018) and Discontinuous Deformation Analysis (Shi, 1992). However, these newer methods often come at a higher computational cost. Since our aim is to develop a method that can rapidly predict stress-strain behavior of granular aggregates for a range of stresses and microstructures (porosity, and mineralogical composition), this makes "classical" DEM a highly suitable method for modeling granular media.
We employ the Particle Flow Code (PFC) (Itasca, 2019) as a DEM approach, which uses rigid disks or spheres to represent particles in 2D or 3D, respectively. Using spherical particles facilitates the contact detection process (i.e., neighbor search, see Munjiza et al. (2011) for details about contact detection algorithms) and the rigid particle assumption reduces the computational time by decreasing the number of degrees of freedom. Particle penetration, sliding, rotation, and separation (i.e., contact interactions) thus form the sources of displacements occurring in PFC models. Particles have contact with adjacent particles or with walls (or facets), i.e. rigid walls typically used to apply boundary conditions. The behavior of particle-particle or particle-wall contacts is described by a particle-interaction model, or going forward, simply referred to as a "contact model". This is essentially a force-displacement law linking normal and shear contact forces to normal and shear displacements across contacts. In a simple contact model (e.g., linear elastic contact model), Hooke's law is applied, in the normal and shear directions, to relate the normal force ( n E F) and shear force ( s E F), through the normal stiffness ( n E k ) and shear stiffness ( s E k) of the contacts, to the relative normal and shear displacements ( n E d and s E d, respectively) of particle centers across the contact. However, force displacement laws can be more complicated by limiting forces with friction, adding dashpots, or using nonlinear elastic interaction laws, for example (Itasca, 2019).

10.1029/2021JB021722 6 of 23
The PFC approach is an explicit approach, using the current state of the model aggregate to determine the next state. Timesteps between states must be small enough to assume that "disturbances" of one particle do not propagate beyond adjacent particles (Cundall & Strack, 1979), while ensuring the model remains conditionally stable. Each timestep calculation cycle includes a sequence of operations, starting with the calculation of the timestep, followed by the determination of the new position and velocity of each particle based on Newton's second law of motion. From the new particle positions, the new contact population is established per particle. Finally, the forces and moments are updated for each contact, based on the selected contact model (i.e., the selected force-displacement law). This cycle of calculations continues from the initial state until the model reaches the desired end state (Itasca, 2019).
Continuum-based numerical modeling approaches typically use macro-mechanical parameters, such as Young's modulus, shear modulus or Poisson's ratio, as input for continuum-based constitutive laws, to describe the bulk behavior of an aggregate. By contrast, PFC does not require the direct input of such macro-mechanical parameters. Instead, it can be populated with micromechanical parameters, relevant for describing the behavior of grain-to-grain contacts, through the selected contact model. These micromechanical parameters are adjusted to simulate behavior in accordance with the macro-mechanical properties of the aggregate (Mehranpour et al., 2018;Potyondy & Cundall, 2004). Therefore, using a physically realistic contact model plays a key role in properly modeling deformation behavior in PFC, such as by including pressure solution (Bernabé & Evans, 2014) or stress corrosion cracking (Potyondy, 2007). In the case of the Slochteren sandstone, the contact model should represent the grain-scale mechanisms controlling compaction, including irreversible/inelastic processes as observed in experiments (e.g., clay rim deformation). Here, we use experimental data obtained on Slochteren sandstone (Pijnenburg et al., 2019a) to constrain the macro-mechanical properties of our simulated aggregates so as to obtain a new contact model describing the observed grain-to-grain interaction processes in terms of appropriate micromechanical parameters.

Aggregate Generation and Boundary Conditions
Aggregate generation is an important step in the PFC modeling procedure (see Potyondy & Cundall, 2004). Based on a given particle size distribution and porosity, spherical particles can be generated inside a model vessel of predefined size and shape (cuboid for this research). To achieve typical sandstone porosity values in the range of 10%-30%, particle overlap (i.e., particle penetration) is inevitable for spherical assemblies, which leads to intergranular forces. To remove these artificial intergranular forces, the relative normal deformation between each pair of overlapping particles after particle generation was numerically removed by setting their contact reference gap ( r E  ) equal to the current contact gap (δ), which is the negative value of the current relative normal deformation (i.e., n E d    ). For each pair of particles, relative normal deformation is measured by subtracting the distance of particle centers (d) from their average diameter (L). The reference gap comes as an extra term into this equation (i.e., n r E d d L     ) to help us set the relative normal deformation equal to zero for the initial state of aggregate. After this step the new contact model was assigned to the contacts and the desired boundary conditions were applied.
For the calibration of the model to experimental data obtained for Slochteren sandstone, three different aggregates were generated with porosities of 13.4%, 21.5%, and 26.4% to simulate the samples used by Pijnenburg et al. (2019a). Since no grain size distributions were given for these samples, we selected grain size distributions measured for samples from similar depths and the same location (Zeerijp-3a well) obtained in a separate study (Hol et al., 2018). The grain size distributions for the three different samples are shown in Figure 3a. Using a wide particle distribution dramatically increases computational time. Therefore, grain sizes below the 20th percentile and above the 80th percentile were excluded from our simulations. Within the included size range (shaded area in Figure 3a), particle size distributions are near-linear, suggesting a uniform grain size distribution. Therefore, it was assumed that particle sizes are uniformly distributed with a minimum diameter (D min ) of 130 , 110, and 150 μm and a maximum diameter (D max ) of 330, 300, and 390 μm for aggregate porosities of 13.4%, 21.5%, and 26.4%, respectively. The simulated aggregates have a height (H) of 6.25 mm and a square cross-section with sides of 2.5 mm length (L, i.e., H:L = 2.5-Figures 3b-3d). The average coordination number varies from 10.52 to 9.55 and 8.97, for porosity values of 13.4%, 21.5% and 26.4%, respectively. 10.1029/2021JB021722 7 of 23 Boundary conditions representing hydrostatic loading/unloading and deviatoric loading/unloading were simulated by displacing the walls bounding the simulated aggregate, taking the walls to be frictionless. Specifically, fixed displacement increments were applied in the vertical direction, while virtual servo-control displaced the lateral walls in a way that maintained the lateral stresses either equal to the axial stress, for hydrostatic loading/unloading, or equal to a constant value for deviatoric loading/unloading. Note that the principal stresses acting at the sample scale are given as i E  , where, i = 1, 2, or 3. In the present study, only hydrostatic ( i c E P   ) and axisymmetric ( 1 compressions are simulated, with c E P being the confining pressure in the terminology of triaxial testing. Moreover, the pore pressure is taken as zero, so that all normal stresses are equivalent to Terzaghi effective stresses and poroelastic effects of fluid pressure are neglected. Sample wall stresses are calculated by dividing the summed particle forces acting normal to each wall by total wall area. Stress in each principal direction is given as the average stress exerted on opposite walls in that direction. Principal strains are calculated from the ratio of the relative displacement imposed between opposite walls to the initial wall separation.  Hol et al. (2018) to represent the three Slochteren sandstone samples studied by Pijnenburg et al. (2019a), which form the basis for our model calibration. Note that the GSDs vary with porosity (13.4%: solid green line; 21.5%: blue dashed line; and 26.4%: red dashed line). Only the GSDs included in the shaded area were used in our simulations (20th-80th percentile). Simulated sandstone aggregates with (b) a porosity of 13.4%, a D min of 130 μm and a D max of 330 μm (4,435 particles and 20,570 contacts), (c) a porosity of 21.5%, a D min of 110 μm and a D max of 300 μm (5,586 particles and 23,845 contacts), and (d) a porosity of 26.4%, a D min of 150 μm and a D max of 390 μm (2,315 particles and 8,957 contacts).

Slochteren Sandstone Contact Model (SSCM)
The linear elastic contact model and the Hertzian contact model, together with a frictional slip law, are commonly used in DEM for describing grain-to-grain interactions in loose sand aggregates or poorly cemented sandstones (Marketos & Bolton, 2009;B. Wang et al., 2008). For well-cemented sandstones, bonded contact models with linear elastic normal and shear properties only, are typically used (Cheung et al., 2013). Since compaction of Slochteren sandstone is significantly influenced by intergranular clay layers, we assume that the particles representing the quartz grains in the model sandstone are unbonded, i.e. we assume negligible clay film cohesion in tension and shear. In addition, our Slochteren Sandstone Contact Model (SSCM) allows for elastic deformation of the sand grains, as well as irreversible consolidation and shear of the intergranular clay films and quartz grain failure, in line with experimental observations (Pijnenburg et al., 2019a(Pijnenburg et al., , 2019b. Note that we define axial strain as e L L a   / 0 , and volumetric strain as e V V v   / 0 . For the very small strains obtained in our models, the logarithmic strain ε ≈ engineering strain e.

Elastic Behavior and Cohesive/Frictional Strength of Contacts
During hydrostatic loading, Slochteren sandstone initially exhibits a linear relation between mean stress p and elastic volumetric strain e v E e, followed by power-law behavior with an exponent of two Figure 2), which cannot be described using a typical Hertzian contact model for which Goddard, 1990;Johnson, 1985). For the linear elastic part of our contact model, the elastic normal stiffness ( el n E k ) is considered as a spring, represented by an elastic cylinder with an initial radius ( ini E r ) and length L, defined as the average diameter of two adjoining particles (see Figure 4a).
When the relative elastic normal displacement ( el n E d ) exceeds a critical value ( c E  ), the elastic normal stiffness increases linearly with increasing normal displacement through the coefficient el E C . Increasing elastic normal stiffness can be interpreted physically as increasing contact area between two particles (see two shaded cylinders in Figure 4a). As the elastic normal stiffness increases linearly with the relative normal displacement, normal force increases quadratically, leading to nonlinear elastic behavior. To model the complex elastic behavior of Slochteren sandstone, the elastic normal stiffness ( el n E k ) in the new contact model is defined as, where, g E E is the Young's modulus of the particles, L is the average diameter of contacting particles  is the critical effective overlap, specifying the transition from linear elastic behavior to nonlinear elastic behavior, ini E r is the initial radius of active contact area, r E R is the initial radius ratio of active contact area (R r r r ini  / ),r is the minimum radius of the and el E C is the elastic coefficient (see Figure 4a).
The linear and nonlinear elastic, contact-normal behavior predicted using this SSCM is schematically shown in Figure 4b. Note that the shear behavior of contacts is represented as a spring (i.e., elastic) and slider (i.e., inelastic) configuration, with the elastic shear stiffness being taken as a fraction of the normal elastic stiffness, through a stiffness ratio of r E k ( k k n el s el / ). Shear force development is limited through a friction coefficient μ, which in effect represents the friction coefficient for sliding on the intergranular clay films. Since in our SSCM, particles are unbonded, the contacts cannot bear any tensile force (see behavior for negative values of el n E d in Figure 4b) or any shear force when the normal force is tensile.

Inelastic Consolidation Behavior of Contacts
Pijnenburg (2019a) Here n E   is a reference normal stress of 5 MPa (Brown et al., 2017) and γ is a consolidation constant, determined to be about 0.05 for illite compaction If we assume that the clay layers within grain contacts are disks with an initial thickness c E T and an area equal to the initial active contact area ( ini E A ), the inelastic normal stiffness ( inel n E k ) can be calculated using is the inelastic thickness change of the intergranular layer. This inelastic normal stiffness acts in series with the elastic normal stiffness, but is only activated when the compressive normal force at the contact becomes greater, or equal to, the highest, previously applied, normal force ( max n E F ). During unloading (incremental relative normal displacement n E d  > 0), or during loading ( n E d  < 0) when the normal force at the contact max n n E F F  , the normal contact stiffness ( n E k ) is equal to the elastic normal stiffness (see Figure 5). In other words, the contact normal stiffness is defined as It is important to note that elastic deformation of the thin clay films is negligible compared to the sand grains, so is not considered explicitly in our SSCM.

Grain Failure Model
Implementing grain failure in PFC is challenging, due to the difficulty of implementing an appropriate fracturing mechanism and accounting for the post-failure behavior. Hertzian elastic contact theory plus a Griffith failure criterion are typically used to account for grain contact failure (Brzesowsky et al., 2014;B. Wang et al., 2008;Zhang, Wong, Yanagidani et al., 1990). However, since most grain contacts are point or extended contacts containing compressible clay films (Waldmann, 2011), Hertzian contact theory is deemed inappropriate for our SSC model. The quadratic elastic behavior of Slochteren sandstone underpins this. Indeed quadratic behavior is often inferred to be associated with sharp (conical or wedge-shaped) grain 10.1029/2021JB021722 10 of 23 contact indentation (Goddard, 1990). Accordingly, we use an indenter model suitable for both sharp and blunt (long/indented) contact points, as previously described by Swain and Lawn (1976). In this case, the length of a stably growing, intragranular median/cone crack length ( E C) is a function of the applied normal force according to where, n E F is the normal force acting on the contact, IC E K is the mode I fracture toughness and κ is a constant representing the fracture coefficient, which depends on the indenter tip angle (i.e., apical angle) and Poisson's ratio of the grain (Lawn & Fuller, 1975;Swain & Lawn, 1976) i.e. a sharper indenter tip angle results in higher κ. In the grain failure component of our SSCM, the indenter particle is chosen randomly among any two particles in contact. This means that the other particle is the indented particle through which the fracture propagates. Given that Equation 6 assumes stable fracture growth, grain failure is taken to occur when two median cracks within a particle intersect, i.e. when two intragranular cracks reach a length equal to or greater than the particle radius, so-called critical cracks.
To date, several different approaches have been used to model post-grain failure using PFC (Cil & Alshibli, 2014;Couroyer et al., 2000;Fu et al., 2017;Laufer, 2015;Marketos & Bolton, 2009;Sun et al., 2018;B. Wang et al., 2008). In our study, the particle removal approach (Couroyer et al., 2000) is chosen because microstructural observations Pijnenburg et al., 2019b) show extensive grain fragmentation upon failure, hence loss of local load supporting capacity provided subsequent strains are not too large, and because this approach is computationally efficient.

Calibration of the New Contact Model
Our grain contact model contains five micromechanical parameters for elastic deformation ( g E E, Young's modulus; r E R, initial radius ratio; el E C , elastic coefficient; c E  , critical effective overlap; r E k, stiffness ratio) and five micromechanical parameters for inelastic deformation (μ, friction coefficient; γ, inelastic clay consolidation coefficient; c E T, initial clay thickness; IC E K , mode I fracture toughness; κ, fracture coefficient). In order to obtain the values and/or interdependencies of these, we calibrated our model against the triaxial tests performed on Slochteren sandstone by Pijnenburg et al. (2019a), using the procedure below.

Calibration Approach
Calibration is a trial and error procedure in which micromechanical parameter values vary iteratively to fit numerical modeling results with corresponding experimental results. Therefore, calibration procedures can be laborious and tedious, if a precise and explicit procedure is not chosen (Mehranpour & Kulatilake, 2017). The calibration of the SSC model is divided into three behavioral components seen in the macroscopically observed mechanical behavior of the Slochteren sandstone: (a) purely nonlinear elastic behavior (unloading at Stages 1 and 2), (b) inelastic behavior prior to macroscopic yielding (i.e., inelastic deformation of intergranular clay films in Stages 1 and 2), and (c) inelastic behavior beyond macroscopic yielding (grain failure and sliding in Stage 3). To obtain the micromechanical parameter values relevant for each type of behavior, we fitted our SSC model to the mechanical data obtained in the stress-cycling experiments performed on Slochteren sandstone with three different porosities by Pijnenburg et al. (2019a). Given the fact that quartz grains are the major component of the Slochteren sandstone (up to 90%) (Waldmann et al., 2014), we assumed all particles are quartz. Therefore, representative Young's modulus g E E , stiffness ratio r E k, and fracture toughness IC E K values for quartz were used, as summarized in Table 1 (Atkinson, 1979;Potyondy & Cundall, 2004;Wong & Wu, 1995).
It should be kept in mind that the experimental data for each porosity was obtained during three different experiments. Though Pijnenburg et al. (2019a) took great care to limit sample variability, by recovering samples for a specific porosity-range from the same depth, slight variations in rock texture, and hence, mechanical behavior, were inevitable. Therefore, during our model calibration, we aimed to simulate the average of the behavior of different samples with the same porosity, rather than attempt to model each experiment individually. This enables us to develop a model able to predict the mechanical behavior of "averaged" Slochteren sandstone of a specific porosities and grain size distribution.  Figure 2.

Elastic Behavior of Grain Contacts
The elastic component of our SSCM was calibrated using the hydrostatic stress-cycling stages of the nine experiments performed on Slochteren sandstone shown in Figure 2. Inelastic contributions to model fitting ( Figure 6) were avoided by setting c E T in Equation 4 equal to zero and IC E K (Equation 6) equal to infinity. The value chosen for contact friction coefficient (μ) did not affect the modeling results for hydrostatic loading, demonstrating that intergranular shear displacements were too low and homogeneous to activate contact slip.
The linear elastic component of the contact behavior was calibrated by setting the critical effective contact overlap c E  to infinity in Equation (1), i.e. by forcing only linear elastic behavior to occur, and fitting the model simulations to the experimental data by varying r E R (contact initial radius ratio). The nonlinear component of elastic behavior was fitted by adjusting the micromechanical parameters c E  and el E C (elastic coefficient). The calibrated values for the elastic micromechanical parameters are summarized in Table 1. Comparing the fit of our predictions to the experimental data ( Figure 6), it can be seen that the errors, taken at the highest attained mean effective stress of each experimental data curve, and averaged for each sample porosity amount to −0.38%, 1.47%, and 0.32%, with the standard deviation values are 4.21%, 3.37%, and 2.06% for sample porosities of 13.4%, 21.5%, and 26.4%, respectively.

Clay Film Deformation
To calibrate the model for the inelastic behavior resulting from clay film deformation, i.e. to obtain γ and c E T in Equation 4, any contributions by grain failure were inhibited by again setting the grain fracture toughness IC E K to infinity. Elasticity was incorporated using the parameters determined above. The model was calibrated against the same nine cyclic hydrostatic loading experiments employed above (thus precluding any effects of intergranular slip/friction), as shown in Figure 7. Comparing our predictions of permanent strain after unloading in each cycle to that measured in the experiments shows that the errors, averaged for the five experimental load cycles presented at each porosity, are −1.24%, −1.07%, and 3.58% for sample porosities of 13.4%, 21.5%, and 26.4%, respectively, with the standard deviation values of 6.38%, 3.44%, and 5.11%, respectively. During calibration, first γ was obtained to match the inelastic deformation ratio between cycles, after which the clay thickness c E T was varied to match the magnitude of inelastic deformation. The calibrated values for γ and c E T are summarized in Table 1. Note that we assumed γ and c E T to be the same for all contacts within an aggregate.

Intergranular Sliding and Grain Failure
The inelastic behavior of Slochteren sandstone, due to intergranular sliding and failure along localized shear planes, is governed by the friction between grains (controlled by μ) and the failure strength of the grains (controlled by IC E K and the fracture coefficient κ in Equation 6). Given the assumption that all grains consist of quartz, K is chosen to be equal to 1.0 MPa.m 0.5 (Atkinson, 1979). Both μ and κ were calibrated against deviatoric loading experiments, performed at a confining pressure of 20 MPa (see Pijnenburg et al., 2019a). This specific confinement was selected as it covers both compactive and dilatant behavior, enabling us to precisely calibrate both μ and κ under conditions where grain failure occurs but the number of failed grains and therefore the number of removed particles is low enough to be sure the model remains valid. The calibrated values for μ and κ are summarized in Table 1, while the model fitting results are shown in Figure 8. The errors of calibrated models to predict the strength of samples are 0.45%, 1.21%, and 1.01% for sample porosities of 13.4%, 21.5%, and 26.4%, respectively.

Inter-relationships Between Calibrated Parameters and Porosity
In the calibration of our SSC model process three parameters were assumed to be constant ( g E E , r E k, and IC E K ) on the basis of literature data, while the remaining seven parameters were calibrated using the experimental data of (Pijnenburg et al., 2019a) (Table 1). As seen from Table 1, the inelastic clay consolidation coefficient (γ) was found to be independent of porosity, as was the initial contact radius ratio ( r E R) within a few percent. The remaining five micromechanical parameters ( el E C , c E  , c E T, κ and μ) display a clear dependence on aggregate porosity.
Furthermore, some degree of interdependency between the fitted parameters emerged from the calibration procedure (see Figure 9), which we believe reflects real physical effects. First, as suggested in Section 3.2.1, the nonlinearity in elastic behavior seen in Figures 2 and 6 is likely related to the roughly conical/indenter shape of grain contacts. In general, higher porosity sandstone will be characterized by sharper contact points, with the result that the elastic coefficient ( el E C ) will be lower. Therefore, contact stiffness will increase less rapid with increasing contact force than in low porosity material where contacts are blunter. Based on the elastic foundation model (Johnson, 1985), the elastic coefficient values emerging from the calibration process imply cone angles of ∼165° for a sample with 30% porosity to ∼175° for a sample with 10% porosity. Furthermore, sharper contacts mean that grains fail at a lower force, which through Equation 6 means that the fracture coefficient (κ) should increase with increasing porosity. As a result, the elastic coefficient ( el E C ) shows not only a clear correlation with porosity ( Figure 9a) but also with the fracture coefficient (κ- Figure 9c).
As reported in petrographic studies of the Slochteren sandstone, clay content also increases with porosity (Waldmann, 2011). The accompanying increase in intergranular clay film thickness ( c E T) means that intergranular rearrangements/displacements are easier at the start of the compaction process when contact stresses are still low, and that larger sample-scale strains and larger critical overlap ( )  In addition, increasing the clay thickness between grains results in reducing the friction coefficient between the grains, as the chance for quartz asperities present in the grain contact to interact with each other decreases. Once a critical clay thickness is achieved (∼5 μm), corresponding to a porosity in the Slochteren sandstone of ∼20%, the friction coefficient stays constant (µ = 0.3) even if the clay thickness increases (see Figures 9e  and 9b), similar to observations made for rock joint surfaces (Indraratna et al., 2008).
It is worth noting that the porosity independency of inelastic clay consolidation coefficient (γ) and initial contact radius ratio ( r E R) is also expected from the microstructure of the Slochteren sandstone. As is mentioned earlier, in the Slochteren sandstone, the clay films consist of illite. Therefore, inelastic clay consolidation coefficient value should remain constant for different samples despite their porosity. The value that we obtain for this parameter is 0.06, which is closely similar to the values obtained from macroscopic compaction experiments on illite clay by Brown et al. (2017). For the initial contact radius ratio, the calibrated values lie in the range 0.16-0.17, meaning an initial contact area of ∼2.5%. This is consistent with the contact geometries observed in the microstructure of the Slochteren sandstone by Waldmann (2011), which indicated predominantly point contacts.

Quality of Calibration
To speed up the calibration procedure for intergranular sliding and grain failure, we modeled monotonic loading. However, the experimental data (Pijnenburg et al., 2019a) were based on incremental cyclic loading tests, i.e. load cycling to progressively higher stress levels. Although, Pijnenburg et al. (2019a) showed that the cyclic and monotonic loading behavior of similar samples are very similar in terms of the enveloping stress-strain behavior, we still need to verify our model under cyclic loading. Therefore, incremental cyclic loading tests were simulated, using the same micromechanical parameter values presented in Table 1. These simulations were performed not only for a confining pressure of 20 MPa, i.e. the pressure that we used for calibration, but also for a confining pressure of 40 MPa to check the accuracy of our model. These additional simulations show the new model can properly reproduce experimental data ( Figure 10). Furthermore, comparison shows that a curve enveloping the cyclic loading modeling results ( Figure 10) is very similar to that of the monotonic loading results (Figure 8).
It should be noted that due to the limitation of the method with respect to modeling widespread grain failure, our model is not valid for the highest differential stresses achieved in the experiments. For example, we did not simulate the post-failure portion of experiments performed at a confining pressure of 40 MPa (or higher). Roughly speaking, the stress level at which the numerical model broke down is the same as the macroscopic yield stress determined from laboratory experiments. Beyond the macroscopic yield point, gradually deformation is increasingly controlled by grain sliding and grain failure. Therefore, although the current model has limited predictive capacity for higher confining pressures and mean stresses beyond yield at higher strains, it can still accurately predict the macroscopic yield point, i.e. the onset of significant grain failure, for these conditions.  and friction coefficient μ on intergranular clay thickness. Exponential or linear least-square fit to the data are indicated for each relationship. Since the porosity within the reservoir roughly varies between 10% and 30%, our correlations, are presented for this range.

Grain Failure Predictions Versus Experimental Microstructures
Our model allows prediction of the failure modes expected as particles start to break at the onset of Stage 3. For the simulations presented in Figure 10, the frequency distributions of critical crack orientations agree favorably with microstructural observations (Figure 11). Critical cracks in each particle are those two cracks which result in particle failure (see Section 3.2.3). In Figure 11, the frequency distributions of two crack angles are presented. First, the angle between the two critical cracks within a single grain (i.e., the angle of failure; Figures 11a-11c) and second, the acute angle between a critical crack and the axial (maximum) loading direction (Figures 11d-11f). As can be seen, by increasing the confining pressure from 20 to 40 MPa more grains fail by fractures coming from adjacent contacts (i.e., more acute angles), rather than opposite ones (Figures 11a-11c). In addition, increasing confinement leads to more fractures originating from the contacts that are perpendicular to the axial loading axis (i.e., parallel to lateral loading direction-see Figures 11d-11f). These results agree well with the experimental observation reported by Pijnenburg and Spiers (2020) for Slochteren sandstone.

Effect of Model Scale on Simulated Behavior
For DEM modeling, model length-scale (size) is important to eliminate any randomness in particle generation and boundary effects, e.g. the bulk of the central grains should no longer experience interference from the model walls (Mehranpour & Kulatilake, 2016). However, larger digital samples come at significant computational cost. We studied the effect of model-scale by increasing the dimensions of our model by a factor Figure 10. Comparison between the experimental data (dashed lines; Pijnenburg et al., 2019a) and the calibrated SSC model (solid lines) for Slochteren sandstone with porosities of 13.4% (green), 21.5% (blue), and 26.4% (red). Note that cyclic hydrostatic loading was performed by incrementing the confining pressure by 5 MPa per cycle, while cyclic deviatoric loading was performed using 10 MPa increments. Figure 11. (a-c) The distribution of the angle of failure and (d-f) the distribution of the critical fracture angle with respect to the loading axis, i.e. the axial stress, at the maximum stress obtained in our simulations (cf. Figure 10).

Figure 12.
Simulated compaction behavior for an aggregate with a porosity of 26.4% subjected to a confining pressure of 20 MPa followed by axial compression. The standard digital samples used in this study had dimensions of 2.5 mm width and depth, and 6.25 mm height (blue lines). The effect of size was studied in the simulation depicted here by increasing the model dimensions by a factor of 1.5 (red lines). Note that there is virtually no effect, implying that the standard dimensions were sufficient to eliminate significant boundary and statistical effects. of 1.5, and by simulating the behavior of an aggregate with 26.4% porosity subject to a confining pressure of 20 MPa (micromechanical property values as shown in Table 1). For an increase in model dimensions by a factor of 1.5, the compaction behavior of the aggregate does not significantly change ( Figure 12).

Comparison of the Model With Independent Uniaxial Compaction Tests on Slochteren Sandstone
We have developed and calibrated a new contact model to simulate the compaction behavior of the reservoir rock constituting the seismogenic Groningen gas field, i.e. the Slochteren sandstone (SSCM). Our new discrete element model is able to simulate all key aspects of behavior seen in triaxial lab experiments on the sandstone (Pijnenburg et al., 2019a), namely combined linear and nonlinear elastic behavior, inelastic deformation due to intergranular clay film compression and shear, as well as the effects of grain failure and sliding at high differential stresses.
To validate our SSC model calibrated as described above, and to assess its more general applicability to the uniaxial strain path expected in the Groningen gas field during depletion, we simulated an independent set of uniaxial compaction experiments on Slochteren sandstone reported by Hol et al. (2018). These uniaxial compaction experiments, were performed at the initial (pre-depletion) reservoir in-situ stress conditions We simulated a total of seven uniaxial compaction tests performed on samples with porosities in the range 10%-30%. Particle size distributions of the model were based on grain size distributions provided by Hol et al. (2018). Overall, the GSD can be divided into three groups of high (>23.5%), intermediate (19%-23.5%), and low porosity (<19%), such that for each porosity range a single GSD can be selected for our simulations (black, dashed lines in Figure 13). From this averaged GSD, the portion between the 20th and 80th percentile was used in the simulations, similar to the model calibrations (Figure 3a).
The aggregates were generated according to the above specifications for porosity and GSD. The micromechanical properties of these aggregates were specified using the calibrations and correlations presented in Figure 9. Since, our model cannot simulate pore pressure depletion, we emulated this by increasing the effective stresses acting on the aggregate, in line with the increase in effective stress observed in the experiments, while maintaining a zero lateral displacement boundary.
Our simulations indicated that substantial grain failure occurs at high effective axial stresses (>60 MPa) in highly porous (>20% porosity) samples. Notably, for samples with 30% porosity, pervasive grain failure was Figure 13. Grain size distributions (GSDs) for Slochteren sandstone, as provided by Hol et al. (2018). Note that the GSDs can be divided into three groups based on their porosity. The black, dashed lines indicate the averaged GSDs, with the GSD in the gray areas (20th to 80th percentile) being used for our simulations, for each porosity range.
predicted at an axial effective stress of about 60 MPa, i.e. following a pore pressure drop to 10 MPa, leading to a break-down of the numerical model beyond this point. The depletion-induced compaction rate (defined as the compaction strain produced per 1 bar or 0.1 MPa reduction in pore pressure or increase in axial effective stress) measured in the laboratory experiments (blue circles) and computed in our simulations (red crosses), after 32 MPa of pore pressure depletion (67 MPa effective axial/vertical stress), are shown in Figure 14. As can be seen, the modeling results predict depletion-induced compaction rates that are slightly higher than the average rates measured in laboratory tests. Nonetheless, both data sets fall in the same range and both exhibit a similar increase in compaction rate with increasing porosity.

Comparison With Field Data and Future Perspectives
The ultimate aim of our modeling effort is to provide a microphysical basis to simulate compaction of the Groningen reservoir and perhaps other, similar reservoir rocks too. Given the sheer size of the giant Groningen gas field, there will be significant lateral and horizontal variations in porosity, grain size distribution and composition related to changes in depositional environment. A discrete element model of the present type can help overcome limitations such as the number of available samples and operational costs, which make experimental investigation of a broad range of samples, covering the whole field, infeasible. Coupling experimental observations and microphysical theory with numerical modeling capability, as attempted in this study using PFC, augmented with our new contact model, enables us to perform an unlimited number of numerical simulations under any boundary conditions and for porosities, clay contents and particle size distributions representing the full range observed within the Slochteren sandstone. These simulations can form a basis for studying of the compaction behavior of the reservoir at local or full gas field scales.

Comparison of Reservoir Compaction Simulations With Field Measurements
To take a first step toward applying our model to the Groningen reservoir and to validating it at this scale, we have compared our model predictions to in-situ compaction data obtained from the Zeerijp-3a well at the center of the field. This well is equipped with Distributed Strain Sensing (DSS) technology, which allowed continuous strain monitoring during the period April-November 2016 (7 months), as described by Cannon and Kole (2017). Along the well, the Slochteren sandstone reservoir layer was segmented into 31, ∼7 meter long, DSS measurement zones or blocks, covering the entire ∼210 m reservoir section. The compaction strains characterizing each block were obtained by averaging the internal strain values taken over the 7-month test period. During this period, the pore pressure could not be measured directly but was estimated, from empirical relations, to have dropped by about 0.15 MPa (1.5 bar) against a background reservoir pore pressure of roughly 8 MPa (Cannon & Kole, 2017). The data obtained are plotted in Figure 14 in terms of depletion-induced strain per bar of pressure drop versus the averaged well-log porosity for each 7 m measurement interval of the reservoir.
As seen from Figure 14, our model predictions for uniaxial strain boundary conditions fall at the lower limit of the field data on depletion-induced compaction per bar pore pressure drop. Most of the laboratory-derived, depletion-induced compaction rate data also fall just below the field rates (30 out of 34 data points). These discrepancies between field data versus lab data and modeling predictions could be related to the short period of in-situ field data collection, which introduces significant uncertainty into both the strain measurements and the pore pressure change estimates. Furthermore, appreciable uncertainties may be present in the average porosity value for each 7 m well section. Since the relation between porosity and Figure 14. Depletion-induced compaction rates (amount of strain per bar pore pressure depletion) as a function of porosity obtained from uniaxial compaction experiments (blue circles; Hol et al., 2018) and numerical simulations (red and orange crosses; this study), as well as from in-situ field data (green squares; Cannon & Kole, 2017). Note that the simulation predictions are performed for two different initial stress conditions: red crosses, initial stress conditions as employed by Hol et al. (2018); and orange crosses, slightly lower initial effective stress conditions, taking into account measurement uncertainty. Simulations undergoing pervasive grain failure prior to reaching the maximum axial stress, leading to breakdown of the model, are indicated with a large red circle.

10.1029/2021JB021722
19 of 23 depletion-induced compaction rate is highly nonlinear, using an average value introduces another level of error, especially when high porosity sandstone layers exist in an interval. The slower rate of pore pressure drop at the field can be another reason to have a higher compaction rate. At a slower rate, time-dependent mechanisms such as stress corrosion cracking and/or pressure solution creep can result in higher permanent deformation (Pijnenburg et al., 2018).
An additional source of uncertainty affecting the experimental and modeling results shown in Figure 14 stems from the uncertainty in the value of the minimum effective horizontal stress ( ,min eff h E  ) within the reservoir before depletion started. In-situ measurements suggest a range of values between 8 MPa (Breckels & van Eekelen, 1982) and 22 MPa (Hol et al., 2018). An analysis of the effect of initial stress conditions on the depletion-induced compaction rate was explored using our model (red crosses and orange stars in Figure 14). Note that we were not able to investigate the lower ,min eff h E  value of 8 MPa, due to premature failure of the aggregate with high porosities prior to achieving the initial, in-situ stress state. Moreover, from the experimental results reported by Pijnenburg et al. (2019b), for a Slochteren sandstone sample with 26.4% porosity, it can be inferred that the strength of this sample is lower than 30 MPa under the confining stress of 8 MPa. This suggests that the 8 MPa estimate of effective minimum horizontal stress may be an underestimate of the true ,min eff h E  . For our simulations, a minimum horizontal stress of at least 10 MPa was needed to stabilize sample behavior prior to simulating pore pressure depletion. Crucially, the stress boundary conditions have an increasingly large impact as aggregate porosity increases from 13% to 26%, leading to an almost 50% increase in predicted depletion-induced compaction rate. This underpins the need for accurate measurement of in-situ stress state in the Groningen reservoir.

Future Perspectives
In this paper, a new contact model has been proposed, which can simulate important features of sandstone compaction behavior based on physically observed grain-scale mechanisms. Aspects captured include nonlinear elastic behavior combined with time-independent permanent deformation up to and beyond macroscopic yield. To the best of our knowledge, this is the first such development in the context of DEM and to understand the properties and stability of the model fully, a more complete sensitivity analysis should be performed in future work. Although the contact model is developed for Slochteren sandstone, it can also be applied to other sandstones displaying similar behavior due to the operation of similar grain-scale mechanisms. For the new contact model, obtaining accurate values for the micromechanical parameters characterizing the grain-scale processes included is crucial. This depends primarily on calibration using experimental data obtained for sandstone samples of known porosity, mineralogy and/or particles size distribution, and on microstructural studies to verify that the micro-scale mechanisms that govern mechanical response are included in the contact model.
It should be noted that DEM also allows the generation of layered sandstone samples to investigate how layers with different properties influence the bulk sample behavior. This also offers a way of investigating how choices in representative elementary volume size may affect the results of large scale numerical modeling of compaction of layered sandstone reservoirs.
Another major benefit of employing DEM is in the application of numerous boundary conditions and stress paths in repeated analyses of identical "digital" samples. In future studies, we aim to apply our new contact model to explore the mechanical behavior of "digital" Slochteren sandstone samples under a wide range of loading paths relevant for the entire history of gas production from the Groningen reservoir, i.e. since depletion started. This investigation will help place better constraints on the initial in-situ stresses in the Groningen reservoir (and in some of the smaller, stratigraphically equivalent sandstone reservoirs nearby) by comparing different model scenarios with measurements of the current in-situ stress state, residual gas pressure and surface subsidence, for example. We will also be able to better predict the future mechanical response of the Groningen and surrounding smaller reservoirs to any continued gas extraction, and, in the case of the nearby smaller fields, to injection of CO 2 or hydrogen for geological storage.
Generating a wide range of "digital" samples and modeling them under numerous boundary conditions will also help us generate a comprehensive data bank. This offers opportunities such as fitting analytical 10.1029/2021JB021722 20 of 23 constitutive models to the simulated data to obtain computationally efficient strain hardening deformation laws and macroscopic yield or failure criteria for the sandstone, in the form of fully general 3-D descriptions.

Conclusions
In this study, we developed a new DEM to explicitly account for mechanisms of nonlinear elasticity, intergranular clay film deformation, and grain breakage, as observed to operate during deformation of Slochteren sandstone of the seismogenic Groningen gas field in the Netherlands. SSCM allows for modeling the compaction behavior of Slochteren sandstone in the porosity range 10%-30%. We employed Particle Flow Code as the DEM platform, using rigid spheres representing the sandstone grains. To our knowledge, this is the first time that observed microphysical mechanisms have been incorporated into DEM. Overall, we can conclude the following: 1. The SSC Model requires five micromechanical parameters to describe the elastic deformation of the Slochteren sandstone ( g E E , Young's modulus; r E R, initial radius ratio; el E C , elastic coefficient; c E  , critical effective overlap; r E k, stiffness ratio) and five micromechanical parameters for inelastic deformation (μ, friction coefficient; γ, inelastic clay consolidation coefficient; c E T, initial clay thickness; IC E K , mode I fracture toughness; κ, fracture coefficient). All parameters, with the exception of γ and r E R, are dependent on porosity. Values for these micromechanical parameters were obtained through calibration against both hydrostatic and conventional triaxial compression experiments performed on Slochteren sandstone, for porosities of 13.4%, 21.5%, and 26.4% 2. Validation of the SSCM against an independent set of pore pressure depletion experiments conducted under uniaxial strain conditions showed good agreement between the experimental data and the simulations. Our simulations suggest significant grain breakage to occur in highly porous samples (>20% porosity) at effective axial stresses of >60 MPa. Predicted depletion-induced compaction rates (i.e., amount of strain per bar of pore pressure depletion) are slightly higher than the average rates measured in the laboratory tests, but fall in the same range and show a similar dependence on porosity 3. Comparison of our simulations, under uniaxial strain boundary conditions, to field data show that predictions fall at the lower limit of the field-measured depletion-induced compaction rate. This discrepancy is inferred to be due to the short period of in-situ field data collection, uncertainties in average porosity over the field measurement interval, and/or uncertainties in the in-situ state of stress 4. In future, the new contact model combined with DEM simulations allows for modeling a wide range of sandstone samples having different mineralogy and petrophysical properties, but similar mechanical behaviors. Furthermore, the SSCM model makes it possible to study the mechanical behavior of a single, "digital", sample under numerous boundary conditions and loading paths. This offers potential to better constrain in-situ stresses. Additionally, with this capability, it is possible to generate a data bank suitable for fitting analytical constitutive models to the simulated data to obtain strain hardening deformation laws and macroscopic yield/failure criteria for sandstones, which can be used for upscaling

Data Availability Statement
The numerical modelling code used for simulations in PFC and the simulation data used for this paper are available through https://doi.org/10.24416/UU01-575EWU. Yoda data publication platform, Utrecht University.