Multi-scale simulation of rock compaction through breakage models with microstructure evolution

Regional subsidence due to fluid depletion includes the interaction among multiple physical processes. Specifically, rock compaction is governed by coupled hydro-mechanical feedbacks involving fluid flow, effective stress change and pore collapse. Although poroelastic models are often used to explain the delay between depletion and subsidence, recent evidence indicates that inelastic effects could alter the rock microstructure, thus exacerbating coupling effects. Here, a constitutive law built within the framework of Breakage Mechanics is proposed to account for the inherent connection between rock microstructure, hydraulic conductivity, and pore compaction. Furthermore, it is embedded into a 1-D hydromechanical coupled finite element analysis (FEA) to explore the effects of micro-structure rearrangement on the development of reservoir compaction. Numerical examples with the proposed model are compared with simulations under constant hydraulic conductivity to illustrate the model capability to capture the non-linear processes of reservoir compaction induced by fluid depletion.


Introduction
Various anthropogenic activities are responsible for the occurrence of surface subsidence and damage to buildings and infrastructure. Among these, rock compaction and land subsidence are major concerns in hydrocarbon reservoirs (Zoback, 2010). Reservoir depletion lowers the fluid pressure, leads to an increase in the stress carried by the rock skeleton, and generates non-negligible compaction. Diagnostic tools are necessary to understand and forecast regional subsidence, but a challenging obstacle in such analyses is the interaction of multiple physical processes operating over a wide range of time and length scales. Poroelastic modeling frameworks are often used to compute either near-field reservoir compaction or far-field subsidence (Du and Olson, 2001). However, poroelastic models restrict the analysis to linear deformation processes, thus ruling out possible inelastic changes in rock microstructure (i.e., particle crushing, particle rearrangement, and pore collapse) caused by compaction and the consequent hydro-mechanical couplings deriving from them (Esna Ashari et al., 2018). In this circum-stance, more advanced constitutive laws are needed to capture the material properties and connect the particle-scale processes to macro-scale fluid flow characteristics. For this purpose, a new constitutive relationship within the framework of Breakage Mechanics is developed and the particular case of the Groningen gas field (Roels, 2001;Hol et al., 2015;van Thienen-Visser and Fokker, 2017) is used to validate the proposed model. Different modeling scenarios (e.g., elastic and inelastic compaction, constant hydraulic conductivity, and concurrent change of compressibility and hydraulic conductivity) are simulated to illustrate the potential impact of inelastic processes on reservoir compaction.

Partially coupled breakage model
Breakage mechanics is a constitutive framework built within the frame of continuum thermodynamics to study the inelastic deformation of granular solids subjected to changes in grain size distribution due to particle crushing. Such approach requires the definition of a Helmholtz free energy ,  which in case of linear elasticity can be expressed as: where ε e v and ε e q are the volumetric and deviatoric elastic strain; B is the breakage index; K and G are the elastic bulk and shear modulus; d m and d M are the minimum and maximum particle diameters; θ is a grading index associated with the initial and ultimate grain size distribution (GSD), p 0 (d) and p u (d) (Einav, 2007a). Here, a Pearson distribution is used to describe mathematically the grading characteristics of a granular rock, eventually using experimental data to constrain the model parameters (Hol et al., 2018).
Generally, brittle granular materials present two coupled dissipative mechanisms: plasticity and breakage. This work simply assumes that coupling exists only in compression, while the shear response is solely governed by frictional plasticity: in which the individual terms are being factorized as δ (Einav, 2007b), and: where E B is the breakage energy. E c is the critical breakage energy. p and q are the mean effective stress and deviatoric stress. ε p v and ε p q are the volumetric and deviatoric plastic strain. M f = 6 sin φ/ (3 − sin φ), q 0 = 6c cos φ/ (3 − sin φ). c is the cohesive shear strength and φ is the effective friction angle. It leads to two separate yielding criteria, the compressive breakage criterion f B and the frictional shear failure f f , here expressed as (Einav, 2007a, b): and it obeys the following flow rules (Das et al., 2011): where B and f are non-negative plastic multipliers.

Simulation of cyclic depletion-inflation paths
Standard plasticity models are usually unable to replicate plastic strains prior to the intersection of a yield surface and neglect any inelastic deformation caused by loadingunloading cycles taking place in the pre-yielding regime. However, ample evidence has proved that plastic strains start to accumulate much earlier than the classic yielding cap (van Thienen-Visser and Fokker, 2017; Pijnenburg et al., 2018;Hol et al., 2018). To circumevent this limitation, hypo-plastic modeling techniques can be incorporated (Einav, 2012). This choice implies the use of flow rules adjusted through multiplicative coefficients taking into account the distance of current stress state from the yield surface: where s B and s f are the hypo-plastic model parameters for breakage mechanics and frictional plasticity, respectively.

