An adaptive pluri-Gaussian simulation model for geological uncertainty quantification

1. Introduction

For the implementation of any assisted history matching (AHM) model in the area of reservoir engineering, one of the most important topics to address is the reservoir geology. Usually, the initial knowledge of the geology is very uncertain and the aim of the AHM exercise is to reduce this uncertainty while preserving key features such as for example facies proportions. In this study, we are interested in the facies distribution and try to reduce its initial uncertainty using additional information given by the production data when available. In order to quantify the geological uncertainty, a geological simulation model has to be firstly defined. From this model, one generates possible facies realizations in the reservoir domain. These facies instances must be consistent with the prior knowledge of the reservoir geology (number of facies types, facies distributions, possible transition between facies types, direction of the facies, geometry of the facies, etc.). The second part consists of conditioning these facies realizations to production data from a Bayesian point of view. The literature on Bayesian history matching methods is vast (Aanonsen et al., 2009, Oliver and Chen, 2011, Chen and Oliver, 2013, Emerick and Reynolds, 2013, Stordal and Lorentzen, 2014, Stordal, 2015). Here, the focus is on improved geological simulation using the truncated plurigaussian simulation (TPS) model. (Astrakova and Oliver, 2014). The TPS method generates facies realizations by combining Gaussian random fields with a truncation map. The truncation map is either user defined based on prior knowledge (Liu and Oliver, 2005) or is defined through an internal construction of the model (Sebacher et al., 2013). The TPS method is a natural generalization of the Gaussian truncation method (Matheron et al., 1987) where a single Gaussian field was truncated using some thresholds defined on the real axis. The TPS method was first introduced by Galli et al. (1994) and further developed by Le Loch et al. (1994) and Le Loch and Galli (1997). In Le Loch et al. (1994), two Gaussian random fields were truncated using a truncation map defined in the bi-dimensional real space by the intersection of some vertical and horizontal lines. These lines split the space in rectangular sub-domains and a facies type is assigned to each sub-domain. Using different geostatistical properties of the Gaussian fields, different realizations of facies distributions were generated, corresponding to different characteristics of the facies types. Le Loch and Galli (1997) investigate the relationship between the variogram of the Gaussian fields and the indicator variogram of the facies for case of stationary proportions (stationary TPS). They also present a conditional simulation technique for the non-stationary case of the proportions of facies (non-stationary TPS). Important prior information about facies fields is related to the possible contacts among different facies types. Xu et al. (2006) addressed this issue through a binary dynamic contact relation matrix in the context of the pluri-Gaussian simulation model. One of the challenges of the TPS methodology is the conditioning of the facies fields to hard data (facies observations) collected at the exploration wells (Lantuejoul, 2002, Armstrong et al., 2011). Practical issues for implementation of the TPS methodology conditioned to hard and soft data were investigated in Emery, 2007, Emery, 2008 and computer programs were proposed for conditioning the facies simulation and validating of the algorithms. Its easy implementation within an inverse modeling process has made this methodology to be used either in synthetic fields (Le Ravalec-Dupin et al., 2004) or even in real fields (Deraisme and Farrow, 2004).

The coupling of the truncated pluri-Gaussian simulation method with an ensemble based history matching method was first introduced by Liu and Oliver (2005), where the facies boundaries were automatically adjusted when the Gaussian random fields were updated with new measurements. The authors used both the ensemble Kalman filter (EnKF) and the randomized maximum likelihood as HM methods and compared the performance of both. In this study, the EnKF gave better results. For the facies simulation model, the truncation map used was built with three lines, whose intersections in a two dimensional plane generated seven regions, each having a facies type assigned. The same truncation map with EnKF as the HM method was used by Agbalaka and Oliver (2008) in a three dimensional model, with three facies types (each pair could be in contact), and using in addition, a distance-based localization scheme for the estimated covariance matrix. After the assimilation of the production data, the Gaussian field values are changed and, by truncation, yield facies fields not necessarily having consistent facies observations. An additional iterative procedure was introduced to correct this misfit. For the facies observation operator, a proxy function was used. This function has only two values: 0 if the facies type is as observed and 1 if the facies type is not as observed. Zhao et al. (2008) integrate the truncated plurigaussian simulation model, introduced by Liu and Oliver, into an iterative EnKF process, assimilating production data in addition to seismic data, for a three dimensional reservoir model with three vertically uncorrelated layers. Agbalaka and Oliver (2009) applied EnKF to facies fields that exhibit non-stationarity in facies proportions whilst simultaneously updating the petrophysical properties of the facies. Sebacher et al. (2013) used a probabilistic operator for the facies observation, and the truncation map of TPS model was introduced based on a maximization criterion. The authors defined the concept of ”probabilities fields” to estimate the binary fields of each facies type. In addition, the truncation map parameters were introduced in the state vector and estimated for a better uncertainty quantification of the facies proportions. Astrakova and Oliver (2014) applied the TPS methodology coupled with the LM-EnRLM with an interior point method suitable of inequality constraints. The interior point formulation for the objective function, used in the LM-EnRLM, makes it possible to maintain the correct facies types at the well locations.

