KINE3D: a New 3d Restoration Method Based on a Mixed Approach Linking Geometry and Geomechanics

Résumé — Une nouvelle méthode de restauration 3D basée sur une approche couplée géométrie-mécanique — Une nouvelle méthode de restauration 3D a été développée en couplant un logiciel de modélisation géométrique avec un code d’élément fini de mécanique. Cette méthode permet à l’utilisa-teur d’imposer le déplacement sur le plan de failles afin d’obtenir la géométrie initiale désirée tout en cal-culant une déformation continue dans les blocs par le programme de mécanique. Le maillage du bloc 3D permet de tenir compte des éventuelles hétérogénéités induites par les couches géologiques et l’évolution latérale des faciès. Les résultats sont discutés dans des cas extensifs et compressifs; en rétro-déformation ainsi qu’en déformation directe. Le cas compressif correspond à la restauration d’un anticlinal faillé constitué de bancs massifs, sable et argile en alternance. Il a aussi été restauré en 2D, via une mise à plat surfacique, et les conclusions des deux restaurations 2 et 3D sont comparées. Le second cas, une modéli-sation directe, est un cas de glissement gravitaire


INTRODUCTION
2D and 3D geological restorations are classically achieved through a geometrical approach.The rheology of the natural material is approximated by deformation modes, namely essentially flexural slip and simple shear.These two modes allow to represent the bulk behaviour of sedimentary rocks in addition to the rigid block rotation (domino style), the flow (for the salt and some shales) and obviously the mass evolution due to compaction and layer parallel shortening.Under the flexural slip hypothesis, the main rock discontinuities, and therefore the main shear zones, are the lithological bed interfaces (Suppe 1983).At the opposite, with the simple shear mode, the material acts as a granular medium (such as unconsolidated sand) and the shear zones are planar and correspond to the shear angle of the material (White et al. 1986).This global approach is the base of restoration tools such as LOCACE (Moretti and Larrère 1989), GEOSEC (Rowan and Kligfield 1989) and 2D/3D MOVE.These methods can easily be incorporated in forward as well as in backward modelling approaches.It has been for instance done in THRUSPACK (Divies 1997).The displacements are defined by giving a known part of the final geometry that could either be the geometry of a layer or the displacement along a fault.Final state in this case refers to the computation: it is the present one or the past one for backward and forward modelling respectively.Usually in forward modeling the target geometry is given at the bottom of the layers, assimilated to the top of the décollement level whereas in backward modeling the target geometry is often given by at the top of the layer, supposedly flat at a given time step.Nevertheless the methods will allow us to implement other choices.Purely geometrical approaches are also used in usual surface unfolding techniques (Bennis et al. 1991, Rouby et al. 2000, Rouby et al. 2002).Such methods are now all based on a parameterization of the surfaces when in the years 90's trial based on triangulated surfaced has been tested (Gratier and Guillier 1993).They happen to be too dependent on the triangulation and search direction of the best fit during the flattening of the triangles.
On the other hand, mechanical approaches allow to compute the deformation given initial conditions: geometry, behaviours, boundary forces and/or displacements with their evolution versus time.In Earth Science, since the boundary forces are usually unknown, the boundary conditions are often imposed displacements and/or thermal conditions.This approach is largely used at large (asthenospheric/ lithospheric) scales to understand the Earth evolution (see for instance Poliakov et al. 1993, or Burov andPoliakov 2001 with the PARAVOZ code), and at a smaller scale to quantify site evolutions (for instance the site of a dump in a valley, the slope stability near a drilling rig, etc).
At the basin scale (between 100 and 5 km length and 2 to 10 km deep) the geomechanical approach is still poorly used for various reasons.First of all, the modeller has to face the major difficulty of postulating a constitutive law for geological period and kilometric scale materials.A second challenge is to determine the initial conditions, especially the preexisting discontinuities that will localize the shear zones.Starting with a homogeneous material, one cannot localize faults at their known positions.On the contrary, introducing weak zones over-specifies the problem, therefore producing results which are only reflecting the hypotheses.Moreover, the friction on the fault planes is a poorly known parameter, laboratory measurements on small scale samples, as well as data deduced from the seismicity at small time scale and large space scale, being not easily extrapolated to the large time scale of a basin evolution.
Nevertheless, the increase of the computer CPU facility and the quality of the seismic image let now the explorationist dream about a full 4D modelling, meaning a complete representation of the evolution of a zone during its geological history with quantification of the deformation versus time.Such an approach would constitute a major step in the frame of basin modelling in complex areas and of the understanding of fractured reservoirs.
Knowledge of the dynamics of the deformation is necessary to predict the fracturation of a deformed reservoir, as the geometry of the drainage area versus time is a required input for the oil maturation/migration software.During the last couple of years, research into deformation in gravity gliding contexts (West African margin, Gulf of Mexico, etc.) have testified that the rheology of the components and the associated facies boundaries played a major role in the stress orientation, and therefore could not be neglected (Calassou andMoretti 2002, Moretti et al. 2003).
Based on this analysis of the problem, we started to develop a 3D-restoration tool based on a mixed approach: geometry and geomechanics.The geometry at various time steps, as well as the major fault offsets, will be defined using the GOCAD geomodeller (Mallet 2001), whereas the restoration within the main "blocks" will be computed by a finiteelement mechanical code.We have chosen Code_Aster, EDF property, which will be briefly described in the following pages before discussing, on examples, the interest of this method.A prototype was developed in 2004 by F. Lepage (F.Lepage et al., 2004), the industrial tools is now developed in the frame of KINE3D, IFP-Earth Decision Sciences joint project.Facing the same problems, other authors have also developed GOCAD plugins based on a dynamic-relaxation solution (Santi et al. 2003;Muron, 2005), or stand-alone algorithms (like DYNEL for example).

