# A New Type of Conditioning of Stationary Fields and Its Application to the Spectral Simulation Approach in Geostatistics

Niyaz S. Ismagilov Ismagilov.NS@gazpromneft-ntc.ru

Mikhail A. Lifshits mikhail@lifshits.org

Andrey A. Yakovlev Yakovlev.AAle@gazpromneft-ntc.ru

**Journal «****Mathematical Geosciences****»**

**Abstract**

This paper considers the problem of conditional stochastic simulation in geostatistics. A new type of conditioning method that generalizes the well-known two- step conditioning approach is introduced. The distinctive feature of the new method is that it solves the modeling problem for multiple stochastic fields in a general setup. A generalized kriging procedure is developed in the paper, using linear combinations of the simulated fields as input data instead of the commonly used separate fields’ values. Although the new conditioning method was developed initially in the framework of a particular approach to geostatistical simulation (the spectral method), it can be applied in many more general settings of conditional simulation of stationary stochastic fields of arbitrary dimension. The workings of the method and its applicability are illustrated with the results of several numerical experiments, including simulation on real oil field data.

**Keywords** Geostatistical simulation · Stationary fields · Conditional simulation · Kriging · Spectral simulation · Horizontal wells

# 1. Introduction

Modeling of rock properties in unsampled locations is a difficult problem in petroleum geology. To address this difficulty, a number of concepts and assumptions are applied. One of the key concepts is that of a random variable used to describe various geolog- ical phenomena, including the subject of primary interest in this paper, which is the spatial distribution of geophysical and petrophysical reservoir properties. The random variable concept is also one of the main pillars of geostatistics, a branch of statis- tics focusing on quantitative descriptions of the spatial or spatiotemporal distribution of natural variables. Originating from the pioneering works of G. Matheron and his research group at the Center of Geostatistics in the Paris School of Mines for appli- cations of the mining industry (Matheron 1962, 1965), geostatistical methods rapidly expanded to the oil industry and are nowadays widely used in petroleum geology for the purposes of reservoir characterization.

The main assumptions of geostatistics are stationarity, ergodicity, and a Gaussian distribution of the modeled geological variables (in some models, the isotropy assump- tion is also used). These assumptions have led to the development of two broad classes of geostatistical modeling methods including estimation and simulation techniques. The aim of the methods from the first class, which are also called deterministic meth- ods, is to provide a single estimation model of a geostatistical variable. Usually the obtained model turns out to be overly smooth with respect to reality. One of the most widely used deterministic methods is kriging. Many variations of kriging and the entire family of kriging-based methods are described in Chilès and Delfiner (2012) and Stein (2012). The aim of simulation methods is to generate multiple equiprobable models, called realizations, which represent spatial variability of the simulated phe- nomena. Due to the random nature of the simulated models, these methods are also called stochastic simulation methods. Among the widely applied methods of stochastic simulation, sequential Gaussian simulation (SGS) (Dubrule 2003; Pyrcz and Deutsch 2014) and the Fourier integral method (FIM) are worth mentioning (Pardo-Iguzquiza and Chica-Olmo 1993; Yao 1998).

Each of these two methods represents a special group of conditional simulation methods. The first group, represented by SGS, performs a direct conditional sim- ulation of stationary fields. The second group, represented by the FIM, is a set of two-step methods that start with generating realizations of unconditional simulations and matching the required statistics, but they do not honor the prescribed constraints, for example, the data measured at the wells. A kriging-based smooth additive correc- tion is then calculated, correcting values at the wells so that they honor the constraints (Chilès and Delfiner 2012).

All methods mentioned above are important in the context of the present work, which has two main aspects. First, it presents a description of a geostatistical approach to conditional simulation that is called the spectral approach (method). This method was introduced in Baykov et al. (2010) and was later developed in Baykov et al. (2012) and Khasanov et al. (2014). The motivation behind this method is to loosen the stationarity assumption of the simulated property. The non-stationarity of the data is a common case in practice, usually observed as a variogram not reaching the sill. This issue has been addressed by many researchers (Cressie 1986; Sampson and Guttorp 1992; Machuca-Mory and Deutsch 2013; Boisvert et al. 2009), and the most common approach is to reduce in some way the non-stationary data problem to a stationary case.

