Modeling Fluid Flow in Faulted Basins.

— Mode´lisation des transferts ﬂuides dans les bassins faille´s — Cet article pre´sente un nouveau simulateur de bassin spe´cialement conc¸u pour mieux prendre en compte l’impact des failles sur l’e´coulement. Il simule les e´coulements, les transferts thermiques et la ge´ne´ration des hydrocarbures. Compare´ aux simulateurs de bassins classiques, il utilise un concept innovant de maillage e´volutif sufﬁsamment ﬂexible pour suivre les de´formations d’un bassin en contexte tectonique complexe. Le maillage suit les se´diments et se de´forme continument pour suivre la se´dimentation, la compaction et les de´formations cine´matiques. Bien qu’il soit organise´ suivant les couches stratigraphiques, le maillage est globalement non-structure´. Il peut contenir des mailles de type varie´ telles que des hexae`dres, des te´trae`dres, des pyramides, des prismes ou des hexae`dres de´ge´ne´re´s. Les failles Abstract — Modeling Fluid Flow in Faulted Basins — This paper presents a basin simulator designed to better take faults into account, either as conduits or as barriers to ﬂuid ﬂow. It computes hydrocarbon generation, ﬂuid ﬂow and heat transfer on the 4D (space and time) geometry obtained by 3D volume restoration. Contrary to classical basin simulators, this calculator does not require a structured mesh based on vertical pillars nor a multi-block structure associated to the fault network. The mesh follows the sediments during the evolution of the basin. It deforms continuously with respect to time to account for sedimentation, erosion, compaction and kinematic displacements. The simulation domain is structured in layers, in order to handle properly the corresponding heter-ogeneities and to follow the sedimentation processes (thickening of the layers). In each layer, the mesh is unstructured: it may include several types of cells such as tetrahedra, hexahedra, pyramid, prism, etc. However, a mesh composed mainly of hexahedra is preferred as they are well suited to the layered structure of the basin. Faults are handled as internal boundaries across which the mesh is non-matching. Different models are proposed for fault behavior such as impervious fault, ﬂow across fault or conductive fault. The calculator is based on a cell centered Finite Volume discretisation, which ensures conservation of physical quantities (mass of ﬂuid, heat) at a discrete level and which accounts properly for heterogeneities. The numerical scheme handles the non matching meshes and guaranties appropriate connection of cells across faults. Results on a synthetic basin demonstrate the capabilities of this new simulator.


NOMENCLATURE
Set of face-to-face contacts C N,i,j Set of node-to-face contacts F A face of the grid g Intensity of gravitỹ g Gravity vector h Height of a 1D vertical segment hs Solid height of a 1D vertical segment k Time superscript or cell index in 1D K m Permeability M ðtÞ The set of cells of the grid N ðtÞ The set of nodes of the grid N i A set of nodes N s The set of cells that have s as a node N d The set of cells that share at least one node with face d n Normal to a surface P Pressure q s Deposition rate of sediment q w Deposition rate of water s A node of the grid SðtÞ The set of faces of the grid Water viscosity q s Solid density q w Water density r Effective stress r l Lithostatic pressurẽ r Lithostatic potential

INTRODUCTION
Thermogenic origin of hydrocarbons began to be understood during the sixties. Shortly after, the first numerical simulators of sedimentary basins and their petroleum systems were proposed (Doligez, 1987). They were originally limited to vertical 1D models. They evolved into 2D models of vertical basin sections during the eighties (Ungerer et al., 1990), and into full three-dimensional models during the mid-nineties (Schneider et al., 2000). All these simulators had in common to be limited to basin architectures only controlled by erosion, sedimentation, vertical compaction and vertical shear deformation resulting from basement deformation. In terms of tectonic deformation, only vertical slip along vertical faults could be taken into account through vertical shear zones.
This hypothesis was too limitative for most sedimentary basins having experienced a complex tectonic history, especially compressive tectonic settings, for two reasons. First, as horizontal displacement was not allowed, these models failed to account for the evolution through time of the relative position of source rocks and traps, because of relative horizontal displacement along decollements, for instance. Secondly, in faulted environments, simulation of fault-controlled fluid flow was oversimplified: the lateral connectivity of the deposits could not be interrupted geometrically at fault contact because of mesh constraints. So, faults were modeled as vertical shear zones, along which one had to artificially adjust rock permeability in order to account for enhanced fault flow or flow barrier effect.
At the same time, fault control on fluid flow and on oil migration began to be studied in detail (Knipe et al., 1998), revealing a complex internal structure of the faults. Associated to fault slip and especially in shaly stratified environments, fault core can be progressively filled with a shaly impervious material (named "fault gouge") which decreases across-fault flow efficiency, whereas a damage zone on both compartment walls can develop, thus enhancing along-fault flow. Simultaneously, petroleum system simulation was extended from exploration purpose to overpressure prediction. The latter requires accurate modeling of fault fluid flow, as overpressure occurrence is at first order controlled by architecture and permeability of poorly permeable bodies and by faults which can either enhance fluid flow or prevent for fluid escape.
Two approaches were proposed during the last decade to simulate petroleum systems in these faulted contexts: Ceres (Schneider et al., 2002), and Tecklink (Baur et al., 2009). Both approaches are limited to vertical shear deformation of each fault-delimited structural block, the structural blocks being allowed to slip against each other. Ceres is limited to 2D cross-sections and Tecklink is based on a very coarse time-discretization named "Paleo stepping" which affects computation accuracy.
All these factors motivated the development of a newgeneration petroleum system simulator as well as a modified modeling workflow which take up the technical challenges of structurally complex basins, in 2D and 3D, without limitations on deformation mode and time discretization, and which benefit from the newest advances of computer science and technology.

Structurally Complex Basins: An Updated Modeling Workflow Involving Structural Restoration
The usual basin modeling and petroleum system simulation workflow is made of the following components: -geomodel construction and gridding in order to provide a discrete representation of the basin architecture at Present; -geomodel backstripping: it is an automatic backwards procedure which aims at restoring the model architecture at all times since the beginning of the basin subsidence. Backstripping procedure iteratively modifies the discretisation grid at Present in order to take into account vertical expansion, creation or removal of grid cells due to backwards unloading, erosion and sedimentation. The basic hypothesis here is that only basement subsidence, deposition, erosion and vertical compaction controls this time-dependant architecture. The backwards movement of any solid particle in the basin is only vertical, from Present to deposition. Backstripping is an automatic procedure which can be considered as a pre-processing of the forward calculation. It can be easily re-run at each new hypothesis on Present-day basin architecture, lithology distribution and compaction parameters; -forward simulation: it is the procedure which calculates the evolution through time of all relevant physical variables such as temperature, porosity, fluid pressure and hydrocarbon saturation. Forward simulation is based on the resolution of partial derivative equations that model the physical processes of heat transfer, mechanical equilibrium, hydrocarbon geochemistry and fluid motion. The discretisation grid is the time-dependant grid calculated by backstripping.
As soon as basin architecture is affected by tectonic movement, the assumption that particle displacement is vertical all along the basin history is no longer valid. The restoration of past basin architecture can no longer be calculated automatically by backstripping, and the restoration of basin past architecture must be performed according a new procedure named "structural restoration".
Structural restoration provides a structural scenario. The latter is a description of the basin architecture as a function of geological time. It is delivered in the form of a discrete number of deformation stages between basin initiation and Present.
Different computer-based approaches are proposed to the structural geologist to build this structural scenario. Most of them are 2D and operate on vertical cross-sections, as 2D-move (Gibbs, 1984), Locace and Kine3D2-XS (Moretti and Larre`re, 1989), Lithotect (Geiser et al., 1988), Dynel (Maerten and Maerten, 2006), and Geosec (Kligfield et al., 1986). They provide the user with tools to interactively apply a complex deformation field to a given cross-section described as a set of horizons and faults. This includes erosion and deposition, compaction, and sliding movement along faults. Although these applications were initially designed for cross-section balancing, a whole scenario can be built by sequentially applying deformation steps to the cross-section of interest. First commercial releases of 3D restoration tools like Kine3D-3 (Moretti et al., 2006) and 3DMove were also recently made available but are still in their infancy.
Running a forward petroleum system simulation on a structural scenario requires that the time-dependant discretisation grid used by forward modeling fits to the deformation scenario delivered by structural restoration. With most restoration applications, the grid must be built subsequent to restoration, according to the requirements of the forward simulator. In the case of Ceres (Schneider, 2003) or in the case of the workflow proposed in this paper which involves Kine3D-3 as a gridbuilder, the structural restoration operates on the same grid as the forward simulator. The grid is built to fit the architecture at Present and is subsequently deformed by the restoration procedure. The time-dependant grid obtained in this way is then used as an input of forward simulation.
Creating the time-dependant mesh that best fits the needs of forward petroleum system simulation remains a technical challenge. "Paleo Stepping" approach presented by Baur et al. (2009) proposes to discretise tectonic episodes in a set of successive static grids representing successive deformation stages. At the breakpoint between two stages, the physical values relevant to forward simulation such as pressure and temperature are instantaneously transferred from one grid to the next one. According to this technique, the grid as well as physical values are not continuous through time, which may appear as a risky simplification of reality.
According to the technology described in this paper, the forward simulator operates on a unique grid that gradually deforms through time, by linear interpolation between deformation stages. Such time-dependant grid is provided by Kine3D-3 for 3D cases, or by the restoration module of Ceres in 2D, as mentioned above.
Considering other restoration tools, we need to create the gradually deforming grid subsequently to restoration. A new computer tool designed to build such a mesh is currently under development, for 2D cross-sections.
After all, the basin modeling and petroleum system simulation workflow presented in this paper to achieve modeling of tectonically complex basins is made of the following components: -structural restoration, in order to restore the past architecture of the basin in the form of a structural scenario; -gridding of the obtained structural scenario, if gridding is not imbedded in the structural restoration; -forward petroleum system simulation.
More information on structural restoration and grid building is given in Section 3. Unlike backstripping, structural restoration is a heavy task involving a lot of interactivity. Consequently, it can not be easily re-run at each new hypothesis on structural architecture, lithology distribution or compaction parameters, during basin model calibration procedure. Conventional calculators have their own technique to overcome backstripping inaccuracy introduced by the backstripping decompaction model which slightly differs from their compaction model. However, in our approach we had to design our forward simulation in order to be even more flexible regarding to past geometry inaccuracy. In this sense, we propose the decoupling of porosity and geometry described further.
In the following sections, we will mainly focus on the technology of forward simulation. Section 1 describes the constraints on time and space discretisation which are required by our forward simulator and Section 2 focuses on forward simulation itself. Section 3 illustrates this procedure with a semi-synthetic example.

