In order to perform calculations with qualitative properties it is necessary to transform them into numerical values beforehand. This can be done using the indicator function: the indicator of a given lithotype is equal to 1 when the observed point belongs to this lithotype and 0 otherwise. There are as many indicators as there are lithotypes involved in the deposit description, which turns the qualitative properties into a multivariable numerical setup. Moreover, as the indicators are numerical variables, it is now possible to consider their spatial characteristics through traditional tools such as simple and cross-variograms. The indicator variogram measures the probability that two points do not belong to the same lithotype (Eq. 1) as a function of their distance:

The indicator cross-variogram measures the probability that two points belong to two different lithotypes (Eq. 2) as a function of their distance:

As the lithotypes constitute a partition of the space, when the indicator of a given lithotype is equal to 1, the indicators of the other lithotypes are equal to 0. This leads to particular relationships relating simple and cross-variograms of the indicators.

For instance, in Fig. 1a, a simulation has been performed using an object based model (Boolean simulation): in a white background, elongated grey objects are overlaid by black circular objects. Fig. 1b represents the simple variogram of each indicator calculated in north–south and east–west directions, and the cross-variograms between any pair of indicators (Fig. 1c). The variogram of the black lithotype indicator is isotropic (same shape in all directions), while the variogram of the grey lithotype indicator is anisotropic: the dark curve (NS) stabilizes at a larger distance than the light curve (EW). The variogram of the white indicator is slightly anisotropic as complementary of the other two. The three cross-variograms (Fig. 1c) exhibit combinations of the different anisotropies. In order to produce a simulation that restitutes all these characteristics, it is important to use a consistent multivariate model that fits all these curves simultaneously.

The truncated Gaussian model is an elegant solution for simulating a set of facies which corresponds to a consistent multivariate model for their indicator variables. All the lithotypes are obtained simultaneously from one normalized underlying GRF: they are simply obtained by thresholding this GRF (Armstrong et al., 2011).

The GRF displayed as a map (Fig. 2a1) is composed of values simulated according to a normal distribution (Fig. 2a2) and is characterized by its spatial structure (isotropic cubic model) summarized by its variogram (Fig. 2a3). In the truncated Gaussian framework, one interval on the GRF corresponds to a probability (hence a proportion) derived from the Gaussian cumulative density function. Conversely for a given proportion, the width of the interval (hence its bounds) depends upon the location of the interval along the cumulative density function. Fig. 2b and Fig. 2c represents two intervals for the black lithotype corresponding to the same proportion (∼18%). These two constructions are summarized in a diagram, the lithotype rule (Fig. 2b1 and Fig. 2c1), where the Gaussian values are conventionally represented along the horizontal axis. The interval corresponding to the black lithotype is represented along this axis. Applying these two different thresholds to the same outcome of the GRF (as shown in Fig. 2a1) and with the same proportion, we obtain two black lithotype outcomes with different spatial structures: compact bodies for the first interval (Fig. 2b2) and elongated ones for the second one (Fig. 2c2).

This difference is also visible on the variograms of the indicators (Fig. 2b3 and Fig. 2c3) (the maximum distances for the variogram computation correspond to the quarter of the field size). An interval located close to the center of the Gaussian distribution leads to shorter range for the indicator variogram (Fig. 2c3) than when the interval concerns the lowest (or the largest) Gaussian values (Fig. 2b3): the lithotype changes at short distance are more frequent in Fig. 2c than in Fig. 2b. Conversely in Fig. 2c, the black lithotype exhibits clearly connected components longer than in Fig. 2b. It is important to underline that the variogram does not give information concerning the connectivity of the lithotypes. This characteristic that could be important for some studies can only be a consequence of the model. Connectivity used as an input constraint for lithotype simulations is still a difficulty. Several authors have proposed solutions for sequential algorithms, for instance in truncated methods (Allard, 1994) adding constraints in the simulation process, or in the multipoint approach (Renard et al., 2011) keeping only the multipoint configurations that satisfy the connectivity constraint.

The truncated methodology can be extended to represent several lithotypes ordered sequentially, with a single GRF. The lithotype rule describes the ordering of the lithotypes along the Gaussian cumulative density function, and dictates the neighboring relationships between lithotypes. In Fig. 3a a GRF with an isotropic cubic model is considered. The lithotype rule (Fig. 3b) indicates that there should be no contact between black and light grey lithotypes, as demonstrated in the lithotype outcome (Fig. 3c).

The simple and cross-variograms of the three lithotypes are shown in Fig. 3d: they do not reproduce the cubic shape of the initial GRF. Note that because the GRF is isotropic each variogram is calculated without paying attention to the direction (omnidirectional). The sills of simple variograms depend on the lithotype proportions and the ranges depend on the Gaussian interval bounds: the intermediate lithotype (dark grey) has a shorter range than the others. The cross-variograms are all negative: as a matter of fact the lithotypes compose a partition of the space (at a given point only one indicator is equal to 1 the others being equal to 0), then the product of two terms non null is always negative (see Eq. 2). Another important remark for designing the lithotype rule is that the cross-variogram between the two lithotypes located at both extremities of the Gaussian axis (black and light grey) remains flat at short distances (equal to 0) as these lithotypes are not in contact.

