Geostatistical Modelling of Cyclic and Rhythmic Facies Architectures

A Pluri-Gaussian method is developed for facies variables in three dimensions to model vertical cyclicity related to facies ordering and rhythmicity. Cyclicity is generally characterised by shallowing-upward or deepening-upward sequences and rhythmicity by the repetition of facies at constant intervals along sequences. Both of these aspects are commonly observed in shallow-marine carbonate successions, especially in the vertical direction. A grid-free spectral simulation approach is developed, with a separable covariance allowing a dampened hole-effect to capture rhythmicity in the vertical direction and a different covariance in the lateral plane along strata, as in space-time models. In addition, facies ordering is created by using a spatial shift between two latent Gaussian functions in the Pluri-Gaussian approach. Rapid conditioning to data is performed via Gibbs sampling and kriging using the screening properties of separable covariances. The resulting facies transiograms can show complex patterns of cyclicity and rhythmicity. Finally, a three dimensional case study of shallow-marine carbonate deposits at outcrop shows the applicability of the modelling method.


Introduction
Spatial distributions of facies in sedimentary rocks are commonly characterised by cyclicity and rhythmicity in the vertical direction and a variety of lateral patterns along stratigraphy. The resulting facies architectures control heterogeneity in hydrocarbon reservoirs and groundwater aquifers. It is, therefore, important to represent them in three-dimensional geostatistical earth models that are used as input to flow simulations and reserves quantification (Pyrcz and Deutsch 2014).
Cyclicity is defined as a characteristic facies ordering in vertical successions (Wilkinson et al. 1997;Burgess 2016). The characterisation of cyclicity needs to be addressed statistically (Wilkinson et al. 1997) in order to apprehend the variability of the resulting facies patterns and to reproduce them in earth models. Facies cycles show preferential transitions between successive facies, such that one facies tends to be observed on top of another facies. This is also called asymmetry (Carle and Fogg 1996), because the transitions between facies differ between the upward and downward directions. For example, shallow-marine carbonate rocks at outcrops (Strasser 1988;Goldhammer et al. 1990) and in subsurface reservoirs (Lindsay et al. 2006) are typically characterised by facies cycles that record upward shallowing (regression) and that consist of subtidal facies overlain by intertidal facies then by supratidal facies. The facies succession that records upward deepening is commonly incomplete or absent, due to non-deposition or erosion, such that supratidal facies are directly overlain by subtidal facies, which mark the base of a new cycle. Such sequences are illustrated in Fig. 1c, d.
Classical geostatistical methods such as sequential indicator simulation (Alabert 1989) and object-based methods (Deutsch and Tran 2002) do not in their traditional form reproduce this asymmetric facies ordering. Multipoint statistics (Strebelle 2002) (a) (b) (c) (d) Fig. 1 a-d Synthetic examples of facies sequences: a cyclic (two cycles), non rhythmic, b cyclic (two cycles) and red facies presents rhythmicity, c perfect cyclicity, non rhythmic d perfect cyclicity and rhythmicity should be able in theory to reproduce cyclicity, but in practice, it is not easy to obtain a three-dimensional training image showing cyclicity. The representation of asymmetry and facies ordering is, however, straightforward in one dimension with Markov chains (Carle and Fogg 1996;Parks et al. 2000;Li 2007;Purkis et al. 2012) or renewal processes (Matheron 1968), which are based on probabilities of transition between facies, but are difficult to generalise to two or three dimensions. Rhythmicity is another important aspect observed in vertical facies successions. It is characterised by the repetition of a facies at a constant interval along a sequence, a feature that has commonly been used to interpret periodic processes of deposition (e.g., via analysis of Fischer plots) (Read and Goldhammer 1988). Note that periodic processes (i.e., repetitive processes when looking at a time series) can either result in rhythmic sequences (when sedimentation rates are similar from cycle to cycle) or non-rhythmic sequences (if sedimentation rates change between cycles). Because the space domain is considered here, we use the term rhythmicity rather than periodicity. Rhythmic stacking of facies cycles has been observed in many shallow-marine carbonate successions (Goldhammer et al. 1993;Egenhoff et al. 1999;Lindsay et al. 2006). This aspect is also shown by the red facies in Fig. 1b, d, in which thicknesses between different beds of the same facies are constant.
Rhythmicity can be quantified by geostatistical tools, such as the variogram and the transiogram, which show oscillations or dampened oscillations called hole-effects (Pyrcz and Deutsch 2014). By looking at the probability density function (pdf) of facies thicknesses, Ma and Jones (2001) show that, as the coefficient of variation of this pdf decreases, the hole-effect becomes more pronounced. This observation is in agreement with the previous remark that rhythmicity is associated with low thickness variability of vertically stacked facies cycles. This also explains why Markov processes cannot create hole-effect transiogram models (Dubrule 2017), as the corresponding thickness pdf is exponential (coefficient of variation equal to one). On the other hand, some renewal processes may be able to create dampened hole-effect transiograms, because they offer the possibility to choose a thickness pdf with a lower coefficient of variation (Matheron 1968). However, the transiograms derived from renewal processes are not always known analytically, and are thus difficult to fit to the observed rhythmicity.
Truncated (Pluri-)Gaussian methods have been successfully used to create facies models (Armstrong et al. 2011), and they have been applied to shallow-marine carbonate reservoirs Amour et al. 2012;Le Blévec et al. 2017). The contacts between facies are defined by the truncation rule applied to a random Gaussian function, which provides control on facies juxtapositions. However, in its traditional form the method does not incorporate cyclicity and rhythmicity. Le Blévec et al. (2017) have extended the Pluri-Gaussian method to the modelling of facies asymmetry in vertical successions. They produce asymmetric transition probabilities between facies by introducing a shift in the correlation of two random Gaussian functions as suggested by Armstrong et al. (2011). This is similar to the approaches of Langlais et al. (2008) or Renard and Beucher (2012), but with more flexibility in the resulting facies transiograms. However, the use of this approach to model cyclicity and rhythmicity has not yet been investigated. Pluri-Gaussian Simulations (PGS) enable the use of hole-effect models (Beucher and Renard 2016) and may lead to hole-effect facies transiograms that could model rhythmicity.
Although cyclicity and rhythmicity are common features of vertical facies successions, they may have a variable expression laterally, depending on the formative depositional processes and controls. Laterally extensive facies in shallow-marine carbonate strata are generally attributed to external (allogenic) controls that operated over an entire carbonate platform or shelf, such as relative sea-level variations (Goldhammer et al. 1990). Facies of limited lateral extent may be attributed to the nucleation, vertical build-up and lateral shifting of tidal flat islands across a carbonate platform or shelf (Pratt and James 1986). This mechanism is internal to the dynamics of the carbonate platform depositional system (autogenic) and may generate both vertical and horizontal asymmetry in the stacking of facies if the tidal-flat-island deposits obey Walther's Law (Burgess et al. 2001;Le Blévec et al. 2016). Cyclic and rhythmic facies successions can also be overprinted by diagenetic facies after deposition; for example, hydrothermal dolomite bodies associated with faults and igneous intrusions are observed to cut across shallow-marine carbonate platform deposits characterised by rhythmic facies cycles (Jacquemyn et al. 2015).
In order to model cross-sections and volumes that exhibit cyclic and rhythmic vertical facies successions but different lateral facies patterns, it is necessary to use different vertical and lateral covariance models. This is possible via the use of separable anisotropic models (Chiles and Delfiner 2012), although such models have rarely been used for facies modelling (Matheron et al. 1988).
The aim of this paper is to extend the PGS approach of Le Blévec et al. (2017) to model facies cyclicity and rhythmicity in the vertical direction, and a range of appropriate lateral facies patterns using space-time (lateral-vertical) separable covariance models. After presenting the main definitions, the three key aspects of the modelling method and their impact on the transiograms are presented: cyclicity, rhythmicity, separability. A new method for simulating the resulting complex facies architectures is then presented, firstly for unconditional simulations and then for simulations conditioned to data. Finally, the method is applied to a case study from the Triassic Latemar carbonate platform (northern Italy), which has been widely interpreted to show cyclicity and rhythmicity.

