Sci. Tech. Energ. Transition
Volume 78, 2023
Characterization and Modeling of the Subsurface in the Context of Ecological Transition
Article Number 22
Number of page(s) 17
Published online 29 August 2023

© The Author(s), published by EDP Sciences, 2023

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [24] 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 [1013]. 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 [1417] and simulators [1820]. 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, 2628] using the prototype code A2, 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 A2 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 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 km2 and spans a time interval of 10My. Our results under both oedometric and non-oedometric conditions are compared with those obtained using A2 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 ϕ the Eulerian porosity of the medium and Vα the velocity field associated with phase α ∈ {s, w} (solid or water). We write ρα the density of phase α and and the effective densities. The medium is assumed to be fully saturated with water.

By convention, we orient the orthonormal basis (x, y, z) anticlockwise so that z points upwards and gravity reads g = (0, 0, −g) with g > 0. The origin of the corresponding axis triplet (x, y, 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 Model equations

We write ϕ the Eulerian porosity of the medium and Vα the velocity field associated with phase α ∈ {s, w} (solid or water). We write ρα the density of phase α and and the effective densities. The medium is assumed to be fully saturated with water.

By convention, we orient the orthonormal basis (x, y, z) anticlockwise so that z points upwards and gravity reads g = (0, 0, −g) with g > 0. The origin of the corresponding axis triplet (x, y, 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 Flow

We give here the equations on water and solid flow.

Conservation of mass . For any phase α ∈ {s, w} we have conservation of mass:(1)where qα is the rate of deposit of phase α due to sedimentation at the top of the basin and Vα is the velocity of phase α.

Darcy’s law . The filtration, or mean percolation, velocity Uw := ϕ(Vw − Vs) is assumed to be small enough to follow Darcy’s law, that is,(2)where ηw is the water mobility, K the permeability tensor and pw 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 pw and temperature T:(3)where T0 and p0 are reference values for the temperature and water pressure, ρs,0 and ρw,0 are reference densities, αw is the water thermal expansion and β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 ηw is simply the inverse of its viscosity μw, which we assumed is temperature-dependent. Specifically, we choose the following law (cf., for example, Eq. (48) of Ref. [28]):(4)where aw, bw and cw 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(5)where A is an anisotropy tensor and K, the intrinsic permeability, follows the Kozeny–Carman law [24]:where S0 is the specific surface area of the porous medium and k1, k2, n1, n2, m1, m2 and ϕ0 are constants related to lithology. To ensure the continuity of K at ϕ0, we impose the following relation on the parameters: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,(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 σ 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(7)where σi and σij are the normal and shear stresses for all i, j ∈ {x, y, z} with j ≠ 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(8)where is the homogenized density and f := (fx, fy, fz) is the vector containing the volumic external forces other than gravity. More explicitly, (8) rewrites as(9)where we make the assumption thatThis assumption allows us to compute σ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 σyz, σxz, and σ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 ϕ, 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 ϕ 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 ϕ, 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 ≤ 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 Ω the basin and by top its top, i.e., ∂Ωtop = {(x, y, z) ∈ Ω|z = s(x, y)}, we impose(10)Here, psup is the pressure stemming from the weight of seawater and atmosphere lying above the basin. On the rest of the boundary Ω, we impose zero-flux conditions, i.e.,where ν is the outward unit normal of Ω.

Regarding stresses, a note from (9) that these boundary conditions only close the problem for the vertical stress σz. Indeed, they do not give us the horizontal stresses σx and σy, nor do they give us the shear stresses σyz, σxz, and σ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 generically write t0 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 boundary conditions given in (10):

3 Compaction

For the sake of presentation, we consider only one sedimentary layer starting from its initial time t0 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, 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 σ′, which is defined aswhere b is Biot’s tensor. In [35], the authors link Biot’s tensor to the compressibility of the medium according towith Ks and K the moduli of compressibility of the solid phase and the skeleton, respectively, and I the identity matrix. Biot’s coefficient b satisfies b ≃ 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, 3840].

The effective stress tensor readswhere for all i ∈ {x, y, z}. The tensor σ is in fact the skeleton’s stress tensor (cf. Remark 2.1).

Note that, by equilibrium with the seawater and the atmosphere, we have pw = psup, that is, on the top boundary top; in particular, on top when b = 1. Equivalently, we have(11)in particular, when b = 1.

3.1.2 Skeleton’s infinitesimal strain

We denote by ε the skeleton’s infinitesimal strain, or deformation, tensor. (In the following, we omit the term “infinitesimal” when referring to ε.) Because the basin is assumed to be elastic, ε coincides in fact with the elastic strain tensor and satisfies |ε| ≪ 1.

We use the following notation:where εi and εij are the normal and shear strains for all i, j ∈ {x, y, z} with i ≠ 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 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 ϕ can be approximated as ϕ ≃ ϕ F, where(12)and F is a continuous, vector-valued function (experimentally fitted), is the diagonal vector of σ and · stands for the Euclidean inner product. Following [40] and adapting the terminology therein to our nonplastic setting, we sometimes refer to 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 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 F only depends on effective stress:(13)We call (13) the compaction law, which is eventually nonlinear. We can retrieve an expression for ϕ F by integrating it between the initial time t0 and a time t > t0:(14)This completes our standard basin model given in Section 2 provided we find an appropriate elastic function 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(15)

3.3.1 Vertical compaction law

Horizontal effects on porosity being negligible in this case, we choose F in (13) to depend exclusively on the vertical stress :(16)where β is a function determined by lithology, so that (13) becomes(17)where the notation ϕS replaces ϕ 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 β [24]:(18)where ϕ1, ϕ2 ⩾ 0 (porosities) and σ1, σ2 > 0 (stresses) are known parameters depending on lithology. The compaction law (17) is thus referred to as Schneider’s law and ϕS as Schneider’s porosity. Let α be a primitive of −β, that is,(19)for some additional parameter ϕr called the residual porosity. At t0, we impose(20)where we refer the reader to (11); this becomes ϕS(t0) = ϕr + ϕ1 + ϕ2 if b = 1. Then, (14) yields(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, ϕ1 and ϕ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 ϕ2 = ϕr = 0 in (18) and (19):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 and ) 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 σ and ε verifywhere G is a continuous function with fourth-order tensor values, called the stiffness tensor, and : stands for the tensor product between 3 × 3 × 3 × 3 and 3 × 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, G only depends on ε and so, by time integration, σ′ is an explicit function of ε; when furthermore the elasticity is linear, G is, in fact, independent of both σ′ and ε and Hooke’s law applies.

Strain-free stiffness tensor. Compaction impacts the form of the stiffness tensor G. Indeed, as shown later (cf. (34)), assuming ϕ ≃ ϕ F with a compaction law of the form (13) requires G to depend on σ′. In fact, in contrast with classical elasticity, we drop the dependence on ε and only keep that on σ′:(22)Equivalently, writing H(σ′) = G(σ′)−1, we have(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(24)where ν ∈ [0, 0.5) is Poisson’s coefficient and E > 0 is Young’s stress-dependent modulus. Equivalently, we rewrite (23) as(25)For simplicity, we refer to both (24) and (25) as Hooke’s law.

As a direct consequence of (25), we have(26)The trace of is the dilatation or relative volume change of the skeleton. Also, from the first two equations in (25), we establish(27)and so, by summing these two equations and adding , we get(28)Initially, recalling (11), we impose(29)which becomes 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(30)Thus, (26) gives(31)Then, the oedometric rates of effective stress and strain tensors read as(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),(33)Then, using (17), (21), and (31), and recalling the assumption that the porosity is well approximated by Schneider’s porosity, we get(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 ), 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 ν.

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 (σ * , ε *) whose components verify(35)and(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 ϕ *:(37)where β is given by Schneider’s function (18). We further require that(38)with ϕ(t0) given in (20).

Approximated tensors. We let σ and ε satisfy the constitutive relationwhere we recall that H is given by (25). We directly get the analog of (26):Also, following the computations in Section 3.4.1, we yield(39)and(40)in agreement with the oedometric forms given in (30).

Approximating assumption. The assumption we wish to make is that the approximated porosity ϕ * is close to the actual porosity ϕ (or at least the empirical one ϕ F given in (13)), that is,(41)so that ϕ can be replaced by ϕ * wherever needed. Thanks to (38), this is equivalent to(42)For this approximation to be justified given (37), we still need to choose an appropriate vertical approximated stress .

Vertical approximated stress. To account for the horizontal effects of compaction, must include the horizontal deformations εx and εy. Moreover, for consistency, it should coincide with in oedometric conditions.

Convex combination. We ask for a convex combination of rates of horizontal stresses to be preserved, that is, for some λ ∈ [0, 1] we let(43)Given (27) and (39), this yieldswhere we recall that E(σ) is Young’s stress-dependent modulus. Then,(44)note that this formulation only holds if ν ≠ 0.

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

Remark 3.1

In oedometric conditions, we note that and, thanks to (36), that ; furthermore, (17) and (37) give ϕ* = ϕ.

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:(46)which, thanks to (28) and (40), would lead toThis relation, as opposed to our choice (45), disregards the relative magnitude of the horizontal-strain rates, but makes sense when ν = 0. Interestingly, thanks to (40) and (46), we have E(σ′) tr(ε) = E(σ * ) tr(ε * ), where tr(ε) and tr(ε * ) 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 σ * by tr(σ′) 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:(47)for some strain-dependent, continuous function ζ such that ζ (0, 0) = 0. For instance, in our case (cf. (45)), we have(48)One could also integrate geometrical constraints related to the basin’s shape into the function ζ. We leave a detailed study of the many possible choices of ζ to for future work.

Porosity. We split the elastic function F in the compaction law (13) according to(49)where Fo follows the oedometric formulation given in (16), that is,(50)Recalling ϕ ≃ ϕ F, this leads to a decomposition of the porosity as(51)where ϕo and ϕno satisfy and and are the oedometric and non-oedometric contributions to the porosity change ; in particular, there holdswhere β 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 aswhere, according to the stress-strain constitutive law (23), we have(52)with H given in (25). From (27) and (32), we choose(53)From (52) and (53), we get(54)

Stress-dependent modulus. From (33), there holdsBy identification with (51), we, therefore, find that ϕo and tr(εo) satisfy(55)Then, (54) and (55) lead to(56)

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

Validity of the model. Let us come back to the question of the validity of the assumption that ϕ*= ϕ (cf. (41) and (42)). When there are no horizontal strains, we already know from Remark 3.1 that our model is such that . and so δ = 0 (cf. (42)). However, we would like to quantify how small and need to be for δ 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 ϕ ≃ ϕ F (cf. (13)), we havewhere 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 ϕ < 1. Recall that ζ is given by (48) in our case. We thus see that if we have a control on how close is to , then the second and third terms in the right-hand side of the last inequality can be bounded by and . 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 A2. 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).

thumbnail 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).

4.2 Code A2

As mentioned in the introduction (Sect. 1), and the references therein, code A2 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 A2 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].

Table 1

Chosen parameters of the model equations.

In addition, the horizontal strain rates, measured in s−1, are of the formwhere, for i ∈ {x, y}, the parameters αi and βi are selected from the list below:

(zero) αx = βx = αy = βy = 0.

(ref) αx = 0.02, βx = 0.0025 and αy = βy = 0.

(dble) αx = 0.04, βx = 0.005 and αy = βy = 0.

(xy) αx = αy = 0.01 and βx = βy = 0.00125.

(neg) αx = − 0.0025, βx = −0.0003125 and αy = β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 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 A2, 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] (A2) 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 A2. 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 A2, 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.

