Comparison of Ordinary and Universal Kriging interpolation techniques on a depth variable (a case of linear spatial trend), case study of the Sandrovac Field

Abstract

Until now, Universal Kriging has not been used for the mapping of geological data in Croatia. However, it is one of the most frequently used methods of Kriging, probably the most adequate in cases when the input data is marked by a common trend. That exact feature is often an attribute of deep geological data, and thereby that of structural maps. Mapped surfaces in a row of examples have a structural trend towards one cardinal direction, or a sequence of geological structures, like anticlinorium, is a part of a structural unit of a higher order such as regional monocline. An example is given of geographical trend recognition in e-log Z' surface spread in Sandrovac Field as well as successful mapping of that marker depth variable by using Universal Kriging.

Headnote

Abstract

Until now, Universal Kriging has not been used for the mapping of geological data in Croatia. However, it is one of the most frequently used methods of Kriging, probably the most adequate in cases when the input data is marked by a common trend. That exact feature is often an attribute of deep geological data, and thereby that of structural maps. Mapped surfaces in a row of examples have a structural trend towards one cardinal direction, or a sequence of geological structures, like anticlinorium, is a part of a structural unit of a higher order such as regional monocline. An example is given of geographical trend recognition in e-log Z' surface spread in Sandrovac Field as well as successful mapping of that marker depth variable by using Universal Kriging.

Key words

Universal Kriging, Ordinary Kriging, trend, structural map, Sandrovac Field

Sazetak

Usporedba obicnoga i univerzalnoga krigiranja na primjeru interpolacije dubinske varijable s lineamim trendom, analiza polja Sandrovac

Univerzalno krigiranje do sada nije koristeno pri kartiranju geoloskih podataka u Hrvatskoj. Medutim, radi se o jednoj od cescih metoda krigiranja, vjerojatno najprimjerenijoj u slucajevima kada su ulazni podatci obiljezeni zajednickim trendom. Upravo takvo obiljezje vrlo je cesto svojstvo dubinskih geoloskih podataka, a time i strukturnih karata. Kartirane plohe u nizu primjera imaju strukturni nagib prema jednoj strani svijeta. Cesto je niz geoloskih struktura, poput antiklinorija, dio strukturne jedinice vecega reda poput regionalne monokline. U radu je dan primjer prepoznavanja zemljopisnoga trenda u pruzanju plohe EK markera Z' u polju Sandrovac te uspjesnoga kartiranja te plohe, odnosno dubinske varijable, univerzalnim krigiranjem.

Kljucne rijeci

univerzalno krigiranje, obicno krigiranje, trend, strukturna karta, polje Sandrovac

(ProQuest: ... denotes formulae omitted.)

1.Introduction

Geostatistics is widely used in geology as an interpolation approach for mapping different geological variables like depth, thickness, porosity or permeability. One of the most used estimators for the interpolation of spatial data are Kriging techniques. They are considered an interpolation method for the estimation of a regionalized variable at selected grid points that predict values from interpolation without bias, and with minimum variance. There are different types of Kriging techniques, such as Ordinary Kriging (abbr. OK), Universal Kriging (abbr. UK), Indicator Kriging, Co-kriging and others. The most commonly used method is OK, but the choice of which kriging to use depends on the characteristics of the data.

Geostatistical interpolation and variogram analysis has been the standard analytical tool in Croatian geology over the last 10-15 years. Even so, UK has not been used yet for the mapping of geological data in Croatia. It is a method used for the estimation of spatial means when data has a strong trend and the trend can be modelled by simple functions. UK should also be used when the dependent variable does not meet the criterion of second order stationarity necessary for kriging. Second order stationarity means that the mean and variance are the same on the entire area and that the correlation between any two observations depends only on their relative position in space. Even so, UK also has its limitations, like a surface area. If we analyze a large surface, like a big country or a continent, it is difficult to follow a trend along the direction of spreading, so in that case it is not advisable to use a UK. However, when describing a structural local phenomenon, this technique should be applied and the results compared to other kriging techniques in order to decide whether to use it or not.

The aim of this paper is to describe the UK method and introduce it as another statistical tool for mapping variables in Croatian geology. Sandrovac Field was selected as a case study. In this area, the OK and UK methods will be compared and described by mapping a depth variable, and their advantages and disadvantages will be discussed.