Geostatistical Quantification with Transiograms
The random function representing a facies is the indicator function I (x). If the facies i is present at a spatial location x, I i (x) = 1 and if not, I i (x) = 0. In the stationary case, the probability of having a facies at a location x is equal to the first statistical moment or proportion The presence of a facies also depends on the surrounding facies, quantified by the covariance function c i j (h) for two facies i and j. With the stationary assumption, it is assumed that the covariance depends only on the vector h separating two locations (Chiles and Delfiner 2012). This paper uses the non-centered indicator covariance from which can be derived the transiogram t i j (h), which is the probability to transition from facies i to facies j along a certain vector h This transiogram can be asymmetrical, which means that it is different in opposite directions and can thus quantify asymmetrical facies successions (Carle and Fogg 1996;Le Blévec et al. 2017). If the transiogram is symmetric, it means there is no polarity in the facies succession and the transitions between facies are equivalent in opposite directions. The transiograms between different facies can be gathered in a transiogram matrix that the simulation method aims to reproduce. For instance, a transition matrix between three facies (1, 2, 3) is The terms on the diagonal are the auto-transiograms and those off-diagonal are the cross-transiograms. With stationary proportions, the tangent at the origin of the autotransiograms T ii (which is negative) is related to the mean length of the facies i along this direction (Carle and Fogg 1996). The tangent at the origin (at right) of the crosstransiogram is called the transition rate Transition rates are of interest when studying asymmetry because they are related to juxtaposition between different facies along the studied direction h. According to Carle and Fogg (1996) (under the stationary assumption) if a facies j tends to follow a facies i in the direction h, rather than preceding it, then Transition rates are also related to embedded transition probabilities r i j (given the presence of a facies i, the probability that it precedes j in the direction h) which can also be gathered in a transition matrix R = ⎡ ⎣ 0 r 12 r 13 r 21 0 r 23 r 31 r 32 0 ⎤ ⎦ .
The diagonal equals zero because embedded Markov chains only record the transitions between different facies. The following important properties of the embedded matrix are used in this paper (Carle and Fogg 1997) In the next section, it is shown how a geostatistical simulation method, the truncated Gaussian simulation, relates to these quantities.

