A constitutive and numerical model for mechanical compaction in sedimentary basins

Computers and Geotechnics 33 (2006) 316–329 www.elsevier.com/locate/compgeo

A constitutive and numerical model for mechanical compaction in sedimentary basins Denise Bernaud, Luc Dormieux *, Samir Maghous ENPC-LMSGC, 6–8, avenue Blaise Pascal, Cite´ Descartes, Champs-sur-Marne, 77455 Marne-La-Valle´e, France Received 30 January 2006; received in revised form 10 May 2006; accepted 15 May 2006 Available online 2 October 2006

Abstract The aim of the paper is to provide new elements concerning the constitutive behavior of sedimentary rocks and the numerical aspects for basin simulators. A comprehensive model for mechanical compaction of sedimentary basins is developed within finite poroplasticity setting. Particular concern is paid to the effects of large porosity changes on the poromechanical properties of the sediment material. A simplified micromechanics-based approach is used to account for the stiffness increase and hardening induced by large plastic strains. A key challenge for numerical assessment of sedimentary basin evolution is to integrate multiple coupled processes in the context of open material systems. To this end, a numerical approach inspired from the ‘deactivation/reactivation’ method used for the simulation of excavation process and lining placement in tunnel engineering, has been developed. Periods of sediments accretion are simulated by progressive activation of the gravity forces within a fictitious closed system. Fundamental components of the constitutive model developed before (hydromechanical coupling, dependence of poroelastic properties on large plasticity, impact of irreversible porosity changes on the hardening rule, evolution of permeability with porosity) are included into our finite element code. Illustrative examples of basin simulation are performed in the one-dimensional case. Various aspects of the constitutive model are investigated. Their influence on the corresponding basin response is analyzed in terms of compaction law, porosity and fluid pressure profiles.  2006 Elsevier Ltd. All rights reserved. Keywords: Mechanical compaction; Sedimentary basin; Finite poroplasticity; Numerical approach

1. Introduction Simulation of sedimentary basins is a complex multidisciplinary problem involving geological, chemical and mechanical aspects. Reconstructing the stress and deformation history of a sedimentary basin is a challenging and important problem in geoscience. The potential applications include petroleum exploration, reserve assessment and production. Many efforts have been done during the past decades as regards the integration of geological data and the simulation of hydrocarbon generation, primary migration, secondary migration, and accumulation processes. However, *

Corresponding author. Fax: +33 1 64153748. E-mail address: [email protected] (L. Dormieux).

0266-352X/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2006.05.004

progress in this field has been hampered by the absence of a comprehensive mechanical description of the geological material, which accounts for the strongly coupled nature of the deformation problem. Key components for mechanical basin modeling are  large irreversible strains setting,  hydromechanical coupling,  interaction of macroscopic poromechanical constitutive response with evolving morphology of the sediment microstructure. Sedimentary basins, such as those in the North Sea, form when waterborne sediments in shallow seas are deposited over periods of tens of millions of years. The deposited material then compacts under its own weight,

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

causing a reduction of porosity and hence the expulsion of pore fluid. Eventually, as depth increases, chemical reactions occur, such as cementation of granular aggregates and pressure solution. Actually, the mechanical response of a sedimentary basin is the consequence of complex processes involving mechanical, geochemical, geophysical and geological aspects. For the mechanical analysis at the macroscopic scale, two phenomena contribute to the compaction of sediments: (1) purely-mechanical compaction which originates mainly from rearrangement of the solid particles during burial and (2) mechano-chemical compaction resulting from dissolution–precipitation mechanisms, generally induced by stress (pression–dissolution). Purely mechanical phenomena prevail in the upper layers, whereas chemical compaction dominates for deeper burial as stress and temperature increase (e.g. [37]). Disregarding chemical aspects, the paper will focus on the modeling of purely mechanical compaction. The basic models of mechanical compaction are still based on phenomenological relationships relating porosity to effective vertical stress. The concept of porosity versus Terzaghi’s effective stress dependence has been early introduced by Hubert and Rubey [29] and later by Smith [41]. These ideas have been widely adopted and implemented in numerical finite element models that have yielded valuable contributions to the understanding of the evolution of sedimentary basins [40,10,1,45,46,11,38,39,20,21,31,47,42, 16,33,44,2,28,43,22,36], to cite a few. Still, a more comprehensive description of the mechanics involved in basin simulation should be achieved within the tensorial formalism of the constitutive model. This is necessary for addressing the boundary conditions encountered in the case of tectonic activity, or for the determination of the horizontal stresses induced by the compaction phenomenon. To deal with, the main difficulty rises from the large porosity changes involved in compaction problem. This requires that the poromechanical constitutive law be formulated in the framework of finite irreversible strains [23,24,1,45,43,13,17,15,4,3,5]. The coupled nature of the deformation problem may be understood as follows. Large strains modify the microstructure, microstructure modifications lead to a change in the poromechanical properties of the sediment material, which in turn affect the basin response. Theoretical formulation accounting only for the effects of the porosity changes on the elastic and plastic properties of the porous material have been proposed in [19,8,6,17]. It is shown in particular that the time dependence of the poroelastic parameters with large plastic strains leads to additional terms in the constitutive equations. Appropriate micromechanics estimates have been proposed to account for the increase of elastic stiffness as the porosity decreases. A large part of this work is devoted to the numerical aspects of basin simulations. It is well known that dealing with open material systems like sedimentary basins in the framework of finite element procedures in a strongly non-

317

linear context is not an easy task. A specific computational technique is developed for this purpose. Schematically, it consists in replacing the real open system material by a fictitious one whose evolution under loading-unloading (i.e., deposition–erosion) is directly connected to that of the basin. The technique is widely inspired of the so-called ‘deactivation–reactivation’ method used for the simulation of excavation process and lining placement in tunnel engineering (see for instance [27,9]). This technique is implemented within a finite element tool incorporating the constitutive model formulated in large poroplasticity. The paper is organized as follows. Section 2 describes the micromechanics-based constitutive model for the sedimentary model within the finite poroplasticity framework. Finite element analysis using the updated Lagrangian approach is then briefly presented. In Section 3, the ‘deactivation/reactivation’ procedure is developed to deal with basins simulations during phases of deposition or erosion. Examples of numerical simulations under oedometric conditions are finally given in Section 4. 2. Constitutive equations During compaction process, the particles of the porous material (i.e., sediment material) are subjected to large volumetric strains: the porosity change may reach 50% as the burial goes on. Appropriate prediction of the basin response should thus take into account both the constitutive nonlinearities and the geometric nonlinearities induced by large strains. 2.1. State equations and complementary relations At the scale of the pores and solid grains, large strains modify the microstructure of the sediment material. In turn, these microstructural changes are responsible for the evolution of elastic and plastic mechanical properties. The sedimentary rock is modelled as a fully saturated poroelastoplastic material undergoing large strains. The anisotropy of the mechanical properties of the sedimentary material induced by compaction is disregarded. The elastic part of the deformation gradient of the skeleton particles is assumed to remain infinitesimal. This means that large strains involved during compaction process are of irreversible (plastic) nature. In the framework of finite poroplasticity, the constitutive behavior comprises two state equations in rate-type formulation together with complementary relations specifying the flow rule [19,8]. Let r denote the Cauchy stress tensor and p the pore pressure. The first state equation _ the pore pressure rate p, _ and the relates the stress rate r, strain rate tensor d defined as the symmetric part of the velocity gradient DJ r0e ¼ r_ 0e þ r0e  X  X  r0e ¼ c : ðd  d p Þ þ c_ : c 1 : r0e    Dt ð1Þ