The spectral method presents an alternative approach to this problem. It is based on Fourier analysis of the well log data and does not assume stationarity along the vertical axis. Variability patterns along the vertical axis are captured as coefficients of the Fourier series expansion. A short analysis of this feature based on synthetic and real data is presented in Ismagilov et al. (2019).

One of the disadvantages of the spectral method is its reliance on the verticality of the wells used in the simulation process. At this point, the second aspect of the paper becomes relevant. A new conditioning method presented here provides the means to condition spectral simulation models on data from horizontal wells. In this setting, the two-step approach mentioned above is used, but the distinctive feature here is that, due to the specificities of the spectral method’s algorithm, conditioning constraints are imposed on multiple stochastic fields simultaneously, and instead of a simple or ordinary kriging usually used for the calculation of a correction function, an original generalized kriging is utilized. In this generalized setting, the input data for the kriging comprise a linear combination of the values of several stochastic fields at sample points. This provides a means for admitting a broader type of conditioning constraints on simulated fields when equality constraints are imposed, not just on each simulated stochastic field separately but also on their linear combinations.

Although the new conditioning method was developed specifically for the purpose of extending the spectral approach described in the next section, it may certainly find other applications due to its general mathematical nature, and it is interesting in its own right.

The article is organized as follows. In Sect. 2 the basic spectral method of geostatis- tical modeling is briefly described. In Sect. 3 the mathematical setting for the method’s improvement is provided. In Sect. 4 and 5 a solution to the stated problem is given. It can be viewed as a generalized kriging procedure. In Sect. 6 the described approach is illustrated by a series of numerical experiments of increasing complexity, including a model of a real oil field. In Sect. 7 the aspects of computational performance are briefly discussed. Finally, the conclusions are stated in Sect. 8.

# 2. The Spectral Method

The term spectral simulation usually refers to a specific type of stationary field simu- lation process in various applications of the theory of stochastic processes. It involves integral representation of stationary processes based on Bochner’s theorem. This method is also applied to the problem of geostatistical simulation and is usually referred to as the Fourier integral method (FIM), which was already mentioned in the introduction.

The spectral approach presented in this paper has a different origin. It was developed in the works of Baykov et al. (2010, 2012) and later elaborated by Khasanov et al. (2014) as an attempt to overcome the disadvantages of the classical methods mentioned above, imposing weaker stationarity and isotropy assumptions on the considered three-dimensional property. For an interesting alternative example of the spectral approach, we refer the reader to Zagayevskiy and Deutsch (2016).

The main idea is based on Fourier analysis of well logs (hence the term spectral) followed by the estimation of the correlation structure of the Fourier expansion coef- ficients and their conditional simulation. This is achieved through the same FIM and kriging conditioning.

Here a rather short description of the spectral approach is presented, primarily in order to introduce the method to a wider audience (the cited works are available only in Russian). For more details, see Baykov et al. (2012). The extended setting considered in this paper relies on details of the spectral method and was developed as an improvement to it.

The expansion of Eq. (1) to (2) can be performed for every point x ∈ X, thus giving a parametrized expansion of the field G

represents random fields on X . The right-hand side of Eq. (2) is a limit of a sum of random variables and can be considered as a Gaussian random variable (Baykov et al. 2012). This leads to the assumption of the Gaussianity of random fields c j(x).

Further, according to the spectral method, random fields c _{j} (x) are assumed to be quasi-stationary (Baykov et al. 2012) for more detailed discussion), which means that

where functions a j(x) and σj(x) are supposed to be deterministic and represent the mean and the variance of the stochastic field c j(x), while ηj(x) is a centered stationary random field with unit variance and some correlation function K j(·). The latter, coupled with the Gaussianity assumption of the stochastic fields c j(x), puts the model into the general geostatistical setup of simulation of Gaussian stationary fields.