Generating a Suitable Simulation Mesh for 3D Deformation Modelling
As previously said, the entire geometric modelling part, and therefor the meshing, is undertaken by the chosen geomodeller, the finite-element geomechanical code being only used to perform the simulation (computation of a "deformed" geometry).This step is crucial for the restoration process and will be so describe briefly here.
Through the use of high-resolution data like seismic cubes or well loggings, geological characterisations often lead to the construction of high-resolution structural models that capture most of the details of the reservoir geometry.From the geometrical point of view, in 3D, these models are a set of discrete representations of the geological objects of dimension 0, 1 and 2, namely faults, horizons, and their borders.These objects and their borders have various topologies: for example, they can be either closed or open.Building such structural models is a crucial research topic and numerous techniques have been proposed so far (one can for example refer to Caumon et al. 2003).In the following, they will be considered as an input.
Theses structural models must be used as constraints for building the 3D mesh dedicated to the finite-element simulation.However, in such models, there is most of the time a poor description of the contacts between surfaces, although cases where faults or fractures are inter-connected and organised in complex networks are quite common.Contacts are just an information, stating that a set of objects are supposed to intersect, and that this intersection is materialised by a set of objects of lower dimension.At this step, contacts have thus no real geometrical or topological meaning, and corresponding structural models can not be considered as valid (Lepage 2003, Prévost et al. 2004).To achieve model consistency, a Soft Frame Model approach can be used, as described in Lepage 2003.This gives us a general frame for describing and remeshing the input data properly, so that it is ready to be used as constraints for the generation of the 3D simulation mesh.
Although the mechanical finite-element code we use could handle various kinds of elements for 3D meshes, only tetrahedral elements have been considered for our studies, because of the flexibility they provide and the relative ease of shape and size control.To keep meshes as coarse as possible and to prevent over-refinements while keeping an overall satisfying shape quality for tetrahedra, boundaries (elements of the constraining structural model) are recovered using a Boundary-Constrained Delaunay mechanism, and the shape of elements is optimised through a modified Delaunay refinement process, combined with an edge-swap and a smoothing technique (Lepage, 2003).It is possible as well to vary the resolution inside the volume of interest.Resulting meshes have been previously successfully tested in the frame of thermal finite-element (Kohl et al. 2003) and flow finite-volume (Lepage 2003, Prévost et al. 2004) simulations, so their suitability for geomechanical applications like 3D restoration would demonstrate that they are a good support for most of the numerical methods solving classical PDEs in Geoscience.
In the particular context of 3D restoration (or forward deformation modelling), fault blocks must be free to move independently from the other ones, so they must be made of their own mesh elements.Unfortunately, the tetrahedralization process described above results in matching elements (apices, edges and faces of tetrahedra) across interfaces.Thus, as a post-processing, nodes are split along main faults and horizons.This can be done efficiently, using topological considerations only (see Caumon et al. 2003 for the analogous 2D problem).This information will also be useful for defining material regions (groups of connected tetrahedra) or setting boundary conditions, where one has to identify nodes lying on a given fault or a given horizon.