2.Study area

The entire Sandrovac Field covers a small surface area of about 38 square km, located 18 km east from the city of Bjelovar. It represents a classical uplifted structure which is geographically placed in the Bjelovar Subdepression (see Figure 1), a separate regional geotectonic unit in the Croatian part of the Pannonian Basin (abbr. CPBS).

The Bjelovar Subdepression is situated in the south-western part of Drava Depression. It is partially isolated due to uplift of Bilogora Mt. during the Pliocene, Pleistocene and Holocene. The Subdepression is also separated from the Sava Depression by Moslavacka gora Mt. On a larger scale, Bjelovar Subdepression is located on the southern margin of the Pannonian Basin System (abbr. PBS). The creation of PBS has been previously explained in many papers (e.g. Royden 1988, Rögl 1996; 1998), as were the specifics of the Bjelovar Subdepression (e.g. Malvic 1998; 2003; 2004; 2011, Pavelic 2001).

The opening of the Bjelovar Subdepression was a result of the activities at the main transcurrent fault system together with other fault systems diagonal or transversal to subdepression (e.g., Malvic 2003; 2011) or to the direction of the Central Drava Fault. In the Early Badenian (16.4-15.0 Ma) sedimentation took place in alluvial fan environments and its rate was lower than the rate of transtension during this time. The 1st transtension sediments were breccias, conglomerates and sandstones. Sedimentation changed to fine-grained clastics and carbonates during the Middle and Late Badenian (15.0-13.0 Ma). At the beginning of the 1st transpression (during Sarmatian-Early Pannonian, 13.0-9.3 Ma), salinity decreased due to a reduced sea area, and brackish environments were formed. The 2nd transtension lasted from the Late Pannonian (9.3-7.1 Ma) to the Early Pontian (7.1-6.3 Ma). It was marked by a deposition of monotonous series of sandstones and marlstones.

The 2nd transpression started in the Late Pontian (6.3-5.6 Ma) and was characterized by different weakly consolidated and unconsolidated clastics with sporadical peat and coal (e.g. Malvic 2011, Mesic Kis and Malvic, 2014).

The Bjelovar Subdepression was not situated on the main direction of transposal of materials, particularly in the transpressional periods for which the sediment income was significantly lower. Therefore, the thickness of the Neogene-Quaternary clastics above the Mesozoic carbonate basement or Paleozoic igneous/metamorphic basement, rarely exceeds 3000 m in contrast to 7000 m in the central Drava Depression zone.

The area of the Bjelovar Subdepression is considered a mature petroleum province. Accordingly, Sandrovac is an oil - gas field that was put into production since 1967. Its exploitation was in accordance with the strategic plan of INA d.d. (Zornjak, 2009).

3.Methods

Dataset was provided from the doctoral thesis and an article published by Malvic (2003, 2011). Figure 2 shows a palaeostructural map of marker Z' from which 18 well data coordinates were taken. E-log marker Z' was selected as a mapping horizon due to its chronostratigraphically importance in the Drava Depression. It is situated approximately in the middle of the 2nd transtension and represents a border between the Pannonian and Pontian (and between the formations of Ivanic Grad and Klostar Ivanic).

Table 1 shows well coordinates and values of the "depth" variable that were used for variogram analysis and as input data for Kriging techniques. The dataset was analyzed using OK and UK. Since a variogram analysis is necessary for map interpolation using Kriging methods, it is also briefly described. The emphasis will be on UK, an interpolation method that has not been used for the mapping of geological data in Croatia.

3.1.Variogram analysis

A variogram (2γ) is one of the basic geostatistical tools that is used to determine spatial dependence. It is often referred to as a semivariogram (γ), which has exactly the same characteristics, except that in the denominator of the equation, the number 2 is eliminated. A variogram is mathematically expressed by Equation 1:

... (1)

Where:

N(h) - number of data pairs at distance "h" (inside searching neighbourhood area),