Truncated Gaussian Simulations
Truncated Gaussian (TGS) and PGS (Armstrong et al. 2011) consist of simulating one or several standardised Gaussian random functions that are then truncated into a facies field. Any pair (Z (x), Z (x + h)) of a Gaussian random function Z (x) is a bi-Gaussian random vector, with Z (x) and Z (x + h) correlated to each other according to the non-centered covariance function The truncation rule determines which facies is present at location x from the value of the random variables Z (x). For instance, the truncation rule with the case of only two facies 1 and 2 and one Gaussian function (TGS) controls the indicator functions where q is the threshold of the truncation rule. It is possible to relate every moment of the facies field mathematically to those of the Gaussian function. According to Eq. (12a) the proportion of facies one (first order moment) is with g(x) the standardised Gaussian pdf. The cross-transiogram (second order moment) between facies one and two is where g ρ(h) is the standardised bi-Gaussian probability density with correlation matrix defined by the covariance ρ(h). For q = 0 the two facies have the same proportion 1/2 and the analytical solution of this bi-Gaussian integral (Lantuéjoul 2002;Le Blévec et al. 2017) gives the following auto-transiogram for the two facies and the cross-transiograms For different thresholds, the transiograms can be derived by numerical integration (Genz 1992) or expansions into Hermite polynomials (Chiles and Delfiner 2012). More than one Gaussian function can also be used. This is the PGS approach, which provides more flexibility thanks to a larger number of possible truncation rules. In this paper, two Gaussian functions are used, with two thresholds defining three facies (1, 2, 3) ( Fig. 2) As for TGS, the corresponding indicator statistical moments can be derived by numerical integration (Genz 1992). Therefore, the facies transiograms can be derived from the parameters of the truncated Gaussian model which are summarised in Table 1. A more detailed discussion on the link between the Pluri-Gaussian parameters and the transiograms is given in Le Blévec et al. (2017). This link is used in this paper to match different transiograms and a key objective is to use a covariance ρ(h) such that the transiograms t i j (h) show rhythmicities and asymmetries.

Understanding Cyclicity and Rhythmicity with Transiograms
Wilkinson et al. (1997) define cyclicity as an apparent ordering between facies. Therefore, this definition directly relates to the transition rates (Eq. 6) or the embedded transitions (Eq. 8) that define the juxtapositions between facies. For instance, the perfect cyclic sequences of Fig. 1c, d, have the upward embedded transition matrix  One can quantify cyclicity with the probability P c to observe a given cycle above a given facies, which can be written as the probability of observing the sequence 2-3-1 above an occurrence of facies one (under the Markov assumption) We see that the notion of cyclicity is related to asymmetry by using the closing relations of embedded transition matrices (Eqs. 10a, 10b) The probability P c increases as r 12 increases and r 21 decreases. Therefore, for a sequence with three facies, asymmetry between two facies results in cyclicity. The cyclicity studied here ( Fig. 1c, d) is such that each facies appears only once per cycle. This is incompatible with symmetric cycles, which are not considered here. Rhythmicity is defined by the repetition of a facies at constant intervals. This is usually observed on experimental transiograms that show dampened hole-effects (Journel and Froidevaux 1982;Johnson and Dreiss 1989;Ma et al. 2009). Rhythmicity cannot be quantified by embedded transition rates as they are independent of facies thicknesses. Rhythmicity can first be understood when studying two facies. If those facies have constant thicknesses, the auto-transiogram varies with a constant wavelength that is equal to the sum of the two facies thicknesses . This is similar with more facies as we can still regard this as the succession of two facies, the one of interest and its complement. This is interesting to consider in combination with cyclicity because the resulting sequences show a constant cycle thickness (Fig. 1d). Thicknesses along sequences are not usually constant but can show low variability which results in non-perfect rhythmicity and dampened hole-effects. The method developed here models cyclic and rhythmic sequences, quantified by transiograms.

The Cyclical Pluri-Gaussian Approach
In this section, the classical PGS is extended to render cyclicity (or asymmetry) (Sect. 3.1), rhythmicity (Sect. 3.2) and separable anisotropy (Sect. 3.3). The first two sections (Sects. 3.1, 3.2) present results in one dimension and the last section (Sect. 3.3) extends them to three dimensions.