NUMERICAL MODELING OF BASIN ARCHITECTURE AND ITS EVOLUTION THROUGH TIME
To describe the evolution of a 3D sedimentary basin in complex settings, we introduce an innovative evolving grid concept able to follow the various geometrical and topological changes supported by the basin through time. Before going into the details, we recall some of the characteristics of the grids introduced previously in the context of basin modeling.

Time and Space Basin Discretisation Based on Vertical Pillar Grids
First introduced for a vertical column, the principle of forward basin simulation is to compute the thickness of the different layers through time, starting with the deposition of the deepest layer, assuming conservation of rock volume and solving fluid flow equation coupled with compaction. In this process, the input geometrical data is the depositional rock thickness of each layer which is in fact unknown, the only true geometrical input being the geometry of the column at Present. Therefore, depositional rock thicknesses (or deposition rates) are first estimated from present day layer thicknesses using a priori porosity/depth curves. This estimation can also account for non deposition and erosion events which can modify the maximum depth reached during a layer history. This is usually known as the "backstripping" process. Forward modeling of layers deposition coupled with fluid flow then leads to layer thicknesses at Present which may be different from the input present day geometry. It can then be necessary to iterate the forward simulation process to calibrate the rock deposition thicknesses in order to satisfy present day geometry.
To describe the underlying space-time grid concept, we introduce a decomposition of the time interval ½ÀT ; 0 into K subintervals, where -T stands for the time where the first layer started to deposit and 0 for present day: corresponds to a stratigraphic event namely the deposition of a layer, a non deposition or an erosion period. The time subintervals indexes give a natural indexation of the layers as time can be related to depth through the deposition time of a layer and therefore allows a simple management of the evolving grid: the grid is always composed of K segments which possibly have a zero height. More precisely, let us consider a given time t. We denote by K(t) the time index such that t 2 ÀT KðtÞÀ1 ; ÀT KðtÞ Â Ã (Fig. 1) If the time subinterval ÀT kÀ1 ; ÀT k Â Ã corresponds to an erosion or a hiatus, the associated grid cell has always a zero thickness i.e. h k ðtÞ ¼ 0 for t 2 ÀT ; 0 ½ . Erosion can also lead to additional zero thickness cells if layers are fully eroded.
At time t, the grid is then entirely defined by the values h k ðtÞ for k ¼ 1; :::; KðtÞ.
Flow and temperature equations are solved on this evolving grid that follows the sediment. To this end, intermediate time-steps can be introduced within each time subinterval where the grid definition deduced from the layer thicknesses still holds. We can already notice that this grid structure requires handling zero length cells. In the framework of finite volume methods, it then becomes necessary to manage an additional mesh connectivity which gives, for a given cell, the upward neighboring cell that has a non-zero height.
The 1D vertical approach has been extended to 2D or 3D basin where no horizontal displacement occurs, by introducing vertical pillar grids for which the horizontal section of the basin (in the (x, y) plane) is meshed using a cartesian grid. The time interval ½ÀT ; 0 is still decomposed into a set of successive time subintervals. They include all the stratigraphic events that have to be represented in the 3D model. Then, a 1D vertical grid (similar to the one considered in 1D) is associated to each node (i, j) of the horizontal cartesian grid. The nodes are now indexed with (i, j, k) and are completely defined by the layer heights of the column. As before, the grid has a constant structure all over the time interval, layers deposition or erosion being handled through zero height edges. The backstripping process is applied to all the columns and gives a first estimation of the depositional height of each layer.
In the 3D space, this leads naturally to a 3D structured grid where each cell is defined by its four vertical edges (Fig. 3, 4). However, zero height edges can lead to cells which are geometrically degenerated hexahedrons, and that can even have a zero volume. To manage these degenerated items, it is necessary to introduce an additional connectivity in the vertical direction (Fig. 3).
This process, which is based on a tight relation between time and space, is very efficient for settings dominated by normal deposition and erosion. In complex settings, where faulting leads to important lateral movements, this approach is no more suitable. Indeed, the evolution of the basin geometry can no longer be managed through a grid whose horizontal structure is fixed in time. Moreover backstripping is not valid anymore. It is replaced with the structural restoration of the basin, which gives the paleo-geometries at given ages. Different approaches have been proposed to handle Time axis Figure 1 Decomposition of the time interval into subintervals ÀT ; 0 these complex geometric evolutions. Schneider (2003) proposed an approach to manage 2D sections of a basin whose geometry also evolves by salt or mud creeping and displacement along faults. Still based on a time interval discretised into stratigraphic events, this approach considered that the fault network decomposes the 2D section into a set of blocks. The blocks are organized with respect to each other through a bottom to top ordering that defines how they are stacked one above the other (Fig. 5). Each block can individually be described as a standard model. It is discretised by a vertical pillar 2D grid which can handle, during a stratigraphic event, deposition, erosion, compaction, but also other continuous deformations which are determined during the restoration stage, as long as they keep the pillars of the grid vertical (translation, vertical shear). Introducing a more general decomposition of the time interval into "paleo steps" instead of "event steps", Baur et al. (2009) describes an innovative technique designed to handle 3D basin in complex geological settings. The corresponding simulator takes the complete geometry of the basin at each paleo-step as an input, and directly jumps Schematic representation of a 2D vertical pillar grid with degenerated items. Vertical pillars are indexed with i while stratigraphic events and layers are indexed with k. Here, stratigraphic event k + 1 corresponds to a hiatus from vertical pillar i + 1 up to the right boundary. For any i' > i, nodes (i', k) and (i', k + 1), represented as black dots, have the same coordinates leading to zero area cells. For instance, cell (i + 1, k) whose nodes are (i + 1, k), (i + 2, k), (i + 2, k + 1), (i + 1, k + 1), is squashed. This modifies the vertical connectivity: upward neighboring cell of cell (i + 1, k À 1) is cell (i + 1, k + 1). Lateral connectivity is also modified but in a much weaker sense as degenerated items only suppress lateral neighboring cells.

Figure 4
Example of a 3D vertical pillar grid (Schneider et al., 2000). from one paleo-geometry to the next one. At a given step, the domain is described as a set of standard blocks separated by faults in order to manage multiple-z value for an individual layer. This decomposition can be different from one paleo-step to the other since a block can be split into sub-blocks due to fault activation. As the complete geometry is an input data, porosity evolution is decoupled from the geometry evolution.

A New Evolving Grid Concept
Going a step forward in considering complex settings, we introduce a new 4D grid concept to describe the evolution of a sedimentary basin with a particular attention to faults. Faults are geometrically described as discontinuity surfaces along which sliding can occur. They can have a major impact on fluid flow paths and therefore on overpressure distribution. Indeed, due to displacement along the fault surface, stratigraphically disconnected layers can become juxtaposed across a fault creating a new flow path or on the contrary breaking an existing one. As overpressure development is the result of a transient process, these flow path modifications should be captured as continuously as possible.
Compared to the previous techniques presented above, the new evolving grid concept has been designed with two major motivations. First, it should be able to follow the deformation of the basin in a continuous way with respect to time as it is an important factor of overpressure development. In particular, sliding along a fault surface should lead to continuous modifications of the flow paths which can be hardly obtained with an approach where the geometry jumps from one configuration to the other at paleo-steps. Second, it should be more flexible concerning the fault network description and the possible structural deformations. For instance, the grid structure should not be too tightly linked to the fault network configuration like in a multi-block approach where the fault network must split the basin into individual blocks.
As in the approach of (Baur et al., 2009), the evolving grid concept relies on a workflow where input geometries are given once for all by the restoration process. Since restoration usually provides a quite rough geometry evolution and compaction description, we propose a decoupling between the porosity evolution and the geometry evolution. We use a formulation of the model which certainly remains an approximation but which tries to limit the impact of the decoupling. Before describing the details of the evolving grid, we present how porosity and geometry are decoupled.

Decoupling Porosity and Geometry Evolutions
To understand the idea underlying the proposed approach to decouple porosity and geometry evolution, we consider the case of a 1D basin mentioned previously that furthermore undergoes only deposition and compaction processes. Since this decoupling mainly involves fluid flow and compaction, we omit thermal transfer and temperature dependencies.
In the following, dependency of the variables with respect to the z coordinate or to time is not specified unless it is necessary for understanding. Example of a 2D multi-block grid used in Ceres (Vilasi et al., 2009). The vertical lines are the vertical pillars of the different blocks.
Fluid flow coupled with compaction is governed by the following set of equations: -mass conservation of solid: q s stands for the rate of sediment deposited at the top of column, / the porosity,ṽ s the solid velocity, q s the solid density; -mass conservation of fluid: q w stands for the rate of water deposited at the top of the basin which is expressed as a function of the deposition porosity: q w ¼ q s / 0 =ð1 À / 0 Þ, q w the water density,ṽ w the water velocity; -Darcy's law: where l w is the water viscosity, K m the permeability, P the pressure,g the gravity; -sedimentary load: where r l is the weight of the above sedimentary column or lithostatic pressure; -compaction law: where porosity / is expressed as a function of effective stress r ¼ r l À P through a lithology dependent law. These equations are completed with constitutive laws (permeability as a function of porosity, fluid density and viscosity as a function of pressure) and appropriate boundary conditions: -on the top boundary which is assumed here to be under the sea level, the pressure and sedimentary load are equal to the weight of the above water column: -on the bottom boundary, there is a no fluid flow condition:ũ w Áñ ¼ 0 To highlight the coupling between the different equations and their link with the geometry, we reformulate the equations using the overpressure oP and lithostatic potentialr unknowns defined asr ¼ r l À P h and oP ¼ P À P h , where the hydrostatic pressure P h satisfies the following equation: and the boundary condition P h ¼ P top on the top boundary.
For the sake of simplicity, we assume here that water density and viscosity are constant but the principle remains valid with varying densities and viscosities. We first rewrite Equations (3) to (5).
Subtracting (6) from (4) gives: Furthermore, Darcy's law (3) simplifies to: Finally, noticing that the effective stress satisfies r ¼r À oP, the compaction law (5) writes: The boundary conditions at the top of the sedimentary column formulated with the unknowns oP and lithostaticr are written as follows: oP ¼ 0 andr ¼ 0.
To reformulate the mass conservation equations, we first discretise them in space following a cell centered finite volume method (time discretisation is not discussed here, as it is not mandatory for the understanding of the approach). We consider a 1D Lagrangian mesh, in which each individual node moves with the velocityṽ s , where discrete unknowns in cell x k are denoted with the subscript k. Considering as main unknowns oP k ,r k , / k and h k for k ¼ 1; . . . ; KðtÞ, the set of discretised equations is then (details of the derivation are presented in Appendix): r KðtÞ ¼ ðq s À q w Þhs KðtÞ g=2 with the additional discretised Darcy's law: In this formulation where the rate of deposited sediment is assumed to be an input data, we notice that the unknown heights h k only appear in Darcy's law discretisation (15) through the transmissibility. This leads us to propose an approach where on one hand, solid thickness becomes an unknown, which is independent of the real thickness by relaxing the constraint "hs k ¼ h k ð1 À / k Þ", and, on the other hand, the real thickness is an input data (obtained by the backstripping stage for instance). This process is then integrated in an outer "solid content loop" that aims at calibrating the rates of deposited sediments to satisfy the relation hs k;t¼0 Â ð1 À / k;t¼0 Þ ¼ h k;t¼0 at present time. We do so to account for the fact that thicknesses at Present are in fact the only true geometrical input data. This outer "solid content loop" can be schematically defined by the following iterative procedure (exponent i refers to iteration i): -initialization: estimate the sedimentation rates q 0 s;k and the geometry evolution h 0 k ðtÞ (with a backstripping stage for instance); define a correction factor of the sedimentation rates: ; Until e i k ¼ 1. At convergence of the iterations, solid and real thicknesses satisfy the constraint: and the computed unknowns satisfy the solid and water mass balances. The only "non-consistent" equation is the discretised Darcy's law where the transmissibility is expressed through the estimated real thickness at a given time. This can be further corrected using the modified approximation: which can be seen as the transmissibility computed with an approximate geometry but for a corrected vertical permeability: For a vertical column and with the permeability correction, this iterative loop is similar to the geometrical loop usually performed in classical contexts to match the present day geometry although it keeps the real thickness as an input data. This formulation can moreover be extended to a multi-dimensional approach keeping an input geometry evolution (obtained by the backstripping stage for instance) and replacing solid and real thicknesses by solid and real volume of each cell. Although less accurate than the standard approach with a true multi-dimensional unknown geometry, this process somehow preserves some of the essential features of the physics of fluid flow in a compacting porous medium. We will follow this approach for complex settings since it allows accounting for compaction even though the evolution of the input geometry is an input data.

A 4D Grid
The new evolving grid concept is based on a deformable mesh that follows the input basin evolution through time provided by restoration. This evolution includes stratigraphic modifications, cell deformations and sliding along faults. The basin is represented by a fully unstructured mesh that does not rely on any preferred direction and that can therefore follow deformations that do not preserve the vertical direction. In the grid, faults are handled as pairs of internal boundaries that can slide one with respect to the other through time. Moreover the evolution of the grid is described as a succession of incremental modifications. This give much more flexibility than an approach where changes can only affect node coordinates of a global grid whose structure includes all the layers at any time. To define precisely the evolving grid, we again introduce a decomposition of the time interval into subintervals: where L, the number of subintervals, needs only to be larger than K the number of stratigraphic intervals. Each subinterval can then include one or both of the following events: -one stratigraphic event (deposition, erosion), -one structural deformation event.
It must be emphasized that a deformation event, seen from the geological point of view, can be further decomposed into smaller time intervals if the corresponding deformation or displacement is very non linear. To remain consistent with stratigraphic evolutions, a deposition or erosion event can then also correspond to several successive time events.
In the following sections, we first give the characteristics of the grid at a given time. Then, we specify how the grid can evolve through time and we finally end with a focus on connectivity across fault surfaces.

Grid at a Given Time
At a given time, the grid is fully unstructured and described as a set of cells that are of one of the following types: hexahedra, pyramid, prism, tetrahedra, degenerated hexahedra (Fig. 6).
A cell is given as a list of nodes and a type which implicitly define its faces. A face can have two neighboring cells (i.e. the two cells that contain the nodes of the given face) and it is then a topologically internal face. It can alternatively have only one neighboring cell and it is then a topologically boundary face.
In the unstructured grid, a layer corresponds to a set of cells which means that a horizon corresponds to a set of faces.
A fault is handled as two internal boundaries, one for each side of the fault (Fig. 7). An internal boundary is understood as a set of faces that are topologically boundary faces but that are not geometrically located at the boundary of the domain. The two boundaries are geometrically close to each other but are not glued to one another via the structure of the mesh. The mesh is non-matching across the fault as the nodes from one side of the fault do not belong to the other side of the fault. Connectivity across a fault is said to be a geometrical property, that is to say that it results from geometrical computations that use nodes coordinates ( Sect. 1.2.2.3). It is different from intrinsic connectivity between cells that share a face, which is a topological property that does not depend on the nodes coordinates.
This definition of faults gives a lot of flexibility in the fault network definition. For instance, a fault surface does not need to be extended to the boundary of the domain; a particular fault face can belong to several fault surfaces allowing for crossing faults. Moreover, there is not any concept of block ordering.

Evolution of the Grid
Starting with an empty grid, grid evolution is defined by successive incremental modifications defined for each time subinterval ÀT lÀ1 ; ÀT l Â Ã that can include: -addition of cells, which corresponds to the deposition of a layer. Each new cell is defined by its type and nodes and by the coordinates of the new nodes at time of deposition; -deletion, modification (cut) of cells, that correspond to an erosion. The deleted cells are given as a set of cells while cells that change of type or of geometry due to erosion have to be redefined; -displacements of nodes given through their coordinates at the end of the time subinterval, that describe: cell thickness increase during deposition, compaction, structural deformation, sliding along the fault surfaces. The two first types of modification are called topological evolutions of the grid as they modify its connectivity. On the other hand, the last one is called a geometrical Figure 6 The different types of cells. Top, from left to right, hexahedra, pyramid, prism and tetrahedra. Bottom, degenerated hexaedra. evolution. In that framework, sliding along fault surfaces corresponds to geometrical evolution. Although many types of cells can be used, we can notice that hexahedra or prism type cells are better suited to follow layer thickening associated to deposition. Given these incremental modifications, the grid is computed for any time t in the considered subinterval ÀT lÀ1 ; ÀT l Â Ã (Fig. 8): -the topological modifications are taken into account instantaneously, right at the beginning of the subinterval i.e. at ÀT lÀ1 ; -the geometrical changes are taken into account in a continuous way using linear interpolation with respect to time between the coordinates of the nodes at time ÀT lÀ1 and time ÀT l . This last point could appear as not enough accurate when large displacements or deformations are considered. To overcome this drawback, additional time subintervals can be added that give intermediate input geometries.

Connectivity Across Faults
To simulate mass and heat transfer on the considered grid, intrinsic connectivity is not sufficient as it does not contain any information on how the two sides of a fault are glued together. This information is for instance necessary to measure how a layer is juxtaposed to other layers across a fault surface. As already mentioned, the grid is topologically non-matching across a fault, there can even be gaps and overlaps due to the decomposition of the fault surface into elementary faces but also due to approximations on how node coordinates change during an event (linear interpolation) resulting in an approximate sliding of the two fault surfaces. At a given time, standard geometrical property (cell volumes, face area, etc.) computation needs therefore to be complemented with a specific computation of connectivity across faults, which is performed knowing the values of the nodes coordinates. The first additional connectivity is called "node-to-face" contact and gives information on how a node of one surface projects itself on the other surface. The second one called "face-to-face" contact gives pairs of faces (one from each surface of the fault) that are geometrically in contact across a fault.
To precise the actual content of the fault contacts, let us consider a particular fault defined as the pair of two sets of faces S 1 and S 2 (Fig. 9a for a schematic representation in 2D). We also introduce N 1 and N 2 the sets of corresponding nodes.
The node-to-face contact (Fig. 9c) of S i on S j ði ¼ 1; 2; j ¼ 1; 2; i 6 ¼ jÞ computes the following set of node, face pairs: Evolution through time of the simple basin considered in Figure 7. From a) to d): grid after the deposition of respectively one, two, three and four layers. From e) to i): snapshots of the grid during the last event when deposition of the fifth layer is simultaneous to sliding along the fault.
Schematic illustrations of contacts across the surfaces of a fault in a 2D section. a) Sketch of a 2D section composed of 3 layers cut by one fault; b) the fault is defined as the pair of surfaces S 1 and S 2 (blocks have been artificially shifted for the sake of clarity); c) node-to-face contact of node n 2 of surface S 2 and face F 1 of surface S 1 ; d) face-to-face contact between faces F 1 of S 1 and F 2 on S 2 .
The face-to-face contact (Fig. 9d) computes the following set of face, face pairs: It also gives the contact area between the faces of each pair.
Although this face-to-face contact is written using formal geometric intersection of faces, it must be understood in a weaker sense as fault surfaces are not as simple as in the schematic illustration shown in Figure 9. We use a projection of both surfaces on a common mean surface to compute the desired intersections. Detailed description of the algorithms used to effectively compute these contact properties is out of the scope of this paper.
Intersecting faults can also be taken into account; actually a face located near the intersection can appear in pair of contact faces belonging to several faults.

