Articles | Volume 382
Pre-conference publication
22 Apr 2020
Pre-conference publication |  | 22 Apr 2020

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

Giuseppe Buscarnera, Yanni Chen, José Lizárraga, and Ruiguo Zhang

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.

1 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 (Zoback2010). 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 Olson2001). 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 circumstance, 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 (Roels2001; Hol et al.2015; van Thienen-Visser and Fokker2017) 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.

Figure 1Simulations of the UPPD/M tests with two end-member porosities (test ID: ZRP-3A_18CV and ZRP-3A_39BV) and the medium porosity (test ID ZRP-3A_99BV).


Table 1Summary of model parameters typical in the Groningen gas field.

Download Print Version | Download XLSX

2 Constitutive modeling

2.1 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:

(1) Ψ = 1 - θ B [ 1 2 K ε v e 2 + 3 2 G ε q e 2 ] , θ = 1 - d m d M p u x x 2 d x d m d M p 0 x x 2 d x

where εve and εqe are the volumetric and deviatoric elastic strain; B is the breakage index; K and G are the elastic bulk and shear modulus; dm and dM are the minimum and maximum particle diameters; θ is a grading index associated with the initial and ultimate grain size distribution (GSD), p0(d) and pu(d) (Einav2007a). 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).

Figure 2Schematic description of model setup: (a) model geometry; (b) initial and boundary conditions.


Figure 3Simulations of a single reservoir compaction during and after depletion: (a) time history of reservoir compaction; (b, c) 25-year isocrones of the nodal pore pressure and displacement; (d) 25-year isocrones of hydraulic conductivity at gauss points.


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:

(2) δ Φ = δ Φ v p * 2 + δ Φ B * 2 + δ Φ q p

in which the individual terms are being factorized as δΦvp*=sin-1ωδΦvp and δΦB*=cos-1δΦB (Einav2007b), and:

(3) δ Φ v p = p E c E B 1 - B 2 δ ε v p , δ Φ q p = M f p + q 0 δ ε q p , δ Φ B = E c 1 - B 2 δ B

where EB is the breakage energy. Ec is the critical breakage energy. p and q are the mean effective stress and deviatoric stress. εvp and εqp are the volumetric and deviatoric plastic strain. Mf=6sinϕ/3-sinϕ, q0=6ccosϕ/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 fB and the frictional shear failure ff, here expressed as (Einav2007a, b):

(4) f B = E B 1 - B 2 E c - 1 0 , f f = q M f p + q 0 - 1 0

and it obeys the following flow rules (Das et al.2011):

(5) d ε v p = 2 Λ B E B 1 - B 2 p E c sin 2 ω , d ε q p = Λ f 1 M f p + q 0 , d B = 2 Λ B 1 - B 2 E c cos 2 ω

where ΛB and Λf are non-negative plastic multipliers.

2.2 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 loading-unloading 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 Fokker2017; Pijnenburg et al.2018; Hol et al.2018). To circumevent this limitation, hypo-plastic modeling techniques can be incorporated (Einav2012). 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:

(6) d ε v p * = f B + 1 s B d ε v p , d ε q p * = f f + 1 s f d ε q p , d B * = f B + 1 s B d B

where sB and sf are the hypo-plastic model parameters for breakage mechanics and frictional plasticity, respectively.

2.3 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 Buscarnera2018). 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:

(7) k k 0 = K r K 0 = n n 0 3 D H D H 0 2 , 1 D H = 1 - B 1 D H 0 + B 1 D Hu

where n and n0 are the current and initial porosity; k and k0 are the current and initial permeability; Kr and K0 are the current and initial hydraulic conductivity resulting from Kr=ρfgk/μf, where ρf and μf are the density and viscosity of reservoir fluid, and g is the gravitational acceleration; DH is the current harmonic mean grain size, with DH0 and DHu being the values of DH corresponding to the initial and ultimate GSD:

(8) 1 D H 0 = d m d M p 0 ( x ) 1 x d x , 1 D Hu = d m d M p u ( x ) 1 x d x
3 Example of application to the Groningen gas field

3.1 Model performance

