A simpli ﬁ ed vertical and horizontal geomechanical model for compaction in sedimentary basins

,


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][3][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][11][12][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][15][16][17] and simulators [18][19][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][27][28] using the prototype code A 2 , which couples the finite-volume 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 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 threedimensional, 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 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 one-dimensional 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 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 2 and spans a time interval of 10My. Our results under both oedometric and non-oedometric 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.

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.

Model equations
We write / the Eulerian porosity of the medium and V a the velocity field associated with phase a 2 fs; wg (solid or water). We write q a the density of phase a and q s ¼ q s ð1 À /Þ and q w ¼ q w ð/Þ 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.

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.

Model equations
We write / the Eulerian porosity of the medium and V a the velocity field associated with phase a 2 fs; wg (solid or water). We write q a the density of phase a and q s ¼ q s ð1 À /Þ and q w ¼ q w ð/Þ 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.

Flow
We give here the equations on water and solid flow.
Conservation of mass. For any phase a 2 fs; wg we have conservation of mass: where q a is the rate of deposit of phase a due to sedimentation at the top of the basin and V a is the velocity of phase a.
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, where g 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: where T 0 and p 0 are reference values for the temperature and water pressure, q s,0 and q w,0 are reference densities, a w is the water thermal expansion and b 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 g w is simply the inverse of its viscosity l w , which we assumed is temperature-dependent. Specifically, we choose the following law (cf., for example, Eq. (48) of Ref. [28]): 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 where A is an anisotropy tensor and K, the intrinsic permeability, follows the Kozeny-Carman law [2][3][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.

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, 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 r 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 r ¼ r x r xy r xz r xy r y r yz where r i and r ij are the normal and shear stresses for all i; j 2 fx; y; zg with j 6 ¼ 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 where q :¼ q w þ q s 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 @r x @x þ @r xy @y þ @r xz @z ¼ f x ; @r xy @x þ @r y @y þ @r yz @z ¼ f y ; where we make the assumption that @r xz @x þ @r yz @y ¼ 0: This assumption allows us to compute r 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 r yz , r xz , and r 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]).

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.

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.

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 X the basin and by oX top its top, i.e., @ top ¼ fðx; y; zÞ 2 jz ¼ sðx; yÞg, we impose Here, p sup is the pressure stemming from the weight of seawater and atmosphere lying above the basin. On the rest of the boundary oX, we impose zero-flux conditions, i.e., where m is the outward unit normal of X.
Regarding stresses, a note from (9) that these boundary conditions only close the problem for the vertical stress r z . Indeed, they do not give us the horizontal stresses r x and r y , nor do they give us the shear stresses r yz , r xz , and r 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.

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

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

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.

Effective stress
Biot's theory [34] for consolidation couples fluid flow with rock deformation and introduces the effective stress tensor r 0 , which is defined as where b is Biot's tensor. In [35], the authors link Biot's tensor to the compressibility of the medium according to with 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][39][40]. The effective stress tensor reads The tensor r 0 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,

Skeleton's infinitesimal strain
We denote by e the skeleton's infinitesimal strain, or deformation, tensor. (In the following, we omit the term "infinitesimal" when referring to e.) Because the basin is assumed to be elastic, e coincides in fact with the elastic strain tensor and satisfies |e| ( 1.
We use the following notation: where e i and e ij are the normal and shear strains for all i; j 2 fx; y; zg with i 6 ¼ j. In agreement with our stress convention, any positive normal strain corresponds to a compression.

General compaction law
Given any differentiable time-dependent function f, we write _ f its time derivative, which we also refer to as the rate of f.

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 and F is a continuous, vector-valued function (experimentally fitted), r 0 n :¼ ðr 0 x ; r 0 y ; r 0 z Þ is the diagonal vector of r 0 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 G /; r 0 n À Á Á r 0 n to the right-hand side of (12) to account for viscous effects [24,40,41].

Porosity-free elastic function
For simplicity, we assume that F only depends on effective stress: 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 : 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.

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

Vertical compaction law
Horizontal effects on porosity being negligible in this case, we choose F in (13) to depend exclusively on the vertical stress r 0 z : where b is a function determined by lithology, so that (13) becomes where the notation / S replaces / F . This law is empirical, and one-dimensional in that it only takes vertical stress into account.

Schneider's law
We choose Schneider's function for b [24]: where / 1 , / 2 P 0 (porosities) and r 1 , r 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 a be a primitive of Àb, that is, for some additional parameter / r called the residual porosity. At t 0 , we impose where we refer the reader to (11); this becomes When positive, the residual porosity helps numerically avoid Schneider's porosity becoming negative as vertical effective stress increases. Schneider's law is empirical and its parameters need to be found by experiment, which makes Schneider's porosity a fitted value as opposed to a direct measurement. The double exponential formulation allows for a better fit at both upper and lower sediment layers; indeed, for layers close to the surface (say, above 1000 m), where the effective stress is weak, the porosity generally shows steep variations, while this is not the case at lower depths, where the stress is stronger but the variations flatter. To accommodate both regimes, the double exponential, with its five parameters, proves to be a good choice (cf. [7,40]); in this respect, / 1 and / 2 represent characteristic porosities at their corresponding stress regimes. This is in contrast with the single-exponential law proposed by Athy [42], which is recovered when / 2 = / r = 0 in (18)  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].

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.

