Simulating Non-Stationary Seismic Facies Distribution in a Prograding Shelf Environment

Résumé — Simulation de la distribution de faciès sismiques non stationnaires dans un environnement de plate-forme progradante — Un réseau dense de profils sismiques du delta sousmarin du Rhône en Méditerranée est utilisé pour construire un modèle stochastique d’un environnement de dépôt de plate-forme progradante, analogue à un réservoir pétrolier. L’altitude des limites supérieures et inférieures des unités sédimentaires est tout d’abord estimée par pointage sur les profils sismiques. La géométrie complète dans l’espace des enveloppes des sept unités stratigraphiques identifiées est ensuite estimée par une méthode géostatistique non stationnaire, après l’inférence de leurs fonctions de covariance généralisée par la méthode des incréments. Un modèle architectural complet en 3-D est alors construit en remplissant l’espace compris entre ces enveloppes par la simulation stochastique de la répartition des faciès internes. Ces faciès sont définis par leur signature sismique, caractérisée par la configuration du signal réfléchi, sa continuité et son amplitude. Le programme géostatistique HERESIM est utilisé en mode non stationnaire pour analyser puis simuler la variabilité verticale et horizontale des faciès sismiques. Cette variabilité verticale est mesurée pour chaque unité par une courbe de proportion verticale, tandis que la variabilité horizontale est analysée par les variogrammes horizontaux de faciès, qui mesurent la fonction d’autocorrélation d’un faciès donné en fonction de la distance. Les autres paramètres du modèle HERESIM sont calés pour reproduire les données sismiques observées sur les profils. Les résultats des simulations stochastiques sont utilisés pour interpréter les mécanismes de dépôt, d’érosion et de déformation tectonique qui ont construit cet environnement de plate-forme.


INTRODUCTION
For the purpose of fluid flow simulation, it is important to have "good" reservoir models representative of the distribution of the petrophysical heterogeneities.When petrophysical variables are linked with lithofacies, the development of reservoir models may require modelling of the lithofacies distributions.This is an area where geostatitical techniques are increasingly employed, in particular in the oil industry, see e.g.Goovaerts (1997) or Le Ravalec-Dupin (2005).These techniques have the ability to produce equiprobable 3D images of the reservoirs based on the input constraints of the model while respecting the existing data, which directly affect reservoir continuity and geometry.Moreover the geostatistical approach allows the association of an uncertainty with the obtained results.
We present here a methodology developed and tested on a case study (the underwater Rhone River delta) at the French Institute of Petroleum (Chihi, 1997) to characterize the architecture and to model in 3D a deltaic sedimentary sequence, which can be considered as an analogue of a hydrocarbon reservoir.The modelling is based only on high-resolution seismic data.Although some superficial borehole data were available (e.g.Gensous and Tesson, 1996;Rabineau et al., 2005Rabineau et al., , 2006)), which could have been used to better calibrate the upper portion of the seismic lines, they do not sample the major stratigraphic unit of interest in this article.We therefore decided to limit the data set to the seismic profiles to test the applicability of this methodology to the initial geophysical survey of a potential reservoir.Very recently, Bassetti et al. (2008) reported on a 100 m long core in this area, with sediment dating, which will be used for comparison in the results of simulation section.
The three main steps of the methodology can be summarized as follows: -constitution of the database.The type and nature of the data will condition the choice of the regionalized variables to define and which ones will be most likely to create the desired reservoir architectures; -the second is the construction of the 3D geometrical model.We first model the spatial variability of the geometry, through its structural analysis, for the different regionalized variables that will be used to model the stratigraphic units.This allows the best estimation of the surfaces of each sedimentary unit with minimum uncertainty.The superposition of the different units to construct the 3D geometric model of the whole sedimentary sequence is straight-forward; -the last step is the simulation of the facies within each unit.It requires specific technical approaches.We have to take into account the data type (only seismic sections), the information content of stratigraphic geometries as portrayed on seismic lines and the geostatistical concepts (matrix of proportions, reference level, etc.) available to provide realizations that reproduce the data at the seismic section locations.In this paper, steps 1 and 2 will be briefly summarized, and can be found in more detail in Chihi (1997), and Chihi et al. (2000Chihi et al. ( , 2007)).Thus step 3 will be the main focus.

