# A simplified vertical and horizontal geomechanical model for compaction in sedimentary basins

Francesco Saverio Patacchini, Marie-Christine Cacas-Stentz, Nicolas Maurand, Walid Saber-Cherif, Françoise Willien

► To cite this version:

Francesco Saverio Patacchini, Marie-Christine Cacas-Stentz, Nicolas Maurand, Walid Saber-Cherif, Françoise Willien. A simplified vertical and horizontal geomechanical model for compaction in sedimentary basins. Science and Technology for Energy Transition, 2023, 78, pp.22. 10.2516/stet/2023019. hal-04387433

HAL Id: hal-04387433

https://ifp.hal.science/hal-04387433v1

Submitted on 11 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
REVIEW ARTICLE

OPEN ACCESS

# A simplified vertical and horizontal geomechanical model for compaction in sedimentary basins

Francesco Saverio Patacchini<fnref n="1" />,*, Marie-Christine Cacas-Stentz<fnref n="1" />, Nicolas Maurand<fnref n="1" />,<fnref n="2" />, Walid Saber-Cherif, and Françoise Willien<fnref n="1" />

Received: 5 January 2023 / Accepted: 18 July 2023

**Abstract.** In the context of mechanical compaction in sedimentary basins, we introduce a simple model including lateral deformations with the goal to improve the results obtained under oedometric conditions (i.e., neglecting horizontal strains) without losing much computational time. The model is based on a modified vertical porosity-stress law where horizontal strains are inserted and on an elastic stress-strain law with stress-dependent Young modulus. Though it is not three-dimensional and does not involve plasticity, we manage to validate the model on a geometrically and lithologically complex test case by comparing our results with those obtained on the same case using a full-dimensional finite-element simulator. We conclude that our model offers a significant improvement in accuracy against an oedometric model, with little loss in computational time, and so provides a useful tool to users who want a quick insight into results before running longer and more accurate simulations.

**Keywords:** Sedimentary basin, Mechanical compaction, Horizontal deformation, Geomechanical model, Finite-volume simulator.

## 1 Introduction

When modeling the evolution of a sedimentary basin, one is often interested in determining the changes in rock porosity and permeability taking place over a geological time interval. Indeed, these quantities control the overpressure distribution in the basin and need to be predicted, for instance, for drilling risk assessment [1]. Since permeability is obtained directly from porosity *via* the Kozeny–Carman law [2–4] or other sediment-specific laws [5, 6], we choose porosity as our main unknown.

Porosity evolution is mainly driven by mechanical and chemical compaction [7]. Usually, mechanical compaction is the dominating phenomenon in the upper layers of the basin, while chemical compaction dominates in the lower layers. The depth transition between mechanically and chemically dominated layers depends on lithology, varying from 0.1 km for carbonates [8, 9] to 1.5 km for sandstone [10–13]. In the present work, we neglect chemical compaction so that our model is best suited, but not limited, to upper layers.

Because mechanical compaction is mostly driven by gravity through burial, it is often described as a vertical,

thus one-dimensional, phenomenon in basin models [14–17] and simulators [18–20]. This approach provides good accuracy in oedometric contexts, where the horizontal strains are negligible. However, since possibly horizontal phenomena having a significant influence on porosity occur [21, 22], e.g., as a result of tectonics, this is not entirely satisfactory.

One solution is proposed in [7], where Schneider’s porosity-stress law [23, 24] is integrated into the Drucker–Prager model for plasticity [25], resulting in a three-dimensional compaction model. This model is tested numerically in [7, 26–28] using the prototype code A<fnref n="2" />, which couples the finite-volume basin simulator ArcTem [11] with the finite-element geomechanical simulator Code Aster [29] following an iterative algorithm [30]. (See [31] for further A<fnref n="2" /> simulations in the context of rock failure.) This approach reaches very good accuracy but requires long simulation times due to its three-dimensional formulation.

Here, rather than derive a three-dimensional model such as in [7], we simply insert the horizontal strains into the vertical porosity-stress law and so keep a one-dimensional compaction model without neglecting the remaining dimensions. More precisely, we define an approximated porosity using a vertical porosity-stress law where the horizontal strains are added to the vertical stress; we then suppose

<footnote n="1">¹IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote n="2">²CEA Cadarache, F-13108 Saint-Paul-Lez-Durance, France</footnote>

<footnote>* Corresponding author: francesco.patacchini@ifpen.fr</footnote>
that the approximated porosity is close to the actual one. For this closeness supposition to be justified, we need to insert the horizontal strains in a physically motivated way. In this paper, we propose one formulation based on an elastic stress-strain constitutive law with stress-dependent Young modulus. We thus suppose that the basin is an elastic medium, as opposed to elastoplastic, and consequently place ourselves in the context of small deformations. We refer to this method as a simplified geomechanical model, where “simplified” contrasts with the more complex three-dimensional, elastoplastic approach of [7]. This simplified model bears a similarity with that investigated in [32], namely, that the vertical porosity-stress law is altered to account for a specific effect: in our case, the three-dimensional mechanical strain; in their case, the anisotropy.

Our main goal is to improve the results obtained with the oedometric assumption, in particular in the more strained lower layers, without losing significant computational time. Of course, as already mentioned, this entails simplifying assumptions, such as purely elastic deformations and a one-dimensional porosity stress law, so that our approach is by no means intended to reach the accuracy of any full-dimensional, elastoplastic model such as presented in [7]. Furthermore, we neglect the effects of the lateral deformations on the temperature; see for example [33] for a complete thermal model. Rather, we expect our model to be useful to users who might need a first, quick idea of the basin’s porosity distribution before deciding whether to pursue more accurate and longer simulations. Our model is integrated entirely into ArcTem.

In Section 2, we recall the standard equations underlying the fluid and solid mechanics of a sedimentary basin; we place ourselves in the setting of a single-phase water flow. This section can be skipped by the experienced geomechanics and can be referred to only for notation. We then discuss in detail mechanical compaction (simply referred to as compaction in the sequel) and present our simplified geomechanical model in Section 3. Then, Section 4 is dedicated to the numerical results obtained with our model coded in ArcTem and applied to a large test case modeling the Vaca Muerta formation of the Neuquén basin in Argentina. This test case, which features high overpressures, covers a surface of about 35,000 km² and spans a time interval of 10My. Our results under both oedometric and non-oedometric conditions are compared with those obtained using A² in [26] with the model derived in [7]. Furthermore, CPU times are given for various strain configurations. Finally, we summarize our results and give an outlook in Section 5.

## 2 Standard porous medium model

We briefly discuss the standard equations that we use to model the flow and mechanical equilibrium in our porous medium, i.e., our sedimentary basin. We also give the corresponding boundary and initial conditions.

### 2.1 Model equations

We write $\phi$ the Eulerian porosity of the medium and $V_\alpha$ the velocity field associated with phase $\alpha \in \{s, w\}$ (solid or

water). We write $\rho_\alpha$ the density of phase $\alpha$ and $\bar{\rho}_s = \rho_s(1 - \phi)$ and $\bar{\rho}_w = \rho_w(\phi)$ the effective densities. The medium is assumed to be fully saturated with water.

By convention, we orient the orthonormal basis $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ anticlockwise so that $z$ points upwards and gravity reads $\mathbf{g} = (0, 0, -g)$ with $g > 0$. The origin of the corresponding axis triplet $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ is situated at the present sea level so that $z = 0$ at the sea surface, and the basin can have either negative or positive $z$-coordinates.

## 2 Standard porous medium model

We briefly discuss the standard equations that we use to model the flow and mechanical equilibrium in our porous medium, i.e., our sedimentary basin. We also give the corresponding boundary and initial conditions.

### 2.1.1 Model equations

We write $\phi$ the Eulerian porosity of the medium and $V_\alpha$ the velocity field associated with phase $\alpha \in \{s, w\}$ (solid or water). We write $\rho_\alpha$ the density of phase $\alpha$ and $\bar{\rho}_s = \rho_s(1 - \phi)$ and $\bar{\rho}_w = \rho_w(\phi)$ the effective densities. The medium is assumed to be fully saturated with water.

By convention, we orient the orthonormal basis $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ anticlockwise so that $z$ points upwards and gravity reads $\mathbf{g} = (0, 0, -g)$ with $g > 0$. The origin of the corresponding axis triplet $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ is situated at the present sea level so that $z = 0$ at the sea surface, and the basin can have either negative or positive $z$-coordinates.

#### 2.1.1.1 Flow

We give here the equations on water and solid flow.

**Conservation of mass.** For any phase $\alpha \in \{s, w\}$ we have conservation of mass:

$$ \frac{\partial \bar{\rho}_{\alpha}}{\partial t} + \operatorname{div}(\bar{\rho}_{\alpha} \mathbf{V}_{\alpha}) = q_{\alpha}, \qquad (1) $$

where $q_\alpha$ is the rate of deposit of phase $\alpha$ due to sedimentation at the top of the basin and $V_\alpha$ is the velocity of phase $\alpha$.

**Darcy's law.** The filtration, or mean percolation, velocity $U_w := \phi(V_w - V_s)$ is assumed to be small enough to follow Darcy's law, that is,

$$ U_w = \eta_w \mathbf{K} (-\nabla p_w + \bar{\rho}_w \mathbf{g}), \qquad (2) $$

where $\eta_w$ is the water mobility, $\mathbf{K}$ the permeability tensor and $p_w$ the water, or pore, pressure.

**Densities.** The solid density is assumed to be constant and the water density to follow a law which is linear in pressure $p_w$ and temperature $T$:

$$ \rho_s = \rho_{s,0} \quad \text{and} \quad \rho_w = \rho_{w,0}(1 + \alpha_w(T - T_0) + \beta_w(p_w - p_0)), \qquad (3) $$