z(") - value at location n,

z(n+h) - value at location n+h.

Calculation of an experimental variogram is necessary input for different geostatistical interpolation techniques, like Kriging. An experimental (semivariogram) curve is shown in Figure 3. On such a curve, the following parameters can be read: nugget - Co, sill (variance) - C, range - a, and distance - h.

A given experimental variogram is then approximated by the theoretical model (see Figure 4). In the interpretation of geological variables, commonly used theoretical models are spherical, exponential and a Gaussian model.

The following equations determine the behavior of a certain variable that was approximated by a selected theoretical model (e.g. Isaaks i Srivastava, 1989):

... (2)

... (3)

... (4)

Where:

y(h) - variogram value at distance "h"

C - constant,

h - distance,

a - range.

3.2.Ordinary Kriging

OK is one of the most commonly used Kriging techniques, a geostatistical interpolation technique that is described by the acronym BLUE - "best linear unbiased estimator". It is the "best" because it aims at minimising the variance of the errors, "linear" because its estimates are weighted linear combinations of the available data, and "unbiased" since it tries to have the mean residual or error equal to 0.

The principle of Kriging is to estimate values of a regionalized variable at a selected location (Z¿), based on the surrounding existing values (Z,). Selected locations are assigned a relevant weighting coefficient (λ,) which represents the influence ofparticular data on the value of the final estimation at the selected grid node. Variogram values express the relationship between the existing (hard) data and the estimation point, or by covariance in case of second order stationarity (e.g. Malvie and Balie, 2009).

OK premise is a constant unknown mean in the local neighborhood of each estimation point. Local variance of the data within the search ellipsoid is used for estimation, which is useful in the case of a small number of input data (15 or 20). Then, the global variance often does not reflect local changes, so deviations of the mean and estimation can be large. In the OK technique, the amount of kriging variance is minimized using a linear external parameter called the Lagrange factor (μ). The limiting factor minimizes error and assessment becomes impartial. The condition when assessing the OK technique is that the sum of all weights is equal to 1 (e.g. Malvíe et al., 2008). The OK equation in matrix form is:

... (5)

Where:

γ- variogram values,

zi....zn - real value at location 1 to n,

x - location where new value is estimated,

μ - Lagrange factor.

3.3.Universal Kriging

The aim of the UK method is to predict Z(x) at an unsampled area as well. It splits the random function into a linear combination of deterministic functions, the smoothly varying and nonstationary trend, that is also called a drift μχ [varepsilon] R, and a random component Y(x): = Z(x) - μ(χ) representing the residual random function (Wackernagel, 2003).

OK assumes a stationary, i.e. constant mean of the underlying real-valued random function Z(x). But in reality, the mean value varies, it is often not constant across the entire study area and the variable is said to be nonstationary. A nonstationary regionalized variable can be considered as having two components (Davis, 1973); drift (average or expected value of the regionalized variable) and a residual (difference between the actual measurements and the drift).

The method of UK assumes that the mean m(x) has a functional dependence on the spatial location and can be approximated by a model with the equation (e.g. Kumar, 2007):

... (6)

Where:

ai - ith coefficient to be estimated from the data,

fi - ith basic function of spatial coordinates that describes the drift,

k - the number of functions used in modelling the drift.

UK is also known as kriging with a trend or kriging in the presence of a drift. Spatial trend or a drift represents any detectable tendency for the values to change as a function of the coordinate variables.

The mean can be a function of the coordinates in linear, quadratic or higher form. The UK considering linear trend was used for depth spatial interpolation. The mathematical model with two variables is (e.g. Kastelec and Kosmelj, 2002):

... (7)

Where:

Z(s) - random variable at location s,

xi(s) - longitude at location s,

x2(s) - latitude at location s,

βο, βι, β2 - coefficients of the linear trend,

S(s) - intrinsic stationary random process with existing variogram 2y(h), where h = Si - Sj is a distance vector between locations st and Sj.

The kriging estimate is as follows (e.g. Kumar, 2007):

... (8)

Where:

Z' (so) - estimated value at so,

Z(Si) - observed values at points si,

N - sample size,

λ - weights chosen for si so as to satisfy following statistical conditions:

(a) unbiasedness

... (9)

(b) minimum variance

... (10)

According to Journel and Huijbregts (1978), the minimisation of Equation (10) that is subject to the constraint of Equation (9), using the Lagrange multiplier (μ), results in the following UK system:

... (11)

Where:

y(si, sj) - semivariogram between two points si and sj,

μι - Lagrange multiplier associated with the Ith unbiased condition.

The optimum weights λι can be obtained by solving these equations simultaneously. For the calculation of λ-s, a variogram which is a measure of spatial continuity of the data, is required (Kastelec and Kosmelj, 2002). The variance of this estimation is given in Equation 12 (e.g. Kumar, 2007):

... (12)

The process of using UK can be described through the following steps:

1) it is necessary to understand why the trend exists based on the nature of our data,