1-2 Setting Boundary Conditions on Displacement
As explained in the introduction, the boundary conditions will be based on displacement, not on forces.Like in the 2D restoration process, this displacement field corresponds to the definition of the geometries of selected markers and/or contacts between two interfaces (horizons and/or faults).We selected a solution that works in both forward and backward models.
These boundary conditions on displacement have two major purposes: on one hand, to manage the contacts between the different fault blocks and their relative movements, on the other hand, to make sure that part of the final computed geometry matches a given "target" known geometry, which is in most of the cases the one of an horizon of the model.In the following, for sake of clarity, these two different kinds of boundary conditions will be treated separately.However, it will be shown how they can be simultaneously taken into account during simulations.

Matching a Given Target Geometry (Numerically Known)
This kind of boundary condition consists of imposing known displacement values (in the X, Y and/or Z directions) on some nodes of the mesh, letting the other nodes free to move and accommodate the induced deformation.As previously said, this aims at matching a given target geometry defining part of the final state.Resulting displacement constraints can thus come from various requirements.For example, in the case of backward deformation modeling (restoration), a simple target geometry could be that of the flattened top horizon of a folded structure.An even simpler one could be that of the plane defined by Z = 0 for the same horizon.The geologist may also wish to maintain vertical a given plan, that would correspond to the 3D extension of the 2D pin-line concept.
Clearly the uniqueness of the simulated final state, has to be guaranteed by a mechanism avoiding rigid body movements, especially rotation around a vertical axis; it could be done by giving the final position of at least 2 selected nodes.

Managing Contacts and Relative Movements of Fault Blocks
Interfaces, i.e. surfaces, separating pairs of adjacent fault blocks can be either horizons or faults.Managing node displacements on horizons is mandatory if one considers them as slipping interfaces.In this work, they will be rather considered as sticking interfaces and in that case, relative block movements on horizons must be forbidden.If mesh elements match across horizons and if nodes lying on them are not split, then this requirement will be automatically satisfied.On the contrary (if mesh elements do not match or if nodes lying on horizons are split), specific displacement constraints are needed (see Fig. 1).This is also obviously the case for faults, which must be treated as slipping interfaces.
Various approaches can be used for describing the relative movements of fault blocks: -Given displacement values (direction and norm) on slipping interfaces like faults, one can compute a constraining displacement for all the nodes of the mesh that lie on these surfaces.The main difficulty is to compute such a consistent displacement field, which must take into account the contacts between the interfaces.Moreover, in the case of faults, displacement values are often seen as a result of the modeling, more than as an input.However, this approach can perfectly be used when considering horizons as sticking interfaces (see Fig. 1); -One can impose displacements only on a restricted set of nodes (for example those that lie on the borders of the horizons affected by a given fault, when setting boundary conditions for a fault), and let the other ones move freely, as shown in Figure 2. Note that this does not guarantee to keep fault blocks in contact.In KINE3D, a classical "master/slave" formulation has been adopted: the nodes lying on the "master" side of a given interface are relatively fixed whereas the nodes lying on its "slave" side are mapped to (projected on) the "master" side.In the above first two approaches, this mapping function associates every slave node to a projected location which is an impact point somewhere on the master side: -In the first approach, impact points are located inside the tetrahedra faces defining the master sides of the considered interfaces (see Fig. 1); -In the second approach, impact points are located inside the tetrahedra edges defining the borders of the master sides of the considered interfaces (see Fig. 2).
Impact points can thus be defined by their barycentric coordinates relatively to at most 3 (first approach) or at most 2 (second approach) nodes of a master side.This can easily be extended to the case of branching faults: whenever some of the nodes of the master side of a given fault are part as well of the slave side of another (branching) fault, linear relationships must be combined until the impact point formulation only depends on nodes that are not part of any slave side.
Note that the location of these impact points clearly depend on the choice of the master and slave sides associated with interfaces.The choice of the mapping function associated Example of a simple model made of three isopach layers.Meshes do not match across interfaces between layers.If horizons are considered as sticking interfaces, specific displacement constraints must be set: every yellow-circled node is associated to a projected location somewhere inside a tetrahedra face lying on the bottom of the above layer, which defines a linear combination between the overall displacements of the yellow-circled node and the ones belonging to the tetrahedra face.
with each interfaces has a very strong impact on the final solution as well, as it is directly used for imposing "hard" displacements during stress/strain simulations.Such functions must thus be defined carefully and have a realistic geological meaning.In the particular case where the considered interfaces are faults and where only the nodes lying on the borders of a given horizon are constrained (second approach), both slave and master nodes form 3D piecewise polygonal lines, so a default choice for the mapping function could be to associate a curvilinear abscissa to every slave node and to associate it with an impact defined by the point of same curvilinear abscissa on the master polygonal line (see Fig. 2).
Imposing displacements only on a restricted set of nodes (same as 2), and forcing the other ones to slip on the fault/horizon being considered would be the more realistic solution, as it guarantees to keep fault blocks in contact.The development of this kind of boundary condition is currently under development in IFP with Code-Aster (Bahloul, Master thesis; 2005).We are also working to increase the user freedom on defining the slip on the fault plane.