The idea of truncation of Gaussian variables has been used in the context of multi-point geostatistical simulation model (MPS, Caers and Zhang (2004)) in Sebacher et al., 2015, Sebacher et al., 2016 where random fields, only marginally Gaussian are defined and truncated using different truncation maps. These random fields were used for parameterization of complex channelized reservoirs generated with MPS in combination with a training image. The random fields do not have a Gaussian field structure, being defined in each grid (marginally) by sampling from standard Gaussian distribution. The truncation maps are defined based on the probability fields of the facies calculated from an ensemble of realization generated with MPS. These random fields are built by random sampling from conditional Gaussian distribution and are not simulated like in standard TPS. Moreover, their variogram structure is inherited from the training image and is, in general, not Gaussian.

Here is assumed that initial probability fields of each facies type can be a priori designed by a group of experts in the early phase of the reservoir exploration, using various data, such as core information, seismic data, well logs data, outcrops information, etc. In addition to the probability fields, extra information could be available for constructing facies maps, for example the possible contacts among the facies types, the expected facies proportions (global), the geometry of the facies and the relative position among the facies types. All this information is part of the expert knowledge of the reservoir and is used to construct a geological simulation model. The main idea is to generate an ensemble of geologically plausible facies fields, for which the probability fields of the facies types estimated from the ensemble are consistent with the probability fields specified by the experts (Armstrong et al., 2011).

In the stationary TPS, a single truncation map is defined for all cells. In this study we develop a geological TPS simulation model with increased flexibility, allowing for the definition of a truncation map for each grid cell individually. This model generalizes the non-stationary TPS (Armstrong et al., 2011) that is based on proportion matrices which are not conditioned to facies observations. The proposed model is based on probabilities that can be conditioned to facies observations. The simulation of the facies fields is performed using spatially correlated uniform random variables having spatial correlation by the means of marginal transformations of random Gaussian fields (each component of the Gaussian field is transformed to a variable with a uniform distribution). Here we use the probability integral transform (Casella and Berger, 2002) that provides a way to link any continuous random variable with the uniform random variable defined on (0,1). For each grid cell, the uniform random variables define the facies type in that grid cell based on a partition (simulation map) of the ${(0,1)}^{q}$ hypercube, where q is a priori chosen based on the number of facies types and connections. In our approach, we define the simulation map (truncation map) as a suitable partition of the ${(0,1)}^{q}$ hypercube contrary to truncated plurigaussian simulation method where the truncation map is defined as a partition of the multi-dimensional real space (with dimension q). The simulation map is adaptive in the sense that its shape depends on the values of the probability fields of the facies provided by the experts.