Modelling Asymmetrical Facies Juxtapositions in Vertical Successions
A method to simulate asymmetrical facies successions is summarised here (see Le Blévec et al. (2017) for a detailed treatment). For simplicity two standardised Gaussian random functions with a Cartesian truncation rule are used (Fig. 2) as in Eq. (17). The correlation between the Gaussian functions Z 1 (x) and Z 2 (x) is based on the linear model of co-regionalisation (Wackernagel 2003) with a shift α between the Gaussian random functions as proposed by Armstrong et al. (2011) with Y 1 (x) and Y 2 (x) uncorrelated standardised Gaussian random functions, ρ 1 (h) and ρ 2 (h), respectively, their covariances and ρ the correlation coefficient between The proportions of the different facies define the thresholds of the truncation rule as in Eq. (13), which are also impacted by the correlation coefficient ρ. The resulting transiograms between facies can be computed by Gaussian integral on the facies domain defined by the truncation rule as in Eq. (14). Le Blévec et al. (2017) demonstrate analytically and numerically that these transiograms are asymmetric as in Eqs. (4) and (7). A sensitivity study is also carried out on the parameters α and ρ showing that the asymmetry can be controlled by varying these two parameters (Le Blévec et al. 2017). Here, the covariance models ρ 1 (h) and ρ 2 (h) used for the Gaussian functions Y 1 (x) and Y 2 (x) are Gaussian with parameters r 1 and r 2 the scale factors of the models. This gives in total four parameters for the method: α, ρ, r 1 and r 2 . An example of a vertical sequence generated with this method and the corresponding facies transiograms is shown in Fig. 3. The simulation method of the latent Gaussian functions is given in Sect. 4 and the transiograms are computed numerically (Genz 1992) based on Eq. (14). For clarity, only four transiograms instead of nine are shown for the three facies ( Fig. 3a-d), because the transiograms of the third facies can be directly derived from the transiograms of the other facies. However, in practice, one can work with all the transiograms.
It is clear in the realisation shown in Fig. 3 that the facies are statistically organised in shallowing-upward cycles as highlighted by the low tangent at the origin of the cross transiogram from intertidal (orange) facies to subtidal (red) facies upwards (Fig. 3c). Indeed, this transition is absent along the vertical section while the opposite transition (from subtidal facies to intertidal facies) occurs six times. However, the succession is not perfectly cyclic (Fig. 3e) because of a low, but non-null probability of the subtidal (red) facies to transition upwards directly to the supratidal (white) facies. A high variation of facies thicknesses is also noted, resulting in a non-rhythmic sequence.
One realisation e with corresponding transiogram matrix model a-d and Gaussian functions Z 1 (x) and

Modelling Rhythmicity in Vertical Successions
Using covariance functions with hole-effects for the latent Gaussian function is expected to produce hole-effect facies transiograms. This is verified in this section in which a new indicator transiogram model representing facies rhythmicity is obtained. In one dimension, the cosine function is a valid covariance model (Chiles and Delfiner 2012) and produces by truncation vertical sequences with a constant cycle thickness. As this is rarely observed, a model in which the oscillations attenuate with distance, called dampened hole-effect (Pyrcz and Deutsch 2014) gives more flexibility in reproducing cycle thicknesses. Ma and Jones (2001) define the Gaussian cosine covariance as the product of two valid one-dimensional covariance functions This covariance gives dampened oscillations controlled by the two parameters r and b. If the scale factor r tends to infinity, the model is the cosine function. According to Eq. (15) the resulting hole-effect transiogram model is As seen in Fig. 4, this transiogram model has the same wavelength as the latent covariance model (Eq. 24) and the attenuation of the oscillations is also similar. Therefore, this transiogram model seems a good candidate to model rhythmicity because of its flexibility.
It is important to remember that Eq. (25) is valid only if the threshold is zero, that is the two facies have equal proportions. When this is not the case, the model t 11 (h) can be numerically computed (Genz 1992). Figure 5 gives the simulation of a vertical sequence with three facies (Fig. 5e) and corresponding hole-effect transiograms (Fig. 5a-d). The covariances of the two Gaussian random functions have the form of Eq. (24) with respective parameters r 1 , b 1 and r 2 , b 2 . Figure 5 clearly shows the effect of rhythmicity on the transiograms and the corresponding realisation. The facies cycles are repeated in the vertical succession (Fig. 5e) with a rhythmicity controlled by the latent Gaussian functions (Fig. 5f). Asymmetry in facies stacking is also added to create a cyclical vertical succession. After developing covariance models in one dimension, it is necessary to expand these models into two and three dimensions while incorporating anisotropy.