GEOLOGY
A grid of shallow high-resolution seismic profiles was established on the continental shelf of the Rhone River delta on the French Mediterranean coast.The grid was created using a 50 joule SIG Minisparker ™ and a resolution of about 1.0 m with a penetration depth of 170 m (200 ms two-way traveltime) below the sea floor.Several transects showed the regional stratigraphic relationships within late Quaternary deposits.The survey consisted of lines covering a surface area of 5226 km 2 on the western portion of the Rhone delta, of which 82.5 km run in an East-West direction (X) and 63.35 km in a North-South direction (Y) (Fig. 1).The data contain 15 seismic reflection profiles, 9 oriented NE-SW and 6 others oriented NW-SE.The profiles are regularly spaced with a mesh of about 6.6 km.
A sequential analysis shows that the sediments within the Rhone delta are organized into a system of stacked seaward prograding wedge-shaped units (seven seismic units labelled in descending order from U1, the upper unit, to U7, the deeper unit) (Fig. 2).Each wedge is composed of seaward prograding clinoforms, and is bounded by a basal downlap surface and an upper toplap surface.The geometry and stratal patterns of the wedges result from a combination of high frequency eustatic cycles and episodic shelf tilting during the terminal phase of the Quaternary (Gensous and Tesson, 1996;Vella and Proveansal, 2000).Tesson et al. (1993), Gensous and Tesson (1996) have proposed that the prograding wedges were constructed during the last fourth-order glacio-eustatic cycles (Würm) of the Pleistocene.The individual wedges present high-frequency parasequences formed during fifthorder sea level changes.The differential subsidence caused the wedge form of the stratigraphic units.The total subsidence allowed for an aggradational pattern within this sequence.Glacio-eustasy combined with subsidence (i.e., relative changes of sea level) and sediment flux variations produced cycles of shoreline migration across the shelf.These cycles of sedimentation were punctuated by erosional bounding surfaces or discontinuities (Tesson et al., 2000).
The wedges U2, U4, U5 and U7 are interpreted as prograding wedge-shaped units of regional extent.Each wedge tapers landward, onlapping onto a mid-shelf unconformity, and thickens seaward, reaching a thickness of up to 50 m at the shelf edge.They accumulated during a period of relative sea level lowering (lowstand system tracts).Internal down-shift surfaces represent successive downward shifts of coastal onlaps.They subdivide a progradational wedge into a number of successive lozenge-shaped subunits and indicate that progradation occurred under conditions of "forced regressions" (Posamentier et al., 1992).The toplap surface of the wedges is composed of two types of surfaces: -a subaerial erosional surface associated with relative sea level fall during progradation (lowstand surface of erosion); -an erosional flooding surface during subsequent relative sea level rise (ravinement surface).
The wedges U3 and U6 are interpreted as intercalated seismic units of limited areal extent.They are the result of transgressive erosion and longshore transport of underlying lowstand wedge deposits during episodes of relative sea level rise that occurred between periods of wedge progradation: transgressive system tracts (Gensous et al., 1993;Gensous and Tesson, 1996).The erosional discontinuity at the top of wedge U2 (surface S1, Fig. 2) marks the maximum lowstand.The overlying postglacial (early to middle Holocene) deposits (unit U1) constitute the transgressive and highstand systems tracts.

GEOMETRIC MODELLING
A geometric model was built by first mapping the boundaries of each sedimentary unit without considering its filling with facies: this mapping is based on the measurements of the topographic elevation Z of the unit boundaries on each seismic profile.We used all the seismic profiles already interpreted (Tesson et al., 1993;Tesson, 1996;Chihi, 1997) and based our work on the elevation of the 8 identified surfaces.We recorded the depths (measured in seismic time in milliseconds) every 15 minutes of the boat travel time on each profile, and then processed these records to transform them (by linear interpolation) into samplings with a constant   spatial step over each profile.A total of 39 748 data points were extracted from the seismic profiles for the 8 surfaces.The rather dense sampling along the profiles is meant to adequately collect the information given by the survey, without taking a quasi continuous sampling, which would make the calculations cumbersome.The geometrical modelling was then made in two main steps: -the geostatistical structural analysis of the "Depth" variables using the generalized covariance; -the estimation of the surface boundaries using kriging.