318

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

where dp denotes the plastic strain rate while X is the rotation (spin) rate tensor which aims at taking the large rotation of the elementary volume into account. This equation involves:  a rotational time derivative of the Biot effective stress tensor r0e ¼ r þ bp1, where b is the Biot coefficient;  a term related to the particulate derivative c_ of the ten sor of drained elastic moduli c. The evolution of the  elastic properties is related to those of the microstructure which cannot be neglected in the domain of large strains. Relationship (1) extends classical rate form formulations [34,13,18] to the case of variable elastic properties. The second state equation relates the pore volume change to the rate p_ of the pore pressure and to the strain rate d. In view of its derivation, it is convenient to normalize the pore volume dXpore in any configuration of the representative elementary volume by the total volume dX0 of the latter in a reference configuration [25]. This yields the so-called lagrangian porosity / = dXpore/dX0 that is related to the classical porosity u (defined as the pore volume fraction in the current configuration) by / = Ju. J is the jacobian of the transformation of the elementary volume, defined as the ratio between its volume in the current configuration and its initial one. The plastic part /p of the lagrangian porosity is defined as the lagrangian porosity in the unloaded configuration of the elementary volume. Similarly, the plastic part Jp of the jacobian is defined as the jacobian in this unloaded configuration. Note that the rates J_ and J_ p are given by J_ ¼ J tr d

and

J_ p ¼ J p tr d p

ð2Þ

With these notations, the rate form of the second state equation reads [8] ! _ /_  /_p M p p_ ¼ M b tr ðd  d Þ þ þ p  M b_ tr ðc 1 : r0e Þ p  M J ð3Þ _ and b_ in (3) M is the Biot modulus. The terms involving M are related to the influence of large plastic strains on the poroelastic properties. The complementary equations prescribe the plastic flow rule. In the theory of poroplasticity, the latter deals with the plastic strain rate and with the rate of the plastic part of the lagrangian porosity (e.g. [12]). Generalizing the concept of plastic potential, we therefore introduce the function g(r, p) which derivatives yield the sought plastic rates d p ¼ v_

og ; or

og /_ p ¼ J p v_ op

ð4Þ

where v_ is a non-negative plastic multiplier. The plastic incompressibility of the solid phase is assumed in the sequel. This implies that the plastic part of the pore volume change is identical to the plastic part

of the total volume change and yields /_ p ¼ J_ p . Owing to (2), it is further obtained that /_ p ¼ J p tr d p . Combining this equation with (4) reveals that the plastic potential depends on r and p through the so-called Terzaghi effective stress r 0 = r + p1 (e.g. [12]) gðr; pÞ ¼ Gðr0 Þ;

d p ¼ v_

oG or0

ð5Þ

The relevancy of the above model depends on the capability to specify: 1. the elastic–plastic coupling characterized by the dependence of tensor c as well as of poroelastic coefficients  b and M on large plastic strains; 2. the influence of large plastic strains on the evolution of the yield surface as well as on the hardening rule. This is the purpose of the next section, where this question is addressed within the framework of a simplified micromechanics-based approach. 2.2. Influence of microstructural changes on poroelastic and hardening properties In the case of plastic incompressibility of the solid matrix, large plastic strains of the porous material are expected to induce significant porosity and pore shape changes. In order to capture the influence of the plastic strains on the poroelastic properties, the idea is to resort to micromechanical estimates of the latter. For the sake of simplicity, the anisotropy induced during the loading process is disregarded in this study. The pore space is therefore entirely characterized by its volume fraction, namely the porosity u. We herein adopt the Hashin–Shtrikman upper bounds which are known to reasonably model the elastic properties of isotropic porous media (see for instance [48]). Accordingly, the bulk K and shear l moduli of the porous medium now appear as functions of the porosity as well as of the elastic properties of the solid phase (that are assumed to be constant) KðuÞ ¼

4k s ls ð1  uÞ ; 3k s u þ 4ls

lðuÞ ¼

ls ð1  uÞð9k s þ 8ls Þ k s ð9 þ 6uÞ þ ls ð8 þ 12uÞ ð6Þ

where ks and ls are bulk and shear moduli of the solid phase. It is recalled that the Biot coefficient and modulus are connected to K through bðuÞ ¼ 1 

KðuÞ ks

and

1 bðuÞ  u ¼ MðuÞ ks

ð7Þ

From a practical point of view, the determination of the poroelastic parameters K(u), l(u), b(u) and M(u) should not be based on a direct measurement of ks and ls, but from the identification of the latter from the inverse analysis of macroscopic tests. More precisely, K(u0), l(u0) are first evaluated from laboratory tests performed on samples

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

of the porous material in its initial state (i.e., / = /0). In a second step, ks and ls are determined from Eq. (6) written for / = /0 and in which K(u0) and l(u0) are known quantities. We now need to relate the porosity change to the plastic strains undergone by the porous material. To do so, we first observe that the condition of plastic incompressibility of the solid constituent reads J p  /p ¼ 1  u0

ð8Þ

where u0 denotes the initial value of the porosity. In the framework of infinitesimal elastic strains, it is possible to neglect the variation of the pore volume and that of the total volume between the loaded and the unloaded configurations of the elementary volume, since these variations are reversible by definition. This justifies the following approximations: Jp  J;

/p  / 1  u0 1  u0 1 J Jp

pc ðuÞ ¼ pco ð10Þ

In view of (10), Eqs. (6) show that the macroscopic stiffness tensor c is a function of the (total or plastic) jacobian :  c ¼ c ðJ p Þ. The same conclusion holds for the poroelastic   coefficients b and M. Neglecting the induced anisotropy therefore amounts to the fact that the elastic–plastic coupling properties is only governed by the plastic volumetric strains. As regards the evolution of the plastic properties of the porous medium, it is first assumed the yield surface f is that of the standard modified Cam–Clay (e.g. [14,35]) 3 f ðr0 ¼ r þ p1; pc Þ ¼ s : s þ M 2cs p0 ðp0 þ pc Þ 2

ð11Þ

where s ¼ r  13 tr r is the deviatoric stress tensor, p0 ¼ 13 tr r0 is the mean effective stress. pc is the consolidation pressure and represents the hardening parameter of the model. The constant Mcs is the slope of the critical state line. The plastic flow rule is associated, i.e., g = f. It is assumed in the sequel that the shape of the yield locus is not affected by large plastic strains. In contrast, the hardening law, that is the influence of large plastic strains on the consolidation pressure, is a crucial feature of the model. In the framework of small strain plasticity, let us recall that the classical hardening law in the Cam–Clay model can be written as a function of the plastic volumetric strain ep in the form p