2) it is required to use a simple form of the trend if possible and avoid extrapolation beyond available data,

3) after the trend selection, it is subtracted from the observed data to obtain the residuals,

4) those residuals are used for variogram computation and spatial estimation,

5) and finally, kriged residuals are added back to the trend (e.g. Isaaks and Srivastava, 1989).

4.Results

Results include maps of variogram surface for variables "depth" and "residuals", experimental and theoretical variogram models, maps interpolated using OK and UK that were also compared using cross-validation. Variogram surface maps, as well as experimental and theoretical variogram models were obtained using Variowin 2.0, a software for spatial data analysis in 2D. Maps were interpolated using the licensed version of the Surfer 8.0 program.

4.1. Ordinary Kriging

A prerequisite for using a kriging method is the construction of experimental variograms. A variogram surface map for the variable "depth" was created (see Figure 5) with the following parameters:

X: lag spacing: 1000 m number of lags: 6

Y: lag spacing: 1000 m number of lags: 6

A variogram surface map calculated for a "depth" variable indicates 2 axes: a primary axis that has a NW-SE trend (135°-315°) and a secondary axis with a NE-SW trend (45°-225°). Primary and secondary axes were defined according to the structural map (see Figure 2).

Figure 6 shows an experimental variogram and a Gaussian theoretical variogram model of Sandrovac Field for a "depth" variable (primary axis), with the following parameters:

nugget 0

range 2280 m

sill 45000

Figure 7 shows an experimental variogram and Gaussian theoretical variogram model of Sandrovac Field for a "depth" variable (secondary axis), with the following parameters:

nugget 2300

range 4000 m

sill 33120

According to the structure of the Bjelovar Subdepression, the ratio of primary and secondary axis is 2:1. Using the variogram data, a map was interpolated by OK using Surfer 8.0(TM) (see Figure 8).

By the method of cross-validation MSE was calculated for a "depth" variable, which is 29269. The root mean square error of the variable "depth" is 171. The most underestimated data is at X = 6423000, Y = 5085000, Z (depth) = 1700 where the estimated value of depth was 1419. The most overestimated data is at X = 6422000, Y = 5088200, Z (depth) = 1040 where the estimated value of depth was 1552.

4.2. Universal Kriging

Spatial correlation between variables can been shown in a scatter plot, if there is such. That is also the starting point in the regression analysis. Figure 9 shows a positive linear trend of depth with longitude and a negative linear trend with latitude. In this case, Gauss-Kruger coordinates are independent variables and depth is a dependent variable. As mentioned before, a variogram for interpolation using UK is constructed using different variables - "residuals".

Since there are two independent variables, multiple regression analysis is used for trend removal. An estimated model has the following Equation 13:

... (13)

Estimation of β coefficients is given by resolving normal equations of least square method:

...

which means that the estimated regression equation is:

... (14)

If real values of independent variables x are included in an estimated regression equation for every i = 1, 2, ..., 18; regression values y of dependent variable y are given. Accordingly, the first regression value y is given by including xi (longitude) and X2 (latitude) variables from Table 1 into Equation 14:

...

Residuals are given by subtracting the real value of dependent variables from the estimated. Specifically, the first residual value is:

...

Variogram surface map for the variable "residuals" was created (see Figure 10) with following parameters:

X: lag spacing: 800 m number of lags: 5

Y: lag spacing: 800 m number of lags: 5

A variogram surface map calculated for the "residual" variable indicates 2 axes: a primary axis that has a NW-SE trend (135°-315°) and a secondary axis with a NE-SW trend (45°-225°), same as for the "depth" values.

Figure 11 shows an experimental variogram and spherical theoretical variogram model of Sandrovac Field for a "residual" variable (primary axis), with the following parameters:

