NOE0415380416 ch052
NOE0415380416 ch052
NOE0415380416 ch052
sedimentary basins
ABSTRACT: The aim of the paper is to provide new elements concerning the constitutive behavior of sed-
imentary rocks and the numerical aspects for basin simulators. The mechanical model is developed within
the framework of finite poroplasticity. Special attention is paid to the effects of large porosity changes on the
poromechanical properties of the sediment. A simplified micromechanics-based approach is used to account for
the stiffness increase and hardening induced by large plastic strains. The numerical method is inspired from the
deactivation-reactivation method used for the simulation of excavation processes. Periods of sediments accretion
are simulated by progressive activation of the gravity forces within a fictitious closed system. The influence of per-
meability and that of the accretion ratio on the compaction law, the porosity and pore pressure profiles are exam-
ined. In particular, the importance to take into account the coupling between porosity changes and permeability
is emphasized. Cam-Clay inspired hardening laws are shown to yield negative porosities at high compaction
levels. These limitation are solved by means of the hardening law derived from the micromechanical reasoning.
339
Copyright 2005 Taylor & Francis Group plc, London, UK
is briefly presented; (3) some examples of numerical elementary volume, defined as the ratio between its
simulations are given. 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 configura-
2 CONSTITUTIVE EQUATIONS tion of the elementary volume. Similarly, the plastic
part J p of the jacobian is defined as the jacobian in
At the scale of the pores and solid grains, large strains this unloaded configuration. Note that the rates J and
J are given by
p
modify the microstructure of the sediment material.
In turn, these microstructural changes are responsi-
ble for the evolution of elastic and plastic mechanical
properties.
The sedimentary rock is modelled as a fully sat- With these notations, the rate form of the second state
urated poroelastoplastic material undergoing large equation reads (Bernaud et al. 2002):
strains. The anisotropy of the mechanical properties
of the sedimentary material induced by compaction is
disregarded. The elastic part of the deformation gra-
dient of the skeleton particles is assumed to remain
infinitesimal. This means that large strains involved
during compaction process are of irreversible (plastic)
nature. M is the Biot modulus. The terms involving M and b
In the framework of finite poroplasticity, the con- in (3) are related to the influence of large plastic strains
stitutive behavior comprises two state equations in on the poroelastic properties.
rate-type formulation together with complementary The complementary equations prescribe the plastic
relations specifying the flow rule (Dormieux and flow rule. Generalizing the concept of plastic poten-
Maghous 2000) (Bernaud et al. 2002). Let denote tial, we therefore introduce the function g(, p) which
the Cauchy stress tensor and p the pore pressure. The derivatives yield the sought plastic rates :
first state equation relates the stress rate , the pore
pressure rate p , and the strain rate tensor d defined as
the symmetric part of the velocity gradient.
rotation (spin) rate tensor which aims at taking the Owing to (2), it is further obtained that = J trd .
p p p
large rotation of the elementary volume into account. Combining this equation with (4) reveals that the plas-
This equation involves: tic potential depends on and through the so-called
a rotational time derivative of the Biot effective Terzaghi effective stress
= + p1 (Bourgeois et al.
stress tensor
e = + b p 1, where b is the Biot 1995):
coefficient; .
a term related to the particulate derivative c of the
tensor of drained elastic moduli c . The evolution
of the elastic properties is related to those of the The relevancy of the above model depends on the
microstructure which cannot be neglected in the capability to specify:
domain of large strains.
1. the elasticity-plasticity coupling characterized by
Relationship (1) extends classical rate form formu- the dependence of tensor c as well as of poroelastic
lations to the case of variable elastic properties.
The second state equation relates the pore volume b coefficients M on large plastic strains;
change to the rate p of the pore pressure and to the 2. the influence of large plastic strains on the evolution
strain rate d. In view of its derivation, it is conve- of the yield surface as well as on the hardening rule.
nient to introduce normalize the pore volume dpore In the case of plastic incompressibility of the solid
the so-called lagrangian porosity = dpore /do that matrix, large plastic strains of the porous material are
is related to the classical porosity (defined as the expected to induce significant porosity and pore shape
pore volume fraction in the current configuration) by changes. In order to capture the influence of the plas-
= J . J is the jacobian of the transformation of the tic strains on the poroelastic properties, the idea is to
340
Copyright 2005 Taylor & Francis Group plc, London, UK
resort to micromechanical estimates of the latter. For A simple way to generalize the classical hardening
the sake of simplicity, the anisotropy induced during law in the Cam-Clay model consists in adopting
the loading process is disregarded in this study. The
pore space is therefore entirely characterized by its vol-
ume fraction, namely the porosity . We herein adopt
the Hashin-Shtrikman upper bounds which are known However, such a hardening law does not enforce the
to reasonably model the elastic properties of isotropic condition J p > 1 o , corresponding to total pore clo-
porous media. Accordingly, the bulk K and shear sure (see (8)). It is therefore expected that a simple
moduli of the porous medium now appear as func- extension of the classical Cam-Clay hardening law in
tions of the porosity as well as of the elastic properties the form (10) might yield negative porosities under
of the solid phase k s and s (that are assumed to be high isotropic compression. That is of course not sat-
constant): isfactory. Instead of (10), we shall hereafter resort to a
micromechanical approach to the hardening law which
overcomes the above difficulty. The idea is to identify
pc (J p ) 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
It is recalled that the Biot coefficient and modulus are medium. Considering the case of a solid phase of the
connected to K through Tresca type, the limit load appears to be a function of
the current porosity (Deud et al. 2004):
341
Copyright 2005 Taylor & Francis Group plc, London, UK
and may be easily incorporated within existing finite fluid mass exchange associated with the pore pres-
element tools. sure dissipation). The system is termed here as closed
For simplicity, the basin undergoing compaction because the corresponding finite element grid com-
is modeled as an horizontal infinite layer. The base- prises a constant number of finite elements and nodes.
ment rock is located at z = 0. Assuming it remains The geometrical domain o is regarded as the initial
horizontal, the upper boundary is defined by the time configuration of the fictitious system material whose
dependent plane equation z = H (t), where H (t) refers constituent is defined as follows: mass density = p f ,
to the total thickness of the basin at time t. Besides, permeability = |k| ko , poroelastic and poroplastic
the sea level is located at the plane z = Lo H (t). In parameters equal to those of the sediment material in
absence of tectonic activity, the gravity forces g ez its reference state. The initial stress and fluid pressure
and the fluid pressure at the top of the basin constitute fields correspond to hydrostatic state. The spatial dis-
the loading parameters in the compaction process. The cretization of o is carried out as follows. The time
accretion rate M d (t) is considered as given. The bound- interval of study is subdivided into n sub-intervals
ary conditions consist in prescribing the values of the [ti1 , ti ], with to = 0 and tn = T . During the time incre-
stress vector T = ez , displacement vector , fluid ment ti = ti ti1 , the dump of sediments would
pressure p and fluid flux q ez correspond (if the material were rigid) to the thickness
increase Hi of the basin. The geometrical domain
z = H (t) (upper surface) o is subdivided into n layers of thickness Hi , each
one being discretized using for instance triangular ele-
ments (6 nodes for the displacement approximation, 3
z = 0 (basin basement) nodes for the pressure approximation).
The issue consists in determining the evolution of
the fictitious closed system material as the sediment
The key idea of the numerical technique consists in is supplied and the layers (F.E. grid) deform with the
transforming the real open material system (the basin) compaction process. As mentioned above, the sedi-
into a fictitious closed system. The evolution in time ment layers are activated as they added by setting their
of the latter must of course correspond to that of the physical parameters to those of the real material in its
real system. The main advantage of this procedure lies reference state.
j j
in the fact that the finite element approach such as that The layer i is defined by Hi1 z Hi , where
j
developed in (Bernaud et al. 2002) is applicable. z = Hi refers to the current (i.e. at time tj ) position of
The principle is inspired from the so-called deacti- the material plane located initially at z = Hi . The geo-
vation/ reactivationmethod (Hanafy and Emery 1980) metrical volume occupied by the (fictitious) porous
(Bernaud et al. 1995). This technique was initially material at time tj is denoted as j :
developed and implemented for numerical simulation
in tunnel engineering (simulation of excavation pro-
cess and lining placement). The idea consists here
in activating sediment layers as they are effectively
deposited by substituting their fluid characteristics If i j, the layer i has already been activated and
with material sediment ones. This process amounts to the mechanical properties are those of the sediment
a stepwise activation of the gravity forces in the upper material in the current state, since they evolved with
layers. compaction from time ti of activation.
More precisely, let us assume that the simulation In contrast, if i > j, the layer i is not activated yet
of the basin formation is carried out during the time and the mechanical properties correspond to those
interval t [0, T ]. We hence introduce the geometrical of the initial fictitious configuration.
domain o defined by 0 z H, where
The local fields at time tj are the solutions to successive
finite poroplastic boundary value problems performed
during time intervals [ti1 , ti ].
342
Copyright 2005 Taylor & Francis Group plc, London, UK
6000 0.8
0.6
4000
porosity
H(t) (m)
0.4
2000 (a1)
(a) (a2)
(b) 0.2 (b1)
(b2)
0 0
0 500 1000 1500
0 2000 4000 6000
t (million years) depth (m)
Figure 1. Compaction curves. (a) ko = 1010 MPa1 m2 s1 ; Figure 2. Porosity profiles. (a1 ) ko = 1010 MPa1 m2 s1 ,
(b) ko = 109 MPa1 m2 s1 . t = T ; (a2 ) ko = 1010 MPa1 m2 s1 , t = 2T ; (b1 ) ko =
109 MPa1 m2 s1 , t = T ; (b2 ) ko = 109 MPa1 m2 s1 ,
Each layer is subdivided into finite elements. The t = 2T .
time increment tj corresponding to the activa-
tion (supply) of layer j is itself subdivided into 100
time sub-increment in order to better capture the
hydromecanical evolution of the basin. (a)
The numerical procedure has been described assum- (b)
80
ing that the layers remain horizontal (i.e. 1D com-
pore pressure (MPa)
343
Copyright 2005 Taylor & Francis Group plc, London, UK
process of basin compaction by means of a simpli- Bernaud, D., P. de Buhan, and S. Maghous (1995). Numerical
fied micro-to-macro framework. The magnitude of simulation of the convergence of a bolt-supported tunnel
porosity changes induced by compaction imposes through a homogenization method. Int. J. Numer. Anal.
to adopt the framework of finite poroplasticity. The Meth. Geomech.19, 267288.
Bernaud, D., V. Deud, L. Dormieux, S. Maghous, and
derived model is able to account for the stiffness D. Schmitt (2002). Evolution of elastic properties in finite
increase and for the hardening induced by large plastic poroplasticity and finite element analysis. Int. J. Numer.
strains through the corresponding variation of the pore Anal. Meth. Geomech. 26(9), 845871.
volume. An important feature of this model is that the Bourgeois, E., P. de Buhan, and L. Dormieux (1995). Deriva-
expression of the hardening law avoids development of tion of an elastoplastic constitutive law for saturated
negative porosities that are encountered with classical porous media at finite strains. Comptes Rendus Acad. Sc.,
models of the Cam-Clay type. Paris 321(IIb), 175182.
From a numerical point of view, a finite element Bourgeois, E. and L. Dormieux (1997). Prise en compte
formulation using the updated Lagrangian scheme is des non-linarits gomtriques dans la modlisation de
la compaction de sdiments. Oil & Gas Science and
implemented in order to analyze the time dependent Technology 52(1), 2334.
large deformation of the poroplastic basin. It is empha- Deud, D., L. Dormieux, S. Maghous, J. Barthlmy, and
sized that the geometrical non linearity induced by D. Bernaud (2004). Compaction process in sedimen-
large strains does not allow to neglect the geometry tary basins: the role of stiffness increasr and hardening
change between the two configurations reached by induced by large plastic strains. Int. J. Numer. Anal. Meth.
the system at t and t + t respectively. The numer- Geomech. 28, 12791303.
ical analysis takes into account the variations of the Dormieux, L. and S. Maghous (2000). Evolution of elastic
mechanical and hydraulic properties of the material properties in finite poroplasticity. Comptes Rendus Acad.
with the porosity change during compaction. A numer- Sc., Paris 328(IIb), 593600.
Hanafy, E. and J. Emery (1980). Advancing face simula-
ical approach to the process of sediment accretion tion of tunnel excavations and lining placement. 13th
has been proposed: its principle consists in activating- Canad. Rock Mech. Symp.: Underground rock engineer-
deactivating the sub-layers. ing, Toronto, Canada.
These theoretical and numerical developments have Hubert, M. and W. Rubey (1959). Role of fluid pressure in
been implemented in a finite element code. Several mechanics of overthrust faulting, i. mechanics of fluid-
numerical basin simulations have been performed in filled porous solids and its application to overthrust
a one-dimensional setting, illustrating the ability of faulting. Bulletin of Geological Society of America 70,
the approach to account for the different couplings 115166.
(hydromechanical coupling, elastic-plastic coupling, Smith, J. (1971). The dynamics of shale compaction and
evolution of fluid pressures. Mathematical Geology 3(3),
hardening-porosity change coupling) involved in such 239263.
a problem. Tuncay, K., A. Park, and P. Ortoleva (2000). Sedimen-
To finish with, let us review some issues that still tary basin deformation: an incremental stress approach.
need to be addressed: 1) the chemical compaction due Tectonophysics 323, 77104.
to pressure solution phenomena should be taken into Wood, D. (1990). Soil behavior and critical state soil mechan-
account in the constitutive equations ; 2) the anisotropy ics. Cambridge University Press.
induced by compaction ; 3) more efficient numeri-
cal techniques is needed in order to deal with larger
degrees of freedom, particularly for three-dimensional
simulations.
REFERENCES
344
Copyright 2005 Taylor & Francis Group plc, London, UK