Variographic Analysis of the "Depth" Variables
A geostatistical analysis of the data was performed for each variable "Depth": S0 to S7, in order to fit the best model of variogram or covariance for estimation of the surfaces.The "Depth" experimental variograms show in general a parabolic behaviour (Fig. 3).The variable "Depth" is a non-stationary random function with an obvious anisotropic drift in both the NE-SW and the NW-SE directions, and resulting in the parabolic form of the variograms.This non-stationary behaviour becomes increasingly marked from the sea bottom to the deepest surfaces.The a priori variance generally increases with increasing depth of the surfaces.The variograms grow, reach the variance for the different variables, and continue to grow beyond the variance.We will try to correctly represent the experimental behaviour of the variogram mostly in the zone of the distances that will be used later to estimate the "Depth" by kriging, i.e. with a moving neighbourhood of 25 km.
These variograms cannot be adjusted by a theoretical model involving stationary or intrinsic random functions.We therefore used the intrinsic random functions of order k (IRF-k) approach with the "method of increments" (Chilès and Delfiner, 1999).This involves three stages: -calculation of the increments of order 1 of the variable Y(x) = Z(x + h) -Z(x); -calculation of the experimental variograms of the increments of order 1 of the variable Z; -determination of the generalized covariance of the variable Z by the relation that exists between the variogram of the increments of order 1 and the generalized covariance (Chilès, 1999): where γ(h) is the variogram of the increments of order 1, K(h) is the generalized covariance and C n m is the number of combinations of m numbers taken n × n.This expression can be generalized for increments of a higher order.
To determine the generalized covariance of the variable, we first selected an arbitrary model of generalized covariance, the program deduces the theoretical variogram of the increments and we compared it to the experimental variogram of the increments.Thus, we adjusted the generalized covariance of the variable by successive trials.The method was used on all the variables (S0-S7).For each variable "Depth", the model of the generalized covariance of the variable is composed of the sum of three structures: a spherical one, a Gaussian, and a generalized covariance of order 3. K(h) = σ 1 + σ 2 + alhl 3 where σ 1 is the model of the spherical covariance, σ 2 is the model of the Gaussian covariance, and alhl 3 is the term of the generalized covariance of order 3. Table 1 summarizes the different model parameters obtained for the variables.Figure 4 shows, as an example, the variogram of the increments of surfaces S1 and S6.

Calculation of the Surfaces by Kriging
The estimation of the surfaces requires three stages: -definition of the estimation grid; -definition of the neighbourhood; -kriging of the surface, using the ISATIS geostatistical software (Bleines et al., 2000).
The estimation grid of the studied domain was defined so as to include all the surfaces: S0-S7 and to minimize to the utmost the extrapolation at the corners of the grid.Its origin along X is 270 km and along Y: 4761 km; cell dimension: dx = dy = 0.5 km and number of cells: nx = 166, ny = 128, which results in 21248 cells covering the entire domain.The selection of the neighbourhood was subjected to a certain number of constraints so as to ensure the continuity of the estimator, take the drift into account and to achieve maximum precision, with the following criteria: -the neighbourhood is divided into eight octants centered on the node of the grid being estimated; -the data points used as a neighbourhood are selected so as to have the same neighbourhood density in all of the octants to ensure neighbourhood continuity; -the number of points in each octant and the maximum distance are large enough to allow the selection of the neighbourhood among a sufficient number of seismic profiles in order to account for the drift; -in order to better account for the drift, a minimum distance "d min " is prescribed between two data points used as a neighbourhood; -in order to take into account the anisotropy, a maximum distance "d max " is used: all the samples located on the maximum distance of the ellipsoid are considered to be at the same distance from the point to estimate.
No sample located outside the maximum distance of the ellipsoid is ever selected (for further details see (Chihi, 1997;Chihi et al., 2000)).These criteria led to the selection of d min = 5 km, d max = 25 km, number of points taken as a neighbourhood in each angular sector = 4 and number of angular sectors = 8.
The variables "Depth" were estimated by kriging in IRF-k using the models determined by the method of increments (Tab. 1) and the moving neighbourhood defined above.Figure 5 shows, for instance, the maps of surfaces S1 and S6.
Figure 6 displays the final geometrical reconstruction with the stacking of the seven units.
The maps and the 3D visualisation are tools which allow a better understanding of the precise succession of sedimentary sequences as determined by the seismic data and thereby clarifying the deposition mechanisms and the successive transgression-regression cycles (Chihi, 1997;Chihi et al., 2007;Tesson, 1998).

FACIES SIMULATION
After the envelopes have been reconstructed, the corresponding units can be filled with internal facies distribution.A complete 3D architectural model is then produced.This filling with facies is done using the HERESIM geostatistical software, which allows the vertical and horizontal facies variability to be analyzed and simulated (Matheron et al., 1987;Rudkiewicz et al., 1989): -the vertical variation of facies is analyzed by means of vertical proportion curves of facies; -the horizontal variation of facies is analyzed by means of horizontal proportion curves and horizontal variograms of facies, which measure the auto-correlation of a given facies as a function of the distance.These parameters are used as input constraints for simulating the sedimentary units using the method of Truncated Gaussian Functions.
The input data used for describing the units, that will be simulated, are derived from the seismic sections.They are generated in order to be adapted to the applied simulation.The seismic sections are considered as important factors in controlling the lateral and vertical heterogeneity, and thus for the building of the internal unit architecture (Doligez et al., 1999(Doligez et al., , 2003)).Therefore the data set has to be analyzed (type of available data, spacing of the data), described and calibrated on a consistent basis.