nugget 0

range 1638 m

sill 2450

Figure 12 shows an experimental variogram and spherical theoretical variogram model of the Sandrovac Field for a "residual" variable (secondary axis), with the following parameters:

nugget 0

range 2280 m

sill 2470

Using the variogram data, a map was interpolated by UK using Surfer 8.0(TM) (see Figure 13).

MSE for "depth" variable is 3473. The root mean square error of the variable "depth" is 59. The most underestimated data is at X = 6419000, Y = 5088000, Z (depth) = 950 where the estimated value of depth was 815. The most overestimated data is at X = 6422000, Y = 5088200, Z (depth) = 1040 where the estimated value of depth was 1127.

5.Discussion

Figure 7 shows that there is no sill intersecting in the experimental and theoretical Gaussian variogram model of the secondary axis for a "depth" variable, which means that we have an extremely large spatial dependence. Increasing the number of lags also did not lead to sill intersecting, so the range value of 4000 m was used. According to a regional structural map of e-log marker Z' (see Figure 2), the primary axis has a direction of 135-315°. Adirectional variogram of a primary axis for depth variable was used for OK interpolation, and the ratio was 2:1 since the structure is 2 times longer than wide.

The resulting OK map (see Figure 8) with a Gaussian variogram model does not respect the rule (even geological mapping axiom) that lines of some values should not pass between two points (as the minimum and maximum), where there isn't some intermediate value, through which the line can pass. That is considered as a basic rule of simple interpolation. The second rule of interpolation is that isolines with greater or less values than the maximum and minimum of the measured points, cannot occur on the map. Both rules are not respected on this map. Therefore, this map cannot be considered good, which can be confirmed with a bigger MSE value.

The Sandrovac Field was already analysed by the OK method in a previous study (Mesic Kis and Malvic, 2014) where a linear variogram model was used because it was proposed as a "default" variogram model in SAGA GIS software, and it was compared with Thiessen polygon method. Such a model is usually considered "handy" because it describes regression dependence with a line and does not use any complex variogram "mathematics". Therefore, such krigings are usually considered better than inverse distance or other simple methods because simple regression can be much more easily applied in plain interpolation algorithms. In fact, such kriging implies that there is a general trend in data that is described by simple regression with one dependent variable.

Figure 9 shows that we have two independent variables (x and y coordinates) and one dependent variable (depth). There is a positive linear trend of depth with longitude and a negative linear trend with latitude. Since UK assumes a general polynomial trend model such as a linear trend model, it was used for this study. The problem with 2 independent variables is solved using multiple linear regression which leads to a regressional equation from which you can calculate the value of the regression variables and residuals. A variogram was calculated for a "residuals" variable. Since we had a linear trend, the option "linear drift" in Surfer was chosen to implement the UK. The resulting UK map (see Figure 13) shows a homocline, where the entire Sandrovac Field sinks towards the SE as well as the structure of the entire subdepression.

The reason why Figures 8 and 13 look so different lies in the fact that there is a linear trend in the data. The UK method identified that trend and made a better interpolation. Previous work (Mesic Kis and Malvic, 2014) shows a similar result of the OK method with the UK method used in this work due to a "default" linear variogram model option in SAGA GIS software.

6.Conclusion

This paper describes which kriging technique is more suitable for depth interpolation in Sandrovac Field with e-log marker Z' structural map as input data. Directional variograms of a depth variable was used for the OK interpolation, and the UK interpolation applied the same method except that residuals were used as variables.

Structural maps usually show a trend in input data since mapped surfaces often have a structural trend towards one cardinal direction, or a sequence of geological structures, like in our case a field structure is part of a structural unit of a regional monocline. Therefore, with such a dataset that is characterised by a linear spatial trend or a drift, it is advisable to use UK in relation to the OK because this technique was developed precisely in order to identify and calculate trends in the data, i.e. to describe the regression and give better mapping results. That can be seen by comparing the maps that were obtained by OK and UK. Consequently, it means that UK could be used for the interpolation of depths on all Croatian hydrocarbon fields in the CPBS, where trends such as those described in this paper are observed.

References

7.References

Davis, J.C. (1973): Statistics and Data Analysis in Geology. Wiley, New York, 550 p.