The problem is to simulate stochastic field G in the form given by Eq. (3), which honors (or is conditioned by) data along these S vertical wells. Modeling includes estimation of the mean a j(·), the variance σj(·), the correlation functions K j(·), and the simulation of the Gaussian stationary fields ηj(x) having correlation functions K j(x) and conditioned by the values ηj(xs) calculated by Eq. (2) using the well log data.

Furthermore, in practice, evaluation of an infinite number of Fourier coefficients is not possible. Therefore, the right-hand side of Eq. (3) is substituted by a finite sum of J elements. The problem of choosing an appropriate J can be a subject for discussion, but it makes sense to choose it such that the remainder of the series has a small L2(H) norm. The value of J used in practice is limited mostly by the available computational power. Theoretical estimation of a reasonable value for J , as well as an appropriate choice of a basis (φj) are independent mathematical problems of their own interest in the field of approximation theory. In practical applications, the value of J depends on the amount of information contained in the analyzed well log section. Higher values of J lead to reproducing more details of the well log data, but overly high values may result in capturing noise from data. In particular, for computations on the real oil field model presented in Sect. 6, the value of J was chosen to render the L2(H) norm of the remainder below the 2% level with respect to the L2(H) norm of the whole series. In terms of signal processing, this choice can be interpreted as a low energy level value of the remainder, and suggests that the contribution of the remainder to the simulation result can be neglected.

The stationary fields ηj(x) conditioned by the values expressed in Eq. (6) still need to be simulated. One of the known ways to do this is given by the so-called two-step method. According to this method, first, an unconditioned version of ηj(x) is simulated, for example, using the spectral representation of the stationary fields (Prigarin 2005). The resulting field is then conditioned by adding a correction term built via the kriging procedure (Chilès and Delfiner 2012).

# 3. Problem Statement

It should be emphasized that the version of the spectral method described above is completely based on the assumption of the verticality of wells. However, nowadays petroleum engineering deals with low permeability and high heterogeneity reservoirs, which are developed by means of non-vertical wells, including slant and horizontal ones. Data acquired from these wells must also be used in geological simulation. It should be noted that involving data from horizontal wells in the process of geostatistical simulation is not a problem for conventional pixel-based methods such as SGS and the FIM, which naturally integrate those data. The main goal of the present work is to additionally develop this feature for use in the spectral method described in the previous section.

It seems reasonable to solve the problem in a context broader than that of nonvertical wells: let us consider the measurement results at some additional points