pc ðep Þ ¼ pco eae

ð12Þ

where a is a material constant, that can be fitted from macroscopic isotropic compression tests. Hence, a simple way to generalize (12) in the domain of large strains consists in replacing ep by Jp  1, since both quantities are equal in the range of small plastic strains

p 1Þ

ð13Þ

However, such a hardening law does not enforce the condition Jp > 1  u0, corresponding to total pore closure (see (10)). It is therefore expected that a simple extension of the classical Cam–Clay hardening law in the form (13) might yield negative porosities under high isotropic compression. That is of course not satisfactory. Instead of (13), we shall hereafter resort to a micromechanical approach to the hardening law which overcomes the above difficulty. The idea is to identify pc(Jp) to the limit load of a hollow sphere subjected to isotropic compression in the domain of large strains. In fact, as early noticed by Hashin, the hollow sphere is the simplest conceptual model for an isotropic porous medium. Considering the case of a solid phase of the Tresca type, the limit load appears to be a function of the current porosity u (pore volume fraction)

ð9Þ

Introducing (9) into (8) yields u1

pc ðJ p Þ ¼ pco eaðJ

319

ln u ln u0

ð14Þ

Combining this result with (10) yields the hardening law pc(Jp) in the form   p 1  u0 pc ðJ p Þ ¼ co ln 1  ð15Þ ln u0 Jp As opposed to (12), we note that the consolidation pressure predicted by (15) tends towards infinity when the pore space vanishes (Jp  J ! 1  u0). This is the way this micromechanics-based hardening law avoids the development of negative porosities. Furthermore, it is possible to show that the dominating hardening mechanism at large strains is associated with large geometry change, and not with the increase of trapped elastic energy as in the framework of small strains analysis [6]. 2.3. Finite element analysis Consider a material system defined by porous medium which occupy a geometrical domain X. Within the framework of finite poroplasticity, the quasi-static boundary value problem is defined on X by  the momentum balance equation: div r þ qg ¼ 0

ð16Þ

q is the total (saturated) porous material density and g is the acceleration of gravity.  the fluid mass balance equation: _ qf / þ J div qf q ¼ 0

ð17Þ

where qf is the fluid density and q is the filtration vector. The latter is connected to the fluid pressure through the Darcy’s law

320

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

q ¼ k  ðrp þ qf gÞ

ð18Þ

k denotes the permeability tensor. In the case of isotropy (i.e., k = k 1), effects of microstructural changes on the evolution of k may be modeled by means of the Kozeny–Carman formula 2

k ¼ k0

u3 ð1  u0 Þ u30 ð1  uÞ2

ð19Þ

k0 = k(u0) is the initial value of the permeability.  the constitutive equations: Within the framework of finite poroplasticity, the latter are given by (1) and (3) together with the plastic flow rules (5). The elastic moduli, the yield surface and the hardening law are described by (6), (7), (11) and (15), respectively. In addition, the initial conditions must be specified as well as the mechanical and hydraulic boundary conditions. The latter are written on the boundary oX of X. It must be emphasized that all the equations defining the boundary value problem refer to the mechanical system in its current configuration, which is a priori unknown. The finite element analysis will be outlined hereafter. Fully details are given in Bernaud et al. [8]. The updated Lagrangian scheme [7] is used in order to analyze the large strain behavior of the structure under consideration. In this approach, all static and kinematic variables are referred to an updated configuration in each step time Dt. Denoting by tx the coordinates of any skeleton particle in the configuration of the mechanical system at time t, it is convenient to reformulate the problem in terms of the displacement between t and t + Dt U ¼ tþDt x  t x

ð20Þ

and of the pressure values difference P at points similar within the skeleton transformation at times t and t + Dt P ¼ pðtþDt xÞ  pðt xÞ

ð21Þ

The weak forms of the equilibrium and fluid mass balance equations, expressed at time t + Dt, are derived. A standard 2D finite element formulation employing 6-nodes triangle elements is used. The components of the displacement U are approximated by a quadratic polynomial function of the 6 nodal values. The approximation of the pressure P is linear with corresponding nodes located at the three summits of the triangle. This procedure detailed in Bernaud et al. [8] yields the following matrix equation: " #" # " # KUU KPU U FU ¼ ð22Þ KTUP KPP P FP

where KIJ and FI are the global stiffness sub-matrices and force sub-vectors, respectively. U and P are the global vectors whose components are node values of the displacement U and the pressure P, respectively. The strategy used for deriving (22) consists in rearranging the discretized form of the equilibrium and fluid mass balance equations in such a manner that all the nonlinear terms in U and P are put into the right-hand member T [FU FP]. It results from this strategy that the sub-matrices KIJ do not depend on U and P. Consequently, the above system of equations is nonlinear because of the dependence of FU and FP on U and P. For this reason, at each time, it is necessary to solve (22) by an iterative algorithm until (22) is satisfied up to a required tolerance. Explicit expression for matrices KIJ and vectors FI together with the description of the iteration procedure are given in Bernaud et al. [8]. 3. Numerical procedure for sedimentary basins simulation The specificity of the numerical simulation of sedimentary basins lies in the fact that the system under consideration is an open system. Consequently, appropriate numerical approaches such as finite element one must implement appropriate procedures to model continuous sediment supply. From a technical point of view, this issue may be addressed in the same way as in Tuncay et al. [43]: the finite element grid accretes with sediment infilling, and a new sediment layer is introduced when the sediment layer at the top of the basin reaches a critical thickness. The numerical approach developed in this paper, which is aimed to simulate the time evolution of the basin boundary, is based upon a technique directly inspired from tunnel engineering and may be easily incorporated within existing finite element tools. 3.1. Problem definition The basin undergoing compaction is modeled as an horizontal infinite layer (Fig. 1). The basement rock is located at z = 0. Assuming it remains horizontal, the upper boundary is defined by the time dependent plane equation z = H(t), where H(t) refers to the total thickness of the basin at time t. Besides, the sea level is located at the plane z = L0. These geometrical assumptions are adopted for simplicity. Nevertheless, more complex 2D basin geometries can also be handled by means of the numerical code. In absence of tectonic activity, the gravity forces gez and the fluid pressure at the top of the basin constitute the loading parameters in the compaction process. If q(z, t) denotes the mass density of the sediment material, the mass of sediments supplied per unit area from the initial time t = 0 is equal to that of the vertical column with unit cross-section Z H ðtÞ M d ðtÞ ¼ qðz; tÞ dz ð23Þ 0

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