Combining Boundary Conditions
The two kinds of boundary conditions described above are different: the first one manages relative node movements, whereas the other one manages absolute node movements.Relative node movements are defined by a projection P mapping a node A to a linear combination of the location of some p other reference nodes A i (i in [0, p-1]), with weights a i (see Fig. 3).Let us note u(A) the overall displacement of node A.
If one wants both kinds of constraints to be honoured, the same linear relationship must exist between the final location of A and the final location of "master" nodes A i , so (Eq.1): Let us consider a horizon H and a fault F. Through the projection P, the node A is mapped to the master side of fault F, somewhere between nodes A 0 and A 1 , such that A + P = a.A 0 + (1 -a) • A 1 .At the end of the simulation, the same linear relationship must be found between the final locations of nodes A, A 0 and A 1 , which leads to the following formulation with overall displacements: Setting displacement constraints for cancelling the effect of a fault.Here only the nodes lying on the fault-horizon contact are considered.
Nodes lying on the "master" side of the fault have relatively fixed locations whereas nodes lying on its "slave" side are mapped to the "master" side.Impact points are located somewhere inside tetrahedra edges.Left: the "master" side has "-" polarity.Right: the "master" side has "+" polarity.
The vector P is known and can perfectly be computed before running any simulation.As a result, one can see that for projected nodes, all boundary conditions can be simultaneously taken into account with a linear combination of the overall displacements of a restricted set of nodes, as shown in Figure 3.
Setting displacement constraints can be thus summarised as follows: -Define a "slave" and a "master" side for each of the interface of the model.-Identify the nodes (from the slave sides of the interfaces) which displacement has to be constrained for managing the relative movements of the fault blocks, and compute their projected location ∑ i a i .A i .This defines the projection vector P. -Set the boundary conditions for the contact management, using the above Equation 1. -Set the boundary conditions defining the target geometry, using formulations of type u x (A) = D x and/or u y (A) = D y and/or u z (A) = D z , where D x , D y and D z are known displacement values.-Whenever a node is concerned by both kinds of displacement constraints, only the one coming from the management of contacts must be taken into account.