Modelling Facies Distributions in Two and Three Dimensions with Separable Anisotropy
In sedimentary deposits, a strong anisotropy is always observed between the vertical direction and the one parallel to stratigraphy. A simple way to represent such anisotropy is to use a separable covariance (Chiles and Delfiner 2012). A separable covariance model can be built from the product of a covariance in the vertical direction ρ v (h z ) and a lateral covariance along stratigraphy ρ l (h x , h y ) Along stratigraphy, a standard geometrical anisotropy (which is also separable in this case of a Gaussian covariance) is here used (Chiles and Delfiner 2012). In an earth modelling software, the lateral direction is usually controlled by the stratigraphic grid. According to Eq. (15), the corresponding transiogram model of two facies, after thresholding a single Gaussian at cutoff 0 is, therefore, Equation (28) shows that even though the covariance of the latent Gaussian function is separable, the resulting transiogram is not separable. Figure 6 displays maps (h y constant) of the latent Gaussian function covariance and the resulting transiogram, which both show a dampened hole-effect in the vertical direction and the absence of a hole-effect along stratigraphy. In intermediate directions, the dampened hole-effect is present, but attenuated. Once again, the behavior of the auto-transiogram in two or three dimensions is very similar to the covariance of the latent variable. This suggests that there is about as much flexibility in the transiogram model as in the covariance model.  (Figs. 7, 8). The transiograms have been computed by numerical integration (Genz 1992) (Eq. 14). The parameters for the two Gaussian functions are respectively r 1 x , r 1 y , r 1 z and r 2 x , r 2 y , r 2 z , and vertical hole-effect parameters b 1 , b 2 . Figure 7 shows the combination of asymmetry, rhythmicity and anisotropy in facies distributions in both the realisations and the transiograms. The facies are ordered according to the cycle: subtidal, intertidal and supratidal facies upwards (Fig. 7e). The vertical thickness of these cycles is almost constant, confirming the rhythmicity. The different realisations can be interpreted as facies bodies that pinch out laterally or split into different cycles and inter-finger. Changes in the different parameters can create more complex transiogram models. By using a hole-effect covariance on the second Gaussian function, but not on the first, only the orange facies auto-transiogram shows the combined effect of the Gaussian covariance and the hole-effect covariance (Fig. 8d). This results in more complex geometries for the orange and white facies that show two combined behaviors. They show thin beds corresponding to the high frequency rhythmicity and thicker beds (Fig. 8e) corresponding to a lower frequency. The red facies has a different transiogram from the other facies (Fig. 8a) and appears to truncate both of them. Bodies of the red facies can be interpreted either to be erosional or to be diagenetic in origin, having formed after deposition. In order to have a better understanding of the geometries generated by this transiogram model (Fig. 8), a three-dimensional simulation with the same transiograms is also presented in Fig. 9. The orange and white facies clearly show vertical rhythmicity with low lateral variability in thickness, consistent with a depositional architecture, while the red (diagenetic) facies truncates them.
Until now Gaussian covariances have been used. However, a more general model could be used in order to control the behavior of the transiogram at the origin. For instance, using the stable covariance model (Chiles and Delfiner 2012) with a geometrical anisotropy along stratigraphy and separable anisotropy in the vertical direction would give transiograms of the form (for one latent Gaussian function with a threshold at zero) The smoothness of facies boundaries decreases with the coefficients β and γ (0 < (β, γ ) ≤ 2).

Conditional Simulation of the Cyclical Pluri-Gaussian Model
A simple method is here presented to simulate the latent Gaussian functions with covariance presented in Sect. 3 and to condition them to the facies observed along the wells.

