Comparison of the Polynomial Regression zonal estimation versus the Inverse Distance Weighting and Ordinary Kriging linear interpolator; case study in the Sava Depression Neogene sandstone, Croatia

1. Introduction

The fluid's reservoir modelling is closely based on understanding spatial variations of numerous geological variables like porosity, permeability, saturation, mineral and lithological content. Better knowledge leads to better and longer recovering. Such data are, often, very scarce and irregularly distributed what makes it difficult to interpret them. Also, the most of them are indirectly calculated what complicates their reliability. It puts dilemma in selection of the interpolation methods, which are practically one of the most important interpretative tools.

Here are analysed two datasets, one with 19 (reservoir "L"), and another with 25 (reservoir "K") data, what make them small and larger [1] datasets, retrospectively. Both include porosity values calculated (from logs and cores) as reservoir's averages in real wells and represents the most important available variable for spatial mapping. Regular geological mapping can be done with any number of data, however dozens of them make maps reliable. In subsurface fluid reservoir's mapping such plentiful of data is rarely available, and 10 or so are datasets that need to be handled, more or less, successful for development and recovery prediction. In such cases methods like Kriging are less useful than mathematically simpler like Inverse Distance Weighting (IDW). Also, interpolation often, in very uncertain datasets (because of statistics or data sources), can be replaced with even more basic approach on the zonal estimation. Here are used both of such approach as test for "average" reservoir in the selected area, looking for answer about the most appropriate approach when can be chosen zonal estimation (Polynomial Regression; PR), simpler interpolation (IDW) or more advanced interpolation (Ordinary Kriging; OK).

Analysed reservoirs are chosen in the Croatian part of the Pannonian Basin System (CPBS), i.e., in the western part of the Sava Depression (Figure 1), where numerous hydrocarbon fields had been discovered, and many of them are still in production. Both ae of the Lower Pontian age and mostly consist of the medium-grained sandstones.

And last, but not the least, mentioned methods are widely used in geological researching, even in the CPBS. However, they are not limited only to geology, just the opposite. For example, application of the PR can be found in the medicine [3,4], biology [5], chemistry [6], psychology [7], geography [8], meteorology [9], thermodynamics [10], mechanical engineering [11,12], astronomy and astrophysics [13], computing [14,15], economics [16,17], civil engineering [18], food technology [19], agriculture [20], ecology [21,22], petroleum industry [23,24,25,26,27,28] and other.

Furthermore, the Inverse Distance Weighting method has found applications in geography [29], meteorology [30], geology [31], electrical engineering [32], and computing [33]. The success of Ordinary Kriging is also not limited to a specific field, as it has been applied in medicine [34], geography [35], geostatistics [36], environmental science [37], and economics [38].

2. Geological Settings and Location

The PBS, macro unit of the largest order, encompasses areas in eight countries, generally in the Central Europe. That system is surrounded by Carpathian, Alps and Dinarides as the regional mountain chains (orogens). The PBS starting to be created during the Ottnangian due to convergence between Euroasian and African Plates, i.e., subduction of the Apulian Plate below the Dinarides. So, the PBS is typical back-arc basin system, covered with numerous marine and lacustric environments, among which the largest was the (part of) Paratethys [39].

The CPBS (Figure 2), at the south-west margin, includes four geological macro units of the 2nd order, namely Mura, Drava, Sava and Slavonia-Srijem Depressions. There are discovered about 40 hydrocarbon fields with numerous reservoirs, mostly in the Neogene (Upper Miocene) sandstones.

In the CPBS, two main lithological macro sequences dominate. Their differentiation has been established based on lithology and geological history [40]. The 1st, younger, includes mostly clastic sediments, with minor biochemical (like Lithotamnium limestones) or magmatic (like basalts), of the Neogene and Quaternary periods. The 2nd, older, in the bottom is significantly older, from Mesozoic or Palaeozoic eras, encompassing magmatic or sediment rocks, often metamorphosised. The 1st sequence, in the Sava Depression, is divided into six lithostratigraphic formation [41], where 3 of them, important for this researching, are presented on the Figure 3.