Mechanical Assumptions
Although rock behaviour could be much more complex, in order to have a reversible computed transformation, an elastic behaviour of the materials has been assumed.For sake of simplicity, we use a hyper-elastic relationship for which the second Piola-Kirchhoff stress tensor S derives from a potential of the Lagrangian strain tensor E. The latter represents the deformation accounting for geometrical non-linearity (i.e.finite rotation).Using the Einstein convention, its components are (Eq.2): where u is the displacement vector and X the coordinate vector in the initial configuration.However, in the following study cases, only the Cauchy stress tensor σ (deduced from S) will be considered.Its three eigenvectors σ 1 , σ 2 and σ 3 , with associated eigenvalues σ 1 , σ 2 and σ 3 , will be such that σ 1 > σ 2 > σ 3 .As a consequence, σ 1 , σ 2 and σ 3 will also be respectively referred to as the maximum (or least tensile), intermediate and minimum (most tensile) principal stress directions; σ 1 , σ 2 and σ 3 will be respectively referred to as the maximum, intermediate and minimum principal stress values.It is to be noted that the geometrical non-linearity of the deformation requires an iterative Newton-Raphson algorithm to find the solution in terms of increment of displacement for a given deformation step.
The prototype of this new software has been developed in IFP in 2004 and tested in various cases studies.In the next paragraph two different cases will be described: a backward restoration in a compressive zone and a forward restoration in extension, more specifically in gravity gliding context.

Geological Setting
In compressive areas, traps are often the hinges of the anticlines.Unfortunately, onshore, in mountain belts, the seismic images may still have a poor resolution.Their processing may be difficult and does not always allow the explorationist to locate precisely even the horizontal position of the top of the structure.In this context, restoration is often the best tool to precise the geometry.In addition, in term of reservoirs, in numerous cases, as for instance in the Andeous front belt, the target is deep and the bulk reservoir porosity rather low.In this case, the prediction of the fracture network and of the resulting secondary porosity could be the key parameter for a field economical development.Usually the fracture network computation is inferred based on the final curvature of the beds.A knowledge of the cumulative strain should clearly help to improve this fracture model and the knowledge of the stress based on the rock rheology would constitute another step (Guiton et al. 2003).
The example studied here is an anticline in the Sub Andean Zone of Bolivia.The reservoirs are the sand deposits of the Palaeozoic, and the sources rocks, interbedded with these reservoirs, are the shale of the same Palaeozoic.The Devonian-Silurian sequence represented in the model is about 1.2 km thick, both sand and shale being quite massive.The compression is late tertiary, from 10 Myr to now, see Moretti et al. 1996or Moretti et al. 2002 for details concerning this area.A study of the fracture network in the low porosity sandstone of the reservoir could be found in Florez-Niño et al., (2005), this study highlights the role of the kinetic of the deformation over the fracture and fault network development and links an estimated shear strain to the fracture density.

Geometry
The 3D geometry of this anticline is based on the synthesis of surface, seismic and well data.It is shown in Figure 4 and 5. Bed length analysis indicates a shortening of about 10% of deformation perpendicularly to the fold axis.The top-bottom relationship, in thickness, has been tested using a 3D version of the kink-method 2D approach (Galera et al., 2003).A2 -Restored geometry, the lines are the isovalues of the parametrization.
A3 -Internal strain: spreaded variation of area.B3 -Internal strain: no variation of area.
B2 -Restored geometry, the fit on the fault lip is not perfect.B1 -Upper horizon, the large points represent the imposed links between the fault lips.

Boundary Conditions
The top horizon of the anticline is flattened (in both 2D and 3D reestorations), imposing a constant vertical coordinate for all the nodes of the mesh lying on the top surface.The horizontal coordinates are let free to accommodate the deformation excepted for two points at the extremities of the right border.The faults of this structure are known to be synfolding and there is no thickness variation across them.In this case, faulting is Tertiary whereas the modeled layers are Paleozoic.Therefore, in 3D, faults could be unfaulted by imposing coincidence of both top and bottom horizons.