FLUID FLOW MODEL
We consider here only one-phase fluid flow (water). A more rigorous model of the considered physical problem would rely on a true poromechanical formulation accounting for 3D compaction, large deformations and displacements with frictional contact on fault surfaces, and complex non linear constitutive laws relating stress and strain. Although recent progress has been made towards this objective (Guilmin, 2012), such an approach remains elusive for real case studies. Regarding the lack of control one has at geological time scales, more approximate models with manageable complexity can meet our needs. We present here a model whose objective is to allow the study of the impact of faults on fluid flow. Its main ingredients are on one hand, a given input evolving geometry associated to a particular formulation of the equations as discussed in the previous section and on the other hand, a particular interface fault model able to simulate the different fault behaviors.
Before putting a special emphasis on faults, we describe the equations used to model fluid flow. Heat transfer and hydrocarbon generation are also taken into account but not described here.

Governing Equations
As for standard settings, we assume that the effective stress only depends on the weight of the above sedimentary column. Moreover, we assume that water table is at the surface of sediments (the porous volume is fully saturated with water). The considered equations are then mass conservation of solid phase, mass conservation of the fluid phase, vertical equilibrium for the weight of sediment, Darcy's law and a porosity effective/stress relationship. To these equations, we add constitutive laws for the fluid and the petrophysical properties and appropriate boundary conditions. In the following, we won't detail the different constitutive laws that are considered since they are parameters of the model. In order to avoid unnecessarily complicated notations, we will even often assume constant properties.
In this section, we limit the derivation of the fluid flow model to a basin which is not cut by faults and that has therefore "standard" external boundaries. We denote by XðtÞ the domain representing the basin at a given time t and assume that we can define its "top boundary" C top ðtÞ such that: its "bottom boundary" C bottom ðtÞ that corresponds to the base of the sediment and its "lateral boundaries" C lat ðtÞ. Usual boundary conditions give the pressure and sediment weight on the top boundary, while no flux boundary conditions are set on the other boundaries (more general conditions could however be considered on the lateral boundary). We introduce the following notations for the evolving mesh defined in Section 1.2.2: -GðtÞ the grid at time t: M ðtÞ is the set of cells, N ðtÞ the set of nodes and SðtÞ the set of faces; -S top ðtÞ, S bottom ðtÞ the sets of boundary faces that correspond respectively to C top ðtÞ; C bottom ðtÞ and C lat ðtÞ; N top ðtÞ the set of boundary nodes corresponding to C top ðtÞ; -jðtÞ a cell; -d a face (the time dependence is omitted here); -djðtÞ the boundary of jðtÞ which is a union of faces; -N d the set of cells that share at least one node with d (Fig. 10); s a node, N s the set of cells that have s as a node, J s the set neighboring nodes including s itself (Fig. 10).
Following the approach described in Section 1.2.1, the equations are formulated using porosity, overpressure oP and lithostaticr potential unknowns. Since Finite Volume methods are well suited for conservation equations and capture efficiently the heterogeneities, we use a cell centered finite volume method to discretise the mass conservation equations associated to Darcy's law and define porosity and overpressure unknowns in each cell. To decouple porosity evolution from the geometry evolution, we consider an additional unknown, namely the solid volume of each cell. The corresponding unknowns defined for each cell are then oP j , / j and Vols j .
For the sake of clarity, we consider here that the water has a constant density and a constant viscosity although density and viscosity variations could be accounted for. The mass conservation equations are first discretised in space by integration over each cell jðtÞ 2 MðtÞ that is assumed to follow the sediments: -discrete mass conservation of solid: -discrete mass conservation of fluid: where the water filtration velocityũ w is given by Darcy's law: K m is the permeability which at least depends on porosity via a constitutive law. For the sake of simplicity, the termrP h À q wg which is equal to zero if the basin is underwater, is not written down in Darcy's law although it could be taken into account. Following the Finite Volume principle, the second term of the left hand side is modified using Green's theorem: where F w;d % R dũ w Áñ d ,ñ d is the normal to face d, I j;d ¼ 1 ifñ d is oriented outward with respect to jðtÞ and I j;d ¼ À1 otherwise.
Permeability is discretised constant in each cell; its cell value usually depends algebraically on the cell unknowns. To approximate the flux over each edge, we use a two-point or a multipoint flux approximation (Aavatsmark et al., 1998) that has been especially designed for approximating Darcy's flux on general grids with discontinuous permeability tensor. It leads to an expression of the form: where T d;L depends on cell permeabilities and geometries. Given overpressure or no flux boundary conditions are classically handled in this approximation and not detailed here.