During the entire Neogene—Quaternary period [43] described for the CPBS two transtensional (Badenian and Pannonian—Early Pontian) and two transpressional (Sarmatian and Late Pontian—Quaternary) regional tectonic phases. Neogene and Quaternary cyclic sediments were divided in megacycles with lithostratigraphic formations and members. Late Miocene (Pannonian and Pontian) cycle includes sedimentary association of the Sava Group (the Ivanić-Grad, Kloštar-Ivanić and Široko Polje Formations) in the Sava and western Drava Depressions. Deposition lasted approx. 5.9 Ma, and it is represented by the sequences of grey coloured sandstones, siltites and marls. The maximum total thicknesses are deposited in the central part of depressions, pinching out toward the margins. These sandstones. siltites and marls are deposited in deep lake environment basin sedimentation interrupted with turbidites, with the main sources of the material from the Eastern Alps. Main depositional mechanisms during Late Pannonian and Early Pontian were deep water turbidites and sequences of hemipelagic marls were described as well, both in the marginal and the central parts of basins. Several facies associations has been described: turbidite channel fill facies association - thick-bedded sandstones and thin-bedded sandstones; turbidite overbank-levee facies association - laminated sandstones, siltstones and marls passing into sandstones; distal turbidite facies association - alternating thin sandstones with siltstones and marls; massive marls facies association - marls with rare intercalations of thin siltstone or sandstone laminae.

The Sava Depression has been filled with sediments since the Early Neogene. However, the analysed reservoir rocks ("K" and "L" in Figure 3) are of the Late Neogene (Lower Pontian) age and they are a part of the Kloštar-Ivanić Formation. Arenitic sandstones prevail at the base of the formation, and more fine-grained sandstones intercalated with marl appear more frequently toward the top, as well as in the overlying the Široko Polje Formation. Sandstones units (thick 20-150 m) are the main reservoir rocks, and intercalated grey to brown marls (thick from 30-150 m) present the main isolator rocks. These Lower Pontian ("Abichi beds") extend across the Sava Depression, comprising of sandstones and marls. Sandstones were deposited from turbidite currents in the largest thicknesses in the central part of the depression, while marls were deposited as turbidites were calm down and settled from suspension. Generally, there are developed several unique lithofacies resulted from Bouma sequence, namely (1) interlaminated marlstone and sandy siltite/shale, (2) laminated marlstone interchanged with clayey-calcitic laminas and laminas with mica flakes and fine quartz grains, (3) fine-grained silty lithic greywacke sandstones to mudstones with quartz grains lithic fragments, bounded in matrix. (4) silty marlstones with kerogen; (5) lithic arenite sandstones with quartz grains and lithic fragments; (6) fine-grained fossiliferous lithic arenite sandstones; (7) coarse-grained lithic arenite sandstone to petromictic breccia/conglomerate; (8) clast-supported petromictic conglomerate with sandy matrix. Locally can be recognised some of them, depending on age and palaeoenvironment, and analysed reservoirs belongs to lithofacies 5 and 6.

The analysed reservoirs "L" and "K" are located (Figure 4) at the margin of the western Sava Depression. They largely participate in the hydrocarbon potential of analysed zone. That part includes about 20 hydrocarbon fields still in production, but several others are depleted. The largest ones are the Stružec Field (16x10⁶ m³ Original Oil In Place, abbr. OOIP) and the Ivanić Field (7x10⁶ m³ OOIP) [43]. The entire western part covers area of about 8000 km², where surface projection of the field's areas takes about 930 km² [44]. Generally little deeper and more to the west are located sediments described as the source rocks, with fine-grained sediments mainly enriched in kerogen. They correspond to previously described lithofacies 1, 2 and sporadically 3. All are deposited in the in deeper lacustrine environment from basin sedimentation or at the distal parts of turbidites.

3. Mathematical Basics of Applied Interpolation Methods

Here are described analytical methods applied in this researching. The Polynomial Regression is used as zonal estimator, and the Inverse Distance Weighting and Ordinary Kriging as linear interpolator. The Mean Square Error is calculated as estimator of mapping numerical error.