Simplified geomechanical model
We now derive the compaction model to handle nonoedometric conditions. As already mentioned, we suppose that the horizontal strains e x and e y (and thus their rates _ e x and _ e y ) are known time-dependent functions.

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 r 0 and e verify where 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 e and so, by time integration, r 0 is an explicit function of e; when furthermore the elasticity is linear, G is, in fact, independent of both r 0 and e 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 r 0 . In fact, in contrast with classical elasticity, we drop the dependence on e and only keep that on r 0 : Equivalently, writing H(r 0 ) = G(r 0 ) À1 , we have 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 where m 2 ½0; 0:5Þ is Poisson's coefficient and E > 0 is Young's stress-dependent modulus. Equivalently, we rewrite (23) as For simplicity, we refer to both (24) and (25) as Hooke's law.
As a direct consequence of (25), we have The trace of _ e is the dilatation or relative volume change of the skeleton. Also, from the first two equations in (25), we establish and so, by summing these two equations and adding _ r 0 z , we get Initially, recalling (11), we impose which becomes r 0 x t 0 ð Þ ¼ r 0 y t 0 ð Þ ¼ 0 when b = 1. Oedometric case. We briefly return to the oedometric case (cf. (15)). Stress and strain tensors. Both equations in (27) and (28) imply Thus, (26) gives Then, the oedometric rates of effective stress and strain tensors read as The Author(s): Science and Technology for Energy Transition 78, 22 (2023) 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), Then, using (17), (21), and (31), and recalling the assumption that the porosity is well approximated by Schneider's porosity, we get 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 r 0 z ), even under oedometric conditions. Note also that the choice of taking Young's modulus as stress-dependent over Poisson's coefficient in (24) and (25) is arbitrary and is only motivated by the fact (34) is a simpler expression than its equivalent for m.

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 (r * , e * ) whose components verify 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 / * : where b is given by Schneider's function (18). We further require that with /(t 0 ) given in (20). Approximated tensors. We let r * and e * satisfy the constitutive relation where 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 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, so that / can be replaced by / * wherever needed. Thanks to (38), this is equivalent to For this approximation to be justified given (37), we still need to choose an appropriate vertical approximated stress r Ã z . Vertical approximated stress. To account for the horizontal effects of compaction, r Ã z must include the horizontal deformations e x and e y . Moreover, for consistency, it should coincide with r 0 z in oedometric conditions. Convex combination. We ask for a convex combination of rates of horizontal stresses to be preserved, that is, for some k 2 ½0; 1 we let Given (27) and (39), this yields where we recall that E(r 0 ) is Young's stress-dependent modulus. Then, note that this formulation only holds if m 6 ¼ 0.
The choice for the parameter k is arbitrary at this point. We impose that only the rate of x-stress should be preserved (i.e., k = 0) if only the rate of x-strain is nonzero and viceversa (i.e., k = 1) if only the rate of y-strain is nonzero. Therefore, we choose k ¼ j_ e y j j_ e x j þ j_ e y j ; and (44) becomes if _ e x _ e y 0: The Author(s): Science and Technology for Energy Transition 78, 22 (2023) Remark 3.1 In oedometric conditions, we note that _ r Ã z ¼ _ r 0 z and, thanks to (36), that r Ã z ¼ r 0 z ; 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: which, thanks to (28) and (40), would lead to This relation, as opposed to our choice (45), disregards the relative magnitude of the horizontal-strain rates, but makes sense when m = 0. Interestingly, thanks to (40) and (46), we have E(r 0 ) tr(e) = E(r * ) tr(e * ), where tr (e) and tr(e * ) 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 r * by tr(r 0 ) 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: for some strain-dependent, continuous function f such that f ð0; 0Þ ¼ 0. For instance, in our case (cf. (45)), we have if _ e x _ e y 0: One could also integrate geometrical constraints related to the basin's shape into the function f. We leave a detailed study of the many possible choices of f to for future work. Porosity. We split the elastic function F in the compaction law (13) according to where F o follows the oedometric formulation given in (16), that is, Recalling / ' / F , this leads to a decomposition of the porosity as where / o and / no satisfy _ / o ¼ F o ðr 0 n Þ Á _ r 0 n and _ / no ¼ F no ðr 0 n Þ Á _ r 0 n and are the oedometric and non-oedometric contributions to the porosity change _ /; in particular, there holds where b 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.

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

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

Code A 2
As mentioned in the introduction (Sect. 1), and the references therein, code A 2 iteratively couples a finite-volume 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 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 2 against our simplified geomechanical model (cf. Sect. 4.4.3 below for CPU times).

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 À1 , are of the form À8 My t;

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

Comparison with [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.
In Figure 4, we reproduce the profile obtained in [26] of the overpressure (i.e., the difference between the water pressure and the hydrostatic pressure) at the bottom of the basin and over a cross-section at the center of the basin. It shows both results when no horizontal deformations and when reference horizontal deformations are applied. Figure 5 presents our results for comparison with Figure 4. Again, in accordance with Figures 2 and 3, we see that we underestimate the overpressure under the strain condition (ref) and get much closer to the A 2 results when (dble) holds instead.
The noted underestimating behavior should not diminish the fact that the results given under (ref) and (dble) still show great improvement compared to the oedometric simulation (zero), in particular for the overpressure, pressure, and porosity in the undrained region of the basin, i.e., in the deep layers.

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

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