Theory of deterministical and stochastical indicator mapping methods and their applications in reservoir characterization, case study of the Upper Miocene reservoir in the Sava Depression

Indicator based geostatistical methods (Indicator Kriging and Sequential Indicator Simulations) are the mostly used methods for facies mapping or modelling. Although it is assumed that in facies variables should be discrete, it is possible to apply these methods on continuous variables as well. If a variable is continuous, a cumulative probability distribution function has to be created. Methodology includes a series of cut-o ﬀ s. On the cumulative probability distribution function, probabilities for all cut-o ﬀ s can be de ﬁ ned. Based on cut-o ﬀ s, all the data can be divided into two groups: (1) the data which has values lower than the cut-o ﬀ and (2) the data which has values higher or equal than the cut-o ﬀ . In this way, the indicator variable takes a value of 1 at all locations where the value is higher or equal than the cut-o ﬀ and 0 oth-erwise. The larger number of cut-o ﬀ s, the more precise the cumulative distribution function is. Indicator Kriging maps show the probability of certain lithofacies appearing in some location. On the other hand, stochastical realizations provide a di ﬀ erent number of solutions for the same input data set. Those solutions can be very similar, but never identical. It is important to emphasize that all obtained solutions are equally probable. Results of Sequential Indicator Simulations are also probability maps. There are several advantages for Indicator based methods. They do not ask for a normal distributed input dataset (e.g., can be implemented for bimodal distribution), and show connectivity of the largest or smallest values. Indicator Kriging and Sequential Indicator Simulations are applied for porosity mapping of the Upper Miocene sandstone reservoir in the Sava Depression. During the Upper Miocene, sands were deposited through turbiditic currents in the deepest part of the sedimentation area, forming turbiditic channels. Such channels can be recognized on the probability map for the cut-o ﬀ 19 % in Indicator Kriging as well as in Sequential Indicator Simulation maps.