Data Analysis
In this case study, the available data are seismic, we do not have well data (cores and logs) that describe all of the seven units, even if some core data are available for the upper units (Gensous and Tesson, 1996;Rabineau et al., 2005Rabineau et al., , 2006)).We therefore have to describe and quantify the internal facies variations in terms of seismic facies.A seismic facies consists of a group of seismic reflections whose character differs from adjacent groups within a sequence (Mitchum et al., 1977).
Generally, when seismic sequences are studied, the configuration and character of reflections within sequences are studied in an approach called "seismic facies analysis".In seismic facies analysis, it is assumed that there are direct, if not unique, relations between the position of a seismic facies within a sequence, the external form and internal configuration of its reflections, and the sedimentary environments represented by the facies.Using this assumption, one may plot seismic facies on a seismic section and transfer them to a map view, to be interpreted in the same way as sedimentary facies information and for further analysis and interpretation (Cross and Lessenger, 1988).Sangree andWidmier (1977, 1979) and Michum et al. (1977) suggested that some combinations of the following attributes might be used to distinguish one seismic facies from another: -reflection configuration and continuity; -amplitude and frequency spectra of reflections; -interval velocity; -geometrical relations of reflections to the unconformities bounding the sequences; -three-dimensional form of the seismic facies.
Two commonly occurring reflections have been recognized (Fig. 2, 7, and 21): -parallel patterns, in which reflections are parallel with respect to underlying and overlying reflections; the clinoforms present a weak slope; -progradational or clinoform patterns, in which reflections are inclined with respect to underlying and overlying reflections; the clinoforms may present both moderate and strong slope; -reflection continuity, reveals the continuity of the sedimentary beds from which depositional processes can be inferred.All the seismic sections on the Rhone delta show a good continuity of reflections, thus this characteristic was not used to distinguish one seismic facies from another; -amplitude of reflections, reveals the spacing of the sedimentary beds.Three commonly occurring reflections have been recognized: low, moderate and high amplitudes (Fig. 7).These parameters were evaluated visually on the seismic sections.Reflection characteristics that were used to distinguish one seismic facies from another are summarized in Table 2. Figure 7 presents the different facies that were identified and classified in terms of seismic facies attributes.
For the purpose of this work, the domain to be simulated was defined in such a way that all the available stratigraphic units (or seismic sections) are within the domain (Fig. 1).
Thirty-four fictitious wells were extracted by analyzing the seismic characters within the defined domain and used as ''real well'' data in order to perform the facies simulation (Fig. 8).
The seismic facies were coded (Tab.2) and digitized into files and stored in the database containing: X and Y location of each fictitious well, the seismic facies codes and their depth (Z: ms).
The structural analysis and the simulations were performed for the seven units, but in this article, only one unit U5 is presented and used to illustrate the different stages of the facies simulation.Interested readers can examine the entire study in Chihi (1997).

Structural Analysis
The geostatistical simulation starts with a structural analysis performed for each unit: -proportion curves are calculated in order to provide information about the average facies distribution in various directions, usually vertical and horizontal; -variograms for the different facies are calculated and analysed.Horizontal variograms of facies measure the auto-correlation of a given facies as a function of the distance.

Vertical Proportion Curves
The vertical proportion curves of facies quantify the organization of each sedimentary unit.Each unit is subdivided vertically into successive levels ("relative" mesh) and the proportion of each facies is computed along lines parallel to the chosen reference level.The result is a graph where the proportion of each facies is displayed on the horizontal axis, each level being parallel to the reference level.
The origin of the vertical axis corresponds to the position of the reference level of the unit under study.The definition of the reference level used as a correlation horizon (Fig. 9) is closely related to the stratigraphic model, the depositional system and the subsidence which prevail in the study area.For each unit, the reference level was chosen so that the horizontal and vertical proportion curves show sequential evolutions, in terms of progradation and aggradation.Inventory of seismic facies.
Figure 8 Seismic facies identification, X and Y location of each fictitious well (see Fig. 1 for location).Facies number (F3, F6, F9) are defined in Table 2.
Unit U5 is dominated by seismic reflection geometries dipping seaward and bounded by horizontal to irregular erosional surfaces (Fig. 2): -the upper boundary is commonly flat and dipping gently seaward (S4), whereas the lower boundary (S5) is generally irregular and locally convex upward (when the progradational unit U5 overlies the intercalated unit U6) (Fig. 2, 10); -the seaward terminations of internal reflections are downlapping at a low angle or appear to be tangential to the lower boundary for a long distance.These reflections are toplapping, at a higher angle on the upper boundary.For theses reasons, the most appropriate reference level, at the scale of unit U5, is the upper boundary (S4).

