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

1 Introduction

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

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

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

One solution is proposed in [7], where Schneider’s porosity-stress law [23, 24] is integrated into the Drucker–Prager model for plasticity [25], resulting in a three-dimensional compaction model. This model is tested numerically in [7, 26–28] using the prototype code A², which couples the finite-volume basin simulator ArcTem [11] with the finite-element geomechanical simulator Code Aster [29] following an iterative algorithm [30]. (See [31] for further A² 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 km² and spans a time interval of 10My. Our results under both oedometric and non-oedometric conditions are compared with those obtained using A² in [26] with the model derived in [7]. Furthermore, CPU times are given for various strain configurations. Finally, we summarize our results and give an outlook in Section 5.

2 Standard porous medium model

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

2.1 Model equations

We write ϕ 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 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₀ and p₀ 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 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₀ is the specific surface area of the porous medium and k₁, k₂, n₁, n₂, m₁, m₂ and ϕ₀ are constants related to lithology. To ensure the continuity of K at ϕ₀, 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 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, 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 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 t₀ 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₀ 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 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 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 t₀ and a time t > t₀:(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 ϕ₁, ϕ₂ ⩾ 0 (porosities) and σ₁, σ₂ > 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₀, we impose(20)where we refer the reader to (11); this becomes ϕ_S(t₀) = ϕ_r + ϕ₁ + ϕ₂ 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, ϕ₁ and ϕ₂ 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 ϕ₂ = ϕ_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(σ′)⁻¹, 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 ϕ(t₀) 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 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 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 A². We refer the reader to [26] for a detailed account of the basin’s history and sediments’ properties; importantly, the basin shows high overpressures, and drained conditions with hydrostatic pressure in the upper layers, in contrast with the lower layers.

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

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

• Spatial discretization: 100 620 finite-volume cells;

• Sediments: 7 groups of different materials;

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

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 A²

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

4.3 Parameters

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

In addition, the horizontal strain rates, measured in s⁻¹, are of the 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 A², 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²) 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². 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², 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 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 A² 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 A2 (data points and A2 results from [26]).

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

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 A² 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)).

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 A² team from [26], for a full-dimensional run with the strain condition (ref) lies between 2.5 h and 6.5 h on 72 processors, depending on the chosen underlying linear solver; in contrast, the times given in Table 2 were obtained on a single processor on a personal laptop computer.

5 Conclusion and outlook

The model proposed in this paper to describe mechanical compaction is based on an alteration of Schneider’s vertical porosity-stress law and an elastic stress-strain constitutive law involving a stress-dependent Young modulus. It coincides with classical modeling in oedometric conditions and provides a computation of porosity in non-oedometric conditions which is simpler, although not as accurate, compared to a full-dimensional, elastoplastic model. On the test case given by the Neuquén basin, our simplified model, coded in the finite-volume simulator ArcTem, shows reasonable accuracy in pressure, porosity, and vertical stress along a drilling well, as well as in overpressure at the bottom and over a vertical cross-section at the center of the basin. Overall, our model tends to underestimate the horizontal stresses, the pressure, and the overpressure. Still, it shows significant improvement compared to an oedometric model, in particular in the lower, undrained compartment of the basin, and it produces very advantageous computational times. Therefore, it may offer a quick way to get pre-validating results before spending the time and computational resources required by a three-dimensional, finite-element simulator.

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

Acknowledgments

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

References

All Tables

All Figures

	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

	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

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

	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

	Fig. 5 Overpressure at bottom and over central cross-section 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