Issue 
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  
DOI  https://doi.org/10.2516/stet/2023019  
Published online  29 August 2023 
Review Article
A simplified vertical and horizontal geomechanical model for compaction in sedimentary basins
^{1}
IFP Energies nouvelles, 1 et 4 avenue de BoisPréau, 92852 RueilMalmaison, France
^{2}
CEA Cadarache, F13108 SaintPaulLezDurance, France
^{*} Corresponding author: francesco.patacchini@ifpen.fr
Received:
5
January
2023
Accepted:
18
July
2023
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 porositystress law where horizontal strains are inserted and on an elastic stressstrain law with stressdependent Young modulus. Though it is not threedimensional 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 fulldimensional finiteelement 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.
Key words: Sedimentary basin / Mechanical compaction / Horizontal deformation / Geomechanical model / Finitevolume simulator
© The Author(s), published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), 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 [2–4] or other sedimentspecific 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 onedimensional, 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 porositystress law [23, 24] is integrated into the Drucker–Prager model for plasticity [25], resulting in a threedimensional compaction model. This model is tested numerically in [7, 26–28] using the prototype code A^{2}, which couples the finitevolume basin simulator ArcTem [11] with the finiteelement geomechanical simulator Code Aster [29] following an iterative algorithm [30]. (See [31] for further A^{2} simulations in the context of rock failure.) This approach reaches very good accuracy but requires long simulation times due to its threedimensional formulation.
Here, rather than derive a threedimensional model such as in [7], we simply insert the horizontal strains into the vertical porositystress law and so keep a onedimensional compaction model without neglecting the remaining dimensions. More precisely, we define an approximated porosity using a vertical porositystress 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 stressstrain constitutive law with stressdependent 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 threedimensional, elastoplastic approach of [7]. This simplified model bears a similarity with that investigated in [32], namely, that the vertical porositystress law is altered to account for a specific effect: in our case, the threedimensional 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 onedimensional porosity stress law, so that our approach is by no means intended to reach the accuracy of any fulldimensional, 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 singlephase 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^{2} and spans a time interval of 10My. Our results under both oedometric and nonoedometric conditions are compared with those obtained using A^{2} 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 zcoordinates.
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 zcoordinates.
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 U_{w} := ϕ(V_{w} − V_{s}) 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 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:(3)where T_{0} and p_{0} 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 singlephase setting, the water’s mobility η_{w} is simply the inverse of its viscosity μ_{w}, which we assumed is temperaturedependent. Specifically, we choose the following law (cf., for example, Eq. (48) of Ref. [28]):(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(5)where A is an anisotropy tensor and K, the intrinsic permeability, follows the Kozeny–Carman law [2–4]: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 ϕ_{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 := (f_{x}, f_{y}, f_{z}) 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 zaxis 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, p_{sup} is the pressure stemming from the weight of seawater and atmosphere lying above the basin. On the rest of the boundary ∂Ω, we impose zeroflux 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 stressstrain 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 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 boundary conditions given in (10):
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 singlelayer model in each layer.
As discussed in the introduction, we wish to integrate horizontal strains in a vertical porositystress law to account for horizontal compaction when oedometric conditions are not verified. We consider that the horizontal strains are known timedependent 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 K_{s} 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, 38–40].
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 p_{w} = p_{sup}, 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 timedependent 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, vectorvalued 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 righthand side of (12) to account for viscous effects [24, 40, 41].
3.2.2 Porosityfree 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 t_{0} and a time t > t_{0}:(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 onedimensional 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 t_{0}, we impose(20)where we refer the reader to (11); this becomes ϕ_{S}(t_{0}) = ϕ_{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 singleexponential 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 nonoedometric conditions. As already mentioned, we suppose that the horizontal strains ε_{x} and ε_{y} (and thus their rates and ) are known timedependent functions.
3.4.1 Stressstrain 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 fourthorder 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.
Strainfree 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 stressstrain 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.
Hooketype law. We choose Hooke’s elastic law with stressdependent 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 stressdependent 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 stressdependent 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 stressstrain 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 ϕ(t_{0}) 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 stressdependent 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 xstress should be preserved (i.e., λ = 0) if only the rate of xstrain is nonzero and viceversa (i.e., λ = 1) if only the rate of ystrain 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 horizontalstrain 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 tracepreserving 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 straindependent, 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 F_{o} 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 nonoedometric 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 stressdependent modulus. To this end, we split oedometric and nonoedometric 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 stressstrain 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)
Stressdependent 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 righthand side of the last inequality can be bounded by and . The first term, however, depends on the nonoedometric 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^{2}. 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 EastWest and 180 km NorthSouth;

Spatial discretization: 100 620 finitevolume cells;

Sediments: 7 groups of different materials;

EastWest strain: 4% from −10My to −8My (2 events) and 2% from −8My to present day (8 events).
Fig. 1 Geometrical and lithological model of the Neuquén basin (from [7]). (a) Representation of the mesh with a fivetime 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 A^{2}
As mentioned in the introduction (Sect. 1), and the references therein, code A^{2} iteratively couples a finitevolume basin simulator (ArcTem) with a three dimensional finiteelement 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 threedimensional 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^{2} 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].
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 ydirections, 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.
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]). 
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 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.
Fig. 3 Results under the strain condition (dble) along W10 and comparison with A^{2} (data points and A^{2} results from [26]). 
Fig. 4 A^{2} overpressure at bottom and over central crosssection under no horizontal strains (top) and reference horizontal strains (bottom) (from [26]). 
Fig. 5 Overpressure at bottom and over central crosssection 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 crosssection 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 ydirections (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)).
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 fulldimensional, 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^{2} team from [26], for a fulldimensional 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.
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 porositystress law and an elastic stressstrain constitutive law involving a stressdependent Young modulus. It coincides with classical modeling in oedometric conditions and provides a computation of porosity in nonoedometric conditions which is simpler, although not as accurate, compared to a fulldimensional, elastoplastic model. On the test case given by the Neuquén basin, our simplified model, coded in the finitevolume 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 crosssection 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 prevalidating results before spending the time and computational resources required by a threedimensional, finiteelement simulator.
Although the model shows promising first results, it would still benefit from further numerical validation by comparing results on different wells, crosssections, 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).
Acknowledgments
We thank the authors of [26] for sharing the A^{2} 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 notforprofit sectors.
References
 Schneider F., Burrus J., Wolf S. (1993) Modelling overpressures by effectivestress/porosity relationships in lowpermeability 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 permeabilityporosity relationship for mudstones, Mar. Pet. Geol. 27, 8, 1692–1697. [CrossRef] [Google Scholar]
 Guy N., Colombo D., Frey J., Cornu T., CacasStentz 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., PegazFiornet 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ínMoreno 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 stressporosity 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 twodimensional modeling of heat transfer, fluid flow, hydrocarbon generation, and migration, AAPG Bull. 74, 3,309–335. [Google Scholar]
 Luo X., Vasseur G., Pouya A., LamoureuxVar 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]
 ObradorsPrats 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 élastoplastique 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., CacasStentz 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., CacasStentz M.C., Cornu T. (2020) Poroelastoplastic 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., CacasStentz 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. http://www.codeaster.org. [Google Scholar]
 Guy N., Enchéry G., Renard G. (2012) Numerical modeling of thermal EOR: Comprehensive coupling of an AMRbased 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., CacasStentz M.C., Cornu T. (2019) An assessment of stress states in passive margin sediments: Iterative hydromechanical 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. https://doi.org/10.1002/essoar.10504323.1. [Google Scholar]
 Brüch A., Maghous S., Ribeiro F.L.B., Dormieux L. (2018) A thermoporomechanical 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 threedimensional 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) Twodimensional finite element analysis of gravitational and lateraldriven 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, MidNorwegian 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 ButterworthHeinemann. [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
All Figures
Fig. 1 Geometrical and lithological model of the Neuquén basin (from [7]). (a) Representation of the mesh with a fivetime 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 
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]). 

In the text 
Fig. 3 Results under the strain condition (dble) along W10 and comparison with A^{2} (data points and A^{2} results from [26]). 

In the text 
Fig. 4 A^{2} overpressure at bottom and over central crosssection under no horizontal strains (top) and reference horizontal strains (bottom) (from [26]). 

In the text 
Fig. 5 Overpressure at bottom and over central crosssection under the strain conditions (zero), (ref), and (dble). 

In the text 
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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.