Because the prior facies probability fields incorporate the facies observations at the well locations (in probabilistic terms), the simulations always provide the correct facies types at the observation points without the need for constrained optimization. Therefore, within a HM process throughout the assimilation period, the facies observations are kept for all ensemble members, contrary to the approach reported in Liu and Oliver, 2005, Agbalaka and Oliver, 2009, Sebacher et al., 2013, Astrakova and Oliver, 2014. The geological simulation model is called the adaptive plurigaussian simulation model (APS) due its ability to adapt the shape of simulation map with position in the reservoir domain. The initial ensemble of facies fields is generated using stationary Gaussian random fields defined on the reservoir domain with a user defined variogram model. This approach differs from all other approaches (involving the use of TPS methodology for stationary or non-stationary case) where the Gaussian fields are initially constrained to facies observations leading to a reduced initial variability and possible changes of the original geostatistical properties (non-stationary Gaussian fields). In addition, this approach differ from Sebacher et al. (2016) because here APS is defined as geological simulation model while in Sebacher et al. (2016) the geological simulation model is MPS and the truncated random fields have been used only as a parameterization.

It will be shown that the APS methodology provides a better prior uncertainty quantification than the stationary TPS methodology when conditioned to facies observations.

In the next section, the APS algorithm is introduced and some of its properties. Section 3 presents a comparison of the proposed method with the stationary TPS algorithm. Section 4 presents the implementation of the proposed method with the adaptive Gaussian mixture filter (AGM). Section 5 contains a case study and the paper ends with conclusions.

2. The adaptive pluri-Gaussian simulation model

Let us assume that a prior description of the subsurface geology for a given reservoir obtained from different sources (seismic surveys, core interpretations, outcrops, etc.) is available and includes the number of the facies types and the possible contacts (transition) the facies types. Several methods have been proposed to estimate probability fields for each facies occurrence by combining seismic data, well log facies analysis and statistical rock physics (Avseth et al., 2001, Massonnat et al., 1999, Ng et al., 2008). However, these probability fields describe only marginal properties (no information of the conditional probabilityfor neighboring grid cells). The dependence structure is frequently described by variogram models, however, it is not easy to combine the two sources of information for non-Gaussian random fields. Furthermore, the probability fields are typically not conditioned to the reservoir production outcome. In the following, a simulation algorithm is described that honors both the dependence structure described in a variogram and the marginal probability fields. The algorithm makes it easier to combine different sources of expert knowledge into simulations without adding much complexity (e.g., no need for constrained optimization or complex conditional simulation algorithms).

2.1. Facies simulation map

To give a clear description of the APS method, the methodology is presented for a reservoir model where three facies types are present (denoted $F_{1}$ , $F_{2}$ and $F_{3}$ ). In Fig. 1 we present an example of three possible probability fields, one for each facies type, that could be given to the user. At each location (grid cell), the sum of values of these probability fields is one, a condition fulfilled by construction (irrespective of the methodology). In each field, the probability is 1 at the well locations (or any other location where a facies observation exists), where the associated facies type occurs. This means that the facies observations are already incorporated into the probability fields. Additional important a priori information is the transition (contact) between any two facies types. For illustration, let us consider the case where any two facies types may be in contact each other.

Denote by $p_{k}$ the probability field associated to facies type k (where $k = 1,2,3$ ), the mapping of which is presented in Fig. 1. For each grid cell j, consider the categorical variable, denoted $ℱ^{j}$ that expresses the distribution of the facies types at that location $ℱ^{j} \sim (\begin{matrix} F_{1} & F_{2} & F_{3} \\ p_{1}^{j} & p_{2}^{j} & p_{3}^{j} \end{matrix}) .$

Here $p_{1}^{j}$ , $p_{2}^{j}, p_{3}^{j}$ , are the probabilities at the grid cell j obtained from the probability fields (marginal probabilities at grid cell j). One way of generating a facies type at grid cell j consists of generating a sample from this distribution. This can be done as follows.

1.
Consider a domain A of measure (or area) 1, in a multidimensional real space (in general A could be the hypercube ${(0,1)}^{q}$ ). For three facies types, A will typically be the unit square ${(0,1)}^{2}$ from $ℝ^{2}$ (Fig. 2).
1. Download : Download high-res image (113KB)
2. Download : Download full-size image
Fig. 2. The facies simulation map layout.
2.
Split the domain into three subdomains where subdomain k has measure (area) $p_{k}^{j}$ . This assigns a facies type to each subdomain with the corresponding probability.
3.
Generate an element from a uniform distribution defined on A. In $ℝ^{2}$ this is given by $α^{j} = [α_{1}^{j}, α_{2}^{j}]$ where $α_{1}^{j}$ and $α_{2}^{j}$ are uniform random variables(possibly dependent) defined on the interval $(0,1)$ . Assign the facies type depending on the subdivision of A to which the generated element belongs.