As a conclusion, the spatial structures of the lithotypes described by their indicator variograms are linked to the spatial structure of the GRF (underlying variogram), but also involves the lithotype rule and the lithotype proportions. Note that some properties are transmitted from the spatial characteristics of the GRF to the lithotype indicator ones, such as anisotropy orientation. Moreover, unlike indicators for which the variogram models are specific, there is no limitation on the spatial characteristics of the GRF and then one can obtain very different multivariate indicator variogram models.

In order to choose the lithotype rule, other tools can be used to complete the analysis of the lithotype relationships: the edge effect functions (Eq. 3). These two-point statistics, function of their distance, give the probability that the first point x belongs to a lithotype A and the second point (x  +  h) belongs to another lithotype B, where h is the distance between the two points:

On these edge effect curves (Fig. 3e), the ordering of the lithotypes are obvious. Starting from the black lithotype (left), the probability to be in light grey is null at small distance while the probability to be in dark grey is equal to 1. The same feature is visible when starting from the light grey lithotype (right) and entering in black or dark grey. When starting from the dark grey lithotype (middle) the probabilities are intermediate. For long distances, in a stationary case, all the curves tend to the probabilities ratio P_B/(1–P_A).

The truncated Gaussian model with a single GRF produces layouts where lithotypes are ordered sequentially and presenting the same anisotropy. This method has been generalized by the use of two GRFs, leading to the pluri-Gaussian model (presented in supplementary figures S1, S2, S3). This improvement gives the opportunity to define different anisotropies for different lithotypes and to allow, or forbid, direct contacts between lithotypes. Then the resulting lithotype simulations go from nested to sequential lithotype arrangements.

In most studies, the available data is composed of wells scattered in a domain. The parameters which define the truncated Gaussian model are inferred from these data using the lithotype indicator functions. The lithotype rule is based on the geologist's knowledge ranging from the choice of the lithotypes (based on granulometry for instance) to depositional or transformation periods. The resulting construction is then validated with the simple and cross-variograms and the edge effect functions as explained in the previous section. Given the lithotype rule, the proportions of the lithotypes provide the intervals (thresholds) on the Gaussian values. Ultimately the simple and cross-variograms of the lithotype indicators define the model of the GRF(s).

The lithotype proportions are important ingredients of the method as the resulting simulations honor them. The proportions may vary in space. They are first computed along depth in order to check the vertical stationarity. The proportions of each lithotype are cumulated, from discretized information, along horizontal slices giving access to vertical cumulated proportions curves (VPC). These computations are performed with respect to a reference surface, i.e. within paleo-horizontal slices. Therefore a prior flattening may be necessary (Mallet, 2004; Perrin et al., 2012; Ramón et al., 2012). This flattening depends upon the analysis of the sequential stratigraphy (Ravenne, 2002) in consistency with the geological regional knowledge, in order to represent correctly the sequence of deposit which provides a typical signature on the proportions (Volpi et al., 1997).

It may also happen that the proportions show lateral variations. Then the lithotype proportions should be estimated on a 3D grid to account for a 3D non-stationarity of the lithotypes.

The estimation of the 3D lithotype proportions is a critical step. This estimation is performed either from lithotype indicators at wells or from vertical proportion curves. It can be constrained by auxiliary variables such as seismic attribute (Doligez et al., 2007) or using sequence stratigraphy concepts to define transition zones between wells (Labourdette et al., 2008). Finally this estimation is a challenging technique as the lithotype proportions are compositional variables (they must add up to 1).

With a single GRF, the resulting simulation outcome corresponds to a deposit where the lithotypes are sequentially organized, for instance when lithotypes are assigned to classes with increasing clay contents.

This approach has been used to reproduce a turbiditic environment (Fig. 4). Vertical sections of an analog of an oil reservoir field have been sampled (Felletti, 2002, 2003). The first step of the study is to split the domain into sedimentary units, which will be simulated independently from each other. The top and bottom surfaces of each unit are mapped in a coherent way and reference grids for the proportion computations are chosen according to their geological environment. Fig. 4b represents the global vertical proportion curves of the different lithotypes: sand and mud with different characteristics (massive, laminated, etc.). In order to honor the lithotype non-stationarity, these proportions are estimated in the whole field from local VPCs. Finally the different units are simulated conditionally to the vertical sections and are merged together in order to construct the final realization of the field (Fig. 4c).

The outcrop presented in Fig. 5a, located in Paradox Basin (Galli et al., 2006; Van Buchem et al., 2000), is the result of a submarine depositional sequence giving mainly horizontal structures and some algal mounds, which correspond to the vertical extensions visible on the geological interpretation (Fig. 5b). In the upper unit, the submarine structures follow the shape of the mounds.