Horizontal Proportion Curves
The horizontal variation of facies is analyzed by means of horizontal proportion curves.The curves of lithofacies proportion were calculated for the whole unit.They showed that the distribution of facies is non-stationary: the spatial distribution of facies varies systematically from one zone to another.For the structural analysis, we will thus have to test non-stationary geostatistical simulation algorithms and to calculate lithofacies proportion matrix to take account of this non-stationarity.

Matrix of Lithofacies Proportions
For each unit, the spatial variation of facies is quantified by means of a proportion matrix.The matrix of proportion determines the way of integrating the data to produce an adequate regionalization (representing the reality) of lithofacies proportions.Two methods (Mouliere et al., 1996) can be used for this calculation: -proportion matrix calculated by areas; -proportion matrix estimated from vertical proportion curves (VPC).
The choice of the method plays a very significant role for the simulation results, i.e. the images of the facies distribution.Indeed according to the studies carried out on the various stratigraphic units, the choice of one method rather than another depends on the succession of the lithotypes within the wells and on the lateral variation of facies.

Proportion Matrix Calculated by Areas
The first step of this approach consists in defining a matrix of lithofacies proportions.This matrix is built on a grid, with an elementary cell size of 2000 m × 2000 m (Fig. 11), coarser than the simulation grid (500 m).Then each unit is subdivided into zones (Fig. 11) where the facies variability is considered stationary.The lateral extension of each zone is based on seismic facies analysis, which is carried out by pattern recognition within the different wells, and guided by the seismic sections.For each zone, the vertical proportion curve is computed with the well data located in the considered zone.Finally for each zone (Zi) we affect the computed VPCi to each cell of the proportion matrix.Figure 10 Thickness map of unit U5.
-the west-central zone characterized by only one facies: WS-MA (sky blue); -the eastern zone characterized by two facies (1) WS-MA (sky blue), (2) MS-MA (pink); -and the southern zone at the centre of the proportion matrix, characterized by three facies: (1) WS-MA (sky blue), ( 2) subchaotic (yellow) and ( 3) MS-MA (pink) representing the lowest proportion.This method gives a proportion matrix where the abrupt variation of facies can be seen from one vertical proportion curve to another (Fig. 12).For example, starting from the north-western zone, the proportions of seismic facies SS-HA (brown) and MS-HA (orange) decrease suddenly to the south and to the east, and are replaced by the MS-MA facies (pink).Then the proportions of MS-MA facies (pink) disappear abruptly to the west-central zone where the proportion is 100% WS-MA facies (sky blue).

Proportion Matrix Estimated from Vertical
Proportion Curves Location map showing the well positions and the vertical proportion curves calculated for each area.Grid discretization adopted for computation of VPC and for facies simulation.
similar geological depositional sequences can be grouped in order to compute a VPC from their lithological description).
For each cell of the proportion matrix, a VPC is estimated by ordinary kriging (Mouliere et al., 1996) using the previous local VPC.This method gives a proportion matrix with vertical proportion curves, where all the facies existing within the stratigraphic unit are present with a certain proportion (Fig. 13) that varies according to the range of the variogram.As the kriging is carried out with a unique neighbourhood, each facies can be found at any place in the area.This method does not completely respect the zonation of the facies observed on the wells and even less on that of the simulation images that we compared with the seismic profiles.To solve this problem, we estimated the proportion matrix by combining the previous two methods: "by area" and "from Vertical Proportion Curves".

Proportion Matrix Estimated by Area and from Vertical Proportion Curves
Several initial local VPC can be computed in the different zones and the estimation of the VPC in each cell of the matrix can be made zone by zone.For each group of wells characterized by the presence of a certain number of facies, we calculated a vertical proportion curve which we duplicated on the zone where the corresponding wells are distributed (Fig. 14).These proportion curves are considered as ponderators when kriging the vertical proportion curves on the whole volume to be simulated.Thus in each cell we have a vertical proportion curve with a dominant proportion reserved for the facies really existing in the well data characterizing the zone; and a very small proportion of facies reserved for the neighbouring zones (since the kriging is performed with a unique neighbourhood) (Fig. 15).The 3D matrix computed in this way respects the general characteristics of the initial proportions, and the zonation of the facies distribution is completely respected.Being able to obtain this results is very important for producing "realistic" simulations.