Breakage permeability variation
Microfracturing (i.e., particle crushing and cement damage) can have a dominant role on the diffusive processes responsible for macroscopic rock permeability (Zhang and Buscarnera, 2018). Here, the approach suggested by Esna Ashari et al. (2018) is followed to model permeability changes due to inelastic compaction. Specifically, an augmented breakage-dependent Kozeny equation is used: where n and n 0 are the current and initial porosity; k and k 0 are the current and initial permeability; K r and K 0 are the current and initial hydraulic conductivity resulting from K r = ρ f gk/µ f , where ρ f and µ f are the density and viscosity of reservoir fluid, and g is the gravitational acceleration; D H is the current harmonic mean grain size, with D H0 and D Hu being the values of D H corresponding to the initial and ultimate GSD:

Model performance
In this work, data relative to the Groningen gas field were used to validate the model performance (Roels, 2001 ). Each of the stage was followed by 3 unloading-reloading cycles during which the magnitude of the variation in p w was 5 MPa. The model formulation described in the preceding sections postulates elastic behavior during unloading and accumulation of permanent (plastic) strain during loading cycles. Measurements available for the different loading/unloading cycles were used to constrain the model parameters summarized in Table 1. Figure 1 shows the model performance for the UPPD/M tests characterized by different values of porosity. Despite some differences between measurements and simulations, the model generally displays a satisfactory performance for all the considered cases.

Simulation of reservoir compaction
To numerically investigate the effects of material nonlinearity on the development of reservoir compaction, the calibrated model is embedded into a 1-D coupled hydro-mechanical finite element (FE) solver (Sloan and Abbo, 1999). As schematically described in Fig. 2, a reservoir column with thickness of 250 m is subjected to a prescribed basal depletion p b . Figure 2b illustrates the initial and boundary conditions of the simulated reservoir. As shown, it involves the fixed axial stress equal to 70 MPa and the initial radial stress equal to 57 MPa, which is typical across the Groningen gas field. Based on the reported field history, the basal depletion is imposed with an initial pore pressure of 35 MPa and a depletion rate of 0.5 MPa yr −1 . Figure 3a compares the predicted time history of reservoir compaction obtained with a breakage dependent hydraulic conductivity to those obtained under constant hydraulic conductivity. The two end-member values of hydraulic conductivity (i.e., K r = 2.0×10 −7 and K r = 9.0×10 −7 m s −1 ) in Fig. 3d are selected for comparison. Figure 3b-d plots the 25-year isochrones of pore pressure, displacement, and hydraulic conductivity during and after basal depletion for the proposed model. The simulations readily show that microstructural changes due to particle crushing exacerbate hydro-mechanical couplings by lowering the hydraulic conductivity. As a consequence, fluid flow proceeds at a lower rate and the development of reservoir deformation is delayed. As shown in Fig. 3a, the compaction response with the evolving K r is initially close to the behavior of constant high K r , but as the occurrence of breakage, it tends to approach the response obtained for the lowest (constant) value of K r .

Conclusion and discussion
Microstructure evolution caused by stress changes is a major cause of the inelastic deformation and hydro-mechanical coupling. In this paper, a constitutive law has been developed within the framework of Breakage Mechanics to link the micro-scale and the macro-scale response of a compacting rock, as well as to estimate the permeability variations caused by breakage. A 1-D hydromechanical coupled FE solver has been implemented to explore the impacts of microstructure evolution on reservoir compaction. Data relative to the Groningen gas field have been used to illustrate the model characteristics and test its performance. The model has been calibrated and validated with reference to data for rock cores subjected to the UPPD/M testing protocol. Finally, numerical simulations characterized by both constant and breakage-dependent hydraulic conductivity have been conducted to emphasize the non-negligible impact of inelastic effects. Such analyses show that inelastic processes at the grain-scale can significantly decelerate the pore pressure diffusive process, eventually manifesting as additional sources of delayed reservoir compaction. The geomechanical models discussed in this paper can be a springboard for further studies addressing the challenging problem of simulating fieldscale ground subsidence triggering by compaction, for example by integrating near-field simulations with analytical solu-tions (Geertsma, 1973;Mehrabian and Abousleiman, 2015) or numerical modeling platforms (Lewis et al., 2003). Author contributions. GB and JL designed research; YC and GB performed research; YC and JL contributed new numerical tools; YC and RZ analyzed data; GB and YC wrote the paper.