To reproduce this complex deposit, where sets of lithotypes present different orientations and textures, two underlying GRFs are needed in the truncated pluri-Gaussian model. Moreover, the shape dependency is obtained by introducing a correlation between the two GRFs. The algal mound lithotypes (in green) and the surrounding marine lithotype (in yellow) are located in the lowest part of the lithotype rule and depend on the two GRFs (Fig. 5c). The vertical shape of the algal mounds and the smoothed contact between those lithotypes and the others are well reproduced on the conditional simulation (Fig. 5d).

When dealing with two qualitative variables, such as the sedimentology and the diagenesis index for example, several approaches can be considered (Doligez et al., 2010, 2011). One option consists in filling the unit with sedimentary lithotypes first using the most appropriate methodology, then simulating the diagenetic index with the corresponding property per sedimentary lithotype. This sequential approach assumes independency of the diagenesis from one lithotype to another one (Fig. 6a).

Another possibility is to use the pluri-Gaussian method with a lithotype rule where the first GRF reproduces the organization of the sedimentary lithotype and the second GRF provides the sequence of the diagenesis indices (Fig. 6b). This approach allows to link the diagenesis index to the sedimentary facies, and to reproduce the gradual transitions for both variables. However, it assumes that each variable (sedimentation and diagenesis) presents an ordered organization as it is based on a single GRF.

A generalization which breaks this limitation is the use of two truncated pluri-Gaussian models, one for each phenomenon: hence the name of bi-pluri-Gaussian model (Renard et al., 2008). Then the layout of each phenomenon may be as complex as needed, i.e. two GRFs for each variable possibly correlated. Moreover, this model can honor a heterotopic data set: as a matter of fact, the diagenetic index is not necessarily measured at each sample where the sedimentary facies is known (Fig. 6c).

The lithotype shapes can also be considered as resulting from a sedimentary process that induces some particular layouts. Some efforts have been carried out to reproduce these shapes using truncated Gaussian constructions, choosing appropriate underlying GRFs and/or lithotype rules. On the following simulations the aim is to reproduce particular shapes; they are not conditioned to any dataset.

As an example, one can obtain patterns similar to ripple marks on the beach (Fig. 7a) by using a GRF characterized by a periodic variogram (Fig. 7b), with the period being related to the wave dimension. The non-stationary density of ripples from top to bottom when water depth decreases is obtained by proportion changes (Fig. 7c).

Mixed scale textures can be reproduced using two underlying GRFs (Fig. 8a) and particular thresholding schemes. In these cases, a single facies can have variable sizes or orientations. In Fig. 8b, the black lithotype is ribbon shaped with large ribbons due to the GRF2 and thin at a smaller scale due to the GRF1. In Fig. 8c, the black lithotype is a mix with ribbon and grains depending on the thresholds.

Although the truncated pluri-Gaussian simulation method can generate complex layouts which cope with most of the sedimentary deposits, the outcomes always present a symmetrical pattern. Specific options have been designed to overcome this symmetry limitation and introduce an orientation in the facies organization.

This is the case when representing a uranium roll-front produced by the circulation of an oxidizing fluid in a fluvial deposit, where the mineralized zone is located at the forefront of the fluid penetration (Fig. 9a) (Langlais et al., 2008; Renard and Beucher, 2012).

A solution consists in using two GRFs. The first one distributes oxidized versus reduced facies, and the second one locates the mineralized areas at the interface between the two phases while following the direction of the oxidizing fluid. It suffices to consider the second GRF equal to the first one shifted by a vector oriented along the flow direction and whose size is equal to the average dimension of the mineralized areas (Fig. 9b). Unfortunately this solution does not guarantee that the mineralized and the oxidized areas are contiguous. Moreover it leads to an equal size of the mineralized area over the whole field which may not be realistic.

The other solution starts from the same first GRF, but the mineralized area is obtained as the shadow cast by the first GRF considered as a relief. The parameters of the light position and a possible truncation of the relief are deduced from the extension of the mineralized areas to be reproduced (Fig. 9c).

At a smaller scale (some centimeters), the facies organization may need to account for characteristics of high frequency transitions. For example, to reproduce ripples from marginal marine siliciclastic sediments (Fig. 10a) in sets of different properties (electric resistivity or rock hardness) a substitution model has been applied (Lantuéjoul, 2002). This model is similar to truncated Gaussian model: it corresponds to the coding of one GRF. However, in this model, the GRF has strictly stationary increments and the coding is a 1D simulation of a Markov chain process whose transition probability from one set to another is calibrated on the data (Fig. 10c). The length of the Markov chain is equal to the range of the simulated Gaussian values. On the resulting simulation, the anisotropy is determined by the GRF anisotropy and the lithotypes alternating depends on the transition probability between them (Fig. 10b).

By analogy, it is possible to construct a generalized truncated Gaussian, starting with a stationary underlying GRF, but using a lithotype rule simulated with a Markov chain. In that case the length of the Markov chain is chosen a priori. This construction allows reproducing other arrangements. For instance with a cubic model for the GRF and a cyclic Markov chain, the simulation result looks like a horizontal section in stromatolithes (Fig. 11).