Simulations
After this geostatistical analysis, the calibrated parameters (proportion matrix, variogram ranges) were used to simulate the distribution of facies inside each sedimentary unit using a non-stationary model.The simulation of the facies on the volume of each unit consists in allocating a facies to each cell, so that the well data are respected.An algorithm based on the truncated Gaussian approach (Matheron et al., 1987;Matheron, 1988) is used, and makes it possible to obtain a series of equiprobable images of the studied units.

Parameters Used in the Simulations
The simulated volume (37000 m × 40000 m × thickness of each stratigraphic unit) was built with the data from the 35 fictitious wells.The grid for the 3D simulation was defined on 74 × 80 × n (a number depending on the thickness of the stratigraphic unit) cells of size: x = 500m, y = 500m, z =1ms.

Results of the Simulations
Several simulations were performed with the proportion matrix calculated by the different methods presented above.For each case, the simulation results (2D sections, such as Fig. 16, 17 and 18) were compared with the observed data (vertical seismic sections): -Figure 16 shows the facies simulation of unit U5 using the proportion matrix calculated by zone (Fig. 12).This method yields simulations where the lateral facies variation is abrupt.We notice for instance on the 2D section that the passage from the seismic facies MS-HA (orange) to the facies MS-MA (pink) is not progressive; -Figure 17 shows the facies simulation of unit U5, using the proportion matrix calculated from VPC (Fig. 13).This method gives simulations where the zonation of the facies distribution is not completely respected: we notice the presence of subchaotic facies (yellow) distributed within the facies WS-MA (sky blue).This facies is introduced because kriging is carried out with a unique neighbourhood;  Location map showing: 1) the fictitious well positions; 2) the vertical proportion curves calculated for each area; 3) the vertical proportion duplicated curves within the different areas.
-Figure 18 shows the facies simulation of unit U5 (Fig. 15), using the matrix proportion calculated by the combined method (by zone and from VPC).This method gives simulations where (1) the zonation of the facies distribution is completely respected: we do not observe the presence of subchaotic facies (yellow) within the facies WS-MA (sky blue), (2) the lateral facies variation is progressive from a seismic facies MS-HA (orange) to MS-MA (pink).
We must stress here that, due to the use of non-stationary simulation algorithms with correctly inferred parameters, we are able to reconstruct "realistic shapes".The simulation results are satisfactory from a geological point of view: they represent realistic facies distributions and reproduce the data of the seismic sections.With these facies simulations, a graphical display of the results will provide the necessary resolution for a more detailed and sophisticated geologic interpretation of the seismic data: -with 2D visualizations, a map of seismic facies is a 2D grid where each cell may correspond to the facies distribution at a given level or to the thickness of a given facies as presented in Figure 19; Cross-section showing the simulation results of unit U5, using the proportion matrix calculated by area.
Figure 17 Cross-section showing the simulation results of unit U5, using the proportion matrix calculated from VPC.
Figure 18 Cross-section showing the simulation results of unit U5, using the proportion matrix calculated "by area and from CPV".
Figure 19 Map of seismic facies WS-HA, thickness (ms) of unit U5.
-with a 3D visualization, carried out with the GOCAD software (ENSG, 2009) (Fig. 20), the display confirms and clarifies more accurately the geometry of each stratigraphic unit and the build-up mechanism "by prograding" of the sedimentary materials at the end of a fluvial system (the Rhone River), which is strongly influenced by tectonic and eustatic factors.
The facies spatial distribution within unit U5 exhibits a complex architecture.In this section, we try to interpret it by analyzing seismic sections, depth and thickness maps and results of seismic facies simulations.The facies spatial distribution is the result of several factors including at least: -nature of the sedimentary material; -cycles of shoreline migration; -subsidence rates; -the inherited morphology of the lower bounding surface of progradation.
The seismic facies common in stratigraphic unit U5 is dominated by oblique-tangential reflections, chaotic reflections are observed locally on the south side of the area (Fig. 20,21).
Detailed inspections of seismic data and 3D visualization of facies simulation reveal that the facies distribution or "organization" may be explained through the description of vertical facies evolution and lateral facies variation.