Mathematically, the problem can be formulated in the following way. For a random field of the form expressed in Eq. (3), given its values expressed in Eqs. (4) and (7), it is necessary to estimate its parameters a(·) and σ (·) and the spectral characteristics of the fields ηj (the ideas and methods regarding these estimations can be found in Rozanov (1967) or Baykov et al. (2012) and then simulate a sample path of G conditioned by the values expressed in Eqs. (4) and (7).

# 4. Generalized Kriging

The problem stated above can be solved by a two-step kriging method of conditional simulation (Baykov et al. 2012; Chilès and Delfiner 2012) widely used in geostatistics. According to this method, in the first step, an independent unconditioned sample path for a random field is generated; in the second step, it is corrected via kriging in order to fit the hard data (external measurements). The essential difference between the present setting and the classical one is that in traditional kriging, hard data are represented by the values of the simulated fields at certain locations, while in the present setting this role is played by linear combinations of these values. The idea behind the approach proposed in this paper is that, from the point of view of the theory of Gaussian random functions, there is no major difference between the two types of conditions—both are linear functionals of the simulated field.

The main technical result is formulated as follows.

*Remark*The theorem can easily be generalized to the case of dependent fields ηj and to the case when the field values in the same linear form Li are picked from different space points, but for the current purposes, the above simpler setting is sufficient. It should also be noted that the assumption of stationarity of the fields is not used.

# 5. Building a Conditional Random Field

*Remark*Although the described algorithm is designed to work with any point measurements and could ingest exact coordinates and values at well log scale along the horizontal wells, this approach is not very practical. In practice, one has to deal with a simulation grid, and it makes sense to pass coordinates of the cells along the horizontal wells as coordinates of point measurements and values of well log upscaled onto the grid cells as the values of the point measurements. The choice of the actual upscaling algorithm depends on the simulated property and remains beyond the scope of this paper.

# 6. Numerical Experiments

In this section, the results of three numerical experiments of increasing complexity are provided in order to illustrate the theory developed in the article and to confirm its practical applicability.

First, in Sect. 6.1, a numerical experiment for an intentionally simplified synthetic one-dimensional model is presented. The idea was to show how the stochastic fields involved in the mathematical model change with the conditioning procedure applied. This can be most easily done for one-dimensional fields (processes), hence the reduc- tion to one dimension. In fact, this experiment can be easily reproduced by any interested reader.

What follows is the three-dimensional experiment of Sect. 6.2, giving an idea of how the method extends to the case of three-dimensional simulations, which is the natural scale of real geological models. Simulation data, again, were chosen to be synthetic, so the necessary validations of the correlation structure of stochastic fields involved in the process of simulation could be performed. This experiment, although demanding more effort, can still be reproduced.

Finally, a real model simulation is described in Sect. 6.3. It demonstrates that the proposed conditioning method can be successfully applied to practical simulation in a real oil field model. A brief analysis of the simulation results shows how the new data from the horizontal wells are introduced into the model.

**6.1 One-Dimensional Example**

**The following figures illustrate the results of numerical computations. The figures with sample paths of the unconditioned processes η1(x) and η2(x) simulated independently are not shown, as they are not very informative and do not look much different from “the nature” in Fig. 1.**

**6.2 Three-Dimensional Example**

Next, we extend the experiment to a synthetic model that is closer to practical cases of simulation via the spectral method. Let X be a rectangular parallelepiped of a shape that is 10, 050×10, 050×20 meters, representing a geological layer, which is gridded into 201 × 201 × 100 cells, where each cell has a size of 50 × 50 × 0.2 meters. These sizes agree with the scale of typical geological models. According to the spectral method, the geostatistical property is represented as a Fourier sum given by Eq. (3). In this example, Legendre polynomials (Pj) are used as the basis functions, and the sum is limited to four terms

As the analyzed properties are three-dimensional, so are the differences, but threedimensional properties are difficult to visualize on a two-dimensional image, so Figs. 7(a)–7(c) demonstrate horizontal sections of the differences at the depth level of the horizontal well. In the figures, darker regions correspond to lower differences and lighter regions correspond to larger differences. The red line is the horizontal well and the orange point is the vertical well.

One observes that the “nature” and the preconditioned version of the second simulation differ rather randomly, including the area near the wells (Fig. 7(a). Comparison of the “nature” and the conditioned version of the second field (Fig. 7(b)) shows that the difference between the two fields is preserved, apart from the region close to the well, where the difference reduces to zero. This fact is supported by Fig. 7(c), which provides the difference between the preconditioned and the conditioned version of the second field. One can observe that the two fields differ only near the wells, with the difference going down when the distance from the wells increases

**6.3 Real Data Model**

In order to demonstrate practical applicability of the method, a case of a real oil field with a number of horizontal and vertical wells is considered. The model is represented by a geological layer with a size of 6,000 × 6,000 × 13.5 meters gridded into 120 × 120 × 50 cells, where each cell has a size of 50 × 50 × 0.27 meters. The oil field model has 36 vertical (S = 36) and 10 horizontal wells drilled into the target layer. The layout of the model is given in Fig. 9. One has to simulate the gamma-ray data based on the measurements taken along the vertical and horizontal wells. The simulations were performed using the spectral method and the conditioning algorithm described in Sect. 5. For the simulations, the number of Fourier coefficients was chosen to be equal to 150 ( J 150), which, together with 418 conditioning points ( P 418) on 10 horizontal wells, leads to the linear system expressed in Eq. (13) consisting of I 5, 818 equations. This system is solved numerically via standard methods provided by linear algebra software packages and represents the most computationally costly part of the simulation.

As in the previous section’s three-dimensional example, two simulation versions were considered: the first one is simulated conditioned only on vertical wells and represents the usual spectral simulation algorithm result, while the second one is conditioned on the vertical as well as on the horizontal wells.

In order to eliminate the effects of randomness inherent to stochastic simulations, the seed for the pseudorandom number generator was fixed for both simulations, so the two simulations can be appropriately compared. A model of the three-dimensional property simulated conditioned by horizontal wells data is shown in Fig. 10. A model conditioned only on vertical wells looks similar and thus is not shown.

To assess the quality of conditioning on input data, both models were sampled along one of the horizontal wells to produce synthetic logs. These synthetic logs are compared to real logs, which had been passed as simulation data (see the graph in Fig. 11). The black line on this graph corresponds to the real well log, the red line corresponds to the first model, and the blue line corresponds to the second model. There is a difference in the sampling of the real log and synthetic logs, which is explained by the fact that the sampling of the real log is determined by the actual device used to make measurements, whereas gridding of the synthetic logs is determined by the geological model and the geometry of the well trajectory intersection with the grid cells.

In Fig. 11, it can be seen that the blue line corresponding to the second model reproduces the well log data with precision that can be achieved on the coarse grid of the geological model. In fact, the blue line reproduces upscaled values of the well log, which in practice substitute the actual well log as mentioned earlier. For this particular case, the arithmetic mean was chosen as the upscaling algorithm.

The red line on the same figure corresponds to the first model and does not match the well data, as one may expect. It is worth emphasizing that the applied conditioning procedure is such that the red line is converted into the blue line, and the fact that the blue line matches the well data means that the procedure indeed conditions the model.

The effect that the new data from horizontal wells has on the model can be seen in Fig. 12, where vertical sections of the simulated property along one of the hori- zontal wells (the red line) are shown. Figure 12(a) corresponds to the first model, and Fig. 12(b) corresponds to the second model. As the two models were simulated with the same seed, the absolute value of difference between them shown in Fig. 12(c) high- lights the region where conditioning on horizontal well affects the model and describes the magnitude of this effect. The lighter color indicates a bigger difference and the darker color corresponds to the regions with smaller or no difference between the two

models. The figure indicates that the main impact is produced near the trajectory of the horizontal well.

Similarly, the impact zone of conditioning on horizontal wells can be observed in Fig. 13, where a map of the cumulative absolute value of the difference between the two models is presented. The darker areas correspond to the areas where the difference between the two models is more significant, and the lighter areas indicate the absence of such a difference. These observations show that the impact zones are located near the trajectories of the horizontal wells and are limited by the range of the covariance functions. This fact agrees with the theoretical derivations presented in Sect. 5.

# 7. Performance of the Spectral Method

Some comments should be added regarding the performance of the new method, compared to the classical pixel-based simulation methods. The answer here is not straightforward. The computational complexity of SGS is defined by the size of the simulation grid. In its canonical formulation, SGS has the computational complexity of O(N ^{4}) floating point operations for a grid with N nodes. This is impractical when the grid size is on the order of millions of nodes. Therefore, the practical implementations of the SGS algorithm use local neighborhood estimations to reduce computational complexity to the order of O(Nν^{3}) floating point operations, where ν is the number of neighbor nodes (Dimitrakopoulos and Luo 2004).

The computational complexity of the “canonical” (that is, based only on vertical wells) spectral method has not been evaluated, but the practical simulations showed that the spectral method simulates twice as fast as SGS simulations carried out on commercial software for models larger than four million cells (for models less than a half million cells, the spectral method performs worse than SGS). This is achieved due to the structure of the algorithm admitting parallelization: each Fourier coefficient could be simulated on an individual node; the stationary field simulation algorithm itself is also parallelizable.

For this reason, a direct complexity comparison between the two methods is not completely appropriate, and one should not infer strong conclusions from numerical experiments. In fact, this comparison, as well as the computational efficiency of the spectral method, are subject to subsequent research.

# 8. Conclusions

The present study demonstrates a methodology of conditional simulation involving data from slant and horizontal wells in the process of geological simulation based on the spectral simulation method. The conditioning method introduced in the article significantly widens the application scope of the spectral method, as the number of oil fields developed by horizontal drilling is increasing rapidly.

The conditioning method is implemented in a full scale software package suitable for practical applications. In particular, it enables simulation of a wide range of geo- logical models, such as the one described in Sect. 6.3. Practical simulations of real data showed that models produced by the software are of high quality. Therefore, the software is now used in geological modeling.

More research should be carried out to optimize the details of the workflow, reducing the computation time and memory consumption for simulation of very large geological models. In particular, the implementation of the ideas of sparse Gaussian regression (Quinonero-Caneda and Rasmussen 2005) for reducing the problem’s dimension is already underway.

It is also likely that the mathematical core of the proposed method, namely, the generalized kriging, may have many different applications beyond the spectral method used here.

Acknowledgements The authors are grateful to both referees for careful reading and their useful advice. Research was partially supported by Native Towns, a social investment program of PJSC Gazprom Neft.

# References

Baykov VA, Bakirov NK, Yakovlev AA (2010) New methods in the theory of geostatistical modeling.

Vestnik UGATU 14(2):209–215 (Russian)

Baykov VA, Bakirov NK, Yakovlev AA (2012) Mathematical geology, vol 1. Introduction to geostatistics.

Ser. Petroleum Engineering Library, IKI, Moscow-Izhevsk

Boisvert JB, Manchuk JG, Deutsch CV (2009) Kriging in the presence of locally varying anisotropy using non-Euclidean distances. Math Geosci 41(5):585–601

Chilès J-P, Delfiner P (2012) Geostatistics: modeling spatial uncertainty, vol 713, 2nd edn. Wiley series in probability and statistics. Wiley, Hoboken

Cressie N (1986) Kriging nonstationary data. J Am Stat Assoc 81(395):625–634

Dimitrakopoulos R, Luo X (2004) Generalized sequential Gaussian simulation on group size ν and screen- effect approximations for large field simulations. Math Geol 36(5):567–591

Dubrule O (2003) Geostatistics for seismic data integration in earth models. EAGE, Houten

Ismagilov N, Popova O, Trushin A (2019) Effectiveness study of the spectral approach to geostatistical simulation. In: SPE annual technical conference and exhibition, 30 September–2 October, Calgary, Canada

Khasanov MM, Belozerov BV, Bochkov AS, Ushmayev OS, Fuks OM (2014) Application of the spectral theory to the analysis and modeling of the rock properties of the reservoir. Oil Ind J 12:60–64 (Russian)

Machuca-Mory DF, Deutsch CV (2013) Non-stationary geostatistical modeling based on distance weighted statistics and distributions. Math Geosci 45(1):31–48

Matheron G (1962) Treaty of applied geostatistics, vol 1, 1st edn. Technip, Paris (French) Matheron G (1965) Les variables régionalisées et leur estimation. Masson et Cie, Paris (French)

Pardo-Iguzquiza E, Chica-Olmo M (1993) The Fourier integral method: an efficient spectral method for simulation of random fields. Math Geol 25(2):177–217

Prigarin SM (2005) Numerical modeling of random processes and fields. ICM&MG Publisher, Novosibirsk (Russian)

Pyrcz MJ, Deutsch CV (2014) Geostatistical reservoir modeling, 2nd edn. Oxford University Press, Oxford Quinonero-Caneda J, Rasmussen CE (2005) A unifying view of sparse approximate Gaussian process regression. J Mach Learn Res 6:1939–1959

Rozanov YA (1967) Stationary random processes. Holden-Day, San Francisco

Sampson PD, Guttorp P (1992) Nonparametric estimation of nonstationary spatial covariance structure. J Am Stat Assoc 87(417):108–119

Stein ML (2012) Interpolation of spatial data: some theory for kriging. Springer, New York

Yao T (1998) Conditional spectral simulation with phase identification. Math Geol 30(3):285–308 Zagayevskiy Y, Deutsch CV (2016) Geostatical grid-free simulation of natural phenomena. Math Geosci 48(8):891–920