q0(t) = q(H(t), t) is the mass density of the sediment material in the reference (i.e., initial) state. In the case of homogeneous material accretion, that is, if q0(t) is time independent, (26) reduces to H ¼ M d ðT Þ=q0 . Actually, H represents the thickness of the basin at time T if the constitutive material were rigid. It is assumed that no additional sediments are deposited and no erosional phenomena take place beyond time T. Consequently, the system evolves as a closed system for t > T (except for the fluid mass exchange associated with the pore pressure dissipation). Actually, the system is termed here as closed because the corresponding finite element grid comprises a constant number of finite elements and nodes. The geometrical domain X0 is regarded as the initial configuration of the fictitious system material whose constituent is defined as follows:

Fig. 1. Schematic geometry of the basin.

_ d ðtÞ is It is assumed in the sequel that the accretion rate M prescribed. The boundary conditions consist in prescribing the values of the stress vector T = r Æ ez, displacement vector n, fluid pressure p and fluid flux q Æ ez

p ¼ qf gðL0  H ðtÞÞ

ð24Þ

 z = 0 (basin basement) n ¼ 0;

q  ez ¼ 0 ðimpermeability conditionÞ

ð25Þ

3.2. Technical numerical aspects The key idea of the numerical technique consists in transforming the real open material system (the basin) into a fictitious closed system. The evolution in time of the latter must of course correspond to that of the real system. The main advantage of this procedure lies in the fact that the finite element approach briefly described in Section 2.3 is applicable. The principle is inspired from the so-called ‘deactivation/reactivation’ method (see for instance [27,9]). This technique was initially developed and implemented for numerical simulation in tunnel engineering (simulation of excavation process and lining placement). The idea consists here in activating sediment layers as they are effectively deposited by substituting their fluid characteristics with material sediment ones. This process amounts to a stepwise activation of the gravity forces in the upper layers. More precisely, let us assume that the simulation of the basin formation is carried out during the time interval t 2 [0, T]. We hence introduce the geometrical domain X0 defined by 0 6 z 6 H, where Z T H¼ M_ d ðtÞ=q0 ðtÞ dt ð26Þ 0

 mass density: qf;  permeability: jkj  k0;  values of the poroelastic and poroplastic parameters actually do not matter. For convenience, they are set to those of the sediment material in its reference state. The initial stress and fluid pressure fields correspond to hydrostatic state

 z = H(t) (upper surface) T ¼ qf gðL0  H ðtÞÞez ;

321

t¼0

r ¼ qf gðL0  zÞ1;

t¼0

p ¼ qf gðL0  zÞ

ð27Þ

The spatial discretization of X0 is carried out as follows. The time interval of study is subdivided into n sub-intervals [ti1, ti], with t0 = 0 and tn = T. During the time increment Dti = ti  ti1, the dump of sediments would correspond (if the material were rigid) to the thickness increase of the basin Z ti M_ d ðtÞ=q0 ðtÞ dt ð28Þ DH i ¼ ti1