Isaaks, E. and Srivastava, R. (1989): An Introduction to Applied Geostatistics, Oxford University Press Inc., New York, 561 p.

Journel, A. G. and Huijbregts, C. J. (1978): Mining Geostatistics, Academic Press, London, UK, 600 p.

Kastelec, D. and Kosmelj, K. (2002): Spatial interpolation of mean yearly precipitation using universal kriging. In: Mrvar, A. and Ferligoj, A. (eds) Developments in Statistics, FDV, Ljubljana, Slovenia, 149-162.

Kumar, V. (2007): Optimal contour mapping of groundwater levels using universal kriging - a case study, Hydrological Sciences Journal, 52 (5), 1038-1050. DOI: 10.1623/hysj.52.5.1038

Malvíe, T. (1998): Strukturni i tektonski odnosi te znacajke ugljikovodika sirega podrucja naftnoga polja Galovac-Pavljani, Magistarski rad (Structural and tectonic relations and hydrocarbons' attributes of the extended area of the oilfield Galovac - Pavljani, UnpublishedMSc Thesis) Sveuciliste u Zagrebu, Rudarsko-geolosko-naftni fakultet, 112 p. (in Croatian)

Malvie, T. (2003): Naftnogeoloski odnosi i vjerojatnost pronalaska novih zaliha ugljikovodika u Bjelovarskoj uleknini, Doktorska disertacij a (Oil-Geological Relations and Probability of Discovering New Hydrocarbon Reserves in the Bjelovar Sag, Doctoral Thesis), Sveuciliste u Zagrebu, Rudarsko-geolosko-naftni fakultet, 123 p. (in Croatian)

Malvie, T. (2004): Regional geological settings and hydrocarbon potential of Bjelovar sag (subdepression), R. Croatia, Nafta, 55, 273-288.

Malvie, T. (2011): Geological maps of Neogen sediments in the Bjelovar Subdepression (northern Croatia), Journal of Maps, 7 (1), 304 - 317.

Malvie, T., Cvetkovie, M. and Balie, D. (2008): Geomatematicki rjecnik (Geomathematical dictionary), HGD, Zagreb, 77 p. (in Croatian)

Malvíe, T. and Balie, D. (2009): Linearity and Lagrange Linear Multiplicator in the Equations of Ordinary Kriging, Nafta, 59 (1), 31-37.

Mesie Kis, I. and Malvíe, T. (2014): Zonal estimation and interpolation as simultaneous approaches in the case of small input data set (Sandrovac Field, Northern Croatia), RGN Zbornik, 29 (1), 9-16.

Novak Zelenika, K. and Malvie, T. (2011): Stochastic simulations of dependent geological variables in sandstone reservoirs of Neogene age: A case study of the Klostar Field, Sava Depression, Geologia Croatica, 64 (2), 173-183.

Pavelie, D. (2001): Tectonostratigraphic model for the North Croatian and North Bosnian sector of the Miocene Pannonian Basin System, Basin Research, 13, 359-376.

Rögl, F. (1996): Stratigraphic Correlation of the Paratethys Oligocene and Miocene, Mitteilungen Ges. Geol. Bergbaustudenten Österreich, 41, 65-73.

Rögl, F. (1998): Palaeographic Consideration for Mediterranean and Paratethys Seaways (Oligocene to Miocene), Ann. Naturhist. Mus. Wien, 99A, 279-310.

Royden, L. (1988): Late Cenozoic Tectonics of the Pannonian Basin System, In Royden, L. H. and Horvath, F., (eds.) The Pannonian Basin (Memoir 45), American Association of Petroleum Geologists, 27-48, 394 p.

Wackernagel, H. (2003): Multivariate Geostatistics: An Introduction with Applications, (3rd ed.) Heidelberg: Springer, Berlin, 388 p.

Zornjak, D. (2009): Poroznost i propusnost lezisnih stijena u Savskoj i Dravskoj depresiji, Diplomski rad (Porosity and permeability of reservoir rocks in the Sava and Drava Depression, Bachelor Thesis), Sveuciliste u Zagrebu, Rudarsko-geolosko-naftni fakultet, Studij naftnog rudarstva, 54 p. (in Croatian)