-sedimentary load
As for vertical compaction, the lithostatic potential is assumed to satisfy the "vertical" equation: with a boundary condition that givesr ¼ 0 on C top ðtÞ, and porosity is assumed to satisfy a porosity/effective stress relationship: In 1D, or on vertical pillar grids, Equation (22) is easily discretised defining cell centered unknowns as in (11). In a fully unstructured grid, the cell centers are no more vertically aligned and a multi-D discretisation scheme is necessary. A natural extension of what is done for 1D vertical grid would be to still consider the discrete lithostatic potential unknowns located at cell centers. Yet, a consistent approximation of Ào=oz ðrÞ would then link unknowns from different neighboring cells. Since this would require working with a connectivity which is not directly stored in the grid structure, we prefer to consider discrete unknowns located at the nodes. For tetrahedral meshes, (Mello et al., 2009) proposed an approach based on the computation of the intersection of the mesh cells with the vertical lines issued from the nodes. Equation (22) could alternatively be discretised using a finite element approximation with appropriate stabilization due to the fact that the differential operator is of first order (Mello et al., 2001). However, to avoid defining shape functions for the degenerated hexahedral items which are non standard, we choose a discretisation scheme specific to this operator, somehow more geometric and designed to be exact for linear function. We denote byr s the approximation ofr at node s. Equation (22) is first discretised for each node and in each cell using its neighboring nodes which gives: where J j;s denotes the set of neighboring nodes of s in the cell j. The coefficients a j;s;s0 are such that this approximation is exact for linear function. This approximation gives then a residual of (22) for each node of each cell that writes: a j;s;s0rs0 À ðq s À q w ÞgVols j To get an equation for each node, these residuals are condensed on nodes using weights a j;s that introduce some "top to bottom" upwinding coherent with the differential operator: X j a j;s R j;s ¼ 0 For each cell, we define the set of "bottom" nodes P j that is to say nodes for which the coefficient a j;s;s is positive. Then, we simply take a j;s ¼ 1= P K j j if s belongs to P j , a j;s ¼ 0 otherwise, which in particular ensures that the diagonal coefficients of the associated linear system is positive. This leads to the following form of the discretised equation for the lithostatic unknown, for each node s 2 N ðtÞ=N top ðtÞ: X s02Js b s;s0rs0 À X j2N s a j;s q s À q w ð Þ gVols j ¼ 0 Boundary conditions give directly the value ofr s for nodes that belong to N top ðtÞ.
Equation (23) is discretised in each cell: where the lithostatic potential at the cell center is an average of the nodes values: The discretised compaction law writes then:

Fault Models
In the considered geometric representation of a basin, faults are handled as two surfaces that can move past each other. Regarding fluid flow, the grid concept introduced in Section 1.2.2 captures the juxtaposition of distinct stratigraphic layers across a fault and its changes through time in a relatively continuous way. This is an important issue as it is often considered that this provides the basic plumbing of the system. However, a fault is not only a geometric discontinuity but rather a volumetric zone of variable thickness and with rock properties that result from the tectonic history of the basin. Fault zone thicknesses remain small compared to basin scale but their properties namely permeabilities can be very different from the surrounding host rocks and can evolve significantly over time. A fault zone can therefore act as a barrier or on the opposite as a conduit to fluid flow, becoming therefore a major feature of the plumbing system. In standard basin model, fault zones are represented by cells adjacent to fault surfaces to which the corresponding fault properties are assigned.
In order not to overestimate the fault zone thickness, local grid refinement should be used in the vicinity of the fault surface. Already difficult in normal setting, this approach is impractical in complex settings since it would require too much meshing effort. We therefore follow an approach compatible with the surface representation of faults in which the main features mentioned above can be taken into account. As for the presentation of the flow model, our aim here is not to discuss the models that can be used to compute the fault zone properties, using shale gouge ratio for instance (Manzocchi et al., 1999), but to propose a fault flow model in which these properties can be included as parameters.
For the derivation of the discretised set of equations in the presence of faults, we assume that the domain is cut by one single fault although the approach applies to more general fault network. Let us consider a particular fault defined by the boundary surface C 12 ðtÞ which corresponds in the mesh to two sets of faces S 1 ðtÞ and S 2 ðtÞ and the two sets of nodes N 1 ðtÞ and N 2 ðtÞ. Going back to the set of equations derived in the previous section, the fault impacts the sedimentary load equation for nodes that belong to the fault surfaces and the mass balance equation for the cells that are adjacent to the fault faces.
In the next sections, we first consider the impact of faults on the sedimentary load computation. Then, we describe the proposed hierarchy of fault flow models starting with the simplest one, which consists in specific boundary conditions applied at fault walls, and ending with the most complex one which considers flow along and across fault zones.