In this work, data relative to the Groningen gas field were used to validate the model performance (Roels2001). The UPPD/M tests (multicycle Uniaxial-strain Pore Pressure Depletion tests) reported by Hol et al. (2018) are used to calibrate the model parameters. The UPPD/M samples were subjected to three depletion stages (Stage #1: pw decreases from 35 to 25 MPa, Stage #2: from 25 to 15 MPa, and Stage #3: from 15 to 1 MPa). Each of the stage was followed by 3 unloading-reloading cycles during which the magnitude of the variation in pw 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.

3.2 Simulation of reservoir compaction

To numerically investigate the effects of material non-linearity on the development of reservoir compaction, the calibrated model is embedded into a 1-D coupled hydro-mechanical finite element (FE) solver (Sloan and Abbo1999). As schematically described in Fig. 2, a reservoir column with thickness of 250 m is subjected to a prescribed basal depletion pb. 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., Kr=2.0×10-7 and Kr=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 Kr is initially close to the behavior of constant high Kr, but as the occurrence of breakage, it tends to approach the response obtained for the lowest (constant) value of Kr.

4 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 field-scale ground subsidence triggering by compaction, for example by integrating near-field simulations with analytical solutions (Geertsma1973; Mehrabian and Abousleiman2015) or numerical modeling platforms (Lewis et al.2003).

Data availability

The results of the laboratory tests used in this paper are available from the Supplement to the work of Hol et al. (2018) (

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.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “TISOLS: the Tenth International Symposium On Land Subsidence – living with subsidence”. It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17–21 May 2021.

Financial support

This research has been supported by the U.S. Department of Energy (grant no. DE-SC0017615).


Das, A., Nguyen, G. D., and Einav, I.: Compaction bands due to grain crushing in porous rocks: a theoretical approach based on breakage mechanics, J. Geophys. Res.-Solid Earth, 116, B8,, 2011. a

Du, J. and Olson, J. E.: A poroelastic reservoir model for predicting subsidence and mapping subsurface pressure fronts, J. Petrol. Sci. Eng., 30, 181–197, 2001. a

Einav, I.: Breakage mechanics – part I: theory, J. Mechan. Phys. Solid., 55, 1274–1297, 2007a. a, b

Einav, I.: Breakage mechanics – Part II: Modelling granular materials, J. Mechan. Phys. Solid., 55, 1298–1320, 2007b. a, b

Einav, I.: The unification of hypo-plastic and elasto-plastic theories, Int. J. Solid. Struct., 49, 1305–1315, 2012. a

Esna Ashari, S., Das, A., and Buscarnera, G.: Model-Based Assessment of the Effect of Surface Area Growth on the Permeability of Granular Rocks, J. Eng. Mechan., 144, 04018023,, 2018. a, b

Geertsma, J.: Land subsidence above compacting oil and gas reservoirs, J. Petrol. Technol., 25, 734–744, 1973. a

Hol, S., van der Linden, A., Zuiderwijk, P., Marcelis, F., and Coorn, A.: Mechanical characterization of Permian reservoir sandstone from the Moddergat-3 well in the Dutch Wadden Area, Report No. SR, 15, 2015. a

Hol, S., van der Linden, A., Bierman, S., Marcelis, F., and Makurat, A.: Rock physical controls on production-induced compaction in the Groningen Field, Sci. Rep., 8, 7156,, 2018. a, b, c, d

Lewis, R. W., Makurat, A., and Pao, W. K.: Fully coupled modeling of seabed subsidence and reservoir compaction of North Sea oil fields, Hydrogeol. J., 11, 142–161, 2003. a

Mehrabian, A. and Abousleiman, Y. N.: Geertsma’s subsidence solution extended to layered stratigraphy, J. Petrol. Sci. Eng., 130, 68–76, 2015. a

Pijnenburg, R., Verberne, B., Hangx, S., and Spiers, C.: Deformation behavior of sandstones from the seismogenic Groningen gas field: Role of inelastic versus elastic mechanisms, J. Geophys. Res.-Solid Earth, 123, 5532–5558, 2018. a

Roels, H.: Groningen field, past, present and future, Neth. J. Geosci., 80, 12–14, 2001.  a, b

Sloan, S. W. and Abbo, A. J.: Biot consolidation analysis with automatic time stepping and error control Part 1: theory and implementation, Int. J. Num. Anal. Method. Geomechan., 23, 467–492, 1999. a

van Thienen-Visser, K. and Fokker, P. A.: The future of subsidence modelling: compaction and subsidence due to gas depletion of the Groningen gas field in the Netherlands, Netherlands J. Geosci., 96, s105–s116, 2017. a, b

Zhang, Y. and Buscarnera, G.: Breakage mechanics for granular materials in surface-reactive environments, J. Mechan. Phys. Solid., 112, 89–108, 2018. a

Zoback, M. D.: Reservoir geomechanics, Cambridge University Press, Cambridge, UK,, 2010. a