Unconditional Simulation
The continuous spectral approach developed by Shinozuka (1971) is convenient to build simulations based on separable covariances. The Fourier transform of the covariance model is normalised and used as a probability density function (pdf) to sample frequency vectors ν k generating the three dimensional Gaussian function with Φ k random phases sampled from a uniform pdf between 0 and 2π and <> the scalar product. The resulting function is Gaussian when N tends to infinity (central limit theorem). Based on the knowledge of x, the individual value of Z (x) at each location x can be simulated, which enables the algorithm to be coded in parallel and to be grid-free. The spectral approach is well suited for a separable covariance as the multidimensional Fourier transform of a separable covariance is the product of the Fourier transforms. The three-dimensional Fourier transform of a separable covariance model which means the frequencies ν z in the vertical direction are sampled independently from ν x and ν y according to the Fourier transform of their respective covariance models.
In order to simulate the covariance model of Eq. (27), the Fourier transforms of the covariance in the vertical direction and the lateral plane must be known (Eq. 31). For the lateral plane, as a Gaussian covariance with geometrical anisotropy is also separable, the two directions h x and h y are also sampled independently The Fourier transform of a Gaussian covariance is a Gaussian pdf ( Algorithm 1 describes how to simulate a facies field with such a covariance model. Figure 10 shows experimental variograms (1−ρ(h)) of 10 realisations of the simulation of a separable Gaussian cosine covariance in two dimensions (Eq. 27). As expected the covariances of the realisations are centered around the covariance model and represent the hole-effect only in the vertical direction. Some fluctuations of the experimental variograms on the realisations are observed (Fig. 10), which are explained in detail by Lantuéjoul (1994) or Emery and Lantuéjoul (2006). The spectral approach can also be used to simulate a stable separable covariance model (Eq. 29). The Fourier transform of a stable covariance is a stable spectral pdf with stability and scale parameters α and r , skewness and location 0 (Chiles and Delfiner 2012). Therefore, to simulate this covariance, one can use Algorithm 1 by replacing the Gaussian distributions with the stable distributions. These simulations are unconditional, and they must now be conditioned to facies data observed at specific locations in the simulated volume.

Conditioning the Gaussian Simulation to Facies Data
Some methods already exist for conditioning PGS to hard data (Emery and Lantuéjoul 2006;Chiles and Delfiner 2012). This section summarises the main steps and proposes improvements in the case of separable covariances. The truncation rule defines intervals of the Gaussian function corresponding to each facies. At data locations x i , the observed facies constrain the Gaussian random functions with the following inequalities according to the truncation diagram (Fig. 2). The conditional simulation is usually performed in three steps (Chiles and Delfiner 2012), (i) unconditional simulation, (ii) Gibbs sampling at data locations and (iii) conditioning by kriging of the mismatch at data locations. Here, the conditioning is broken down in more steps for optimisation purposes: 1. Unconditional uncorrelated simulations Y u 1 (x) and Y u 2 (x) are performed over the domain as in Sect. 4.1. 2. Joint local conditional simulation of Z 1 (x i ) and Z 2 (x i ) is carried out only at data locations x i using Gibbs sampling such that Eq. (34) is respected (Freulon and de Fouquet 1993). 3. Back-transform of correlated Z 1 (x i ) and Z 2 (x i ) into uncorrelated Y 1 (x i ) and Y 2 (x i ) at data location according to Eq. (22). 4. Two separate simple kriging R 1 (x) and R 2 (x) based on the mismatch at data locations are performed (Chiles and Delfiner 2012) with λ 1 i (x) and λ 2 i (x) the kriging weights (Chiles and Delfiner 2012). The conditional uncorrelated Gaussian random functions Y 1 and Y 2 at every location x are finally obtained by 5. Transform conditioned Y 1 and Y 2 into conditioned and correlated Z 1 and Z 2 according to Eq. (22).
The Gibbs sampling at data locations is not described here. Any algorithm can be chosen (Emery et al. 2014) and applied with the desired covariance model to simulate the values of Z (x i ) at the data locations. Because Z u 1 and Z u 2 are co-simulated, the covariance matrix used for the Gibbs sampling is made of the cross-covariances between Z u 1 and Z u 2 , given in more details in Le Blévec et al. (2017). As this covariance matrix can be singular, it is advised to perform a Gibbs sampling that does not involve its inversion (Dubrule 1983;Emery et al. 2014).
The kriging of the difference (Eqs. 36a, 36b) can be optimised when using a separable covariance such as the model of Eq. (26) and if the wells are strictly vertical. This is important because kriging over the whole domain is computationally expensive and can make simulation prohibitively slow. The simple kriging is optimised thanks to the well known screening properties of separable covariances (Matheron et al. 1988;Chiles and Delfiner 2012). Performing kriging in the uncorrelated space of the Y i allows firstly a reduction of the size of the covariance matrix to invert (and the stability of the inversion) for finding the weights λ 1 i and λ 2 i . Secondly, it ensures that the covariance is separable (and thus its screening property), which is not the case for the covariance of Z 2 .
With a separable covariance model such as that of Eq. (26), the weights λ i (x) associated with data located on different lateral planes from that of the estimated location x are equal to zero (Chiles and Delfiner 2012). This means that the estimation at a given location only depends on the data at the same horizontal level, and the number of kriging weights, therefore, equals the number of wells intersecting this horizontal level, as shown in Fig. 11. Assuming all wells are vertical and have the same length, the number of weights for every kriged point is, therefore, simply the number of wells. If this is not the case, for instance because some vertical wells do not penetrate a particular level, it is convenient to extend them artificially by an unconditional simulation with Gibbs sampling at step (ii) (Fig. 11), so that the geometrical configuration of the data points remains the same at all levels. Therefore, the weights are the same for every horizontal plane and the dual form of two-dimensional kriging may be used, in which the data covariance matrix is inverted only once (Dubrule 1983;Chiles and Delfiner 2012). This enables rapid and efficient kriging. An example of conditional simulation with this method is given in the following section.

Case Study: The Latemar Carbonate Platform, Northern Italy
The Triassic Latemar carbonate platform is superbly exposed in northern Italy and has been chosen as a case study as it provides well-documented examples of facies cyclicity with rhythmicity and asymmetry in shallow-marine peritidal strata (Goldhammer et al. 1990). Rhythmicity has been quantified along a vertical sequence through the entire platform succession using spectral analysis, in order to understand the potential expression of orbital forcing (Milankovitch cycles) on platform growth (Hinnov and Goldhammer 1991;Preto et al. 2001). However, few studies describe sequences in different parts of the platform with the aim of performing facies modelling. Here, the dataset of Peterhänsel and Egenhoff (2008) (upper cyclic Facies stratigraphic interval) is used to evaluate and model asymmetrical (upward-shallowing) facies cycles and facies rhythmicity near the topographically high platform margin (Cimon Latemar) and in deeper lagoonal deposits of the platform interior (Cimon Forcellone). Peterhänsel and Egenhoff (2008) describe five microfacies based on thin sections: peloidal Fig. 11 The screening effect for simple kriging applied on two wells W 1 and W 2 to estimate three locations on three different horizontal planes with one extended well (red) packstone to wackestone (facies one), algae peloidal packstone (facies two), fenestral packstone to wackestone (facies three), packstone to grainstone (facies four), diagenetically overprinted grainstone to packstone (facies five). To illustrate our geostatistical method, these microfacies are re-grouped into three main environments of deposition: subtidal (facies one), intertidal (facies two and three), supratidal storm deposits (facies four and five). The next paragraph provides a quantitative analysis of the variations of these three facies.

Qualitative and Quantitative Study of the Case-Study Dataset
As illustrated by Fig. 12, the eight vertical logs show a high number of facies transitions. The asymmetry is clear as the subtidal facies (red) tends to be on top of the supratidal (white) facies. However, complete upward-shallowing facies cycles, containing subtidal, intertidal and supratidal deposits, occur only 24 times, while there are 56 incomplete cycles in the eight logs, which means that the sequences have some cyclic features. The subtidal facies appears to show regular spacing between beds within some wells (Fig. 12,in logs N8 and N16), which would suggest a rhythmicity of this facies. This is not the case for the intertidal and supratidal facies, which show very different spacings between the beds (Fig. 12).  (Peterhänsel and Egenhoff 2008) subtidal (red), intertidal (yellow) and supratidal (white) All this information can be verified in the experimental transiograms computed on the logs (Fig. 13, grey points). A dampened hole-effect is observable on the autotransiogram of the subtidal deposits (Fig. 13a). The tangent at the origin of the crosstransiogram of intertidal deposits overlain by subtidal deposits (Fig. 13) is low, showing that this transition is rare. More precisely, T 12 = 2.06 and T 21 = 0.19, which means according to Eq. (7) that the intertidal facies is four times more likely to overlie the subtidal facies than to underlie it.

Inference of the Co-regionalisation Model
As seen in Fig. 13a, the subtidal deposit auto-transiogram shows a dampened holeeffect. In order to model this hole-effect properly, it is easier to derive it from only one covariance rather than a combination of two covariances. Therefore, the truncation rule of Fig. 2 where the subtidal facies is defined by only one Gaussian function, is chosen. The use of more complex truncation rules is discussed in Sect. 6.1.
The parameters of the Pluri-Gaussian model are determined using the procedure described in Le Blévec et al. (2017). Only the auto-and cross-transiograms of two of the three facies are necessary to determine the parameters of the model, as the transiograms for the third facies are derived from those of the two other facies. A trial and error procedure is here chosen for determining the parameters, as this gives the possibility to incorporate conceptual knowledge. For instance, as it is known that rhythmicity and asymmetry occur, it is important to incorporate these features during the model construction, which would be difficult with an automatic fitting procedure.
The covariance of the first Gaussian function ρ 1 (h) is used to describe the transiogram of the subtidal facies. The Gaussian cosine model (Eq. 24) gives an appropriate  Table 2 , r 1 y , r 1 z ) (500 m, 500 m, 0.8 m) Scale factors of the first Gaussian covariance (r 2 x , r 2 y , r 2 z ) (500 m, 500 m, 0.4 m) Scale factors of the second Gaussian covariance Vertical frequencies of the two covariances fit although the oscillations observed on the logs seem less pronounced (Fig. 13a). This fit could be improved by using a more complex covariance model, as discussed in Sect. 6.2. Then, the correlation coefficient ρ and the shift α can be chosen to fit the transition rates between the subtidal and intertidal facies. Here the asymmetry (Eq. 7) is clear, because the probability of intertidal facies overlying subtidal facies is four times higher that of the opposite. Therefore, the shift α is equal to the size of one vertical cell and the correlation coefficient ρ is rather high (− 0.6), which allows the model to match these transition rates. The shift is strictly vertical because no information on a possible lateral asymmetry is available (this topic is further discussed in Sect. 6.3). Finally, the covariance of the second Gaussian function ρ 2 (h) determines the transiogram of intertidal facies. A rhythmic covariance is not necessary here, as the transiogram does not show a clear hole-effect. There is not enough information to determine the lateral scale factors of the covariance directly from the data. Therefore, they are derived from visual comparison of the realisations with the outcrop facies panel interpreted between the vertical logs by Peterhänsel and Egenhoff (2008) and are chosen to be isotropic (in the lateral plane). Finally, the thresholds q 1 and q 2 of the Gaussian functions are determined according to the proportions of the facies and the correlation coefficient ρ between the two Gaussian random functions (Armstrong et al. 2011). The parameters are summarised in Table 2. Vertically, a Gaussian cosine covariance is chosen, but laterally an exponential covariance (Chiles and Delfiner 2012) is preferred. A Gaussian or cubic covariance in the lateral direction may lead to a singular kriging system for the conditioning step as the lateral scale factor is large (Sect. 4.2). A stable covariance with a power close to two could also be chosen to obtain a smoother model than that obtained with the exponential.