Vertical Facies Evolution
The oblique-tangential reflections dip seaward from top to bottom in this unit.They dip at a high angle at the point of toplap and at a low angle at the point of downlap.The reflection amplitude also tends to diminish from the point of toplap to the point of downlap.The oblique-tangential reflections are observed to be truncated in the upper part (Fig. 2, 8).The specific sedimentological facies provided by a main source, the Rhone River mouth, are clays grading upward to massive marine sands organized into bottomsets and foresets.The topset bed reflections are absent.They indicate erosion of the unit during relative sea level fall (subaerial erosion) and subsequent relative sea level rise (erosional shoreface processes).The specific organization showing a coarsening-up evolution is pronounced or maintained by the forced regressive process.We observe in Figure 21 that each downward-shift surface is immediately followed by clinoforms with steeper foresets than those immediately subjacent to downward-shift surfaces.It is important to note that the coarsening-up evolution is revealed by the Vertical Proportion Curves (Fig. 11).The progressive evolution from the Weak Slope-Moderate Amplitude facies (sky blue) to Strong Slope-High Amplitude one (brown) is well reproduced by the simulation (Fig. 20).The hypothesis of clays and sand sedimentological facies was advanced by Gensous and Tesson (1996).It is based on seismic and core data analyses of unit U2 (unit f in Gensous' and Tesson's terminology) and on a comparison with other studies of progradational deltaic deposit analogs (Gensous et al., 1993;Gensous, 1994;Gensous and Tesson, 1996;Tesson et al., 2000).

Lateral Facies Variation
Oblique-tangential reflections are present and are spreading out over the entire studied area.However, it is important to note that clinoform slopes (at the point of toplap) are greater in the northern region than in the eastern and western ones.Moreover, in the eastern and western areas, the clinoforms show weak seaward dipping angles and extend a significant distance subparallel to the lower bounding surface.The change of the clinoform slopes at the point of toplap depends on the amount of erosion.Erosion is more or less pronounced and cannot be determined with precision.Gensous et al. (1993) deduced from core data analysis that erosion remained moderate and that only the massive beach sands, deposited above the fair weather wave base, are absent.The isopach map (Fig. 10) describes the areal distribution and the rate of preservation.The stratigraphic unit appears well preserved and continuous, with an axis approximately parallel to the shoreline, in the north and the south.The maximum thickness is in the north-eastern and the southeastern part of the studied area; whereas in the central part, preservation seems weaker.The slope and the extension of the tangential clinoforms on the lower bounding surface are influenced by the sedimentological facies (clays) and by the slope and the topography of the inherited depositional surface (sea bottom).The sedimentary material tends to fill the available space on the depositional surface.A regional seaward tilting has provided the necessary accommodation space for deposits to accumulate.The depth map of unit U5 (Fig. 5) shows that the intensity of tilting increases progressively toward the shelf edge as well as eastward.This may contribute to determine the staking pattern of prograding wedges and the facies configuration: -in the northern part of the study area, unit U5 appears well preserved (Fig. 10, 20) and continuous.The amount of erosion is not too significant so the clinoforms are preserved and present a strong slope at the point of toplap.
Tangential reflections are present but are reduced in extent.This seismic facies may be interpreted as sand-rich; -in the western part of the study area, preservation seems weaker (Fig. 10).This suggests significant erosion of the upper parts of unit U5.Only the tangential patterns of clinoforms are preserved.An important point to note is that the clinoforms, which show weak seaward dipping angles, extend for a significant distance subparallel to the lower boundary.These long tangential segments are interpreted as clay-rich; -in the eastern and south-eastern part of the study area, unit U5 is well preserved (Fig. 10, 20) and continuous.
The lower boundary dips seaward at a higher angle than in the western part of the studied area (Fig. 5).Oblique tangential clinoforms are better preserved; they present a stronger slope at the point of toplap.The tangential portions of clinoforms extend over a significant distance subparallel to the lower boundary.They show moderate seaward dipping angles (Fig. 20, 21).Locally, tangential portions of clinoforms show weak slopes where the underlying unit is better preserved and the inherited morphology of the lower bounding surface of progradation dips gently seaward (Fig. 21); -the southern part of the study area is locally characterised by the presence of chaotic facies (Fig. 20, 21) which suggest sedimentation at the front of a prograding delta affected by mass wasting and a substrate instability area.(Tesson, 1996;Tesson et al., 2000).These results can be compared with the findings by Rabineau et al. (2005Rabineau et al. ( , 2006) ) who used very high resolution seismic data and shallow cores, and by the new deep core (100 m) presented very recently by Bassetti et al. (2008).Although their study site is located further south-west and is related to the Bourcart and Hérault canyons, and not to the Rhone delta, profile P_1052 of Rabineau (Fig. 6, 2005or Fig. 5, 2006) is approximately a continuation of our Profile 10 (Fig. 2).It is possible to relate our Unit U5 to Sequence S2 and Unit U60 of Rabineau et al. (2006); in Bassetti et al. (2008), the extent of our Unit U5 is very limited and corresponds to their unit U57 limited by discontinuities D45 and D30 (~5 m) (Rabineau's and Bassetti's notations are identical).Their findings are that the seismic facies are entirely consistent with the observed core data.Each major sequence S corresponds on average to a 100 ky glacioeustatic cycle, separated by discontinuities D which are polygenic surfaces.Bassetti et al. (2008) recognized, for each sequence, three major facies, from top to bottom: -massive fine to medium-grained sands, with low lamination angles and possibly Swaley cross-stratification, interpreted as a proximal foreshore or upper shoreface; -well-sorted fine sands with Hummocky cross-stratification interpreted as an intermediate lower shoreface, with thin inter-bedded mud; -bioturbated mud with interstratified thin sand beds (storm deposits), corresponding to the distal upper offshore, with a clear erosional base.These three major facies seem to correspond well to our strong slope, moderate slope and weak slope seismic facies, respectively, and therefore to provide a realistic sedimentological basis for our seismic facies classification.