Sedimentary Load
For the computation of the lithostatic potential, we assume that the geometric surface representation of the fault is accurate enough and therefore that the lithostatic potential is continuous across the fault.
On the discrete point of view, the fault introduces additionalr unknowns for the nodes of N 1 ðtÞ and N 2 ðtÞ. We distinguish among these nodes the "bottom" ones (subsets N i;bottom ðtÞ; i ¼ 1; 2) for which discretisation of (22) gives an equation of the form (24). A node is a bottom one if the vertical ray issued from itself towards the top cut one of its neighboring cell. For the other nodes of N 1 ðtÞ and N 2 ðtÞ, we discretise the continuity condition using the node-to-face contact presented in Section 1.2.2.3.
Let us for instance consider a node s of: and d the face of S 1 ðtÞ such that the projection s 1 of s on S 1 ðtÞ belongs to d. The continuity condition at node s gives:r s ¼rðs 1 Þ where the value ofrðs 1 Þ is obtained through linear (or bilinear if the face d is a quadrangle) interpolation of the values ofr at the nodes of d (Fig. 11). For each node of N i;top ðtÞ, i ¼ 1; 2, this procedure gives an equation which can be written: where ðs; dÞ 2 C N;i;j is defined by (16). We can notice that this approach does not make any assumption such as "bottom to top" ordering of the fault surfaces.

Fault as Flow Boundary Conditions
The simplest fault flow model considers that fault surfaces are boundaries where specific flow boundary conditions can be defined. For instance, a fault acting as a perfect barrier can be modeled imposing a no flow condition on the fault surfaces:ũ w Áñ ¼ 0 on C 12 ðtÞ ð 28Þ On the opposite, a fault which acts as a very permeable conduit and which is connected to the surface of the sedimentary stack, can be modeled with a given overpressure boundary condition: These boundary conditions are discretised as standard boundary conditions.

Fault as a Membrane
Another approach for the fault flow model is to assume that a fault acts as a membrane which can retard or impede flow across its surface. In its simplest form, this model assumes that the fault acts only as a juxtaposition surface. In the continuous fluid mass conservation equation, the fault is then an interface across which overpressure and fluid mass flux are continuous. In the discrete model, these continuity conditions are taken into account when approximating the fluxes across faces that belong to S 1 ðtÞ and S 2 ðtÞ, and can be simply seen as the standard flux approximation across non-matching meshes.
In the present work, we consider only a two-point flux approximation for these fluxes. Although more accurate finite volume scheme can be considered, it captures permeable flow paths across the fault interface.
Let us consider a face d 1 that belongs to S 1 ðtÞ for instance and whose adjacent cell (as defined in the connectivity of the grid) is j 1 . We denote by Á d 1 ;2 the set of faces belonging to S 2 ðtÞ that are in contact with d 1 as defined in Section 1.2.2.3: where C F;1;2 is defined in (17). The flux across d 1 is decomposed into a sum of fluxes between j 1 and cells adjacent to the faces of Á d 1 ;2 .
F d 1 ;d 2 is then given by the two point flux approximation between cells j 1 and j 2 : where the transmissibility T j 1 ;j 2 is proportional to the contact area between the two faces (introduced in Sect. 1.2.2.3) and to the harmonic average of the cells permeabilities. Figure 11 Schematic fault of Figure 9, interface condition for the lithostatic potential.
A similar decomposition is performed for fluxes across faces of S 2 ðtÞ. The symmetry of the face-to-face contact and the flux continuity condition associated to the perfect membrane fault model leads to: This flux approximation is the standard flux approximation used in reservoir simulator for computing flux across fault surfaces (Manzocchi et al., 1999). The important feature here lies in the set Á d 1 ;2 which evolves through time and is therefore computed for each time. Following further the approach used in reservoir simulation, we can account for specific fault zone properties in the transmissibility computation. Indeed, fault properties are discretised on the fault surfaces by assigning a fault thickness and a fault permeability to each fault face that belongs to S 1 ðtÞ or S 2 ðtÞ. The transmissibility between the two cells j 1 and j 2 is then computed as if these two cells were separated by two fictitious fault cells (Fig. 12). This approach which requires fault properties defined on the fault faces is consistent with sliding along the fault surface since each surface moves past the other one keeping its properties description. Although out of the scope of this paper, we can point out that these fault properties (permeability and thickness) can be the result of specific models, which account for the neighboring cells properties and for the history of the fault movements.