The polynomial regression, as estimation (here also zonal) method, had been in the focus of this paper as tools that could be useful companion in application of other linear interpolators like the IDW or OK. It could be useful when interpolation methods do not provide reliable solutions to the problem of constructing representative spatial distributions, for example in cases like low number of data, clustering or numerous outliers. Moreover, polygonal methods are considered in geostatistics (beside cell-declustering, (e.g., [46], [47]) as one of the declustering methods, where polygon of influence (known as Thiessen or Voronoi polygon, e.g., [48], [49]) includes all data points that are closer to the sample compared to any other measurements. As result, spatially separated data points will have larger polygons than clustered (grouped) ones. Such polygonal declustering, characteristics for all zonal methods have been tested on two datasets presented in this researching.

3.1. The Polynomial Regression

The polynomial regression is kind of improvement of simple linear regression, in case that the relationship between the data is non-linear. In such vase the polynomial regression can better fit data, adding some polynomial terms into linear regression equation and modelling relationship between the dependent (Y) and independent (X) variables with n^th degree polynomial function. The model is called a quadratic is such polynomial is 2, a cubic for 3 (e.g., Equation 1). The degree of function order can be set on any value (sometimes this order is considered as hyperparameter as well), but selection must be done carefully because each polynomial function can be easily underfitted or overfitted (even if least square method is used for minimizing the error). There is some kind of the "rule of the thumb" that in the subsurface mapping the fitting function could not to be of the 5^th or higher order. So, right polynomial degree would need to lead to reasonable mapping structures (where the main ones are confirmed with other mapping method), but also with smaller cross-validation results (when several models are run).

The Polynomial Regression (PR) method describes relation between independent variable x and dependent y, using polynomial of the n^th order suitable for approximation of their relation [28]. The PR can be calculated using Equitation 1, e.g. [26], and is used for non-linear relations, what can help in looking for connection for more complex datasets.

(1)

where is:

y - dependent variable;

x - independent variable;

n - the n^th order of regression;

ε_i - error for i^th data;

β - values adjusted during the calculation.

The Equation (1) represents general polynomial regression term, which always includes model parameters and hyperparameter at the end. Here is marked the "error" as the correction parameter using for fitting the model during iterations. Such parameter is used to reach optimum the model parameters. In fact, any parameter that models the function shape is called hyperparameter. It could be the ratio between training and validation datasets etc. Here such parameter is called "error", however, in the applied algorithms, the more important hyperparameter is the order of the fitting function, named as the Maximum Total Order (MTO). Here it is tested with several values between 1 to 10, optimising the regression model through iterative model. The goal was to find the best value of hyperparameter, leading to the optimal map.

3.2. The Inverse Distance Weighting

The Inverse Distance Weighting (IDW) is relatively simple linear interpolation where unknown value is calculated from known data inside searching radius, and assumption that closer data stronger participate in interpolation of unknown value. The interpolation of the unknown is defined by Equation (2), e.g., [50]:

(2)

where is:

Z_IU - interpolated value;

d_i - distance to i^th location;

z_i - known value at the i^th location;

p - distance exponent.

It is obvious that interpolation is highly depended on distance exponent (p) and its selection is led by obtaining the logical visual shapes on the map as well as minimal numerical error. In mapping it is very often set on value 2 [51], especially subsurface geological mapping in the Northern Croatia [e.g., 1,42,52,53].

3.3. The Ordinary Kriging

General linear equitation is given in (3), e.g., [54], where higher value of the weighting coefficient points on known value closer to location where value is estimated.

(3)

where is:

z_k - value estimated at location „k“;

λ_i - weighting coefficient at location „i“;

z_i - measured value at location „i“.

The Ordinary Kriging (OK) equation is upgraded Simple Kriging (SK) with adding of the Lagrange factor (µ), aiming to minimizing the Kriging variance. Surely, both techniques minimizing variance, but the OK is more successful if is rightly applied. Also, all other Kriging techniques (like Universal, Disjunctive, Indicator...) has some adding "factors" with the same purpose, designed for the specific analyses. The SK is only algorithm where sum of weightings is not standardised at 1. Also, in the SK the mean is known, but global, and in the OK the local one is used. Eventually, in all Kriging techniques (except of Indicator) the normal distribution of input dataset is highly recommended for more reliable results. It is closely connected with requirements of the stationarity of the 2nd order where is implied that expectation is independent on number and locations of data, and covariance is dependent only on distances among data (variogram). Only the Indicator Kriging assumed the stationarity of the 3rd order that implied independency of expectation and variogram existence (intrinsic hypothesis). The Lagrange function is represented with Equation 4:

(4)

where is: L(x,µ) -Lagrange function;

f(x), h(x) - functions for the "x" values;

µ - Lagrange factor.

Now, the Ordinary Kriging matrix equation can be presented with Equation 5 [e.g., 55]:

(5)

where is:

γ - variogram values;

Z₁....Z_n – measured values at locations "1...n";

Z – locations where new value is estimated;

μ – Lagrange factor.

The mentioned property of the Ordinary Kriging is that technique works with local mean (not global), what can be pretty good feature for mapping geological subsurface data, at least in the CPBS [56]. It is pure deterministic methods as well as Inverse Distance and Polynomial Regression, for example. All are "numerical and subjective", and sure, the stationarity (on some order) can be criteria for choosing the right Kriging technique if researcher knows how to observe stationarity.

3.4. The Cross-Validation

The cross-validation (CV) is set of techniques applied for estimation error calculation. Here is selected the mean square error (MSE) technique as one of the most applied in comparison of several interpolations for the same dataset. As measure of residuals variance, it usually means that the lower MSE would lead to conclusion about better interpolation. There are also some other measures, like mean absolute error (MAE, the average of the absolute difference between measured and estimated values, i.e., residuals average), root mean square error (RMSE, keeping the same unit as input dataset, allowing easier interpretation) etc.

The MSE is iterative procedure where one, randomly selected, measured values is ignored and at the same location its value is again estimated (with selected interpolation method) from the rest of measurements. The difference between those two values is error in that point. The procedure is repealed for all other measured values and eventually the MSE is calculated using Equation 6, e.g., [57]:

(6)

where is:

MSE - mean square error (of the estimation);

N - number of data;

Y_mjereni,i - measured value for loaction "i";

Y_{procijenjeni,i} - estimated value for location "i".

However, the interpretation is not so straightforward, because the maps also can include different isolines shapes, where some forms (structures) could be almost impossible and could eliminate interpolation method although the accompanied MSE is lower than in other methods.

4. Discussion and Results Overview

This researching has been focused on zonal estimation as useful addon to classical interpolation methods in the hydrocarbon sandstone reservoirs. The goal has been set up for the CPBS, where numerous interpolations of mentioned and similar reservoirs had been published. The PR is selected as zonal estimation. The IDW and OK are used interpolations. The datasets are divided into small (19 data, Table 1) and large (25 data, Table 2) although both are of similar size.

Both sets are limited in their representation, but also solely publicly available for analysed zones. Given values are mostly derived indirectly from e-logs, and it is why they are same in the numerous wells. In reality, they are very similar because logged in the same lithofacies. Values represent mean porosity along entire reservoir in the corresponding well, what is often used in numerous mappings. However, some practitioners used porosities pondered with another variables, for example, thickness. Such option is not correct to apply if the researcher does not know (and described) depositional model, i.e., how clastics volume and sizes depends on (well) position in such environment and subsequent compaction. As we did not have such detail information, the pure mean porosity has been mapped.

Listed data had been interpolated or estimated using IDW, OK and PR. Table 3 presents the cross-validation values obtained by changing the value of the maximum total order (MTO). The selected MTO (as previously mentioned as one hyperparameter related on the order of a fitting function) values are 1, 2, 3, 4, 5, and 10. They determine the maximum sum of the exponents of the independent (X) and (dependent) Y variables, i.e., the degree of polynomial regression, ensuring that the sum of the exponents does not exceed the MTO value [58]. For example, for "n" data, the maximum degree of polynomial function and MTO is "n-1". Like it was mentioned for order of the polynomial function, choosing an excessively high MTO value may lead to over-fitting which causes situations where the data or predictions have significant errors or variations, which makes them unreliable [59]. In contrast, a low value of MTO may result in under-fitting of the model leading to failure in representing the relationship between the variables [60]. For all cases the MSE is calculated and maps generated. Regarding the PR there were several options to estimate data changing the option named as "Max Total Order (MTO)". In this analysis such option is varied with values 1, 2, 3, 4, 5 and 10, looking the lowest MSE (Table 3).

When MTO reached 2, the MSE staying constant (and overall, the variations are too small that it is matter). So, this parameter can be considered of low importance, when sill of 2 is reached, for used datasets regarding their abundance and locations. The next step was calculation of the same error using the IDW and OK interpolations in the same reservoirs (Table 4).

In both reservoirs the MSE values are considerably lower for interpolations than for estimator (Table 3 and Table 4). In the reservoir "K" the MSE values are 0.000119 (IDW) and 0.05401 (PR). In the reservoir "L" the MSE are 0.000676 (OK) and 0.040117 (PR). In the reservoir with more data those values are absolutely lower, what indicated on statistically representative datasets, i.e., fact that increasing of values did not include new "outlier" values and confirmed that both datasets are representative. Also, the lowest MSE for OK proved that 25 data were enough abundant for calculation of reliable spatial model (variogram) based (probably, because it is not tested) normal distribution of porosity (as theoretical assumption for sandstone reservoirs).

The maps made for both datasets using the PR zonal estimator are shown on Figure 5a (reservoir "K") and Figure 5b (reservoir "L"). The main difference is, like as in the MSE, in map associate with MTO=1 and maps with MTO>1. As zonal, the maps show sharp borders between different porosity zones and all follow the NW-SE/WNW-ESE strike, what corresponds with the strike of sandstone depositional environment. For reservoir "K" (Figure 5a) the zone values (colours) are almost the same for the MTO=2, 3 and 10, what follows the stabilisation of the MSE on the MTO=2 and higher. For the reservoir "L" almost the same statement is valid, also with additional features reflected in the "parabolic" borders among zones. Also, the PR is sensitive on higher differences between close data, what is visible in well pairs J-168 and J-25, i.e., L-27 and L-131a. In any case when the data are classified mostly in several groups of the same values (Table 1 and Table 2) the zonal estimator will group them pretty clear and simple.

The interpolated maps are given on Figure 6a (reservoir "K", the IDW) and Figure 6b (reservoir "L", the OK). As expected, the interpolators are based on transitional zones between isoporosity lines on the selected equidistance. So, the spatial interpretation and prediction is much easier than in the zonal representation, however, can be tricky when numerous points are of the same values like in the presented reservoirs (Table 1 and Table 2). The main difference is existing of the maximums and minimums, which on the zonal maps were amalgamated into single zone. Back to MSE it would be more appropriate spatial solution than zones.

Application of the OK asked for definition of spatial model, what was not easy task due to relatively low amount of data. It is why the omnidirectional variogam model had been used (for directional there would not be enough data per sectors of the searching ellipsoid). Some theoretical variogram models (Figure 7) had been tested like linear (range about 900 m) and spherical (range about 350 m) with relatively small changes in visual results and MSE. So, we prefer to use models with lower ranges. Also, the Kriged map used larger equidistance (0.02) than the IDW (0.01) just to be sure that different spatial model will not significantly change the shapes of the main "structures". It was not influenced the MSE calculation.

5. Conclusions

The main achievements can be summarised as follows:

1. Zonal estimators and interpolators are two different approaches for subsurface sandstone mapping, both with their own properties and advantages.

2. Larger datasets will favour interpolators, smaller estimators.

3. Here is proven that both datasets can be interpolated with the IDW (19 data, "K", 0.000119) and even using the OK (25 data, "L", 0.000676), and in both cases the MSE will be significantly lower than in the estimation with the PR (0.05401 in "K" and 0.040117 in "L").

4. However, both datasets are very specific meaning that many values are the same, what practically decreased number of different data from 19 or 25 to only several classes.

5. That is result of way how data are calculated, where they are indirectly averaged from mixture of well logs and cores in the sandstone reservoir intervals across the field.

6. In such case pure interpolation will not be so useful, because most transitional isoporosity lines are not confirmed with measured value, i.e., hard data, but are artificially calculated from interpolation algorithm between hard points (e.g., Figure 6aas area among J-25, J-168, J- 173 wells).

7. Alternatively, it could be improved with simultaneously running zonal estimator algorithm, where mentioned value's classes could be easily recognised inside zones, revealing their area and strike (Figure 5a, b).

8. Zonal maps can greatly improve development of presented sandstone reservoirs or their counterparts in other regions of the CPBS if depositional model for mapped structure is described in details.