CONCLUSION
In this paper we have addressed two problems encountered when reconstructing stochastic complete 3D models of reservoir architecture: (1) dealing with non-stationary phenomena in geostatistics and (2) using as input data only seismic measurements.
-geometric modelling: when the variable has a significant drift, it may be impossible to calculate the variogram directly.We have shown that one method that can be used is the Intrinsic Random Function of k order: (IRF-k), and the structural function that replaces the variogram is then the generalized covariance; -the method of increments used here to infer this generalized covariance is more general than the traditional universal kriging method, because universal kriging can only be used for variables with a stationary behaviour over a certain subdomain of the studied zone, which is not always the case; -stochastic simulations of seismic facies offer the added advantage of integrating non-stationary facies distributions: lateral and vertical facies variability, through the 3D proportion matrix.This tool enabled us to use the HERESIM software non-stationary simulation option and to provide realizations that are able to reproduce the data of the seismic sections with the general non-stationarity of the facies distribution; -several methods are available to compute the 3D proportion matrix.The choice of the most appropriate computation method depends on the facies variability and even more on the complexity of the facies architecture: -The computation method by area can give good results when the transition between zones is abrupt; -The computation method from the Vertical Proportion Curve seems more valid when the transition between zones is progressive.However, since the kriging is carried out with a unique neighbourhood, some cells are given a "facies value" from the neighbouring area.This facies will be spread (with a small proportion) in different areas of the studied domain all over the section; -These different methods may be mixed to fill the 3D block of proportions.The heterogeneities are then more accurately simulated, the lateral facies transitions are progressive and the zonation of facies distribution is then completely respected.
Figure 1Location map of seismic data.

Figure 2
Figure 2 Seismic section 10 (see Fig. 1 for location) illustrating the superposition of the stratigraphic units and the internal seismic reflections.(The units are noted U1 to U7 and the bounding surfaces S0 to S7).
Figure 3 Variograms of variable "Depth": S0 to S7. Horizontal axis: distance in km; vertical axis: variogram γ of "Depth" in milliseconds seismic time squared; D1: NE-SW direction; and D2: NW-SE direction.Horizontal dashed line is the variance of field.

*
The anisotropy gives the multiplying coefficient of the range corresponding to direction D1 (NE-SW) or D2 (NW-SE); for example, (1,1) means that the range is the same in both directions; (1.3,1) means that the range is multiplied by 1.3 in D1, and by 1 in D2.

Figure 4
Figure 4 Variograms of the increments of surfaces S1 a) and S6 b).
Figure 5 Depth maps of surface S1 a) and surface S6 b).
Figure 12 illustrates the proportion matrix of unit U5 computed by zones.It is possible to observe the non-stationary distribution of the facies with the appearance of four distinct zones on the proportion matrix: -the north-western zone characterized by four seismic facies (Tab.2, Fig. 7) (1) WS-MA (sky blue) representing the highest proportion, (2) SS-HA (brown) and (3) MS-HA (orange) in equal proportions and (4) MS-MA (pink) with the lowest proportion;

Figure 9
Figure 9 Definition of the reference level.a) Geometrical model of correlation in parallel, b) geometrical model of correlation in proportional.

H
Figure 12Proportion matrix of unit U5 computed by areas.
Figure 13Proportion matrix of unit U5 computed from vertical proportion curves.

H
Figure 15Proportion matrix of unit U5 computed by area and from vertical proportion curves.
Figure 20 3D visualisation of the simulated facies within stratigraphic unit U5.

Figure 21 Tangential
Figure 21Tangential reflection slopes are parallel to the inherited surface, (see Fig.1for location).

TABLE 1
Model of Variable Depth for Each Surface