If $α^{j}$ is in the subdomain k, the facies type $F_{k}$ is assigned to grid cell j. The decomposition of the unit square ${(0,1)}^{2}$ represents the facies simulation map of the grid cell j. Each grid cell can have its own simulation map, including facies probabilities (represented by proportions occupied in the map), by varying the partitioning of the domain into its subdomains. At a location where the facies is observed, its simulation map consists of the square ${(0,1)}^{2}$ , occupied by the observed facies type. In addition, there might be some regions of the reservoir domain where only two facies types occur (in terms of probabilities only two probabilities are greater than zero). In this case, the simulation map contains only two subdomains, each with an area equal to the associated probability of the facies type that occupies that subdomain.

In the following section, is shown how one can easily define a simulation map for each grid cell and how the uniform variables are chosen so that the simulated facies fields have a continuous structure.

2.2. Simulating spatial correlated fields honoring marginal probabilities

There are two situations where the facies simulations will lack continuity. First, the simulation maps of the neighboring grid cells are too different from each other (have different topological structure. Secondly, the samples $α^{j}$ and $α^{j + 1}$ in the neighboring grid cells j and $j + 1$ are independent. The solution for the first situation is to design a layout for all simulation maps. The layout represents a generic decomposition of the domain A, defining in addition, a parameterization of the curves that delimitate the subdomains. For the case presented in the previous section, the simulation map layout is designed as in Fig. 2, from which it can be seen that the parameterization is given by the equations of the horizontal and vertical lines. From this figure, it can also be seen that the simulation map depends on two values, $p_{2} / (1 - p_{1})$ situated on the $O α_{1}$ axis and $p_{1}$ situated on the $O α_{2}$ axis. Consequently, the simulation map of the grid cell j is built using the values $p_{1} = p_{1}^{j}, p_{2} = p_{2}^{j}$ in the layout.

To overcome the second situation, the uniform random fields $α_{ℓ}, ℓ = 1,2$ need to have a spatial correlation. The easiest way is to link these uniform random fields with the Gaussian random fields. The theorem that tell us how to connect a uniform random variable to a Gaussian variable is the probability integral transform (Casella and Berger, 2002). The probability integral transform states that, if X is a continuous random variable with the cumulative distribution function $F_{X}$ , then the random variable $F_{X} (X)$ has a uniform distribution on $(0,1)$ . The cumulative distribution function links any continuous random variable with the uniformly distributed random variable on $(0,1)$ . Therefore, theoretically, the uniform random fields can be defined as piece-wise projections of any spatially correlated fields (defined on the reservoir domain), for which the distribution is continuous, through their cumulative distribution function. In this study, the uniform random fields are defined as piece-wise projections of Gaussian fields because of three important reasons; first, generating Gaussian random fields is easy, second, the correlation structure can be described using only two point statistics (variogram, covariance function etc) and third, the ensemble based history matching methods works well with the Gaussian assumption of the model parameters.

Consequently, if we assume that Y is a (stationary) Gaussian random field defined on the reservoir domain with the marginal $Y^{j}$ having the distribution $N (m; σ)$ , then one can define the uniform random field α in grid cell j as $α^{j} = Φ (Y^{j})$ (i.e. piece-wise projection of the Gaussian field Y with the cdffunction $Φ$ ). Here, $Φ$ is the cumulative distribution function (CDF) of the $N (m; σ)$ distributed variables. Then, $α^{j}$ has a uniform distribution with support $(0,1)$ .

The APS algorithm (for three facies types) is summarized as follows.

1.
Given the prior probability fields $p_{1}$ , $p_{2}$ and $p_{3}$ and information of facies connections, create the layout of the simulation map ( $A_{1}, A_{2}, A_{3}$ where $\cup A_{k} = {(0,1)}^{2}$ , Fig. 2) (we will give later on examples with different simulation maps).
2.
Generate samples from two (stationary) Gaussian random fields, $Y_{1}$ and $Y_{2}$ , with predefined dependence structure.
3.
Transform the Gaussian random fields to uniform random fields, $α_{1}$ and $α_{2}$ , using the integral transform. For each j from 1 to $N_{grid}$ , $α_{1}^{j} = Φ_{1} (Y_{1}^{j})$ and $α_{2}^{j} = Φ_{2} (Y_{2}^{j})$ , where $Φ_{1}$ and $Φ_{2}$ are the marginal Gaussian CDF's of $Y_{1}$ and $Y_{2}$ .
4.
For each j from 1 to $N_{grid}$ , built the simulation map of the grid cell j from the layout and set the facies type in grid cell j, to $F_{k}$ if $α^{j} = [α_{1}^{j}, α_{2}^{j}] \in A_{k}^{j}$ .

For any ensemble based history matching algorithm, the Gaussian random fields can be used as parameters that are updated during the assimilation. Regardless of method, due to the construction of a simulation map in each grid cell, facies observations will be honored during assimilation contrary to implementation of HM methods with TPS.

In Fig. 1 are shown the a priori facies probability fields and in Fig. 3 the facies probability fields estimated from an ensemble of 120 realizations generated with the APS. Here one can observe the main result of the APS; in case of generating an ensemble of sufficient samples of two stationary Gaussian fields, the facies probability fields calculated from the ensemble are consistent with the given facies probability fields.

The APS model can be generalized to an arbitrary number of facies types and to any type of contact between facies types. To do so, a reliable decomposition of the ${(0,1)}^{2}$ square (layout) needs to be defined. If the connectivity is very complex the decomposition may be constructed in higher dimensions, using for example the domain ${(0,1)}^{3}$ and three Gaussian fields for the simulation of the facies fields.

2.3. Impact of the simulation map

In the pluri-Gaussian context, the geological acceptability of the simulated facies fields is achieved by a reliable choice of the geostatistical properties (e.g., variogram) of the Gaussian fields and of the simulation map layout. Usually, the geostatistical properties of the Gaussian fields are decided on the basis of the prior topological characteristics available for the facies types: the geometry of the facies types, the number of facies (rock bodies) that occur in the domain and the orientation of the facies. In the following, is presented an example in which small changes of the simulation map layout, yield different topological characteristics for the simulated facies fields, when keeping the geostatistical properties of the Gaussian fields unchanged. The contact among facies types can be handled by a reliable definition of the simulation map. However, the relative positions among facies types cannot be handled with any simulation map. In Fig. 4(top) are shown four possible realizations of the facies distributions, for reservoirs characterized by the existence of three different rock types, where any two could be in contact with each other. The green facies type is characterized by long anisotropy of which angle with the horizontal direction is obtuse. In contrast, the red facies type is characterized by a large amount of small bodies within the reservoir domain. The blue facies type is the geological environment in which the other two facies types occurred due the depositional system. These geological realizations were generated using the same values of two stationary Gaussian random fields, the first is anisotropicwith long range correlation of 50 grid cells, short range correlation of 5 grid cells and with the principal direction of $120^{\circ}$ with the horizontal. The second Gaussian random field was generated isotropically with a small correlation range of 5 grid cells. The marginal distribution of each Gaussian field is $N (0; 1)$ and the covariance type is also Gaussian. The parameterization and the layout of the simulation maps are presented at the bottom of Fig. 4 and of course, these maps are not unique representations of the contact relationships. The APS model was applied, using uniform probability fields of 0.4 for the facies type one (blue) and type two (green) and respectively 0.2 for the facies type three (red). Even though there are small differences among the simulation map layouts, the facies fields provided exhibit different topological properties. In particular, the relative position of the red facies type with respect to the other facies types. In the first simulation, bodies of the red facies type can be inside both the green and blue facies types. In the second example, bodies of the red facies type may be found within the blue facies type, whereas the green facies type intersects the red facies type only on its border. In the third simulation, the deposits of red facies type are encountered mainly between the borders of the other two facies types, eroding the blue facies type. In the fourth simulation, the relative position of the red facies type with respect to the other facies types exhibits similar characteristics, remaining mainly on the edges of the borders. In this formation, the red facies type has eroded both other facies types. The only difference in these simulations is the layout of the simulation map (the values of the Gaussian fields were the same), Consequently, the design of the simulation map layout and the choice of the geostatistical properties of the Gaussian fields must be done by a group of experts, in a multidisciplinary (geology, mathematics, reservoir engineering, etc.) manner in such way that the simulation results reflect as closely as possible the prior information. As previously discussed, the number of facies is arbitrary although the algorithm is explained for three facies types. Fig. 5(c) presents some random realizations of four facies types with initial probability fields (Fig. 5(b)) and the simulation map layout (Fig. 5(a)). Facies type one is plotted in dark blue, type two in light blue, type three in orange and type four in dark red. Note that facies types one and three are not connected Sebacher et al. (2015).

3. APS versus stationary TPS

In this section is shown why it is important and helpful to have prior probability fields of facies occurrence when one wants to use the plurigaussian technique for simulation of the facies fields in the presence of hard data (facies observations).

The stationary pluri-Gaussian simulation (TPS) model has two main inputs. The first input consists of Gaussian random fields defined on a discrete reservoir grid. The second input is a single truncation map defined on a multidimensional real space. The truncation map is a partition of the multidimensional real space and is defined by intersection of some curves that divide the space into regions, with each assigned a facies type. As shown below, the TPS methodology can be implemented using a procedure similar to the one used for APS. Let us consider for simplicity the case with two Gaussian fields $Y_{1}$ and $Y_{2}$ following, at each location, a normal distribution with CDF $Φ_{1}$ and $Φ_{2}$ . Then, the function(1) $F : ℝ^{2} \to {(0,1)}^{2}, F (x_{1}, x_{2}) = (Φ_{1} (x_{1}), Φ_{2} (x_{2}))$ is a bijection and has increasing components. Consequently, it defines a bijection between the truncation maps designed in the space $ℝ^{2}$ and the simulation maps defined in the space1 ${(0,1)}^{2}$ . The truncation map parameters are calculated from the expected facies proportions (Armstrong et al., 2011). It is easier to define the simulation map in the space ${(0,1)}^{2}$ , simply because is more convenient to use subdomains with predefined areas instead of solving integral equations. Consequently, the APS can be viewed as the TPS model conditioned to probability fields of the facies type (soft data) or as a method that incorporates the probability fields in the truncated plurigaussian context. In addition, the APS with uniform probability fields is the stationary TPS. The first non-stationary version of the TPS was introduced by Le Loch and Galli (1997)and is described in detail in Armstrong et al. (2011), but there TPS is conditioned to non-stationary proportion curves for facies in different areas of the reservoir domain. The proportion curves of facies (or 3D matrix of proportion curves for the 3D case) are not conditioned to facies observations, but the probability fields (or cubes for the 3D case) of facies could be conditioned to facies observations.

In the stationary TPS, where facies observations are available at some locations, the Gaussian field values are generated such that the simulated facies fields satisfy the observations (interval conditioning). This implies a change in the mean and variance functions of the Gaussian fields compared with the case where the facies observations are not present (i.e. the Gaussian fields are no longer stationary). We will show that this modification produces a bias for the probability field of each facies type calculated from an ensemble of realizations.

Consider a rectangular domain with $50 \times 100$ grid cells with three facies types present. One facies type exhibits a long correlation from west to east (denoted facies type two) and a short correlation from north to south. Another facies type (denoted facies type three) is characterized by bodies of rocks that occur most likely on the edge of the other facies and the last facies type is the geological environment (denoted facies type one). The TPS methodology is applied for creating realistic facies distributions, using two Gaussian fields. Initially, both Gaussian fields are stationary with each marginal having a standard normal distribution and having a Gaussian covariance type. The first Gaussian field is anisotropic, modeled with a long correlation range of 50 grid cells and a short correlation range of 5 grid cells and with principal direction being the horizontal direction. The second Gaussian field is isotropic with the correlation range of 5 grid cells. The truncation map used for generating the facies fields is shown in Fig. 6, in which the line parameters are calculated from the expected facies proportions $p_{1} = 0.425$ , $p_{2} = 0.425$ , $p_{3} = 0.15$ . This truncation map is obtained by applying the inverse of function F (Eq. (1)) to the fourth simulation map in Fig. 4 (bottom right). Consequently, the thresholds $α_{1}, α_{2}$ and β (Fig. 6) calculated from the inverse of function F are $α_{1} = Φ_{1}^{- 1} (p_{1} + \frac{p_{3}}{2} - \frac{\sqrt{p_{3}}}{4}),$ $α_{2} = Φ_{1}^{- 1} (p_{1} + \frac{p_{3}}{2} + \frac{\sqrt{p_{3}}}{4}),$ $β = Φ_{2}^{- 1} (2 \sqrt{p_{3}}) .$

Let us consider thirteen locations at which the facies types are observed. The positions and the facies observations are presented in Table 1. The approach from Zhao et al. (2008) was used in order to generate conditional facies fields (to observations). First, at observation grid cells, pairs of Gaussian values are generated so that those pairs yield correct facies observations in accord with the truncation map. A Gibbs sampler can also be used as in Armstrong et al. (2011), considering correlations between Gaussian field values at the well locations. Thus, because the long range correlation of the first Gaussian field is 50 cells and the shorter distance, in the horizontal direction, between two wells with facies type two observed is 40 cells (Fig. 7, center) independent values for the Gaussian fields at the well locations were generated. Of course, if many facies observations are situated within the correlation influence of one or both Gaussian fields then the values of Gaussian fields (at those positions) are correlated, and a Gibbs sampler (Armstrong et al., 2011) should be used. Secondly, the simulated Gaussian values are used in a conditional sequential Gaussian simulation process in order to populate values to the remaining grid cells. With this method, we generate an ensemble of 120 samples of pairs of Gaussian fields. For all the experiments, the Gaussian fields were generated with the sequential Gaussian simulation method implemented in SGeMS (the Stanford Geostatistical Software, Remy 2005). These samples are truncated with the truncation map from Fig. 6 in order to generate 120 different facies fields with the correct facies values at observed locations. The conditioning procedure of the Gaussian fields destroys their stationarity resulting in the mean function of the Gaussian fields are not equal to 0 at each location, and consequently the simulated facies fields are, possible, no longer conditioned to the expected facies proportions.

Table 1. The facies observations.

x coordinate	5	5	5	25	25	50	50	50	75	75	95	95	95
y coordinate	5	25	45	15	30	5	25	45	15	30	5	25	45
Facies observations	2	2	1	3	2	2	3	1	1	2	2	1	3

In Fig. 7 are shown the probability fields of the facies types obtained from the ensemble of facies fields generated by the stationary TPS method and conditioned to the facies observations presented in Table 1. These probability fields are biased in the sense that the geostatistical properties of the Gaussian fields significantly influence their spatial distribution. For instance, the top of probability field of facies type two shows a long correlation in the horizontal direction. This is due to the position of three observations of facies type two within the correlation range of the first Gaussian field. The distance between the observation grids is 40 grid cells whereas the long correlation range of the Gaussian fields is 50 grid cells. This causes the bias in the probability fields, even though the correlation influence is negligible. The bias increases as the number of facies observations within the correlation range increases.

The APS methodology is applied using stationary Gaussian fields having the same geostatistical set up as before. In order to do that, prior probability fields of the facies occurrence that incorporate the facies observations are needed. As this is synthetic example, probability fields have to be created avoiding the strong bias observed from the ensemble obtained using stationary TPS. Using an empirical approach a set of ”prior” probability fields of facies is generated that incorporate facies observations (top of Fig. 8). These fields are generated with a short correlation around the observed locations together with the property that the mean of each probability field is equal to the expected proportion of the associated facies type. This property ensures that the simulated facies fields are conditioned to the expected facies proportions (Appendix A). The layout of the simulation map used in the APS is the last map presented in Fig. 4 (bottom right). An ensemble of 120 samples of pairs of stationary Gaussian fields are generated and use the APS, the simulation map layout and the probability fields to generate the ensemble of facies fields. The probability fields for each facies type are shown at the bottom of Fig. 8, from which one can clearly see that the prior probability fields are very well preserved and the geostatistical properties of the Gaussian fields have little influence.