Conductive Fault
This fault model intends to account not only for flow across the fault zone but also for flow along the fault zone (i.e. parallel to the fault surface) without introducing a specific volumetric grid. The approach is based on a double interface fault model which has been described in details in (Tunc et al., 2012). Its principle is to reduce the 3D flow model in the fault zone to two interface models, one for each side of the fault, where the fault thickness becomes a parameter. Investigated for a domain whose grid does not match across the fault but which does not evolve through time in (Tunc et al., 2012), the method finds an additional interest here as it is naturally designed to cope with fault surfaces that move one past the other through time.
The underlying idea can be schematically understood by introducing a fictitious volumetric mesh of the fault zone obtained by assigning a thickness to each fault face, as shown in Figure 13.
In the fault zone, we assume that fluid flow is governed by mass conservation of fluid coupled with Darcy's Law. We moreover assume that compaction is negligible which leads to a stationary diffusion equation for the overpressure unknown. Since the formal writing of the continuous "double interface" model is rather technical, especially for non planar fault surfaces, we step directly into the discretisation stage. Indeed, once discretised with a cell centered finite volume scheme (Tunc et al., 2012), the "double-interface" model leads to additional overpressure unknowns and mass balance equations in each fault face.
Each fault face d of S i ðtÞ, i ¼ 1; 2, is characterized by: -a thickness d f ;d , -a "along the fault" permeability K fl;d , -a "across the fault" permeability K fa;d . Let us denote oP d the overpressure unknown associated to the face d.    As for the fluxes on internal faces, we can use a twopoint or a multi-point flux approximation (Tunc et al., 2012). For the sake of simplicity, we present here only the two-point flux approximation for which we give the main idea rather than explaining the details of the transmissibility computations: -matrix/fault flux approximation: where T j i ;d is the fault-matrix transmissibility which depends on the permeability in j i , the coordinates of j i nodes, the "across fault'' permeability of d: K fa;d and the half fault thickness ð1=2Þd f ;d . To ensure conservation, this flux approximation is also used in the flux balance for cell j i ; -along fault flux approximation: where d c is the neighboring face of d in S i ðtÞ across the edge c; T d;d c is the "along fault" transmissibility which depends on the "along fault" permeabilities, the fault thicknesses of faces d and d c . It accounts for non planar fault surfaces introducing suitable normal vectors.
Although not detailed here, the along fault flux approximation considers also the boundary conditions at the boundaries of the fault surfaces; -across fault flux approximation: the flux between the fault surface S i ðtÞ and the other fault surface S j ðtÞ is decomposed into a sum of fluxes between the face d and the faces of S j ðtÞ which are in contact. As in previous section (30), we define the set Á d;j of faces that are in contact with d. The across fault flux then writes: where T d;d j is the "across fault" transmissibility which depends on the "across fault" permeabilities, the half fault thicknesses of the two faces and the contact area between them.
This model is particularly well adapted to fault movements since the fictitious fault cells follow the fault surfaces through time. As for the across fault flow model, the fault properties can be the results of appropriate models. It is moreover possible to distinguish a damage and a gouge zone within the fault zone, where the gouge zone properties are accounted for when computing the across fault flux permeability in a way similar to what is done for the "fault as a membrane model" (Sect. 2.2.3)

Solution of the Discrete Equations
We recall that this formulation is inserted in the approach of Section 1.2.1 where the evolution of the grid through time is a given input. An iteration of the "geometrical loop" corresponds to the computation of the main discrete unknowns on the time interval ÀT ; 0 ½ , which are obtained as the solution of the governing equations completed with initial conditions and further discretised in time. We use a linearly implicit Euler time discretisation that is not detailed here. In the sequel, we give an outline of how the discrete equations are solved at each time step.
For each time step t nþ1 ¼ t n þ Át (the superscript n and n þ 1 represent the time level), the grid M nþ1 ¼ M ðt nþ1 Þ at time t nþ1 is updated as explained in Section 1.2.2.2. Then time discretisation of the governing equations leads to a system of algebraic equations for the unknowns at time t nþ1 that are the solid volume, the porosity and the overpressure of each cell (and eventually the overpressure of each fault face) and the lithostatic potential of each node. The proposed formulation induces a sequential resolution that includes the following stages: -solid mass balance Equation (18) gives directly the solid volume Vols nþ1 j knowing Vols n j and the input solid deposition rate, for each cell j 2 M nþ1 ; -once the solid volumes are computed, the lithostatic potential of the nodes are computed solving an algebraic linear system that gathers Equations (24, 27) and the boundary conditions; -once the lithostatic potential are computed, we compute cell porosities, cell overpressures and eventually face overpressures (if the "conductive fault" model is chosen) that are solution of Equations (20), and fault equations that depend on the chosen fault behavior. Equation (32) allows to eliminate the porosity unknowns which gives finally a linear system in the overpressure unknowns. The simulator which implements the computation of the unknowns through time has been developed around Arcane Platform (Grospellier and Lelandais, 2009) and is parallelized using mesh partitioning in order to handle large grids (of the order of one million cells). A specific practitioner has been designed for evolving grids but its description is out the scope of this paper.

ILLUSTRATION ON A SEMI-SYNTHETIC EXAMPLE
This part illustrates the results obtained on a semisynthetic example based on an analogue experiment. It focuses on the impact of fault behavior for fluid flow. This experiment is inspired by geological structures with crustal extension, which cause the plate to split apart and center blocks to drop down relative to the flanking blocks. This feature is the illustration of a rift valley such as the East African rift where still active Arabian and African plates are splitting (Mougenot et al., 1986) or the subsurface structures of the North Sea area dominated by grabens formed during the Triassic and Jurassic periods (Glennie, 1998).

Analogue Experiment Description
Analogue experiment aims at reproducing real geological case at laboratory time and space scales (Ramberg, 1981). In this example, a medical scanner performing computerised X-ray tomography is used (Colletta et al., 1991). The principal advantage of this technique is to produce a real time 4D evolution. They are no unknown on the geometry through time. The model studied is built with an alternation of sand and pounce powder on a silicone basement in a 25 9 70 cm box (Fig. 14). These materials simulate the rheology of geological formation. Total thickness at deposition is 4 cm. The number of layers is 6. Sedimentation rate is 4 cm/h. The deformation is applied in only one-direction. At the end of the experiment, the model has a total length of 32 cm, so an extension of 28%. Normal faults delimit completely the downthrown blocks. Fault dipping is less than 45 degrees.