Results of the 2D Unfolding
The surfacic restoration is based on parameterization of the surface.The unfolding is achived by imposing a depth value and looking for the best result, imposing a preservation of both lengths and angles (Mallet, 2001).If the surface is unfoldable, there is one solution (one imposes also the final position of 2 points to avoid rotation as in 3D).If the surface is not perfectly unfoldable, the research of the best solution will result on imposing some light internal deformations.Two strategies are possible within the frame of the Gocad research algorithms: one may use one parameterization for the full horizon or one by patch.When using different patches the link is given by the position of two points on the lips of the faults.The figure 4 shows the results for both methods applied at the top of the structure which is an Upper Devonian horizon.The large points on the center of the faults represent the location where the boundary conditions have been imposed and the lines on the surfaces represent the isovalues of the parameterization.By using one parameterization (4A), the fit is perfect on the faults, by construction.By using one parameterization by patch (4B), the fault lips fit only on the two points that has been imposed.Both cases point out a geometrical imperfection on the small block in the eastern front of the anticline since there is no way to have a perfect fit between the lips without imposing internal deformation.The one-parameterization approach has the advantage of being fully automatic and to compute a continuous unfolded surface, an internal deformation is imposed to reach this fit that could, or not, be considered as negligible by the geologist.At the opposite, the second approach (4B) has the advantage to highlight the geometrical problem inherited from the initial geometry.This problem may, in this case, be easy solved by changing slightly the geometry of the faults in the flank of the anticline.Nevertheless, at this stage, the explorationist may consider that his interpretation is good enough to implement the well or to answer to any other question he has to face.

Results of the 3D Restoration
The next step, could be the 3D restoration, the results obtained with KINE3D is shown in Figure 5, with satisfying fitting of the blocks.The volume change is less than 10 -3 .The change in top length perpendicularly to the fold axis is about 2 km.It is comparable with the results of a 2D balancing with the flexural slip hypothesis (Moretti et al. 2002).This illustrates that length conservation along the compression is quite satisfied by our methodology.
Despite the strong approximation of hyper elasticity, the stress field resulting from the simulation can be qualitatively interesting in analysing the bulk distribution of the deformation.To this end, we consider the opposite stresses associated to forward deformation.Figure 6 shows the maximum (least tensile) principal stress magnitude σ 1 .It clearly shows the outerarc extension (in blue) and inner-arc compression (in red) in the hinge of the anticline.Figure 6 shows as well the magnitude and orientation (with arrows) of the maximum (least tensile) and minimum (most tensile) principal stresses σ 1 and σ 3 , on left and right pictures respectively.In the left flank and hinge, the extension is oriented along the bedding and the compression perpendicular to it.This illustrates the dominant influence of unfolding in this area.To the opposite, the right flank is dominated by unfaulting deformation.It is particularly obvious for the nearest block from the hinge, where compression from bottom to the top results from misfit in the top geometry.This latter point was already highlighted in a previous work with surface restoration.Therefor the light geometrical problem at the contact between the horizon and the faults in the eastern flank of the anticline that became apparent in the 2D restoration, results now to a wrong stress and strain tensors.

Geological Setting
During the last couple of year, exploration became very active in turbiditic systems of the deep offshore, especially in the south Atlantic margins (Brazil and West Africa) and within the Gulf of Mexico.The reservoirs as well as the tectonic context of these zones are quite particular.The reservoirs are usually constituted of sandy channels interbedded in a shaly/silty matrix.These channels are now well imaged by the high-resolution 3D seismic image even if some of the features imaged by the seismic data are still debated.One of the open questions is the influence of the post-depositional deformation versus the sedimentological features.In the central West African margin (Lower Congo basin and neighbors), the explored channels Figure 7 Forward deformation of a channel embedded into a shaly matrix.Left: initial mesh with a detail showing the mesh within and around the channel.Right: as a boundary condition, a radial extension of 10% is imposed.are post-extension.Nevertheless, they are clearly affected by faults that influence the continuity of the reservoirs.For more than 100 Myr, there is no tectonic activity in this zone but gravity gliding of the margin occurs due to the increasing sedimentary load and to the tertiary uplift of the African continent (Duval et al. 1992).In this context the deformation is clearly multi-directional and the data show extension both parallel and perpendicular to the margin as emphased by the huge rotations deduced from surface unfolding (Rouby et al. 2000(Rouby et al. , 2002)).Analogue models, designed to represent the spreading of a margin, have highlighted that the development of the fault pattern and the fracturation on the channels are strongly dependent on the sand/shale interfaces (Panien et al. 2001, Moretti et al. 2003).We discussed in these papers how the borders of the channels evolved as faults during the gliding and how the stress re-orientation within the channel leads to the development of a fracture network perpendicular to the channel border.Seismic and field data, as analogue and 2D numerical models, have confirmed this evolution (Moretti et al. 2003).
In Moretti et al. 2003, a 2D-elastic geomechanical model with triangular elements and an infinitesimal strain of 1% (radial extension) was used to study the stress orientation in and around the channels versus the Young's coefficient and Poisson's ratio contrasts, with a 2D elastic relaxation solution.Thanks to the link between GOCAD and Code_Aster, a 3D geometry can now be considered for this model, leading to a more realistic description of the channel-matrix interaction.