where $T_0$ and $p_0$ are reference values for the temperature and water pressure, $\rho_{s,0}$ and $\rho_{w,0}$ are reference densities, $\alpha_w$ is the water thermal expansion and $\beta_w$ is the water
compressibility. These laws are justified as long as the solid phase is incompressible and thermally unexpansible and the water phase is weakly compressible and weakly thermally expansible.

**Water mobility.** Because of our single-phase setting, the water's mobility $\eta_w$ is simply the inverse of its viscosity $\mu_w$, which we assumed is temperature-dependent. Specifically, we choose the following law (cf., for example, Eq. (48) of Ref. [28]):

$$ \eta_w = \frac{1}{\mu_w} = a_w \left( T + \sqrt{b_w + T^2} - c_w \right), \quad (4) $$

where $a_w$, $b_w$ and $c_w$ are fitting parameters. If the basin is supposed to be isothermal, then the water mobility in (4) is constant and the water density in (3) only depends on temperature.

**Permeability.** The permeability tensor $K$ is given by

$$ K = K(\phi) A, \quad (5) $$

where $A$ is an anisotropy tensor and $K$, the intrinsic permeability, follows the Kozeny–Carman law [2–4]:

$$ K(\phi) = \begin{cases} \frac{k_1 \phi^{n_1}}{s_0^2 (1 - \phi)^{m_1}} & \text{if } \phi \le \phi_0, \\ \frac{k_2 \phi^{n_2}}{s_0^2 (1 - \phi)^{m_2}} & \text{if } \phi > \phi_0, \end{cases} $$

where $S_0$ is the specific surface area of the porous medium and $k_1, k_2, n_1, n_2, m_1, m_2$ and $\phi_0$ are constants related to lithology. To ensure the continuity of $K$ at $\phi_0$, we impose the following relation on the parameters:

$$ \frac{\phi_0^{n_2-n_1}}{(1-\phi_0)^{m_2-m_1}} = \frac{k_1}{k_2}. $$

The Kozeny–Carman law is applicable to laminar flows in sediments with spherical structures [28], which we may reasonably assume in our illustrative setting. The tensor $A$ in (5) is a diagonal matrix, where each component represents the weight of its direction in the permeability tensor; isotropy of the medium implies that all the components of $A$ are equal, which, without loss of generality, means that $A$ is the identity matrix.

## 2.1.2 Temperature and mechanical equilibrium

We now discuss the equations for temperature distribution and mechanical equilibrium in the basin.

**Temperature.** We suppose that the temperature of the basin obeys a very simple vertical model, namely,

$$ \frac{\partial T}{\partial z} = G, \quad (6) $$

where $G$ is a known temperature gradient, sometimes referred to as a geothermal gradient, which may depend on position and time. We refer the reader to [33] for an example of geomechanical model including the effects of large strains on the temperature, which we do not consider here.

**Mechanical equilibrium.** We denote by $\sigma$ the Cauchy stress tensor. This tensor includes both water and solid stresses and is therefore sometimes referred to as a total stress tensor. We write

$$ \boldsymbol{\sigma} = \begin{pmatrix} \sigma_x & \sigma_{xy} & \sigma_{xz} \\ \sigma_{xy} & \sigma_y & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_z \end{pmatrix}, \quad (7) $$

where $\sigma_i$ and $\sigma_{ij}$ are the normal and shear stresses for all $i, j \in \{x, y, z\}$ with $j \neq i$. We use the sign convention according to which normal stress is positive when there is compression in its direction.

By Cauchy's momentum equation, mechanical equilibrium is given by

$$ \mathrm{div} \, \boldsymbol{\sigma} = \bar{\rho} g + \boldsymbol{f}, \quad (8) $$

where $\bar{\rho} := \bar{\rho}_w + \bar{\rho}_s$ is the homogenized density and $\boldsymbol{f} := (f_x, f_y, f_z)$ is the vector containing the volumic external forces other than gravity. More explicitly, (8) rewrites as

$$ \begin{gathered} \frac{\partial \sigma_x}{\partial x} + \frac{\partial \sigma_{xy}}{\partial y} + \frac{\partial \sigma_{xz}}{\partial z} = f_x; \quad \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_y}{\partial y} + \frac{\partial \sigma_{yz}}{\partial z} = f_y; \\ \frac{\partial \sigma_z}{\partial z} = -\bar{\rho}g + f_z, \end{gathered} \quad (9) $$

where we make the assumption that

$$ \frac{\partial \sigma_{xz}}{\partial x} + \frac{\partial \sigma_{yz}}{\partial y} = 0. $$

This assumption allows us to compute $\sigma_z$ from (9) by simply imposing a Dirichlet condition on it at the basin's top (cf. Sect. 2.2.1).

**Remark 2.1** We assume that the shear stresses within the water phase and between the water and solid phases are negligible with respect to those within the skeleton, i.e., the effect of water viscosity is not considered in the computation of the stresses (although it is in the computation of the flow via (2) and (4)). Thus, in (7), the terms $\sigma_{yz}$, $\sigma_{xz}$, and $\sigma_{xy}$ are in fact the skeleton's shear stresses. This is not true of the normal stresses since the water pressure needs to be taken into account via Biot's diagonal corrective term (cf. Sect. 3.1.1 and the seminal references [34, 35]).

### 2.1.3 Porosity

At this stage, we do not have an expression for the porosity $\phi$, although it intervenes in many of the model equations, namely, (1), (2), (5), and (8). In a standard model for porous media, one could assume $\phi$ to be constant, which would close our model already. However, because of compaction, this is not satisfactory here. We derive a model for compaction, and thus $\phi$, in Section 3.4.

## 2.2 Boundary and initial conditions

We suppose that the sedimentary top is given by the graph of a function $S$, that is, for all $(x, y)$ the value $S(x, y)$ gives
the vertical position of the top right above point $(x, y, z)$, for any $z \le S(x, y)$. Recall that the z-axis is directed upwards.

### 2.2.1 Boundary conditions

Integrating the third equation in (9) and (6) tells us to impose a total vertical stress and a temperature at the sedimentary top. Denoting by $\Omega$ the basin and by $\partial\Omega_{top}$ its top, i.e., $\partial\Omega_{top} = \{(x, y, z) \in \Omega | z = s(x, y)\}$, we impose