Simulation Results
The cell sizes of the grid are 0.1 m vertically and 10 m laterally in both north and east directions, which gives approximately 800,000 cells. The simulation takes approximately one minute to run with a standard Intel processor i7, but it could be much faster by performing the unconditional simulation with the spectral method in parallel, which is the longest computational step. The conditioning by kriging every surface independently, as described in Sect. 4.2, is almost instantaneous. One realisation of the field is shown in Fig. 14.  Fig. 14 One realisation of the Latemar facies field from Peterhänsel and Egenhoff (2008) data (Fig. 12), and two cross-sections obtained with the parameters summarised in Table 2 The visual aspect of the simulation is granular due to the exponential covariance model used in the lateral plane. Visually, the subtidal (red) facies tends to lie on top of the supratidal (white) facies and rhythmicity is confirmed by the regular thickness between two subtidal (red) facies bodies. The transiograms of the wells and the transiograms computed on one realisation are compared in Fig. 13. The sills derived from simulation (black diamonds in Fig. 13) are accurate, which means that the facies proportions match those in the wells. The tangent at the origin of the auto-transiograms is appropriately fitted, which means the average thicknesses of different facies bodies honours the well data. The tangent at the origin of the cross-transiograms matches the transition rates of the wells (T 12 = 2.1 and T 21 = 0.19), which means the asymmetry and cyclicity are respected. The hole-effect in the realisation transiogram is less pronounced than the one used for the theoretical model, but is closer to the one observed at the wells (grey circles in Fig. 13). This might be due to the conditioning of the realisation to the wells, where the experimental transiograms show less pronounced oscillations (black lines in Fig. 13).