Geometry
The geometry of the channel is exactly the one defined in the analogue model presented in Figure 9.The aspect ratio, width and thickness, have been deduced from the channel visible in seismic data, presented in Calassou and Moretti 2003.The Young's modulus E and Poisson's ratio ν for the surrounding shales are respectively set to 85 MPa and 0.20.For the more competent sandy channel, E is set to 100 MPa and ν to 0.20 to analyze the influence of a contrast in Young's modulus.

Boundary Conditions
The study case corresponds to a forward model.A quarter of cylinder represents the part of the margin that is studied, in order to allow a direct comparison with the analogue models.In these models, a wall maintained the silicone in its initial position.This wall was taken off, triggering the gliding of the silicone under its own weight.The silicone is viscous, i.e. with rate-dependent deformation.For rate-independent numerical models with elastic material, we have extended the border by a given value.
A uniform radial extension of 10% has been imposed on the free border of the model as shown by the comparison of the two bottom meshes: blue for the initial and red for the final (see Fig. 7).The two borders of symmetry are vertical and the top and bottom surfaces are imposed to remain flat, thus enforcing plane strain conditions.This last boundary condition does not correspond to natural cases since a subsidence is expected during gravity gliding.Nevertheless, in the numerical experiment presented in Figure 6, there is no gravity and no compaction, therefore, the vertical displacement induced by the gliding is out of the scope of the simulation.

Results
A simple analytical calculation in the context of small deformation (see Annex 2) shows that the radial extension of a homogeneous linear elastic material (i.e.without channel) would result in horizontal equal most tensile and intermediate principal stress values.This means that any two perpendicular horizontal vectors are solutions for the direction of these principal stresses.The numerical solution in large deformation leads to similar conclusions.
Figure 8 shows the computed minimum (most tensile) and intermediate principal stresses σ 3 and σ 2 with scaled arrows.Due to plane strain conditions, the maximum (least tensile) principal stress σ 1 is vertical.This figure also displays a representation of the computed stress field with ellipsoids.Red ellipses are normal to the intermediate principal stress direction σ 2 (thus are defined by σ 1 and σ 3 ) and white ellipses are normal to the minimum principal stress direction σ 3 (thus are defined by σ 1 and σ 2 ).
The results show that there is a complete re-orientation of the principal stresses.Inside the channel, the "least tensile" horizontal stress σ 2 is perpendicular to the edge of the channel and becomes parallel to this edge outside the channel at the sand/shale border.This result is similar to the conclusion of both 2D numerical modeling and analogue experiments (Moretti et al. 2003), as shown in Figure 9.This model suggests the development inside the channel of a fault network perpendicular to the channel boundary, a feature which has been observed in 3D seismic data (Calassou and Moretti 2003).

CONCLUSIONS
The link between Code_Aster and GOCAD that has been presented here let us be very optimistic on the possibility to achieve 3D restoration in complex geological cases and to compute a qualitatively consistent stress field related to the deformation.The possibility to take into account variations in the rock behavior within the main blocks is particularly interesting.Today, the method can still be improved by enlarging the characteristic of the target geometry.The developed methodology is an alternative to the work described in Muron and Mallet 2003, based on energy minimization criteria and the dynamic relaxation approach (Muron 2005;Mueller et al. 2005).
Figure 9 2D analogue experiment showing the fracturation of a channel under gravity gliding.The silicone glides under its own weight, which leads to the development inside the channel of a fault network perpendicular to the channel boundary.
The comparison between the surfacic and the volumetric restoration results may be interpreted on different ways.Restoration has two main goals: -insure the correctness of the interpreted geometry; -compute the past geometries to quantify the deformation and/or the migration pathway evolution at the basin or reservoir scales.
It is obvious that to reach the second goal the volumetric restoration constitute a major step.At the opposite, the incoherence of the geometries may be often shown, as in the compressive case (fig 4 to 6), based on a surface restoration if the hypotheses of this unfolding are constrained enough.In cross-section, a similar analyze could be done when passing from length balancing to full 2D flexural slip restoration.Within KINE3D, various tools, from length balancing to 3D volumetric restoration including cross-section balancing and surface unfolding will be available within the same environment.It will allow the users to solve their problems step by step: complete a 3D meshing on a faulted area to realize that the lengths of the horizons are incoherent is a little bit time consuming.At the opposite full 3D volumetric restoration taking into account the rheologies and the reorientations of the stress induce by the facies borders is a requirement for the next generation of fractured reservoir modelers and basin modeling tools.