Introduction
Geostatistics, as a part of geomathematics, began its development in the second half of the last century. Although it represents a relatively new approach, it is spread in almost every part of geosciences and today, mapping is almost impossible without it.
Geostatistical methods can be divided into deterministical and stochastical. In deterministical methods, all the conditions which can infl uence estimation are completely known. Deterministically, systems must not have randomness of any kind in variables' description and the results can be unambiguously described by the completely known fi nite conditions. It is clear that there is only one geological underground, but since the description of the underground is based on well data (point data) it is not possible to be absolutely sure if the solution obtained with geostatistical methods is absolutely correct. This is why all geostatistical methods contain some uncertainty, they estimate underground conditions and they can provide only the most probable solution or the solution which is the closest to the truth. Some methods give only one solution, i.e. for the same input data and by applying the same method, the same result is obtained. Those methods are called deterministical.
On the other hand, stochastical realizations provide a different number of solutions for the same input data set. Those solutions can be very similar, but never identical. It is important to emphasize that all obtained solutions or results are equally probable. In stochastical processes, a number of realizations can be any number we want. It is obvious that more realizations will cover more uncertainty, i.e. the more realizations there are, the lower the uncertainty is. The important question is how many realizations is enough to cover uncertainty? In practice, 100 realizations is the mostly used number. It is considered that 100 realizations should cover 95% of all possible solutions, so it gives a high enough certainty for the reservoir characterizations. Although there are 100 realizations, for practical reasons, only 5 (P5, P25, P50, P75 and P95) or just 3 realizations (P5, P50 and P95) are The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2017©, , pp. 45-53, DOI: 10.1177©, /rgn.2017 most often presented. Realization P5 is the most rigorous and it shows a solution from which 5% of all solutions have smaller values. On the other hand, realization P95 is the most optimistic and it shows a solution from which 95% of all solutions have smaller values. It is also possible to show the 1st, 50th and 100th realization. In this case, selection criterion is statistical and they were chosen by their occurrence as solutions. Many published books and papers describe the theory of geostatistical methods (e.g., Hohn, 1988;Isaaks and Srivastava, 1989;Hand et al., 1994;Damayanti & Hicks, 1996;Deutsch & Journel, 1998;Sahin & Al-Salem, 2001;Malvić, 2008;Geiger, 2012).
A lot of hydrocarbon reservoirs in the Croatian part of the Pannonian Basin System were interpreted by deterministical and stochastical geostatistical methods, e.g. This paper is focussed on indicator based geostatistical methods (deterministical and stochastical), where only third-order stationarity is assumed. That means that variograms are only representative statistical features. Those methods are Indicator Kriging and Sequential Indicator Simulations. Both of them were applied for porosity mapping of the Upper Miocene reservoir in the Sava Depression. Lithologically, the reservoir is represented with sandstone and marl interlayering. In the Croatian part of the Pannonian Basin, the Upper Miocene was the time with considerable environmental changes. Sedimentation took place through cyclic turbiditic fl ows (Vrbanac, 1996;Rogl, 1996Rogl, , 1998. Sandstones were deposited through turbiditic currents in the deepest part of the sedimentation area, forming turbiditic channels. The main source of material was the Eastern Alps. Detritus was transported through turbiditic currents, redeposited several times before it was fi nally deposited. In the calm period, when turbiditic currents were not active, marl, as typical calm deep water sediment, was deposited (Šimon, 1980;Malvić et al., 2005). Sandstone morphology follows the depositional currents direction (e.g., Saftić et al., 2003). In general, the direction is northwest/southeast.
The Indicator Kriging technique is mostly used for facies mapping or facies modelling, which means that the mapping or modelling variable should be discrete or categorical. But, it is also possible to apply this method on continuous variables like porosity, permeability, thickness, etc. In those cases, the obtained results will be probability maps. If porosity is used as a variable for lithofacies defi nition, Indicator Kriging maps would show the probability of certain lithofacies appearing in some location. Methodology includes a series of cut-off values for indicator transformation of a continuous variable into discrete values 0 and 1, and then, after the variogram analysis for each cut-off, indicator probability maps can be performed. The most important thing during mapping is the number of cut-off values. It actually depends on different parameters, i.e. the mapping variable, reservoir heterogeneity, how detailed the results should be, etc. So, it could be said that the number of cut-offs actually depends on personal experience. Still, there is a general rule that for more than 15 input data, the recommended number of cut-off values is 5 to 11, to satisfy statistical criteria of variance calculation, histograms and the stationarity condition (Novak Zelenika et al., 2010).

Principles of Indicator Kriging
The Indicator Kriging method is used for two categorical or indicator variables mapping, represented with the values 0 or 1. It is usually used to map two different lithotypes, i.e. lithofacies with this method. Lithofacies can be defi ned with porosity cut-offs. It could be assumed that higher porosities point to lithofacies with more sandy components like pure sandstones or marly sandstones. On the contrary, lower porosity values can represent sandy marls and marls. The fi rst step in Indicator Kriging mapping for a continuous variable is the indicator transformation of the input data. If a particular facies can be recognized by a ν cut-off , and if the presence of this facies is assigned with the number 1 and its absence with 0, the indicator variable (Equation 1, Figure  1) can be introduced (e.g., Novak Zelenika et al., 2010): (1) where: I(x) -indicator variable, z(x) -measured value, n cut-off -cut-off value.
analysis for each cut-off. The variograms for every cutoff have to be standardized, estimating the results in continuous intervals of [0, 1]. The variogram models which will be used in Indicator Kriging algorithm must obey the following criteria (e.g., Deutsch & Journel, 1998): • The theoretical (model) function must be the same, • The sill must be identical (standardized variogram), • The nugget effect must be the same, • Only the range can change for different indicator variables. Indicator Kriging is a nonlinear method applied on linear Kriging equations, which means that the weighting coeffi cient λ i must be calculated with the Simple or Ordinary Kriging equations. The obtained values are transformed into an indicator and the probabilities are calculated with the following Equation 3: ( 3) where: P*(A, ν l ) -variable probability for the certain cut-off at the A location F*(ν l ) -unbiased estimate of frequency for the cutoff n l ; λ i (n l ) -weighting coeffi cient for cut-off n l . These steps have to be repeated for every cut-off. So, the results of Indicator Kriging method are not just discrete values 0 and 1, but rather a continuous interval [0, 1], i.e. the probability that the estimated value at the certain location is smaller than the cut-off (Novak Zelenika et al., 2010).

Principles of Sequential Indicator Simulations
Sequential Indicator Simulations are a stochastical geostatistical method. Input data for variogram analysis and mapping are original values, the distribution of the original values and cut-offs. Mapping is based on indicator variables as a result of indicator transformation. Sequential Indicator Simulations are conditional and the estimation is performed only at the locations without any measurement. There are several steps in Sequential Indicator Simulations (Juang et al., 2004): (1) defi ning cut-offs and indicator transformation of input data, (2) variogram analysis for every cut-off, (3) cells for estimation are randomly chosen, (4) for the randomly chosen location, the probability that the mapped value is lower than the cut-off is estimated from the Cumulative Distribution Function (abbr. CDF). Sequential Indicator Simulations estimate values close to average. The consequence is that the differences between realizations are not high. This is because the More than one cut-off can be applied in the Indicator Kriging mapping. For l cut-off, Equation 2 for the indicator transformation is: (2) where: z(x) -mapped variable on the x location, I(x, n l ) -indicator variable on the x location for the l cut-off, n lcut-off -l cut-off. Original input data in this way is transformed into l values (0 or 1). The next step is to perform variogram The Mining-Geology-Petroleum Engineering Bulletin and the authors ©, 2017, pp. 45-53, DOI: 10.1177/rgn.2017 variance of the indicator variable is defi ned as (Equation 4): (4) where: -Cumulative Distribution Function of a continuous variable z k , defi ned as (Equation 5): ( 5) where: P -probability, z -total number of data, z k -number of the data to the k value.
Results of the Sequential Indicator Simulation are probability maps which show the probability that the mapped variable is higher than a certain cut-off.

Results of the mapping by Indicator Kriging and Sequential Indicator Simulations
As was described in Novak Zelenika et al. (2010), the input dataset was porosity data measured in 25 wells equally distributed on the fi eld (see Figure 2). Those were average values of porosity in the Upper Miocene reservoir (the Sava Depression).
During the Upper Miocene, sandstones were deposited through turbiditic currents (e.g., Šimon, 1980;Malvić et al., 2005;Malvić, 2012). Porosity from wells varies between 13.8 and 23.3 % in the analysed reservoir. It was assumed that the lowest value represents marly sand and the highest value clean sandstone, deposited in the deepest and central part of the depositional channel. The presumption was that porosity values can indicate different lithofacies (i.e. from marly sandstone to pure sandstone deposited in the channel centre), and therefore the porosity dataset had been transformed in the 6 indicator datasets or classes, based on the following cut-offs: 14, 18, 19, 20, 22 and 24%. The cumulative probability distribution curve is obligatory input for indicator mapping (see Figure 3).
The results of the Indicator Kriging and Sequential Indicator Simulations are the probability maps (see Figures 4 and 5). They are very similar but opposite. Indicator Kriging maps show the probability that a mapped variable is lower than the certain cut-off while Sequential Indicator Simulations show the probability that mapped values are higher than the cut-off. In both cases, the legend is coloured. The blue colour means that there is no probability or very low probability that the cell value is lower (in Indicator Kriging maps) or higher (in Sequential Indicator Simulations) than the selected cutoff (p is close to 0). The red colour shows that cell value is certainly lower (in Indicator Kriging maps) or higher (in Sequential Indicator Simulations) than the selected cut-off (p is close to 1). On the probability maps (see Figures 4 and 5) for the cut-off 19 %, a depositional channel can be assumed.

Advantages of Indicator Kriging and Sequential Indicator Simulations
The question that arises is what are the advantages of Indicator Kriging and Sequential Indicator Simulations? To answer this question, it is actually worth it to repeat the procedure itself. The Indicator Kriging and the Sequential Indicator Simulations are the methods which can be applied to categorical and continuous variables. Considering a variable (rock type, permeability, thickness, porosity, etc.), its indicator transformation can be interpreted as follows (Equation 6): If a variable is continuous, a cumulative probability distribution function has to be created. Since there is a fi nite number of data, which can be changed by adding additional points, this cumulative distribution function is conditional. It is conditioned by the available spatial locations. The total range of the input values has to be subdivided into several classes using cut-offs. On the cumulative probability distribution function, probabilities for all cut-offs can be defi ned. Based on cut-offs, all the data can be divided into two groups: (1) the data which has values lower than the cut-off and (2) the data which has values higher than the cut-off. In this way, the indicator variable is defi ned and it takes a value of 1 at all locations where the value is higher than the cut-off and 0 otherwise. The conclusion is that the larger number of cut-offs, the more precise the cumulative distribution function is. Indicator Kriging and Sequential Indicator simulations provide a least square estimation of the cumulative distribution function at a particular cut-off. The procedure is repeated for a series of cut-off values. Prob-abilities are calculated for a particular grid point. Since they describe a cumulative distribution function, it could be said that using Indicator Kriging methods, cumulative distribution function at a grid point can be estimated. This estimation is repeated for each grid point. The results of the Indicator Kriging and the Sequential Indicator simulations are the probability maps which strongly depend on the correct number of cut-offs. If a variable is categorical, then Indicator Kriging and Sequential In- dicator Simulations give probabilities of its appearance at a particular location.
Shortly, the advantages of Indicator Kriging and Sequential Indicator Simulations are: (1) they do not need normality for the input dataset, (2) they can be implemented in case of bimodal distribution, (3) they can show connectivity of the largest or the smallest values. It is visible in Figures 4 and 5, where connectivity of the values can be described as the most probable location of the depositional channel.

Conclusion
Indicator based geostatistical methods (Indicator Kriging and Sequential Indicator Simulations) are geostatistical methods where variograms are only representative statistical features. Both of them can be applied to categorical and continuous variables. If a variable is continuous, more than one cut-off has to be introduced and a cumulative probability distribution function for every cut-off has to be created. A series of cut-off values are needed for indicator transformation of a continuous variable into discrete values 0 and 1. After the variogram analysis for each cut-off, indicator probability maps can be performed.
Indicator Kriging and Sequential Indicator Simulations were applied for porosity mapping of the Upper Miocene reservoir in the Sava Depression. The input dataset were average porosity data measured in 25 wells. Porosity from wells varies between 13.8 and 23.3 %. Those values actually can represent different lithofacies from sandy marls to pure channel sandstones, through all transitional lithofacies deposited in the channel edge. The input dataset was divided into 6 classes, based on 6 cut-off values. The results of the Indicator Kriging and Sequential Indicator Simulations are the probability maps, which strongly depend on the correct number of cut-offs. Indicator Kriging maps show the probability that a mapped variable is lower than a certain cut-off while Sequential Indicator Simulations show the probability that mapped values are higher than a cut-off (see Figures 4 and 5). Indicator Kriging maps and Sequential Indicator Simulations showed that for certain cutoffs, a depositional channel can be assumed (see Figures 4 and 5, cut-off value 19 %).