$$ \left\{ \begin{array}{l} \sigma_z(x, y, z) = p_{sup} \\ T(x, y, z) = T_{sup} \end{array} \right. \quad \text{for all } (x, y, z) \in \partial\Omega_{top} \text{ at all times.} \qquad (10) $$

Here, $p_{sup}$ is the pressure stemming from the weight of seawater and atmosphere lying above the basin. On the rest of the boundary $\partial\Omega$, we impose zero-flux conditions, i.e.,

$$ V_w \cdot v = 0 \quad \text{on } \partial\Omega/\partial\Omega_{top} \text{ at all times,} $$

where $v$ is the outward unit normal of $\Omega$.

Regarding stresses, a note from (9) that these boundary conditions only close the problem for the vertical stress $\sigma_z$. Indeed, they do not give us the horizontal stresses $\sigma_x$ and $\sigma_y$, nor do they give us the shear stresses $\sigma_{yz}$, $\sigma_{xz}$ and $\sigma_{xy}$. Nevertheless, the horizontal stresses can be recovered from the vertical ones whenever the horizontal strains are known and a constitutive stress-strain law is imposed (cf. Sect. 3.4.1). As to the shear stresses, they are not of interest to us in this paper since we suppose that they do not impact porosity (cf. Sect. 3.2.1); we, therefore, do not worry about closing the problem for them.

### 2.2.2 Initial conditions

The history of the basin is split into a sequence of geological episodes, called *events*, during which either a new layer of sediments is deposited or an old layer is eroded. We gener- ically write $t_0$ the initial time of any *deposit* event.

Because each simulation starts at the beginning of a deposit event and layers are added to the top only, the initial conditions are already determined by the top bound- ary conditions given in (10):

$$ \left\{ \begin{array}{l} \sigma_z(t_0) = p_{sup}(t_0) \\ T(t_0) = T_{sup}(t_0) \end{array} \right. \quad \text{in new layer.} $$

## 3 Compaction

For the sake of presentation, we consider only one sedimentary layer starting from its initial time $t_0$ of deposition and eventually ending when fully eroded. Modeling the complete basin can then be achieved by applying the single-layer model in each layer.

As discussed in the introduction, we wish to integrate horizontal strains in a vertical porosity-stress law to account

for horizontal compaction when oedometric conditions are not verified. We consider that the horizontal strains are known time-dependent parameters. The resulting model is what we refer to as a simplified geomechanical model.

### 3.1 Relevant stress and strain tensors

We treat our basin as an *isotropic* (in particular, $\mathbf{A}$ is the identity matrix in (5)), *poroelastic* material, thus undergoing small deformations. Note that, in this setting, using the Lagrangian porosity, rather than the Eulerian one, as done here, would also be relevant since the pore pressure, stresses and Lagrangian porosity are thermodynamically linked (cf. Eq. (4.2) of Ref. [36]). Furthermore, the stress and strain quantities of relevance are the *effective stress tensor* and the *skeleton's infinitesimal strain tensor*.

#### 3.1.1 Effective stress

Biot's theory [34] for consolidation couples fluid flow with rock deformation and introduces the effective stress tensor $\boldsymbol{\sigma}'$, which is defined as

$$ \boldsymbol{\sigma}' = \boldsymbol{\sigma} - p_w \boldsymbol{b}, $$

where $\boldsymbol{b}$ is Biot's tensor. In [35], the authors link Biot's tensor to the compressibility of the medium according to

$$ \boldsymbol{b} = \left( 1 - \frac{K}{K_s} \right) \boldsymbol{I} =: b \boldsymbol{I}, $$

with $K_s$ and $K$ the moduli of compressibility of the solid phase and the skeleton, respectively, and $\boldsymbol{I}$ the identity matrix. Biot's coefficient $b$ satisfies $b \simeq 1$ for a weakly compressible solid phase and a highly compressible skeleton, in which case we recover Terzaghi's theory for soil deformation [37]. For explicit expressions of $b$ depending on porosity, we refer the reader to [33, 38–40].

The effective stress tensor reads

$$ \boldsymbol{\sigma}' = \begin{pmatrix} \sigma_x' & \sigma_{xy} & \sigma_{xz} \\ \sigma_{xy} & \sigma_y' & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_z' \end{pmatrix}, $$

where $\sigma_x' = \sigma_x - p_w b$ for all $i \in \{x, y, z\}$. The tensor $\boldsymbol{\sigma}'$ is in fact the skeleton's stress tensor (cf. Remark 2.1).

Note that, by equilibrium with the seawater and the atmosphere, we have $p_w = p_{sup}$, that is, $\sigma_z' = \sigma_z - p_w b = (1 - b)p_{sup}$ on the top boundary $\partial\Omega_{top}$; in particular, $\sigma_z' = 0$ on $\partial\Omega_{top}$ when $b = 1$. Equivalently, we have

$$ \sigma_z'(t_0) = (1 - b)p_{sup}(t_0) \qquad (11) $$

in particular, $\sigma_z'(t_0) = 0$ when $b = 1$.

#### 3.1.2 Skeleton's infinitesimal strain

We denote by $\varepsilon$ the skeleton's infinitesimal strain, or deformation, tensor. (In the following, we omit the term “infinitesimal” when referring to $\varepsilon$.) Because the basin is assumed to be elastic, $\varepsilon$ coincides in fact with the elastic strain tensor and satisfies $|\varepsilon| \ll 1$.
We use the following notation:

$$ \epsilon = \begin{pmatrix} \epsilon_x & \epsilon_{xy} & \epsilon_{xz} \\ \epsilon_{xy} & \epsilon_y & \epsilon_{yz} \\ \epsilon_{xz} & \epsilon_{yz} & \epsilon_z \end{pmatrix}, $$

where $\epsilon_i$ and $\epsilon_{ij}$ are the normal and shear strains for all $i, j \in \{x, y, z\}$ with $i \neq j$. In agreement with our stress convention, any positive normal strain corresponds to a compression.

## 3.2 General compaction law

Given any differentiable time-dependent function $f$, we write $\dot{f}$ its time derivative, which we also refer to as the rate of $f$.

### 3.2.1 Generic formula

Compaction is characterized by a change in porosity under a change in compression. Supposing that the porosity rate does not depend on shear stresses (cf. Remark 2.1), we may assume that the actual porosity $\phi$ can be approximated as $\phi \simeq \phi_F$, where

$$ \dot{\phi}_F = \mathbf{F}(\phi, \boldsymbol{\sigma}'_{\mathrm{n}}) \cdot \dot{\boldsymbol{\sigma}}'_{\mathrm{n}} \qquad (12) $$

and $\mathbf{F}$ is a continuous, vector-valued function (experimentally fitted), $\boldsymbol{\sigma}'_{\mathrm{n}} := (\sigma'_{x}, \sigma'_{y}, \sigma'_{z})$ is the diagonal vector of $\boldsymbol{\sigma}'$ and $\cdot$ stands for the Euclidean inner product. Following [40] and adapting the terminology therein to our nonplastic setting, we sometimes refer to $\mathbf{F}$ as the elastic function, which, in particular, needs to keep the porosity between 0 and 1. If we were to treat chemical compaction in addition to mechanical compaction, then we would need to add a term of the form $G(\phi, \boldsymbol{\sigma}'_{\mathrm{n}}) \cdot \boldsymbol{\sigma}'_{\mathrm{n}}$ to the right-hand side of (12) to account for viscous effects [24, 40, 41].

### 3.2.2 Porosity-free elastic function

For simplicity, we assume that $\mathbf{F}$ only depends on effective stress:

$$ \dot{\phi}_F = F(\boldsymbol{\sigma}'_{\mathrm{n}}) \cdot \dot{\boldsymbol{\sigma}}'_{\mathrm{n}}. \qquad (13) $$

We call (13) the compaction law, which is eventually nonlinear. We can retrieve an expression for $\phi_F$ by integrating it between the initial time $t_0$ and a time $t > t_0$:

$$ \phi_F(t) = \phi_F(t_0) + \int_{t_0}^{t} \mathbf{F}(\boldsymbol{\sigma}'_{\mathrm{n}}(\tau)) \cdot \dot{\boldsymbol{\sigma}}'_{\mathrm{n}}(\tau) d\tau. \qquad (14) $$

This completes our standard basin model given in Section 2 provided we find an appropriate elastic function $\mathbf{F}$. We first discuss the oedometric case and then generalize it to give our simplified geomechanical model.

## 3.3 Oedometric model

Assume in this section that we are within the oedometric hypothesis, that is, the strains satisfy

$$ \epsilon_x = \epsilon_y = \epsilon_{yz} = \epsilon_{xz} = \epsilon_{xy} = 0. \qquad (15) $$

### 3.3.1 Vertical compaction law

Horizontal effects on porosity being negligible in this case, we choose $\mathbf{F}$ in (13) to depend exclusively on the vertical stress $\boldsymbol{\sigma}'_z$:

$$ \mathbf{F}(\boldsymbol{\sigma}'_{\mathrm{n}}) = (0, 0, -\beta(\boldsymbol{\sigma}'_{z})), \qquad (16) $$

where $\beta$ is a function determined by lithology, so that (13) becomes

$$ \dot{\phi}_S = -\beta(\boldsymbol{\sigma}'_z)\boldsymbol{\sigma}'_z, \qquad (17) $$

where the notation $\phi_S$ replaces $\phi_F$. This law is empirical, and one-dimensional in that it only takes vertical stress into account.

### 3.3.2 Schneider's law

We choose Schneider's function for $\beta$ [24]:

$$ \beta(s) = \frac{\phi_1}{\sigma_1} e^{-s/\sigma_1} + \frac{\phi_2}{\sigma_2} e^{-s/\sigma_2}, \qquad (18) $$

where $\phi_1, \phi_2 \ge 0$ (porosities) and $\sigma_1, \sigma_2 > 0$ (stresses) are known parameters depending on lithology. The compaction law (17) is thus referred to as Schneider's law and $\phi_S$ as Schneider's porosity. Let $\alpha$ be a primitive of $-\beta$, that is,

$$ \alpha(s) = \phi_r + \phi_1 e^{-s/\sigma_1} + \phi_2 e^{-s/\sigma_2}, \qquad (19) $$

for some additional parameter $\phi_r$ called the residual porosity. At $t_0$, we impose

$$ \begin{aligned} \phi_S(t_0) &= \alpha(\boldsymbol{\sigma}'_z(t_0)) \\ &= \phi_r + \phi_1 e^{-(1-b)\phi_{\mathrm{sup}}(t_0)/\sigma_1} + \phi_2 e^{-(1-b)\phi_{\mathrm{sup}}(t_0)/\sigma_2}, \end{aligned} \qquad (20) $$

where we refer the reader to (11); this becomes $\phi_S(t_0) = \phi_r + \phi_1 + \phi_2$ if $b=1$. Then, (14) yields

$$ \phi_S = \alpha(\boldsymbol{\sigma}'_z). \qquad (21) $$

When positive, the residual porosity helps numerically avoid Schneider's porosity becoming negative as vertical effective stress increases.

Schneider's law is empirical and its parameters need to be found by experiment, which makes Schneider's porosity a fitted value as opposed to a direct measurement. The double exponential formulation allows for a better fit at both upper and lower sediment layers; indeed, for layers close to the surface (say, above 1000 m), where the effective stress is weak, the porosity generally shows steep variations, while this is not the case at lower depths, where the stress is stronger but the variations flatter. To accommodate both regimes, the double exponential, with its five parameters, proves to be a good choice (cf. [7, 40]); in this respect, $\phi_1$ and $\phi_2$ represent characteristic porosities at their corresponding stress regimes. This is in contrast with the single-exponential law proposed by Athy [42], which is recovered when $\phi_2 = \phi_r = 0$ in (18) and (19):

$$ \beta(s) = \frac{\phi_1}{\sigma_1} e^{-s/\sigma_1} \quad \text{and} \quad \alpha(s) = \phi_1 e^{-s/\sigma_1}. $$
We refer the reader to [43] for an equivalent reformulation of Athy's law on permeability rather than porosity. Also, it is a fact that Athy's law can be derived as a solution to a partial differential equation when compaction happens fast [44].

### 3.3.3 Erosion

Schneider's law is adequate to describe compaction in oedometric conditions. It is less so if erosion is involved, i.e., if decompaction, eventually followed by recompaction, occurs. During erosion, the law in (18) can be adjusted to include decompaction and recompaction as elastic phenomena [24, 40]. We do not discuss this issue further and simply assume that (18) still holds for erosion.

## 3.4 Simplified geomechanical model

We now derive the compaction model to handle non-oedometric conditions. As already mentioned, we suppose that the horizontal strains $ε_x$ and $ε_y$ (and thus their rates $\dot{ε}_x$ and $\dot{ε}_y$) are known time-dependent functions.

### 3.4.1 Stress-strain constitutive law

We start by discussing the relationship between effective stress and strain (in fact, between their rates).

**Generic formula.** We suppose that $\boldsymbol{\sigma}'$ and $\boldsymbol{\varepsilon}$ verify

$$ \dot{\boldsymbol{\sigma}}' = \boldsymbol{G}(\boldsymbol{\sigma}', \boldsymbol{\varepsilon}) : \dot{\boldsymbol{\varepsilon}}, $$

where $\boldsymbol{G}$ is a continuous function with fourth-order tensor values, called the *stiffness tensor*, and : stands for the tensor product between $3 \times 3 \times 3 \times 3$ and $3 \times 3$ tensors; we refer the reader to [45] for a practical way to compute such a product using Voigt's notation. This relation would need to be corrected according to finite strain theory using the spin tensor if we were to consider plastic strains [33, 39, 46]. Note that in classical infinitesimal elastic theory, $\boldsymbol{G}$ only depends on $\boldsymbol{\varepsilon}$ and so, by time integration, $\boldsymbol{\sigma}'$ is an explicit function of $\boldsymbol{\varepsilon}$; when furthermore the elasticity is linear, $\boldsymbol{G}$ is, in fact, independent of both $\boldsymbol{\sigma}'$ and $\boldsymbol{\varepsilon}$ and Hooke's law applies.

**Strain-free stiffness tensor.** Compaction impacts the form of the stiffness tensor $\boldsymbol{G}$. Indeed, as shown later (cf. (34)), assuming $\phi \simeq \phi_F$ with a compaction law of the form (13) requires $\boldsymbol{G}$ to depend on $\boldsymbol{\sigma}'$. In fact, in contrast with classical elasticity, we drop the dependence on $\boldsymbol{\varepsilon}$ and only keep that on $\boldsymbol{\sigma}'$:

$$ \dot{\boldsymbol{\sigma}}' = \boldsymbol{G}(\boldsymbol{\sigma}') : \dot{\boldsymbol{\varepsilon}}. \qquad (22) $$

Equivalently, writing $\boldsymbol{H}(\boldsymbol{\sigma}') = \boldsymbol{G}(\boldsymbol{\sigma}')^{-1}$, we have

$$ \dot{\boldsymbol{\varepsilon}} = \boldsymbol{H}(\boldsymbol{\sigma}') : \dot{\boldsymbol{\sigma}}'. \qquad (23) $$

We call (22) and (23) the *stress-strain constitutive law*. Note that the stress dependence in the stiffness tensor does not alter the mechanical equilibrium given in (8), since its resolution is carried out before the constitutive law is applied.

**Hooke-type law.** We choose Hooke's elastic law with stress-dependent Young modulus for the stiffness tensor, that is, we rewrite (22) as

$$ \dot{\boldsymbol{\sigma}}' = \frac{\boldsymbol{E}(\boldsymbol{\sigma}')}{1 + \nu} \dot{\boldsymbol{\varepsilon}} + \frac{\nu \boldsymbol{E}(\boldsymbol{\sigma}')}{(1 + \nu)(1 - 2\nu)} \mathrm{tr}(\dot{\boldsymbol{\varepsilon}}) \boldsymbol{I}, \qquad (24) $$

where $\nu \in [0, 0.5)$ is Poisson's coefficient and $E > 0$ is Young's stress-dependent modulus. Equivalently, we rewrite (23) as

$$ \dot{\boldsymbol{\varepsilon}} = \frac{1 + \nu}{E(\boldsymbol{\sigma}')^2} \dot{\boldsymbol{\sigma}}' - \frac{\nu}{E(\boldsymbol{\sigma}')^2} \mathrm{tr}(\dot{\boldsymbol{\sigma}}') \boldsymbol{I}, \qquad (25) $$

For simplicity, we refer to both (24) and (25) as Hooke's law.

As a direct consequence of (25), we have

$$ \mathrm{tr}(\dot{\boldsymbol{\varepsilon}}) = \frac{1 - 2\nu}{E(\boldsymbol{\sigma}')^2} \mathrm{tr}(\dot{\boldsymbol{\sigma}}'). \qquad (26) $$

The trace of $\dot{\boldsymbol{\varepsilon}}$ is the dilatation or relative volume change of the skeleton. Also, from the first two equations in (25), we establish

$$ \dot{\boldsymbol{\sigma}}'_{x} = \frac{\nu}{1 - \nu} \dot{\boldsymbol{\sigma}}'_{z} + \frac{\boldsymbol{E}(\boldsymbol{\sigma}')}{1 - \nu^2} (\dot{\boldsymbol{\varepsilon}}_{x} + \nu \dot{\boldsymbol{\varepsilon}}_{y}), \qquad (27) $$

$$ \dot{\boldsymbol{\sigma}}'_{y} = \frac{\nu}{1 - \nu} \dot{\boldsymbol{\sigma}}'_{z} + \frac{\boldsymbol{E}(\boldsymbol{\sigma}')}{1 - \nu^2} (\dot{\boldsymbol{\varepsilon}}_{y} + \nu \dot{\boldsymbol{\varepsilon}}_{x}), $$

and so, by summing these two equations and adding $\dot{\boldsymbol{\sigma}}'_{z}$, we get

$$ \mathrm{tr}(\dot{\boldsymbol{\sigma}}') = \frac{1 + \nu}{1 - \nu} \dot{\boldsymbol{\sigma}}'_{z} + \frac{\boldsymbol{E}(\boldsymbol{\sigma}')}{1 - \nu} (\dot{\boldsymbol{\varepsilon}}_{x} + \dot{\boldsymbol{\varepsilon}}_{y}). \qquad (28) $$

Initially, recalling (11), we impose

$$ \sigma'_{x}(t_{0}) = \sigma'_{y}(t_{0}) = \frac{\nu}{1 - \nu} \sigma'_{z}(t_{0}) = \frac{\nu}{1 - \nu} (1 - b) p_{\mathrm{sup}}(t_{0}), \qquad (29) $$

which becomes $\sigma'_{x}(t_{0}) = \sigma'_{y}(t_{0}) = 0$ when $b = 1$.

**Oedometric case.** We briefly return to the oedometric case (cf. (15)). *Stress and strain tensors.* Both equations in (27) and (28) imply

$$ \dot{\boldsymbol{\sigma}}'_{x} = \dot{\boldsymbol{\sigma}}'_{y} = \frac{\nu}{1 - \nu} \boldsymbol{\sigma}'_{z} \quad \text{and} \quad \mathrm{tr}(\dot{\boldsymbol{\sigma}}') = \frac{1 + \nu}{1 + \nu} \boldsymbol{\sigma}'_{z}. \qquad (30) $$

Thus, (26) gives

$$ \dot{\boldsymbol{\varepsilon}}_{z} = \mathrm{tr}(\dot{\boldsymbol{\varepsilon}}) = \frac{(1 + \nu)(1 - 2\nu)}{(1 - \nu)E(\boldsymbol{\sigma}')^2} \dot{\boldsymbol{\sigma}}'_{z}. \qquad (31) $$

Then, the oedometric rates of effective stress and strain tensors read as

$$ \dot{\boldsymbol{\sigma}}' = \begin{pmatrix} \frac{\nu}{1 - \nu} & 0 & 0 \\ 0 & \frac{\nu}{1 - \nu} & 0 \\ 0 & 0 & 1 \end{pmatrix} \dot{\boldsymbol{\sigma}}'_{z} \quad \text{and} \quad \dot{\boldsymbol{\varepsilon}} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \frac{(1 + \nu)(1 - 2\nu)}{(1 - \nu)E(\boldsymbol{\sigma}')^2} \end{pmatrix} \dot{\boldsymbol{\sigma}}'_{z}. \qquad (32) $$
**Young's modulus.** To find an expression for Young's modulus, we use Schneider's law. Let us first note that, by solid incompressibility (cf. [1] for example),

$$ \dot{\phi} = -(1 - \phi) \text{tr}(\dot{\epsilon}) = -(1 - \phi) \dot{\epsilon}_z. \quad (33) $$

Then, using (17), (21), and (31), and recalling the assumption that the porosity is well approximated by Schneider's porosity, we get

$$ E(\boldsymbol{\sigma}') = \frac{(1 + v)(1 - 2v)}{1 - v} \frac{1 - \alpha(\sigma'_z)}{\beta(\sigma'_z)}. \quad (34) $$

Thus, we see a posteriori, given the form of Schneider's function (18) and its primitive (19), that $E$ must indeed depend on effective stress (specifically, on the vertical effective stress $\sigma'_z$), even under oedometric conditions.

Note also that the choice of taking Young's modulus as stress-dependent over Poisson's coefficient in (24) and (25) is arbitrary and is only motivated by the fact (34) is a simpler expression than its equivalent for $v$.

### 3.4.2 Approximated problem

To determine the porosity, we introduce an *approximated* problem for the strains for which the oedometric hypothesis (15) holds and therefore Schneider's vertical law (17) and (18) is applicable.

More precisely, we take a stress-strain pair $(\boldsymbol{\sigma}^*, \varepsilon^*)$ whose components verify

$$ \varepsilon_x^* = \varepsilon_y^* = \varepsilon_{yz}^* = \varepsilon_{xz}^* = \varepsilon_{xy}^* = 0 \quad (35) $$

and

$$ \sigma_i^*(t_0) = \sigma_i'(t_0) \quad \text{for all } i \in \{x, y, z\}, \quad (36) $$

where we refer the reader to (11) and (29).

**Approximated porosity.** Because the oedometric hypothesis (35) is satisfied by the approximated strain tensor, Section 3.3 motivates the use of the following compaction law to define an approximated porosity $\phi^*$:

$$ \dot{\phi}^* = -\beta(\sigma_z^*)\dot{\sigma}_z^*, \quad (37) $$

where $\beta$ is given by Schneider's function (18). We further require that

$$ \phi^*(t_0) = \phi(t_0), \quad (38) $$

with $\phi(t_0)$ given in (20).

**Approximated tensors.** We let $\boldsymbol{\sigma}^*$ and $\varepsilon^*$ satisfy the constitutive relation

$$ \dot{\varepsilon}^* = \mathbf{H}(\boldsymbol{\sigma}^*) : \dot{\boldsymbol{\sigma}}^*, $$

where we recall that $\mathbf{H}$ is given by (25). We directly get the analog of (26):

$$ \text{tr}(\dot{\varepsilon}^*) = \frac{1 - 2v}{E(\boldsymbol{\sigma}^*)} \text{tr}(\dot{\boldsymbol{\sigma}}^*). $$

Also, following the computations in Section 3.4.1, we yield

$$ \dot{\sigma}_x^* = \dot{\sigma}_y^* = \frac{v}{1 - v} \dot{\sigma}_z^*, \quad (39) $$

and

$$ \text{tr}(\dot{\boldsymbol{\sigma}}^*) = \frac{1 + v}{1 - v} \dot{\sigma}_z^*, \quad (40) $$

in agreement with the oedometric forms given in (30).

**Approximating assumption.** The assumption we wish to make is that the approximated porosity $\phi^*$ is close to the actual porosity $\phi$ (or at least the empirical one $\phi_F$ given in (13)), that is,

$$ |\phi^* - \phi| \ll 1, \quad (41) $$

so that $\phi$ can be replaced by $\phi^*$ wherever needed. Thanks to (38), this is equivalent to

$$ \delta := |\phi^* - \dot{\phi}| \ll 1. \quad (42) $$

For this approximation to be justified given (37), we still need to choose an appropriate vertical approximated stress $\sigma_z^*$.

**Vertical approximated stress.** To account for the horizontal effects of compaction, $\sigma_z^*$ must include the horizontal deformations $\varepsilon_x$ and $\varepsilon_y$. Moreover, for consistency, it should coincide with $\sigma'_z$ in oedometric conditions.

**Convex combination.** We ask for a convex combination of rates of horizontal stresses to be preserved, that is, for some $\lambda \in [0, 1]$ we let

$$ (1 - \lambda)\dot{\sigma}_x^* + \lambda\dot{\sigma}_y^* = (1 - \lambda)\dot{\sigma}_x' + \lambda\dot{\sigma}_y', \quad (43) $$

Given (27) and (39), this yields

$$ \frac{v}{1 - v} \dot{\sigma}_z^* = \frac{v}{1 - v} \dot{\sigma}_z' + \frac{E(\boldsymbol{\sigma}')}{1 - v^2} \left( (1 - \lambda) + v\lambda \right) \dot{\varepsilon}_x + (\lambda + v(1 - \lambda)) \dot{\varepsilon}_y, $$

where we recall that $E(\boldsymbol{\sigma}')$ is Young's stress-dependent modulus. Then,

$$ \dot{\sigma}_z^* = \dot{\sigma}_z' + \frac{E(\boldsymbol{\sigma}')}{v(1+v)} (\dot{\varepsilon}_x + \dot{\varepsilon}_y - (1-v)(\dot{\lambda}\dot{\varepsilon}_x + (1-\lambda)\dot{\varepsilon}_y)), \quad (44) $$

note that this formulation only holds if $v \neq 0$.

The choice for the parameter $\lambda$ is arbitrary at this point. We impose that only the rate of $x$-stress should be preserved (i.e., $\lambda = 0$) if only the rate of $x$-strain is nonzero and vice-versa (i.e., $\lambda = 1$) if only the rate of $y$-strain is nonzero. Therefore, we choose

$$ \lambda = \frac{|\dot{\varepsilon}_y|}{|\dot{\varepsilon}_x| + |\dot{\varepsilon}_y|}, $$

and (44) becomes

$$ \dot{\sigma}_z^* = \dot{\sigma}_z' + \frac{E(\boldsymbol{\sigma}')}{v(1+v)} \begin{cases} \left( \dot{\varepsilon}_x + \dot{\varepsilon}_y - \frac{2(1-v)\dot{\varepsilon}_x \dot{\varepsilon}_y}{\dot{\varepsilon}_x + \dot{\varepsilon}_y} \right) & \text{if } \dot{\varepsilon}_x \dot{\varepsilon}_y > 0, \\ (\dot{\varepsilon}_x + \dot{\varepsilon}_y) & \text{if } \dot{\varepsilon}_x \dot{\varepsilon}_y \le 0. \end{cases} \quad (45) $$
**Remark 3.1** In oedometric conditions, we note that $\dot{\sigma}_z^* = \dot{\sigma}_z'$ and, thanks to (36), that $\sigma_z^* = \sigma_z'$; furthermore, (17) and (37) give $\phi^* = \phi$.

**Other possible choices.** There are various viable conditions to impose other than (43). For example, we could preserve the stress trace, i.e., the first invariant:

$$ \mathrm{tr}(\dot{\boldsymbol{\sigma}}^*) = \mathrm{tr}(\dot{\boldsymbol{\sigma}}'), \quad (46) $$

which, thanks to (28) and (40), would lead to

$$ \dot{\sigma}_z^* = \dot{\sigma}_z' = \frac{E(\boldsymbol{\sigma}')}{1 + v} (\dot{\epsilon}_x + \dot{\epsilon}_y). $$

This relation, as opposed to our choice (45), disregards the relative magnitude of the horizontal-strain rates, but makes sense when $v = 0$. Interestingly, thanks to (40) and (46), we have $E(\boldsymbol{\sigma}') \mathrm{tr}(\varepsilon) = E(\boldsymbol{\sigma}^*) \mathrm{tr}(\varepsilon^*)$, where $\mathrm{tr}(\varepsilon)$ and $\mathrm{tr}(\varepsilon^*)$ are the relative volume changes for the actual and approximated problems, respectively. Note that this trace-preserving approach, coupled with the approximated compaction law (37), yields a model equivalent to replacing $\boldsymbol{\sigma}^*$ by $\mathrm{tr}(\boldsymbol{\sigma}')$ in (37) and appropriately rescaling the parameters of Schneider's law by a factor depending on Poisson's coefficient. This approach should be explored in the future since the trace invariance between the original and the approximated problem is a physically desirable property.

In view of the above discussion, we note that a general expression for the rate of vertical approximated stress is the following additive adjustment of the rate of vertical effective stress:

$$ \dot{\boldsymbol{\sigma}}_z^* = \dot{\boldsymbol{\sigma}}_z' + \zeta(\dot{\epsilon}_x, \dot{\epsilon}_y) E(\boldsymbol{\sigma}'), \quad (47) $$

for some strain-dependent, continuous function $\zeta$ such that $\zeta(0, 0) = 0$. For instance, in our case (cf. (45)), we have

$$ \zeta(\dot{\epsilon}_x, \dot{\epsilon}_y) = \frac{1}{v(1+v)} \begin{cases} \left( \dot{\epsilon}_x + \dot{\epsilon}_y - \frac{2(1-v)\dot{\epsilon}_x \dot{\epsilon}_y}{\dot{\epsilon}_x + \dot{\epsilon}_y} \right), & \text{if } \dot{\epsilon}_x \dot{\epsilon}_y > 0, \\ (\dot{\epsilon}_x + \dot{\epsilon}_y), & \text{if } \dot{\epsilon}_x \dot{\epsilon}_y \le 0. \end{cases} \quad (48) $$

One could also integrate geometrical constraints related to the basin's shape into the function $\zeta$. We leave a detailed study of the many possible choices of $\zeta$ to for future work.

**Porosity.** We split the elastic function $\mathbf{F}$ in the compaction law (13) according to

$$ \mathbf{F} = \mathbf{F}_o + \mathbf{F}_{no}, \quad (49) $$

where $\mathbf{F}_o$ follows the oedometric formulation given in (16), that is,

$$ \mathbf{F}_o(\boldsymbol{\sigma}_n') = (0, 0, -\beta(\boldsymbol{\sigma}_z')). \quad (50) $$

Recalling $\phi \approx \phi_F$, this leads to a decomposition of the porosity as

$$ \phi = \phi_o + \phi_{no}, \quad (51) $$

where $\phi_o$ and $\phi_{no}$ satisfy $\dot{\phi}_o = \mathbf{F}_o(\boldsymbol{\sigma}_n') \cdot \dot{\boldsymbol{\sigma}}_n'$ and $\phi_{no} = \mathbf{F}_{no}(\boldsymbol{\sigma}_n') \cdot \dot{\boldsymbol{\sigma}}_n'$ and are the oedometric and non-oedometric contributions to the porosity change $\dot{\phi}$; in particular, there holds

$$ \dot{\phi}_o = -\beta(\boldsymbol{\sigma}_z') \boldsymbol{\sigma}_z', $$

where $\beta$ is Schneider's function (18).

**Young's modulus.** To close our approximated problem, and therefore our simplified geomechanical model given by (37)–(41)–(45), we need to determine an expression for Young's stress-dependent modulus. To this end, we split oedometric and non-oedometric contributions in porosity, effective stress, and strain.

**Effective stress and strain.** Similarly as for the porosity, we split the effective stress and strain tensors as

$$ \boldsymbol{\sigma}' = \boldsymbol{\sigma}_o' + \boldsymbol{\sigma}_{no}' \quad \text{and} \quad \varepsilon = \varepsilon_o + \varepsilon_{no}, $$

where, according to the stress-strain constitutive law (23), we have

$$ \dot{\varepsilon}_o = \mathbf{H}(\boldsymbol{\sigma}') : \dot{\boldsymbol{\sigma}}_o' \quad \text{and} \quad \dot{\varepsilon}_{no} = \mathbf{H}(\boldsymbol{\sigma}') : \dot{\boldsymbol{\sigma}}_{no}' \quad (52) $$

with $\mathbf{H}$ given in (25). From (27) and (32), we choose

$$ \dot{\boldsymbol{\sigma}}_o' = \begin{pmatrix} \frac{v}{1-v} & 0 & 0 \\ 0 & \frac{v}{1-v} & 0 \\ 0 & 0 & 1 \end{pmatrix} \dot{\boldsymbol{\sigma}}_z'. \quad (53) $$

From (52) and (53), we get

$$ \mathrm{tr}(\dot{\varepsilon}_o) = \frac{(1 + v)(1 - 2v)}{(1 - v)E(\boldsymbol{\sigma}')}, \quad (54) $$

**Stress-dependent modulus.** From (33), there holds

$$ \dot{\phi} = -(1 - \phi)\mathrm{tr}(\dot{\epsilon}) = -(1 - \phi)\mathrm{tr}(\dot{\varepsilon}_o) - (1 - \phi)\mathrm{tr}(\dot{\varepsilon}_{no}). $$

By identification with (51), we, therefore, find that $\phi_o$ and $\mathrm{tr}(\varepsilon_o)$ satisfy

$$ \dot{\phi}_o = -(1 - \phi)\mathrm{tr}(\dot{\varepsilon}_o). \quad (55) $$

Then, (54) and (55) lead to

$$ E(\boldsymbol{\sigma}') = \frac{(1 + v)(1 - 2v)}{1 - v} \frac{1 - \phi}{\beta(\boldsymbol{\sigma}_z')}. \quad (56) $$

Our simplified geomechanical model for compaction is therefore finally given by (37)–(41)–(45)–(56), summarized by

$$ \left\{ \begin{array}{l} \phi = \phi^*, \\ \dot{\phi}^* = -\beta(\boldsymbol{\sigma}_z^*) \dot{\boldsymbol{\sigma}}_z^*, \\ \dot{\boldsymbol{\sigma}}_z^* = \dot{\boldsymbol{\sigma}}_z' + \frac{E(\boldsymbol{\sigma}')}{v(1+v)} \begin{cases} \left( \dot{\epsilon}_x + \dot{\epsilon}_y - \frac{2(1-v)\dot{\epsilon}_x \dot{\epsilon}_y}{\dot{\epsilon}_x + \dot{\epsilon}_y} \right) & \text{if } \dot{\epsilon}_x \dot{\epsilon}_y > 0, \\ (\dot{\epsilon}_x + \dot{\epsilon}_y) & \text{if } \dot{\epsilon}_x \dot{\epsilon}_y \le 0, \end{cases} \\ E(\boldsymbol{\sigma}') = \frac{(1 + v)(1 - 2v)}{1 - v} \frac{1 - \phi}{\beta(\boldsymbol{\sigma}_z)}. \end{array} \right. \quad (57) $$
**Validity of the model.** Let us come back to the question of the validity of the assumption that $\phi^* = \phi$ (cf. (41) and (42)). When there are no horizontal strains, we already know from Remark 3.1 that our model is such that $\dot{\phi}^* = \dot{\phi}$ and so $\delta = 0$ (cf. (42)). However, we would like to quantify how small $\dot{\epsilon}_x$ and $\dot{\epsilon}_y$ need to be for $\delta$ to be reasonably small and thus our model to be valid. Answering this question requires the use of advanced analytical tools, which is out of the scope of this paper. Instead, we validate our model numerically on a complex test case in Section 4.

Still, to give an idea of how to answer the question analytically, we give a formal computation in this direction. Using again $\phi \simeq \phi_F$ (cf. (13)), we have

$$ \delta = |\mathbf{F}(\boldsymbol{\sigma}'_n) \cdot \dot{\boldsymbol{\sigma}}'_n + \beta(\sigma_z^*) \dot{\sigma}_z^*| \le |(\mathbf{F}_{\text{no}}(\boldsymbol{\sigma}'_n) + (0,0,\beta(\sigma_z^*) - \beta(\sigma_z'))) \cdot \dot{\boldsymbol{\sigma}}'_n| + \beta(\sigma_z^*) |\dot{\sigma}_z^* - \dot{\sigma}_z'| \le |\mathbf{F}_{\text{no}}(\boldsymbol{\sigma}'_n) \cdot \dot{\boldsymbol{\sigma}}'_n| + |\beta(\sigma_z^*) - \beta(\sigma_z')||\dot{\sigma}_z'| + \beta(\sigma_z^*) E(\boldsymbol{\sigma}') |\zeta(\dot{\epsilon}_x, \dot{\epsilon}_y)| \le |\mathbf{F}_{\text{no}}(\boldsymbol{\sigma}'_n) \cdot \dot{\boldsymbol{\sigma}}'_n| + |\beta(\sigma_z^*) - \beta(\sigma_z')||\dot{\sigma}_z'| + \frac{(1+\nu)(1-2\nu)}{1-\nu} \frac{\beta(\sigma_z^*)}{\beta(\sigma_z')} |\zeta(\dot{\epsilon}_x, \dot{\epsilon}_y)|, $$

where the first line comes from (13) and (37), the second line from (49) and (50), the third line from (47), and the fourth line from (56), and the fact that $\phi < 1$. Recall that $\zeta$ is given by (48) in our case. We thus see that if we have a control on how close $\beta(\sigma_z^*)$ is to $\beta(\sigma_z')$, then the second and third terms in the right-hand side of the last inequality can be bounded by $|\dot{\sigma}_z'|, |\dot{\epsilon}_x|$ and $|\dot{\epsilon}_y|$. The first term, however, depends on the non-oedometric contribution to the porosity, on which we do not have any information. Hence, for the moment, this analysis is inconclusive.

# 4 Simulations

We now discuss the numerical results obtained with our simplified geomechanical model (57) for compaction coded entirely in ArcTem.

## 4.1 Neuquén basin

We test our model on the Vaca Muerta formation of the Neuquén basin in Argentina. Its geological history involves many episodes of tectonic and sedimentary deformation, which makes it an attractive site to test compaction models.

We subdivide the basin's history into 30 events starting at -200My and ending at the present day. These events include the deposition of sediments such as sandstone and carbonate over the first 20 events and their erosion at the top of the basin in the last 10 events. We use the same geometrical, lithological, and mechanical parameters as in [26] and validate our results by comparison with those obtained there using the code A². We refer the reader to [26] for a detailed account of the basin's history and sediments' properties; importantly, the basin shows high overpressures, and drained conditions with hydrostatic pressure in the upper layers, in contrast with the lower layers.

The main characteristics of the basin, summarized in Figure 1, are as follows:

* Spatial dimensions: 195 km East-West and 180 km North-South;

* Spatial discretization: 100 620 finite-volume cells;

* Sediments: 7 groups of different materials;

* East-West strain: 4% from -10My to -8My (2 events) and 2% from -8My to present day (8 events).

## 4.2 Code A²

As mentioned in the introduction (Sect. 1), and the references therein, code A² iteratively couples a finite-volume basin simulator (ArcTem) with a three dimensional finite-element mechanical simulator (Code Aster). More specifically, during predefined time events, both the basin and mechanical codes carry out their own, independent, time discretizations as well as porosity and pressure evaluations, the former using Schneider's law and the latter using Prager-Drucker plasticity model. At the end of each event, these porosities are compared: if, within a given tolerance, they are equal, then the simulation continues to the next event; if they are too far apart, porosity and pressure are corrected within the basin code using an additive term on the effective vertical stress and the time event is recomputed until convergence, i.e., until an agreement between the basin and mechanical codes. Code Aster has the advantage of being more precise than ArcTem because of its three-dimensional treatment of compaction; however, it is costlier and requires a fixed mesh. Furthermore, ArcTem can accommodate degenerate hexagonal cells, whereas Code Aster needs to reshape such cells and thus requires additional unknowns, which increase the computational cost of running the code A² against our simplified geomechanical model (cf. Sect. 4.4.3 below for CPU times).

## 4.3 Parameters

We give in Table 1 the parameters chosen for our simulations and involved in the model equations given through Sections 2 and 3. As already mentioned, these parameters coincide with those used in [26].

In addition, the horizontal strain rates, measured in s⁻¹, are of the form

$$ f(x) = \begin{cases} 0, & t < -10 \text{ My}, \\ \alpha_i, & -10 \text{ My} \le t < -8 \text{ My}, \\ \beta_i, & -8 \text{ My} \le t, \end{cases} $$

where, for $i \in \{x, y\}$, the parameters $\alpha_i$ and $\beta_i$ are selected from the list below:

(zero) $\alpha_x = \beta_x = \alpha_y = \beta_y = 0$.

(ref) $\alpha_x = 0.02$, $\beta_x = 0.0025$ and $\alpha_y = \beta_y = 0$.

(dble) $\alpha_x = 0.04$, $\beta_x = 0.005$ and $\alpha_y = \beta_y = 0$.

(xy) $\alpha_x = \alpha_y = 0.01$ and $\beta_x = \beta_y = 0.00125$.

(neg) $\alpha_x = -0.0025$, $\beta_x = -0.0003125$ and $\alpha_y = \beta_y = 0$.

Choice (zero) refers to oedometric conditions, where lateral deformations are neglected. Choice (ref) corresponds to that in [26] and to experimentally observed data (cf. Sect. 4.1), and thus gives us a reference point for our simulations. The other three choices, namely, (dble), (xy), and (neg), are variations of the reference case, where
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometrical and lithological model of the Neuquén basin (from [7]). (a) Representation of the mesh with a five-time vertical exaggeration (left); several longitudinal and transversal sections of the basin colored by sediment groups (right) and (b) lithology and burial evolution of each sedimentary group during the Mesozoic and Cenozoic eras at Well 10 (W10).</figcaption></figure>

we either double the horizontal strains, spread them evenly along the *x*- and *y*-directions, or apply negative ones, i.e., apply for an extension.

## 4.4 Results

We first compare the results of our simulations with those given in [26] using $A^2$, where we make use of (**zero**), (**ref**), and (**dble**). Then, we show additional results where the horizontal strains follow (**xy**) and (**neg**). All the results are given at the present day, i.e., at the end of the simulations. Finally, we give the CPU time for each simulation.

### 4.4.1 Comparison with [26]

Figure 2 compares our results with [26] ($A^2$) under the reference strains (**ref**) along Well 10 (W10); there, we display the horizontal total stresses, the water pressure, and

the porosity. We also show experimental data points for the pressure and the porosity under this same condition (**ref**) and give the respective results when no horizontal strains are imposed, i.e., when the oedometric condition (**zero**) holds. We see that our porosity stays very close in the lower layers to that obtained with $A^2$. In the upper layers, the porosities diverge slightly but both stay within the data cloud. The results are not as positive for the horizontal stresses and the pressure; indeed, although the trends are very similar to those with $A^2$, the values gradually separate as depth increases. Still, there is a significant improvement when switching from (**zero**) to (**ref**), in particular for the porosity, which shows that our simplified model does account for horizontal effects and is an improvement compared to an oedometric model, as desired, albeit underestimating the horizontal stresses and the pressure.

To get an indication of how our model behaves with deformation magnitude and how the underestimated
<table><caption>**Table 1.** Chosen parameters of the model equations.</caption><thead><tr><th></th><th>Parameter</th><th>Value (or interval of)</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="6">Densities</td><td>ρ<sub>s,0</sub></td><td>[2645, 2710]</td><td>kg·m<sup>-3</sup></td></tr><tr><td>ρ<sub>w,0</sub></td><td>1161.1</td><td>kg·m<sup>-3</sup></td></tr><tr><td>α<sub>w</sub></td><td>3.9 × 10<sup>-4</sup></td><td>°C<sup>-1</sup></td></tr><tr><td>β<sub>w</sub></td><td>4.5 × 10<sup>-4</sup></td><td>MPa<sup>-1</sup></td></tr><tr><td>T<sub>0</sub></td><td>15</td><td>°C</td></tr><tr><td>P<sub>0</sub></td><td>0.1</td><td>MPa</td></tr><tr><td rowspan="3">Water mobility</td><td>a<sub>w</sub></td><td>21.5</td><td>Pa·s·°C<sup>-1</sup></td></tr><tr><td>b<sub>w</sub></td><td>8078</td><td>°C<sup>2</sup></td></tr><tr><td>c<sub>w</sub></td><td>1200</td><td>°C</td></tr><tr><td rowspan="7">Permeability</td><td>S<sub>0</sub></td><td>[1.7 × 10<sup>12</sup>, 5 × 10<sup>7</sup>]</td><td>m<sup>2</sup></td></tr><tr><td>φ<sub>0</sub></td><td>0.1</td><td>%</td></tr><tr><td>k<sub>1</sub></td><td>20</td><td>—</td></tr><tr><td>k<sub>2</sub></td><td>0.2</td><td>—</td></tr><tr><td>n<sub>1</sub></td><td>5</td><td>—</td></tr><tr><td>n<sub>2</sub></td><td>3</td><td>—</td></tr><tr><td>m<sub>1</sub>, m<sub>2</sub></td><td>2</td><td>—</td></tr><tr><td>Geothermal gradient</td><td>G</td><td>[0.046, 0.068]</td><td>°C·m<sup>-1</sup></td></tr><tr><td>External force</td><td>f</td><td>(0, 0, 0)</td><td>kg·m<sup>-2</sup>·s<sup>-2</sup></td></tr><tr><td rowspan="2">Boundary conditions</td><td>p<sub>sup</sub></td><td>0.1</td><td>MPa</td></tr><tr><td>T<sub>sup</sub></td><td>[18, 26]</td><td>°C</td></tr><tr><td>Biot's coefficient</td><td>b</td><td>1</td><td>—</td></tr><tr><td rowspan="5">Schneider's law</td><td>φ<sub>r</sub></td><td>[0.01, 0.04]</td><td>%</td></tr><tr><td>φ<sub>1</sub></td><td>[0.29, 0.56]</td><td>%</td></tr><tr><td>φ<sub>2</sub></td><td>0</td><td>%</td></tr><tr><td>σ<sub>1</sub></td><td>[1.6, 40]</td><td>MPa</td></tr><tr><td>σ<sub>2</sub></td><td>10</td><td>MPa</td></tr><tr><td>Poisson's coefficient</td><td>ν</td><td>0.24</td><td>—</td></tr></tbody></table>

horizontal stresses and pressure can be impacted, we present in Figure 3 our results along W10 when the strain condition (**dble**) holds, i.e., when the horizontal strains are doubled with respect to the reference ones (**ref**). The results are still showing underestimated values of the horizontal stresses and the pressure, although less significantly and in deeper layers (compare Figs. 2 and 3). There is also a compelling improvement in the porosity, which is now closer to the $A^2$ result and more centered within the data points. The consistent misfit in lateral stresses may be due to the particularly high overpressure encountered in the Vaca Muerta formation, which may reach about 35 MPa (cf. Figs. 4 and 5 below). The way this overpressure dissipates laterally through the rock certainly affects lateral compaction and thus lateral stress variations [26]; this would need to be further investigated.

In Figure 4, we reproduce the profile obtained in [26] of the overpressure (i.e., the difference between the water pressure and the hydrostatic pressure) at the bottom of the basin and over a cross-section at the center of the basin. It shows both results when no horizontal deformations and when reference horizontal deformations are applied. Figure 5

presents our results for comparison with Figure 4. Again, in accordance with Figures 2 and 3, we see that we underestimate the overpressure under the strain condition (**ref**) and get much closer to the $A^2$ results when (**dble**) holds instead.

The noted underestimating behavior should not diminish the fact that the results given under (**ref**) and (**dble**) still show great improvement compared to the oedometric simulation (**zero**), in particular for the overpressure, pressure, and porosity in the undrained region of the basin, i.e., in the deep layers.

### 4.4.2 Additional results

To illustrate the flexibility of our simplified model, Figure 6 shows the horizontal stresses, the pressure, and the porosity along W10 when the horizontal strains are spread evenly in the *x*- and *y*-directions (cf. (xy)) and when negative horizontal strains are applied (cf. (neg)), i.e., when there is a horizontal extension rather than a compression. As expected, condition (xy) yields higher horizontal stresses and pressure and lower porosity than (**zero**), whereas
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison with $A^2$ under the strain condition (**ref**) and results under (**zero**) along W10 (data points and $A^2$ results from [26]).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Results under the strain condition (dble) along W10 and comparison with A² (data points and A² results from [26]).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** A² overpressure at bottom and over central cross-section under no horizontal strains (top) and reference horizontal strains (bottom) (from [26]).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Overpressure at bottom and over central cross-section under the strain conditions (**zero**), (**ref**), and (**dble**).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Results along W10 under the strain conditions (zero), (xy), and (neg).</figcaption></figure>
<table><caption>**Table 2.** CPU time for each simulation run.</caption><thead><tr><td></td><td>Time [min]</td></tr></thead><tbody><tr><td>(zero)</td><td>55.8</td></tr><tr><td>(ref)</td><td>80.8</td></tr><tr><td>(dble)</td><td>90.2</td></tr><tr><td>(xy)</td><td>77.2</td></tr><tr><td>(neg)</td><td>69.7</td></tr></tbody></table>

condition (neg) yields similar horizontal stresses and pressure and higher porosity. Note that, although the total horizontal strain in (xy) is the same as in (ref), the results are different (compare Figs. 2 and 6); this illustrates that our simplified model takes the direction of the strains into account when they are positive (cf. (45)).

When running simulations with (neg), we identified a limitation of our approach. Namely, when the horizontal extension is large, the pressure locally becomes very small while the porosity becomes very large, which prevents the simulation from converging. We were not able to impose any extension stronger than that given in (neg). This nonconvergence was not observed in any of the mentioned compression settings, where the deformations could actually be significantly increased, albeit with the drawback of higher numerical cost, which we discuss just below.

### 4.4.3 CPU times

Table 2 lists the CPU times required to run our simulations. We note that the difference in CPU time between a purely vertical model (i.e., an oedometric model, where no horizontal deformations are taken into account) and our simplified model goes from 24.9% to 61.6%, depending on the magnitude of the horizontal strains. This is a very attractive feature given that for a full-dimensional, finite element approach, as used in [26], the computational time loss in comparison to an oedometric model is much higher. The computational time, as communicated to us by the A² team from [26], for a full-dimensional run with the strain condition (ref) lies between 2.5 h and 6.5 h on 72 processors, depending on the chosen underlying linear solver; in contrast, the times given in Table 2 were obtained on a single processor on a personal laptop computer.

## 5 Conclusion and outlook

The model proposed in this paper to describe mechanical compaction is based on an alteration of Schneider's vertical porosity-stress law and an elastic stress-strain constitutive law involving a stress-dependent Young modulus. It coincides with classical modeling in oedometric conditions and provides a computation of porosity in non-oedometric conditions which is simpler, although not as accurate, compared to a full-dimensional, elastoplastic model. On the test case given by the Neuquén basin, our simplified model, coded in the finite-volume simulator ArcTem, shows reasonable accuracy in pressure, porosity, and vertical stress

along a drilling well, as well as in overpressure at the bottom and over a vertical cross-section at the center of the basin. Overall, our model tends to underestimate the horizontal stresses, the pressure, and the overpressure. Still, it shows significant improvement compared to an oedometric model, in particular in the lower, undrained compartment of the basin, and it produces very advantageous computational times. Therefore, it may offer a quick way to get pre-validating results before spending the time and computational resources required by a three-dimensional, finite-element simulator.

Although the model shows promising first results, it would still benefit from further numerical validation by comparing results on different wells, cross-sections, and test cases altogether. From the modeling and analytical point of view, other choices of approximated stresses need to be investigated, i.e., other functions $\zeta$ in (47) need to be tested; also, the analytical validation of the model requires additional work to show that the approximated porosity is indeed close to the actual one for a given range of horizontal stresses and strains (cf. the computation at the end of Sect. 3.4.2).

*Acknowledgments.* We thank the authors of [26] for sharing the A² data that made the validation of our model possible. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

## References

1. Schneider F., Burrus J., Wolf S. (1993) Modelling overpressures by effectivestress/porosity relationships in low-permeability rocks: empirical artifice or physical reality, *Norwegian Pet. Soc. Spec. Publ.* **3**, 333–341.

2. Carman P.C. (1937) Fluid flow through granular beds, *Chem. Eng. Res. Des.* **75**, S32–S48.

3. Carman P.C. (1956) *Flow of gases through porous media*, Academic Press, New York.

4. Kozeny J. (1927) Über kapillare Leitung des Wassers im Boden, *Sitzungsber Akad. Wiss.*, *Wien* **136**, 2a, 271–306.

5. Doyen P.M. (1988) Permeability, conductivity, and pore geometry of sandstone, *J. Geophys. Res.* **93**, B7, 7729.

6. Yang Y., Aplin A.C. (2010) A permeability-porosity relationship for mudstones, *Mar. Pet. Geol.* **27**, 8, 1692–1697.

7. Guy N., Colombo D., Frey J., Cornu T., Cacas-Stentz M.C. (2019) Coupled modeling of sedimentary basin and geomechanics: a modified Drucker–Prager cap model to describe rock compaction in tectonic context, *Rock Mech. Rock Eng.* **52**, 10, 3627–3643.

8. Beall A., Fisher A. (1969) *Sedimentology. Initial Reports Deep Sea Drilling Project*, vol. **1**, pp. 521–593.

9. Dunnington H.V. (1967) Aspects of diagenesis and shape change in stylolitic limestone reservoirs, in: *7th World Petroleum Congress (WPC)*.

10. Bloch S., McGowen J.H., Brizzolara D.W. (1990) Porosity prediction, prior to drilling, in sandstones of the Kekiktuk formation (Mississippian), North slope of Alaska, *AAPG Bull.* **74**, 9, 1371–1385.

11. Faille I., Thibaut M., Cacas M.C., Havé P., Willien F., Wolf S., Agelas L., Pegaz-Fiornet S. (2014) Modeling fluid flow in faulted basins, *Oil Gas Sci. Technol.* **69**, 4, 529–553.
12 Fuchtbauer H. (1967) Influence of different types of diagenesis on sandstone porosity, in: *7th World Petroleum Congress (WPC)*.

13 Schmidt V., McDonald D.A. (1979) *The role of secondary porosity in the course of sandstone diagenesis*, SEPM Society for Sedimentary Geology.

14 Dasgupta T., Mukherjee S. (2020) *Compaction of sediments and different compaction models*, Springer International Publishing, pp. 1–8.

15 Giles M.R., Indrelid S.L., James D.M.D. (1998) Compaction—the great unknown in basin modelling, *Geol. Soc. Spec. Publ.* **141**, 1, 15–43.

16 Marín-Moreno H., Minshull T.A., Edwards R.A. (2013) A disequilibrium compaction model constrained by seismic data and application to overpressure generation in the Eastern Black Sea basin, *Basin Res.* **25**, 3, 331–347.

17 Okiongbo K.S. (2011) Effective stress-porosity relationship above and within the oil window in the north sea basin, *Res. J. Appl. Sci. Eng. Technol.* **3**, 1, 32–38.

18 Doligez B. (1987) Migration of hydrocarbons in sedimentary basins, in: 2nd IFP Exploration Research Conference.

19 Schneider F., Wolf S., Faille I., Pot D. (2000) A 3d basin model for hydrocarbon potential evaluation: application to Congo offshore, *Oil Gas Sci. Technol.* **55**, 1, 3–13.

20 Ungerer P., Burrus J., Doligez B., Chenet P.Y., Bessis F. (1990) Basin evaluation by integrated two-dimensional modeling of heat transfer, fluid flow, hydrocarbon generation, and migration, *AAPG Bull.* **74**, 3,309–335.

21 Luo X., Vasseur G., Pouya A., Lamoureux-Var V., Poliakov A. (1998) Elastoplastic deformation of porous media applied to the modelling of compaction at basin scale, *Mar. Pet. Geol.* **15**, 2, 145–162.

22 Obradors-Prats J., Rouainia M., Aplin A.C., Crook A.J.L. (2017) Assessing the implications of tectonic compaction on pore pressure using a coupled geomechanical approach, *Mar. Pet. Geol.* **79**, 31–43.

23 Schneider F. (1993) Modèle de compaction élasto-plastique en simulation de bassins, *Oil Gas Sci. Technol.* **48**, 1, 3–14.

24 Schneider F., Potdevin J.L., Wolf S., Faille I. (1994) Modèle de compaction élastoplastique et viscoplastique pour simulation de bassins sédimentaires, *Oil Gas Sci. Technol.* **49**, 2, 1–8.

25 DiMaggio F.L., Sandler I.S. (1971) *The effect of strain rate on the constitutive equation of rocks, 2801T*, Defense Nuclear Agency.

26 Berthelon J., Brüch A., Colombo D., Frey J., Traby R., Bouziat A., Cacas-Stentz M.C., Cornu T. (2021) Impact of tectonic shortening on fluid overpressure in petroleum system modelling: insights from the Neuquén basin, *Argentina, Mar. Pet. Geol.* **127**, 104933.

27 Brüch A., Colombo D., Frey J., Berthelon J., Cacas-Stentz M.-C., Cornu T. (2020) Poro-elastoplastic constitutive law for coupling 3d geomechanics to classical petroleum system modeling, in: *54th U.S. Rock Mechanics/Geomechanics Symposium*.

28 Brüch A., Colombo D., Frey J., Berthelon J., Cacas-Stentz M.C., Cornu T., Gout C. (2021) Coupling 3d geomechanics

to classical sedimentary basin modeling: From gravitational compaction to tectonics, *Geomech. Energy Environ.* **28**, 100259.

29 EDF R&D (2017) *Code Aster*. http://www.code-aster.org.

30 Guy N., Enchéry G., Renard G. (2012) Numerical modeling of thermal EOR: Comprehensive coupling of an AMR-based model of thermal fluid flow and geomechanics, *Oil Gas Sci. Technol.* **67**, 6, 1019–1027.

31 Bouziat A., Guy N., Frey J., Colombo D., Colin P., Cacas-Stentz M.-C., Cornu T. (2019) An assessment of stress states in passive margin sediments: Iterative hydro-mechanical simulations on basin models and implications for rock failure predictions, *Geosci. J.* **9**, 11, 469.

32 Lyakhovsky V., Shalev E., Panteleev I., Mubassarova V. (2020) Directional compaction, *Earth Space Sci*. https://doi.org/10.1002/essoar.10504323.1.

33 Brüch A., Maghous S., Ribeiro F.L.B., Dormieux L. (2018) A thermo-poromechanical constitutive and numerical model for deformation in sedimentary basins, *J. Pet. Sci. Eng.* **160**, 313–326.

34 Biot M.A. (1941) General theory of three-dimensional consolidation, *J. Appl. Phys.* **12**, 2, 155–164.

35 Biot M.A., Willis D.G. (1957) The elastic coefficients of the theory of consolidation, *J. Appl. Mech.* **24**, 4, 594–601.

36 Coussy O. (2004) *Poromechanics*, Wiley.

37 Terzaghi K. (1923) Die Berechnung der Durchlassigkeitsziffer des Tones aus dem Verlauf der Hidrodynamischen Spannungserscheinungen, *Sitzungsberichte der Akademie der Wissenschaften in Wien* **132**, 125–138.

38 Boutéca M., Sarda J.-P., Laurent J. (1991) Rock positive (cf. (45)).mechanics contribution to the determination of fluid flow properties, in: *Second European Core Analysis Symposium (Eurocas II)*.

39 Maghous S., Brüch A., Bernaud D., Dormieux L., Braun A.L. (2014) Two-dimensional finite element analysis of gravitational and lateral-driven deformation in sedimentary basins, *Int. J. Numer. Anal. Methods Geomech.* **38**, 7, 725–746.

40 Schneider F., Potdevin J.L., Wolf S., Faille I. (1996) Mechanical and chemical compaction model for sedimentary basin simulators, *Tectonophysics* **263**, 1–4, 307–317.

41 Schneider F., Hay S. (2001) *Compaction model for quartzose sandstones – application to the Garn Formation, Haltenbanken, Mid-Norwegian continental shelf – part I – theory*, European Association of Geoscientists & Engineers.

42 Athy L.F. (1930) Density, porosity, and compaction of sedimentary rocks, *AAPG Bull.* **14**, 1, 1–24.

43 Gutierrez M., Wangen M. (2005) Modeling of compaction and overpressuring in sedimentary basins, *Mar. Pet. Geol.* **22**, 3, 351–363.

44 Fowler A.C., Yang X.-S. (1998) Fast and slow compaction in sedimentary basins, *SIAM J. Appl. Math* **59**, 1, 365–385.

45 Zienkiewicz O.C., Taylor R.L., Zhu J.Z. (2005) *The finite element method: its basis and fundamentals*, Elsevier Butterworth-Heinemann.

46 Dormieux L., Maghous S. (1999) Poroelasticity and poroplasticity at large strains, *Oil Gas Sci. Technol.* **54**, 6, 773–784.