Radial extension of a homogeneous elastic cylinder
Let us consider a homogeneous linear elastic and isotropic cylinder which initial geometry is defined by a radius R i and a thickness h.This cylinder is submitted to a radial displacement-controlled extension such that the final radius is R f > R i .This extension is uniform along the thickness which is forced to remain unchanged with the deformation, such that plane strain conditions are satisfied.Let us use the cylindrical coordinates (r, θ, z) associated to the orthonormal basis (e r , e θ , e z ), which vectors are along the radius, along the circumference and along the cylinder revolution axis respectively.The displacement field is : Note that it only depends linearly on the radial coordinate.For take of simplicity we write the deformation tensor in its linear form ε ε, with components: Introducing the radial displacement field, we obtain: The Cauchy stress tensor σ σ can then be determined with the linear elastic relationship, σ σ = λtr (ε ε) II + 2 με ε, where λ and μ are the Lamé coefficients, II is the second order identity tensor and tr (.) is the trace operator.We therefore obtain:

Figure 1
Figure 1 Figure42D unfolding, with or without fault lip imposed contacts.

Figure 6
Figure 6 Stress field resulting from the restoration presented in Figure 5. Top: maximum (least tensile) principal stress value.Blue (negative values) means extension and red (positive values) means compression.Bottom left: closer view displaying the least tensile principal stress direction.Bottom right: closer view displaying the most tensile principal stress direction.It is rather parallel to the bedding.

Figure 5 3D
Figure 5 3D restoration of a faulted anticline in the Sub-Andean Zone of Bolivia.Left: initial geometry.Right: restored geometry.

Figure 8
Figure 8Stress field resulting from the forward deformation presented in Figure7(maximum -least tensile -principal stress direction is vertical).Top left: intermediate principal stress direction (perpendicular to the edge of the channel inside of it).Top right: most tensile principal stress direction (parallel to the edge of the channel inside of it).Bottom: views of the stress field displayed in principal components with ellipsoids.Red ellipses are normal to the intermediate principal stress direction, and white ellipses are normal to the most tensile principal stress direction.

theFinalCode_Aster
West African Margin Using 3D Restoration.Journal of Structural Geology, 24, 783-796.Rowan, M. and Kligfield, R. (1989) Cross-Section Restoration and Balancing as an Aid to Seismic Interpretation in Extensional Terrains.American Association of Petroleum Geologists Bulletin, 73, 955-966.Santi, M., Campos, J.-L. and Martha L. (2003) 3D Geological Restoration Using a Finite-Element Approach.Paper presented at the 23 rd Gocad Meeting, Nancy (France).Suppe, J. (1983) Geometry and Kinematics of Fault-Bend Folding.American Journal of Science, 283, 684-721.White, N., Jackson, J. A. and Mc Kenzie, D. (1986) The Relationship Between the Geometry of Normal Faults and That of Sedimentary Layers in Their Hanging Walls.Journal of Structural Geology, 8, 879-909.Code_Aster has been developed and used internally by Électricité de France (EDF) for more than 15 years.Since 1997 and every two years, EDF releases a fully documented free version of the source code of Code_Aster, available online (http://www.code-aster.org)and managed by a group of twenty permanent developers.This version is updated every six months.A "development" version can also be downloaded on the website (current version is now Code_Aster v. 7.3).