Discussion
In order to illustrate the method, standard parameters have been chosen so far for both the synthetic examples and the case study. However, some parameters such as the truncation rule, the covariance model, or the shift can be changed to adjust to different geological environments.

Choosing the Truncation Rule
The same truncation rule has been applied through the paper (Fig. 2) with the subtidal facies defined by the first Gaussian function and the two others by both Gaussian functions. This has implications for the geometries of individual facies bodies and for facies relationships. The facies defined by the first Gaussian function erodes the two other facies (e.g., Fig. 8e), which means that bodies of these facies can have very different geometries from those of the other facies. This behavior can be reduced by increasing the correlation ρ between the two Gaussian functions Z 1 (x) and Z 2 (x).
The truncation rule affects not only the facies transiograms, but also higher order statistical moments, which can have an impact on connectivity (Beucher and Renard 2016). However, these moments are not known analytically and so are difficult to use for defining the truncation rule. It is recommended that the earth modeller tries different truncation rules and inspects the visual aspect of simulations so that it matches with his or her conceptual knowledge.

Adapting the Truncation Rule for More Facies
In this paper, only three facies have been used to illustrate the method. Two methods for generalising the truncation rule to more facies are found in the literature: either the truncation is made more complex (Galli et al. 2006) or the number of Gaussian functions is increased (Maleki et al. 2016). The first method is probably too limited to represent a large number of transition rates between facies, while the second should be able to model all transition rates (but the number of parameters would be very high). The choice between the two methods should depend on the case study and further work on this topic is required. It should also be noted that some methods that create automatic truncation rules have been developed Astrakova et al. 2015). It would be interesting to generalise these methods by incorporating the shift between the Gaussian functions in order to match asymmetric transition probabilities.

Elaborating More Complex Hole-Effect Models
The vertical hole-effect model used in this paper is made of two parameters r z and b (Eq. 24), which provides some flexibility to match observed rhythmicity. However, the case study shows that the observed transiograms can be even more complex (Fig. 13) and two parameters might not be sufficient to represent them. The covariance model could be modified to incorporate more than one structure (Chiles and Delfiner 2012). For instance a Gaussian covariance or a cosine covariance can be added to the Gaussian cosine model (Eq. 24).

Walther's Law
The method developed in this paper models cyclicity only in the vertical direction, which is consistent with observations of most outcrop and subsurface data. However, Walther's Law suggests that the transitions between facies should be equivalent laterally and vertically if no erosion is present (Middleton 1973). This means that the facies ordering could be similar and the transition rates proportional as in Markov chains methods (Doveton 1994;Purkis et al. 2012). Thus, asymmetry could also be observed laterally. Since the shift is actually a three-dimensional vector, it is possible to model such patterns with the presented method by defining a non-vertical vectorial shift between the Gaussian functions (Eq. 22), such that the asymmetry is also lateral (Le . This choice will depend on the depositional environments, processes, controls, and scale to be modelled. These aspects are typically interpreted with reference to an underlying conceptual model, such as those for allocyclically and autocyclically generated facies cycles in peritidal carbonate strata (Pratt and James 1986;Goldhammer et al. 1990). The facies architectures to be modelled are also scale dependent. Environments of deposition generally have large lateral extents (1-10 km), such that few lateral transitions between them are observed at reservoir (1-10 km) and inter-well (< 1 km) scales (Sena and John 2013), which limits the expression of lateral ordering of depositional environments. At smaller scales, the lateral transitions between lithofacies within depositional environments (or facies associations) may be different from the vertical transitions (Hönig and John 2015) because of erosion or lateral changes in palaeotopography. The resulting lithofacies distributions may be highly variable, potentially reflecting a facies migration that is well-ordered and obeys Walther's Law as one end member (Obermaier et al. 2015) or more complex and less ordered facies mosaics as the opposite end member (Wilkinson et al. 1997). The choice of appropriate conceptual model at the scale of depositional environments (facies association) or lithofacies must be made by the earth modeler in collaboration with the geologist, and then used to govern the selection of parameters of the model.

Conclusion
While cyclicity and rhythmicity are commonly observed in facies architectures, few existing geostatistical algorithms can model both patterns in an efficient manner. By addressing this issue, the method developed here is useful for modelling carbonate or shallow-marine clastic reservoirs that contain such cyclical facies successions. Broadly speaking, cyclicity and rhythmicity are quantified by facies transiograms that are computed from data (e.g., vertical facies successions) and fitted with an advanced truncated Pluri-Gaussian model for performing three dimensional simulations.
The model used for the latent Gaussian functions is the linear model of coregionalisation with a spatial shift, which creates the asymmetric cycles. The covariance of the Gaussian functions presents a dampened hole-effect, which captures the rhythmicity. As this hole-effect is generally observed only in the vertical direction, a separable covariance model, which is the product of a lateral and a vertical covariance is used so that no rhythmicity is modelled laterally along the stratigraphy. The space-time separable covariance is simulated readily by the continuous spectral method. The numerical properties of separable covariance allows rapid and efficient conditioning to data via kriging of every horizontal surface independently. The procedure has been applied successfully to model a carbonate platform environment that shows cyclicity and rhythmicity in facies architecture.