The geometrical domain X0 is subdivided into n layers of thickness DHi, i 2 [1, n]. The finite element grid of X0 is therefore generated by discretizing each layer (Fig. 2) using for instance triangular elements (6 nodes for the

Fig. 2. Initial configuration of the fictitious closed system and principle of the F.E. discretization.

322

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

displacement approximation, 3 nodes for the pressure approximation). The layer i in the initial configuration, of thickness DHi, is defined by Hi1 6 z 6 Hi = Hi1 + DHi. The issue consists in determining the evolution of the fictitious closed system material as the sediment is supplied and the layers (F.E. grid) deform with the compaction process. As mentioned above, the sediment layers are activated as they added by setting their physical parameters to those of the real material in its reference state. The layer i is defined by H ji1 6 z 6 H ji , where z ¼ H ji refers to the current (i.e., at time tj) position of the material plane located initially at z = Hi. The geometrical volume occupied by the (fictitious) porous material at time tj is denoted as Xj [ j  Xj ¼ H i1 6 z 6 H ji ð29Þ i¼1;n

 If i 6 j, the layer i has already been activated and the mechanical properties are those of the sediment material in the current state, since they evolved with compaction from time ti of activation.  In contrast, if i > j, the layer i is not activated yet and the mechanical properties correspond to those of the initial fictitious configuration. The local fields at time tj are the solutions to successive finite poroplastic boundary value problems performed during time intervals [ti1, ti]. The validity of this technique requires that the boundary conditions (24), where H ðtÞ ¼ H jj , acting on the top of the real basin (i.e., activated zone) are satisfied. In order to meet these boundary conditions, the poromechanical characteristics of the fictitious material (i.e., in the non-activated zone) must be chosen carefully. The value of the fictitious material permeability is taken very high, i.e., jkj  k0. This implies that the fluid pressure distribution within the non-activated layers corresponds to the hydrostatic regime. Besides, the mass density of the fictitious material being taken equal to qf, the global equilibrium

of the non-activated zone shows that stress vector acting on the top of the activated zone, i.e., at z ¼ H jj , is exactly r  ez ¼ qf gðL0  H jj Þ1. In practice, the value of the Biot coefficient of the fictitious material, taken equal to that of the sediment material in its reference state, is almost equal to unity. Consequently, the whole sub-domain of non-activated layers undergoes a rigid body displacement. We now examine the determination of the hydromechanical state at time tj+1 = tj + Dtj. The layer (j + 1) is activated when starting the calculations for t 2 [tj, tj+1]. The initial conditions for this purpose are (Fig. 3):  i 6 j: the stress, fluid pressure and hydromechanical properties result from previous calculations, actually, those which prevail at time tj.  i = j + 1: the stress and fluid pressure are still given by (27). The hydromechanical properties are those of the sediment material that is deposited at time tj+1 (i.e., reference state): q0(tj+1), u0(tj+1), c ðu0 ðtjþ1 ÞÞ, k(u0(tj+1)), . . .   i > j + 1: the layer i is not activated and the corresponding mechanical properties are those of the initial configuration: qf, jkj  k0, . . . It is observed that conditions (24) and (25) on the boundary of Xj still hold provided H(t) is replaced with H jn . 3.3. Comments  The geometrical domain X0 to be discretized has a bounded extension in the x- and z-directions. Appropriate boundary conditions imposed on the lateral sides depend on the problem nature (one-dimensional compaction, tectonic loading, . . .).  Each layer is subdivided into finite elements. The time increment Dtj corresponding to the activation (supply) of layer j is itself subdivided into time sub-increment in order to better capture the hydromechanical evolution of the basin.  The numerical procedure has been described assuming that the layers remain horizontal (i.e., 1D compaction). However, its applicability goes far beyond this particular situation, and the layers may deform according to the loading process.  Periods of non-deposition and erosional events are simulated in the same way. For example, the simulation of a phase of erosion consists in progressive deactivation of layers already deposited by setting their poromechanical properties to those of the non-activated layers (i.e., qf, jkj  k0, . . .).

3.4. Validation of the numerical technique

Fig. 3. Configuration of the system when starting calculations for t 2 [tj, tj+1].

To check the relevance of our numerical procedure, we compare finite element predictions with the closed form solution obtained for an academic situation of a

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

1 f ðr; pc Þ ¼  tr r  pc 3

6000

4000

H(t) (m)

one-dimensional simulation of the compaction process as detailed in Deude´ et al. [17]. Effects of large strains on the changing elastic and plastic mechanical properties are disregarded. The problem is treated in drained conditions (no fluid overpressure). It must be kept in mind that the purpose here is to verify the numerical procedure described in Section 3.2, and not the constitutive model. The plasticity is modeled by a simplified cap model with associated flow rule, where the straight-line corresponding to positive hardening is defined by ð30Þ

dp m ¼  pc dJ

ð31Þ

is taken constant (i.e., linear variations of pc with respect to the plastic volumetric strains). Denoting by k and l the elastic Lame´ coefficients of the sediment material, expression of the compaction law writes [17] 8  h i g > < kþ2l 1  exp  M ðtÞ t 6 te d q0 g kþ2l  h i H ðtÞ ¼ e > : H e þ bK 1  exp  bg ðM d ðtÞ  M d ðte ÞÞ t P te q0 g ð32Þ

0

pco k þ 2l K g

ð34Þ

40 t (million years)

60

h i 8 q0 g > ðk þ 2lÞ ln 1 þ ðz  H ðtÞÞ > kþ2l > > > > > e < H ðtÞ  H 6 z 6 H ðtÞ rv ¼ h i > > > b ln Ke þ qb0 g ðz þ H e  H ðtÞÞ þ ðk þ 2l  bÞ ln Ke > > > > : 0 6 z 6 H ðtÞ  H e ð36Þ The horizontal stress writes

H ðtÞ  H e 6 z 6 H ðtÞ ð37Þ

0 6 z 6 H ðtÞ  H e

0

ð33Þ

−10 σ v (MPa)

te is the time of occurrence of the plasticity in the basin and it is defined by the relation

20

The expression of the vertical stress rv is

where He is the thickness of the superficial elastic crust k þ 2l ð1  epco =K Þ q0 g

0

Fig. 4. Numerical and analytical results for the compaction law.

8 h i q0 g > < k ln 1 þ kþ2l ðz  H ðtÞÞ h i rv ¼ > : ðb  2lÞ ln Ke þ qb0 g ðz þ H e  H ðtÞÞ þ ðk þ 2l  bÞ ln Ke

M d ðte Þ ¼

Analytical solution Numerical solution

2000

pc is identical to the consolidation pressure of the Cam– Clay model. In addition, the hardening modulus m defined as

He ¼

323

−20 Analytical solution Numerical solution

−30

the parameters b and Ke are −40 2

K ; b ¼ k þ 2l  mþK

qg Ke ¼ 1  0 H e k þ 2l

ð35Þ

0

2000

4000

6000

depth (m) Fig. 5. Vertical stress profile along the basin after t = 60 million years.

324

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

The model parameters are taken as q0 = 0.6468103 kg/ m , Young modulus E = 103 MPa, Poisson’s ratio m = 0.33, hardening modulus m = 100 MPa, pco = 1 MPa. The rate of sediment accretion is constant and equal to 2.05109 kg/s per unit area, which corresponds to 100 m of sediment thickness per million years in unloaded conditions. Numerical calculations have been performed to evaluate the mechanical state of the basin after t = 60 million years. As shown in Figs. 4 and 5, the match between the analytical and numerical results turns to be excellent in terms of predicted stress profile and compaction law (maximum discrepancy less than 4%). 3

4. One-dimensional illustrative basin simulations The basin deformation model described in the previous sections has been implemented in a two-dimensional finite element code. The results presented below illustrate some aspects of basin response in the case of one-dimensional compaction process (i.e., oedometric conditions). In order to emphasize the role of the coupling between large strains and the mechanical properties during sedimentation, three models are considered.  Model M1: standard uncoupled model. This model corresponds to constant elastic moduli and permeability. The hardening law, derived from the evolution of the consolidation pressure, is given by (13).  Model M2: partially coupled model. The effects of large plastic strains on elastic and plastic mechanical properties are accounted for via the porosity changes through (6),(7) and (15). In contrast, these effects on the permeability value k are disregarded, i.e., k is kept constant throughout the simulation.  Model M3: fully coupled model. This model is identical to model M2 except for the evolution of the permeability parameter which is modeled by Kozeny–Carman formula (19).

In view of (b) and (c), the time T (i.e., end of material supply) needed to reach this virtual thickness H would be _d T ¼ q0 H=M

Furthermore, the following material parameters are fixed: initial material density q0 = 1.37103 kg/m3 and initial porosity u0 = 0.72 (data taken from Hamilton [26]), fluid density qf = 103 kg/m3, initial Young modulus E0 = 103 MPa, initial Poisson’s ratio m0 = 0.33, initial Biot coefficient b0 = 0.9715, initial Biot modulus M0 = 1.392105 MPa, coefficient slope of the critical state line Mcs = 1.2, permeability k0 = 1010 MPa1 m2 s1, initial consolidation pressure pco = 1.5 MPa. _ d is taken equal to 2.05109 kg/s The accretion rate M per unit area (i.e., 100 m/My). According to (38), this corresponds to T . 60 million years. Finally, the sea level is located at (see Fig. 1) z = L0 = 8000 m. 4.1. Uncoupled model M1 The model M1 is the simplest one. It is recalled that the poroelastic and hydraulic parameters are kept constant and equal to their initial values (defined above) during compaction. The hardening parameter which is involved in Eq. (13) is taken as a = 3.882. In this section, the simulations are performed for two values of the initial consolidation pressure: pco = 1.15 and 1.5 MPa. The numerical prediction for the compaction law is displayed in Fig. 6. The increasing part of the curves corresponds to the period of sediment deposition t 2 [0, T]. The post-peak phase t > T, during which the sediment supply is interrupted, is the phase of overpressures dissipation. The latter is accompanied by a strong compaction of the basin leading to a significant reduction of its thickness. Once the overpressures have dissipated, the hydrostatic

The simulations carried out in the sequel do not refer to a real case, but only to academic situations. Actually, the purpose here is to illustrate some of the basin phenomena supported by the theoretical–numerical model. For comparison purposes and sake of simplicity, some common characteristics are considered for all models:

6000

4000

H(t) (m)

(a) The yield surface is that of the standard modified Cam–Clay (11). (b) The sediment material supplied at the top of the basin is assumed to have constant mechanical and hydraulic properties at the geologic time scale. _ d is constant at the geologic time (c) The accretion rate M scale (t 2 [0, T]). (d) The basin is made up of an amount of material supply which would correspond to a vertical column of thickness H ¼ 6000 m in absence of compaction (see Fig. 2), or in other words if the sediment material were rigid.

ð38Þ

2000 pco =1.5 MPa pco =1.15 MPa

0

0

T

200

400

600

t (million years) Fig. 6. Numerical prediction of the compaction law for model M1, T . 60 million years.

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

function with a maximum value at the bottom. For t > T, the basin behaves as a closed material system and the dissipation of excess pore pressure proceeds until the hydrostatic regime is recovered for t > t1. Fig. 8 represent the evolution of the porosity along the basin for typical instants. These figures indicate that the porosity strictly decreases with depth. As regards the evolution in time, the porosity significantly decreases due to the expulsion of pore fluid, which is associated with dissipation of excess pore pressures. The simulation for pco = 1.5 MPa predicts at the basin bottom that the ratio u(z = 0, t1)/u(z = 0, T) is of the order of 1/5. Still, porosity does not tend towards zero as asymptotic value. Development of negative porosities at the bottom of the basin, observed here for the lower value of the initial consolidation pressure pco, is attributed to the hardening rule (13) associated with the Cam–Clay model. More generally, for any value of pco, there exists a basin depth beyond

a

pore pressure (MPa)

80

60

0.2

−0.2

t=T t = 1.75 T t = t∞

0

2000

4000

6000

depth (m) pco = 1.15 MPa 0.8

0.6

porosity ϕ

t = 0.5 T t=T t = 2T t = 3.25T t = t∞ (hydrostatic profile)

0.4

0

b 100

0.8

0.6

porosity ϕ

regime is recovered after t = t1, which is of the order of 10 · T for both pco = 1.15 and pco = 1.5 MPa. This last phase, corresponding to the horizontal part of the curves, constitutes the long term basin response. We shall refer to this phase as the stabilized basin configuration. The overall basin deformation may be quantified by the compaction ratio at time t defined as 1  H ðtÞ=H. The magnitudes of this ratio at times T and t1 are: 8.5% and 50% for pco = 1.15 MPa, 7.5% and 42% for pco = 1.5 MPa. It is observed that the large strains analysis is therefore clearly necessary. The prediction of the pore overpressures induced by sediment deposit is of great practical importance in oil drilling operations. Actually, if sedimentation is rapid or the permeability is low, t1 is much greater than T and compaction turns out to be slow, then high pore pressures may subsist and affect drilling operations. In contrast, very permeable sediments or slow sedimentation rate correspond to the situation t1/T is around unity. Compaction is thus said to be fast since the rapid pore water expulsion allows the pore water pressure to return to the hydrostatic value. In so far as mechanical compaction is the only mechanism considered in the present analysis, the above discussion may qualitatively be achieved by examining the dimensionless s qf Þg s parameter k0 ðq _ d =q0 , q being the solid phase density (see M for instance [20] or [30] who introduced a similar parameter involving the ratio between the diffusion coefficient and _ d =q0 ). Values lower (resp. larger) than unity characterize M slow (resp. fast) compaction. The situation considered here manifestly falls within the category of moderate slow compactions. Instead of this global parameter, a local one (i.e., function of depth) may also be introduced to characterize overpressure build-up [46] giving indications about seals localization. The fluid pressure profiles predicted numerically for model M1 with pco = 1.5 MPa, are reported in Fig. 7. During sediments deposition (i.e., t < T), p is a time increasing

325

0.4

t=T t = 1.75 T t = t∞

0.2 40

0 20

0

2000

4000

6000

depth (m) Fig. 7. Evolution of the fluid pressure in the basin (model M1, pco = 1.5 MPa).

0

2000

4000

6000

depth (m) pco = 1.15 MPa Fig. 8. Porosity profiles corresponding to model M1 at different instants.

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

which negative porosity will appear. This drawback should obviously be avoided in basin simulators by means of appropriate modification of the hardening law. In short, the uncoupled model is clearly not satisfactory since it may predict negative porosities within the basin. Besides, it has been shown in Deude´ et al. [17], that disregarding the elastic–plastic coupling leads to erroneous prediction in terms of horizontal stresses. 4.2. Comparison of coupled models M2 and M3

a 100 t = 0.5 T t=T t=2T t = 3.25 T t = t∞,2

80

pore pressure (MPa)

326

60

40

20

0

2000

4000

6000

depth (m)

b 100

t = 0.5 T t=T t=2T t = 3.25 T

80

pore pressure (MPa)

With respect to the previous model, the evolution of (poro)elastic and plastic mechanical properties with porosity changes are now taken into account through (6), (7) and (15). The value of the permeability is maintained constant during the simulation in model M2, whereas model M3 accounts for the decrease of the latter associated with the reduction of porosity (Kozeny–Carman formula (19)). In fact, the comparison of the results derived from model M1 to those of models M2 and M3 turns to be uneasy. This is due to the fact that these models resort to strongly different hardening rules. Indeed, when pore closure is reached (i.e., Jp ! 1  u0), pc tends towards infinity in models M2 and M3, while it takes a finite value in model M1. The dataset used for the simulations is that fixed at the beginning of Section 4. The compaction laws predicted by models M2 and M3 are depicted in Fig. 9. During the phase of sediment deposit (t 6 T), the two models predict quasi identical evolutions of the basin thickness. The differences arise in the consolidation phase (t > T). As expected, the compaction level at any given time t > T is smaller if the permeability decreases with porosity (model M3), since the dissipation of excess pore pressure is slower. This is illustrated in Fig. 10, which show that the pressure profiles do not change at the same time scale. The asymptotic configuration of the basin is reached after t1,2/T . 9 (resp. t1,3/T > 100) for model M2 (resp. model M3). It is

60

40

20

0

2000

4000

6000

depth (m) Fig. 10. Evolution of pore pressure along the basin for models M2 and M3. (a) Model M2 and (b) Model M3.

observed that the long term configuration is almost the same for the two models, in terms of basin thickness and porosity profiles as showed in Fig. 11. Actually, there

6000 0.8

model M2 model M3

0.6

porosity ϕ

H(t) (m)

4000

2000

Model M2, t = T Model M2, t = t∞,2 Model M2, t = T Model M2, t = t∞,3

0.4

0.2

0

0T

500

1000

1500

2000

t (million years) Fig. 9. Compaction laws for coupled models M2 and M3, T . 60 million years.

0

0

2000

4000

depth (m) Fig. 11. Porosity profiles for models M2 and M3.

6000

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

327

100

6000

k o =10 −10 MPa−1 × m2 × s−1, t = T k o =10 −10 MPa−1 × m2 × s−1, t = T pore pressure (MPa)

80

H(t) (m)

4000

2000

60

40

k o =10 −10 MPa−1 × m2 × s−1 k o =10 −9 MPa−1 × m2 × s−1 0

20 0 T

500

1000

1500

0

2000

4000

Fig. 12. Compaction curves (model M3), initial permeability influence, T . 60 million years.

Fig. 14. Pore pressure profiles at t = T (model M3), initial permeability influence.

are significant differences between the predictions of models M2 and M3 as long as the overpressures have not dissipated. In contrast, as regards the long term response of the basin, the evolution of the permeability with porosity may be disregarded.

Within the one-dimensional setting, complementary results derived from simulations implementing the fully coupled model M3 are presented in this section. First, results obtained with a higher initial permeability k0 = 109 MPa1 m2 s1 are represented in Figs. 12–14. These figures illustrate the significant influence of the initial permeability on the compaction curve and on the porosity and pressure profiles. The higher the permeability, the lower the overpressure and the denser the saturated porous rock is.

·

0

2000

4000

6000

depth (m) Fig. 13. Porosity profiles at t = T and t = 2T (model M3), initial permeability influence.

/2

2000

1000

0

0 T ref

200

400

600

800

t (million years) Fig. 15. Compaction curves (model M3), accretion rate influence, Tref . 60 million years.

0.6

0

· M dref

3000

0.6

0.2

·

·

M d = 2 M dref Md =

0.8

=10 −10 MPa−1 × m2 × s−1, t = T =10 −10 MPa−1 × m2 × s−1, t =2T =10 −9 MPa−1 × m2 × s−1, t = T =10 −9 MPa−1 × m2 × s−1, t=2T

·

4000

0.8

0.4

·

M d = M dref

5000

porosity ϕ

porosity ϕ

6000

H(t) (m)

4.3. Further simulations of model M3

ko ko ko ko

6000

depth (m)

t (million years)

0.4 ·

·

M d = M dref, t = T ref / 2 ·

0

·

M d = 2 M dref, t = T ref · · M d = M dref / 2 , t = 2 T ref

0.2

0

1000

2000

3000

4000

5000

depth (m) Fig. 16. Porosity profiles at the end of sediments deposition (model M3), accretion rate influence.

328

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329 100 ·

·

·

· M dref, · M dref

M d = 2 M dref, t = T ref / 2 Md =

pore pressure (MPa)

80

·

Md =

t = T ref

/ 2 , t = 2 T ref

60

40

20

0

1000

2000

3000

4000

5000

depth (m) Fig. 17. Pore pressure profiles at the end of sediment deposit (model M3), accretion rate influence.

The initial permeability is now k0 = 109 MPa1 m2 s1. _ ref ¼ 2:05109 kg=s per unit area, numerical Assuming M d _d ¼M _ ref =2 and calculations are performed with M d _ d ¼ 2M _ ref in order to assess the influence of the accretion M d rate on the global response of the basin. Besides, in order to ensure the same value for H ¼ 6000, the duration of the sediment deposit corresponds to 2Tref and Tref/2, respectively. Results are given in Figs. 15–17. _ d affects the basin response It appears that the value of M essentially during the period of sediment deposit t 6 T. The total amount of sediment supply H being the same for the three simulations, this result may be explained by the fact _ d , the permeability, does not intervene that the parameter M explicitly in the equations governing the hydromechanical evolution of the basin during the consolidation phase t > T. 5. Conclusions A constitutive model has been formulated for sedimentary material. It aims at incorporating some of the fundamentals coupled phenomena involved in the process of basin compaction by means of a simplified micro-to-macro framework. The magnitude of porosity changes induced by compaction imposes to adopt the framework of finite poroplasticity. The derived model is able to account for the stiffness increase and for the hardening induced by large plastic strains through the corresponding variation of the pore volume. An important feature of this model is that the expression of the hardening law avoids development of negative porosities that are encountered with classical models of the Cam–Clay type. From a numerical point of view, a finite element formulation using the updated Lagrangian scheme is implemented in order to analyze the time dependent large deformation of the poroplastic basin. It is emphasized that the geometrical nonlinearity induced by large strains does not allow to neglect the geometry change between the

two configurations reached by the system at t and t + Dt, respectively. The numerical analysis takes into account the variations of the mechanical and hydraulic properties of the material with the porosity change during compaction. A numerical approach to the process of sediment accretion has been proposed: its principle consists in replacing the open material system (i.e., the basin) into a fictitious closed one, and to model the phases of sediments deposit-erosion by activating–deactivating the sub-layers. These theoretical and numerical developments have been implemented in a finite element code. Several numerical basin simulations have been performed in a one-dimensional setting, illustrating the ability of the approach to account for the different couplings (hydromechanical coupling, elastic–plastic coupling, hardening-porosity change coupling) involved in such a problem. To finish with, let us review some issues that still need to be addressed:  The validation of the model by comparison of numerical results with available geological data remains to be done.  Chemical compaction due to pressure solution phenomena should be taken into account in the constitutive equations. At the macroscopic scale, this phenomenon is often addressed through additional viscous terms. A micromechanics-based approach is suitable for this purpose [32].  The anisotropy induced by compaction has not been considered here. In fact, not only porosity changes but also pore-shape evolution affect the elastic and plastic properties of the material. An improved micromechanics-based model taking into account the effects of the microstructural changes in a more comprehensive way is the object of current research.  There is a need for more efficient numerical techniques in order to deal with larger degrees of freedom, particularly for three-dimensional simulations. Conceptually, the numerical tool under its current form can be adopted to simulate basins undergoing three-dimensional evolutions of the type induced by a tectonic loading. However, the increase of the number of degrees of freedom then requires to develop more efficient numerical strategies including remeshing procedures and parallel computing. References [1] Audet DM, Fowler AC. A mathematical model for compaction in sedimentary basins. Geophys J Int 1992;110:577–90. [2] Audet DM. Compaction and overpressure in Pleistocene sediments on the Louisiana Shelf, Gulf of Mexico. Mar Petrol Geol 1996;13:467–74. [3] Barnichon JD. Finite element modelling in structural and petroleum geology. Doctoral thesis, Universite´ de Lie`ge, 1998. [4] Barnichon JD, Charlier R. Finite element modelling of the competition between shear bands in the early stages of thrusting: strain localisation analysis and constitutive law influence. In: Buchanan PG,

D. Bernaud et al. / Computers and Geotechnics 33 (2006) 316–329

[5]

[6]

[7] [8]

[9]

[10]

[11]

[12]

[13]

[14] [15]

[16]

[17]

[18] [19] [20] [21] [22]

[23]

[24]

[25]

Nieuwland DA, editors. Modern developments in structural interpretation, validation and modelling, vol. 99. Geol. Soc. London Spec. Publ.. p. 235–50. Barnichon JD, Havenith H, Hoffer B, Charlier R, Jongmans D, Duchesne JC. The deformation of the Egersund-Ogna anorthosite massif, south Norway: finite-element modelling of diapirism. Tectonophysics 1999;303:109–30. Barthe´le´my JF, Dormieux L, Maghous S. Micromechanical approach to the modelling of compaction at large strains. Comput Geotech 2003;30:321–38. Bathe KJ. Finite element procedures in engineering analysis. Englewood Cliffs: Prentice-Hall; 1996. Bernaud D, Deude´ V, Dormieux L, Maghous S, Schmitt DP. Evolution of elastic properties in finite poroplasticity and finite element analysis. Int J Numer Anal Meth Geomech 2002;26(9):845–71. Bernaud D, de Buhan P, Maghous S. Numerical simulation of the convergence of a bolt-supported tunnel through a homogenization method. Int J Numer Anal Meth Geomech 1995;19: 267–88. Bethke MC, Corbet TF. Linear and nonlinear solutions for onedimensional compaction flow in sedimentary basins. Water Resource Res 1988;24(3):461–7. Birchwood RA, Turcotte DL. A unified approach to geopressuring, low-permeability zone formation, and secondary porosity generation in sedimentary basins. J Geophys Res 1994;99(B10): 20.051–8. Bourgeois E, de Buhan P, Dormieux L. Derivation of an elastoplastic constitutive law for saturated porous media at finite strains. CR Acad Sci, Paris 1995;321(IIb):175–82. Bourgeois E, Dormieux L. Prise en compte des non-line´arite´s ge´ome´triques dans la mode´lisation de la compaction de se´diments. Oil Gas Sci Technol – Rev l’IFP 1997;52(1):23–34. Charlez PA. Rock mechanics, vol. 2, Petroleum applications. Technip; 1997. Charlier R. Approche unifie´e de quelques proble`mes non line´aires de me´canique des milieux continus par la me´thode des e´le´ments finis. Doctoral thesis, Universite´ de Lie`ge, 1987. Connolly JAD, Podladchikov YY. Temperature-dependent viscoelastic compaction and compartmentalization in sedimentary basins. Tectonophysics 2000;324:137–68. Deude´ D, Dormieux L, Maghous S, Barthe´le´my JF, Bernaud D. Compaction process in sedimentary basins: the role of stiffness increase and hardening induced by large plastic strains. Int J Numer Anal Meth Geomech 2004;28:1279–303. Dormieux L, Maghous S. Poroelasticity and poroplasticity at large strains. Oil Gas Sci Technol – Rev l’IFP 1999;54(6):773–84. Dormieux L, Maghous S. Evolution of elastic properties in finite poroplasticity. CR Acad Sci, Paris 2000;328(IIb):593–600. Fowler AC, Yang XS. Fast and slow compaction in sedimentary basins. SIAM J Appl Math 1998;59(1):365–85. Fowler AC, Yang XS. Pressure solution and viscous compaction in sedimentary basins. Geophys Res 1999;104(B6):12.989–97. Garven G. A hydrogeologic model for the formation of the giant oil sands deposits of the western Canada sedimentary basin. Am J Sci 1989;289:106–66. Gibson RE, England GL, Hussey MJL. The theory of onedimensional consolidation of saturated clays, I. finite non-linear consolidation of thin homogeneous layers. Geotechnique 1967;17: 261–73. Gibson RE, Schiffman RI, Cargill KW. The theory of onedimensional consolidation of saturated clays, II. finite non-linear consolidation of thick homogeneous layers. Can Geotech J 1981;18: 280–93. Gue´guen Y, Dormieux L, Boute´ca M. Fundamentals of poromechanics. In: Gue´gyen Y, Boute´ca M, editors. Mechanics of fluidsaturated rocks. Elsevier; 2004.

329

[26] Hamilton EL. Thickness and consolidation of deep-sea sediments. Bull Geol Soc America 1959;70:1399–424. [27] Hanafy EA, Emery JJ. Advancing face simulation of tunnel excavations and lining placement. In: 13th Canadian rock mechanics symposium: underground rock engineering, 28–29 May 1980, Toronto, Canada, 1980. p. 119–25. [28] Hermanrud C. Basin modelling technique: an overview. In: Dore´ AG et al., editors. Basin modelling: advances and applications. Norwegian petroleum society special publication 3. Elsevier; 1993. [29] Hubert MK, Rubey WW. Role of fluid pressure in mechanics of overthrust faulting, I. mechanics of fluid-filled porous solids and its application to overthrust faulting. Bull Geol Soc Am 1959;70: 115–66. [30] Ku¨mpel HJ, Boute´ca M, Dormieux L. Fluid transport in deforming rocks. In: Gue´gyen Y, Boute´ca M, editors. Mechanics of fluidsaturated rocks. Elsevier; 2004. [31] L’heureux I, Fowler AD. A simple model of flow patterns in overpressured sedimentary basins with heat transport and fracturing. J Geophys Res 2000;105(B10):23.741–52. [32] Lemarchand E, Dormieux L, Ulm FJ. A micromechanical analysis of dissolution processes in porous materials. In: Proceedings of the IUTAM symposium on physicochemical and electromechanical interactions in porous media, 18–23 May 2003, Kerkrade, The Netherlands. [33] Lerche I. Basin analysis: quantitative methods, vol. 1. New York: Academic Press; 1990. pp. 208–245. [34] Meroi EA, Schrefler B, Zienkiewicz O. Large strain static and dynamic semisaturated soil behaviour. Int J Numer Anal Meth Geomech 1995;19:81–106. [35] Muir Wood DM. Soil behavior and critical state soil mechanics. Cambridge University Press; 1990. [36] Person M, Raffensperger J, Ge S, Garven G. Basin-scale hydrogeologic modeling. Rev Geophys 1996;34(1):61–87. [37] Schmidt V, Mc Donald DA. The role of secondary porosity in the course of sandstone diagenesis. In: Scholle PA, Schluger PR, editors. Aspects of diagenesis, Society of economic paleontologists and mineralogists special publication 26, 1979. p. 175–207. [38] Schneider F, Potdevin JL, Wolf S, Faille I. Mode`le de compaction e´lastoplastique et viscoplastique pour simulateurs de bassins se´dimentaires. Oil Gas Sci Technol – Rev l’IFP 1994;49(2): 141–8. [39] Schneider F, Potdevin JL, Wolf S, Faille I. Mechanical and chemical compaction model for sedimentary basin simulators. Tectonophysics 1996;263:307–17. [40] Sharp JM, Domenico PA. Energy transport in thick sequences of compacting sediment. Geol Soc Am Bull 1976;87:390–400. [41] Smith JE. The dynamics of shale compaction and evolution of fluid pressures. Math Geol 1971;3(3):239–63. [42] Suetnova E, Vasseur G. 1-D Modelling rock compaction in sedimentary basins using a visco-elastic rheology. Earth Planet Sci Lett 2000;178(9):373–83. [43] Tuncay K, Park A, Ortoleva P. Sedimentary basin deformation: an incremental stress approach. Tectonophysics 2000;323: 77–104. [44] Ungerer P, Burrus J, Doligez PY, Bessis F. Basin evaluation by integrated two-dimensional modeling of heat transfer, fluid flow, hydrocarbon generation and migration. Bull Am Assoc Petrol Geol 1990;74(3):309–35. [45] Wangen M. Pressure and temperature evolution in sedimentary basins. Geophys J Int 1992;110:601–13. [46] Wangen M. A quantitative comparison of some mechanisms generating overpressure in sedimentary basins. Tectonophysics 2001;334:211–34. [47] Yang XS. Modeling mineral reaction in compacting sedimentary basins. Geophys Res Lett 2000;27(9):1307–10. [48] Zaoui A. Continuum micromechanics: survey. J Eng Mech 2002;128(8):808–16.