Geomodel Building and Gridding
First step is the interpretation of the scanned data and the building of the 3D surfacic model. Horizons and faults are picked on the 3D final scanned block imported into a geomodeler (Fig. 15). Around 30 faults were picked, but only 8 were selected. Only faults with significant throws were selected. Faults in X or Y are not considered, because the restoration prototype cannot handle them at the stage of the publication. The fault network description is simplified to allow the realization of the whole basin workflow that means the coupling between an evolving geometry and fluid flow through geological time steps. Figure 15a illustrates the interpretation and simplification done on the data. The model is gridded with an unstructured mesh based mainly on hexahedra. The grid is based on 6 720 cells at the last step of deformation (Fig. 15c). The created mesh conforms to faults and compartmentalizes the model into units which follows the geometry evolution due to deformation or fault activation. Each unit could be delimited by faults or not. The model is not divided into separated meshed blocks but is represented by an unique 4D grid. Each layer is representative of an homogeneous media made of either sand or shale. The mesh follows the stratigraphic deposition.
A quality control in thickness preservation was carried out by comparing the interpreted and built model with the scannered data block. Some parts were re-built and re-meshed. Thickness variations through time is checked and preserved as much as possible.

Structural Restoration
The next step is the building of the kinematics or the different stages of deformation to explain the geometry at present day. This step is solved with a volumetric restoration approach (Durand-Riard et al., 2013;Moretti et al., 2006;Maerten et al., 2001). It is a sequential approach describing the deformation of the model for each stratigraphic layer. At each time step, the target geometry is the top horizon of the structure. Each top surface is restored at its respective depth observed on past scanned steps. Each restoration step corresponds to a paleomodel built at deposition time step. At each time step, the restored layer is removed. Faults are surfaces along which sliding occurs. Some boundary conditions are defined on the model, for example, the left border is fixed during the restoration scenario to ensure coherency of the unfolding in the direction of extension and limit the number of degrees of freedom. The 9 successive steps of restoration are presented in Figure 16, they have been obtained with Kine3D-3 (Moretti et al., 2006).
The result of this step is a time-dependant unstructured grid suitable as input for the forward simulation. This 4D grid reproduces properly the different steps of deformation of the model and the sliding along the faults. In Figure 16, pillars preserve the main direction of the dipping surfaces. The mesh is non-matching along the faults and stratigraphic discontinuities across the faults are modeled. Porosity reduction due to compaction processes is not taken into account in this semi-synthetic example, because the rheology of materials is not varying in this experiment. There is no correction needed for calibrating the change in porosity due to geometry evolution and sediment rates

Forward Simulation and Impact of Faults on Fluid Flow
For the last part of the workflow, which corresponds to the coupling between the deformed geometry through time provided by the volumetric restoration and the fluid flow simulation, particular attention is paid to the impact of faults on fluid flow paths and therefore on overpressure distribution. In this part, each coordinate of the grid has been multiplied by a factor of 200 000 in the z-direction, 50 000 in the x-direction, and 100 000 in the y-direction to produce an example with more realistic dimensions (Fig. 13a). The lithofacies distribution is such as sand is given a reservoir lithological properties and pounce powder some seal properties. Intrinsic permeability follows a Kozeny-Carman's law. Permeability contrasts between sand and pounce powder is 10 4 .
To illustrate the fault hydrodynamics, two different scenarios are considered.
In the first one, faults are considered as membranes. Geology adjacent to fault surfaces is the main driving factor for existing fluid flow distribution. The fault property emphasizes the barrier or conduit behavior across its surface. In the second one, faults are considered as conductive faults. Fluid migration is across the fault zone but also along the dip direction of fault surfaces. In both scenarios, no temporal property for fault permeability is given. Geometric and lithological properties for faults are summarized in Figure 17 and Table 1. Figure 18 illustrates the fault behavior of a membrane. Faults 2, 3, 4 act as a barrier and other fault surfaces as a conduit. For all fault surfaces, stratigraphic contacts across the fault surfaces are similar, with an alternation of sand and shale. In Figure 18b   Scenario 1: fault as a membrane 10 m K = 10 À2 mD K = 0 mD Scenario 2: conductive fault 10 m K = 10 À2 mD K = 0 mD Figure 19 compares the fluid flow direction and velocity in function of the parameterization of the faults. Fluid moves essentially by lithologic contrasts in the units. Due to faults 2, 3 and 4 which are impermeable, fluid is migrating along fault 2 and the drainage is important in the space between faults 2 and 3. Fluid velocity and rate are more important in Figure 19a compared to Figure 19b everywhere. The plumbing system is less active across the faults in Figure 19b, but very active along the dip direction of fault surfaces as illustrated in Figure 19c. In the scenario 2, fluid velocity is more important in the zone along the faults than across the fault zones (Fig. 19c). This behavior is due to the fault permeability which is of the same order of magnitude as the permeability of the sandy layer and is emphasized by the open boundary condition at the top of the faults. The scenario 2 is able to capture across and along fluid movement within the volumetric fault zone.
Figures 20 and 21 demonstrate the impact of the fault model on overpressure distribution through time. The fault parameterization affects fluid flow patterns and is a controlling factor for pressure calibration. In Figure 20, overpressure increases continuously through time to reach 10 MPa in the deepest part of the model. At the opposite, in the scenario 2, overpressure is less important and reaches 5 MPa in the deepest part of the model. The overpressure is mainly concentrated in the shaly layers. Figure 21 shows the overpressure evolution through time in two specific cells, for the two scenarios. We can see that overpressure evolution depends heavily on the juxtapositions that occur across the fault, although the conductive property of the fault for scenario 2 tends to smooth this behavior.

Conclusion on this Example
With this analog experiment, the workflow based on the new developed simulator has been evaluated. This example is a demonstration of the complexity in combining structural modeling, fault analysis and petroleum system analysis. This example illustrates the potentiality of the model to allow for different fault parameterizations that cover a large range of behaviors, going from a perfect barrier to a conduit, using a single evolving grid. Step 1/8 Step 3/8 Step 6/8

CONCLUSIONS
In this paper, we have described the basis of a new forward basin simulator designed to better take faults into accounts. We have introduced an innovative 4D grid concept which is a real step forward compared to previous approaches. It is indeed flexible enough to follow continuously basin deformation in complex settings including displacement along fault surfaces. Besides, the simulator contains a hierarchy of fault flow models. The most complex one is the conductive fault model. It allows accounting for fluid flow along the fault zone coupled with relative displacement of the fault walls, although the fault is only represented as interfaces in the 3D grid. Results on a semi-synthetic test case illustrate the capabilities of the new simulator. Although not shown here, the simulator accounts for thermal transfer and hydrocarbon generation and it is currently being extended to hydrocarbon migration. a) The two cells, close to fault 7, for which pressure evolution is plotted through time. The bottom one is a shaly cell which the upper one is a sandy cell. b) Overpressure evolution for the two selected cells, the blue curves correspond to scenario 1 while the red dashed curve correspond to scenario 2.