thumbnail Fig. 2

Comparison with A2 under the strain condition (ref) and results under (zero) along W10 (data points and A2 results from [26]).

To get an indication of how our model behaves with deformation magnitude and how the underestimated 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 A2 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.

thumbnail Fig. 3

Results under the strain condition (dble) along W10 and comparison with A2 (data points and A2 results from [26]).

thumbnail Fig. 4

A2 overpressure at bottom and over central cross-section under no horizontal strains (top) and reference horizontal strains (bottom) (from [26]).

thumbnail Fig. 5

Overpressure at bottom and over central cross-section under the strain conditions (zero), (ref), and (dble).

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 A2 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 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)).

thumbnail Fig. 6

Results along W10 under the strain conditions (zero), (xy), and (neg).

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 A2 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.

Table 2

CPU time for each simulation run.

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 ζ 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).


We thank the authors of [26] for sharing the A2 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.


  • 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. [Google Scholar]
  • Carman P.C. (1937) Fluid flow through granular beds, Chem. Eng. Res. Des. 75, S32–S48. [Google Scholar]
  • Carman P.C. (1956) Flow of gases through porous media, Academic Press, New York. [Google Scholar]
  • Kozeny J. (1927) Über kapillare Leitung des Wassers im Boden, Sitzungsber Akad. Wiss., Wien 136, 2a, 271–306. [Google Scholar]
  • Doyen P.M. (1988) Permeability, conductivity, and pore geometry of sandstone, J. Geophys. Res. 93, B7, 7729. [CrossRef] [Google Scholar]
  • Yang Y., Aplin A.C. (2010) A permeability-porosity relationship for mudstones, Mar. Pet. Geol. 27, 8, 1692–1697. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Beall A., Fisher A. (1969) Sedimentology. Initial Reports Deep Sea Drilling Project, vol. 1, pp. 521–593. [Google Scholar]
  • Dunnington H.V. (1967) Aspects of diagenesis and shape change in stylolitic limestone reservoirs, in: 7th World Petroleum Congress (WPC). [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Fuchtbauer H. (1967) Influence of different types of diagenesis on sandstone porosity, in: 7th World Petroleum Congress (WPC). [Google Scholar]
  • Schmidt V., McDonald D.A. (1979) The role of secondary porosity in the course of sandstone diagenesis, SEPM Society for Sedimentary Geology. [CrossRef] [Google Scholar]
  • Dasgupta T., Mukherjee S. (2020) Compaction of sediments and different compaction models, Springer International Publishing, pp. 1–8. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • Doligez B. (1987) Migration of hydrocarbons in sedimentary basins, in: 2nd IFP Exploration Research Conference. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Schneider F. (1993) Modèle de compaction élasto-plastique en simulation de bassins, Oil Gas Sci. Technol. 48, 1, 3–14. [Google Scholar]
  • 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. [Google Scholar]
  • DiMaggio F.L., Sandler I.S. (1971) The effect of strain rate on the constitutive equation of rocks, 2801T, Defense Nuclear Agency. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • EDF R&D (2017) Code Aster. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Lyakhovsky V., Shalev E., Panteleev I., Mubassarova V. (2020) Directional compaction, Earth Space Sci. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Biot M.A. (1941) General theory of three-dimensional consolidation, J. Appl. Phys. 12, 2, 155–164. [Google Scholar]
  • Biot M.A., Willis D.G. (1957) The elastic coefficients of the theory of consolidation, J. Appl. Mech. 24, 4, 594–601. [CrossRef] [MathSciNet] [Google Scholar]
  • Coussy O. (2004) Poromechanics, Wiley. [Google Scholar]
  • Terzaghi K. (1923) Die Berechnung der Durchlassigkeitsziffer des Tones aus dem Verlauf der Hidrodynamichen Spannungserscheinungen, Sitzungsberichte der Akademie der Wissenschaften in Wien 132, 125–138. [Google Scholar]
  • Boutéca M., Sarda J.-P., Laurent J. (1991) Rock mechanics contribution to the determination of fluid flow properties, in: Second European Core Analysis Symposium (Eurocas II). [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • Athy L.F. (1930) Density, porosity, and compaction of sedimentary rocks, AAPG Bull. 14, 1, 1–24. [Google Scholar]
  • Gutierrez M., Wangen M. (2005) Modeling of compaction and overpressuring in sedimentary basins, Mar. Pet. Geol. 22, 3, 351–363. [CrossRef] [Google Scholar]
  • Fowler A.C., Yang X.-S. (1998) Fast and slow compaction in sedimentary basins, SIAM J. Appl. Math 59, 1, 365–385. [Google Scholar]
  • Zienkiewicz O.C., Taylor R.L., Zhu J.Z. (2005) The finite element method: its basis and fundamentals, Elsevier Butterworth-Heinemann. [Google Scholar]
  • Dormieux L., Maghous S. (1999) Poroelasticity and poroplasticity at large strains, Oil Gas Sci. Technol. 54, 6, 773–784. [Google Scholar]

All Tables

Table 1

Chosen parameters of the model equations.

Table 2

CPU time for each simulation run.

All Figures

thumbnail 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).

In the text
thumbnail Fig. 2

Comparison with A2 under the strain condition (ref) and results under (zero) along W10 (data points and A2 results from [26]).

In the text
thumbnail Fig. 3

Results under the strain condition (dble) along W10 and comparison with A2 (data points and A2 results from [26]).

In the text
thumbnail Fig. 4

A2 overpressure at bottom and over central cross-section under no horizontal strains (top) and reference horizontal strains (bottom) (from [26]).

In the text
thumbnail Fig. 5

Overpressure at bottom and over central cross-section under the strain conditions (zero), (ref), and (dble).

In the text
thumbnail Fig. 6

Results along W10 under the strain conditions (zero), (xy), and (neg).

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.