Application of the Modified Shepard’s Method (MSM): A Case Study with the Interpolation of Neogene Reservoir Variables in Northern Croatia

Interpolation is a procedure that depends on the spatial and/or statistical properties of the analysed variable(s). It is a particularly challenging task for small datasets, such as in those with less than 20 points of data. This problem is common in subsurface geological mapping, i.e., in cases where the data is taken solely from wells. Successful solutions of such mapping problems depend on interpolation methods designed primarily for small datasets and the datasets themselves. Here, we compare two methods, Inverse Distance Weighting and the Modified Shepard’s Method, and apply them to three variables (porosity, permeability, and thickness) measured in the Neogene sandstone hydrocarbon reservoirs (northern Croatia). The results show that cross-validation itself will not provide appropriate map selection, but, in combination with geometrical features, it can help experts eliminate the solutions with low-probable structures/shapes. The Golden Software licensed program Surfer 15 was used for the interpolations in this study.


Introduction
Geostatistics is a branch of applied statistics that applies theories of deterministic estimation and stochastic processes. These process are described as statistical inferences to different spatial phenomena in numerous branches of science. The most important geostatistical term applied in geology and petroleum geology is the regionalized variable [1], which indicates a measured or estimated magnitude whose behaviour is described in part by a random variable and in part by a deterministic variable. Such variables each show the properties of spatial continuity and, as such, are modelled in space. A qualitative description of the properties of a regionalized variable leads to a more accurate assessment of the regionalized variable.
The interpolation methods used in the Croatian part of the Pannonian Basin System (CPBS) are described in [2], while the final specific interpolation method that we applied is presented in [3]. The IDW method is widely used in the CPBS for mapping geological variables, and its application is described in [4][5][6][7][8].
The modified Sheppard's method (MSM) has been applied all over the world. The following papers were selected as appropriate examples of this method. In [9], the authors applied the MSM to a digital terrain model region of the Chen-Yu-Lan River (Taiwan). The authors in [10] applied the MSM to multicomponent aerosol-cloud parameterization (global climate modelling), and the Stats 2020, 3 FOR PEER REVIEW 3 of 16  A typical geological section is derived from the well data, and clearly outlines the sandstone member (blue borders) with discovered hydrocarbon reservoirs (see Figure 3). The section also reveals a partial structural inversion within the formation, where monocline of deeper members is transformed into the syncline in the upper part. Structurally, in most cases, tectonic blocks are considered to represent separate structural traps, where hydrocarbons are accumulated in the highest parts. Therefore, these tectonic blocks represent separate hydrodynamic units with similar reservoir characteristics. A structural map of the analysed Late Miocene (Lower Pontian) reservoir "K" is shown in Figure 2.  A typical geological section is derived from the well data, and clearly outlines the sandstone member (blue borders) with discovered hydrocarbon reservoirs (see Figure 3). The section also reveals a partial structural inversion within the formation, where monocline of deeper members is transformed into the syncline in the upper part. A typical geological section is derived from the well data, and clearly outlines the sandstone member (blue borders) with discovered hydrocarbon reservoirs (see Figure 3). The section also reveals a partial structural inversion within the formation, where monocline of deeper members is transformed into the syncline in the upper part. The mapped sandstone "K" could be described as one of two typical sandstone types in the northern Croatian Neogene. The younger sandstones, well-sorted, medium-grained, quartz-mica dominated, were deposited in the lacustric, mostly brackish environment of Pannonian Lake. The detritus was transported by sequential turbidites during the Late Miocene (Upper Pannonian and, especially, Lower Pontian) (e.g., [19]). The dominant, primary porosity is intra-granular. The second, older type comprises poorly sorted and coarse-grained sandstones, containing some larger (ruditic) clasts, which are bound together with sparry calcite cement. They contain various lithoclasts derived from older sedimentary and metamorphic rocks: dolomites, quartzites, gneisses and various schists, as well as individual quartz grains which originated from metamorphic rocks. Such calcareous sandstones (e.g., [20,21]) were deposited during the Middle Badenian, in the shallow marine backreef environments of the Paratethys Sea, with a significant input of siliciclastic material from the nearby land (see Figure 4). The dominant, primary porosity is intra-skeletal (see Figure 5).  The mapped sandstone "K" could be described as one of two typical sandstone types in the northern Croatian Neogene. The younger sandstones, well-sorted, medium-grained, quartz-mica dominated, were deposited in the lacustric, mostly brackish environment of Pannonian Lake. The detritus was transported by sequential turbidites during the Late Miocene (Upper Pannonian and, especially, Lower Pontian) (e.g., [19]). The dominant, primary porosity is intra-granular. The second, older type comprises poorly sorted and coarse-grained sandstones, containing some larger (ruditic) clasts, which are bound together with sparry calcite cement. They contain various lithoclasts derived from older sedimentary and metamorphic rocks: dolomites, quartzites, gneisses and various schists, as well as individual quartz grains which originated from metamorphic rocks. Such calcareous sandstones (e.g., [20,21]) were deposited during the Middle Badenian, in the shallow marine back-reef environments of the Paratethys Sea, with a significant input of siliciclastic material from the nearby land (see Figure 4). The dominant, primary porosity is intra-skeletal (see Figure 5).  The mapped sandstone "K" could be described as one of two typical sandstone types in the northern Croatian Neogene. The younger sandstones, well-sorted, medium-grained, quartz-mica dominated, were deposited in the lacustric, mostly brackish environment of Pannonian Lake. The detritus was transported by sequential turbidites during the Late Miocene (Upper Pannonian and, especially, Lower Pontian) (e.g., [19]). The dominant, primary porosity is intra-granular. The second, older type comprises poorly sorted and coarse-grained sandstones, containing some larger (ruditic) clasts, which are bound together with sparry calcite cement. They contain various lithoclasts derived from older sedimentary and metamorphic rocks: dolomites, quartzites, gneisses and various schists, as well as individual quartz grains which originated from metamorphic rocks. Such calcareous sandstones (e.g., [20,21]) were deposited during the Middle Badenian, in the shallow marine backreef environments of the Paratethys Sea, with a significant input of siliciclastic material from the nearby land (see Figure 4). The dominant, primary porosity is intra-skeletal (see Figure 5).   The object of this research, the "K" reservoir, currently produces hydrocarbons using a secondary method (water injection). This is a very maturely explored and developed area, which allows significant improvements and upgrades of geological models in order to increase hydrocarbon recovery. This is especially reflected in the application of statistical and geomathematical analyses of numerous reservoir variables such as porosity and permeability (primary variables) and injected water volumes (secondary variable). Furthermore, most of the existing reservoirs maps (structural, porosity, etc.) should be redrawn in much more detail, using ever-advancing mapping algorithms.
Mathematical algorithms, such as variograms and Kriging, have been the main tools in analysing spatial dependence and the possible anisotropy of the selected variables in the CPBS. However, for a small input dataset, like the injected volume of water (three points of data) and permeability (18 points of data), simpler interpolation methods such as: Inverse Distance, Nearest and Natural Neighbourhood are more effective, e.g., 15 points of data were sufficient to apply the Inverse Distance and Nearest Neighbourhood methods, and even for reservoir characterisation in cases of very low numbers of datasets (e.g., [22]). However, the problem of visually expressed local values ("bull'seye," "butterflies") permanently exists and can lead to misinterpretation. That is the reason for testing the modified Sheppard's methods in this study, as a replacement for using the IDW method, leading to similar cross-validation results and, visually, a lower number of non-acceptable geological shapes.
The interpolation of porosity, permeability and thickness of the Neogene sandstone reservoir "K" in northern Croatia, is valid for all similar reservoirs in the CPBS. This also means that the applied methods, interpretations and conclusions could be applied in the Middle Miocene (Badenian) as well as in the Upper Miocene (Pannonian/Pontian) sandstones in this area, even though their origin is partially different.

Basics of the Applied Interpolation Methods
Interpolation is an estimation of the values of variables in places where no values have been measured. When interpolating, one (primary) or multiple (primary + secondary) variables can be used. The interpolation methods applied in reservoir "K" are Inverse Distance Weighting (IDW) and the Modified Shepard's Method (MSM). Only primary variables (porosity, permeability, thickness) were used.
This was a bi-dimensional interpolation, standardised for hydrocarbon reservoirs. "Standardisation" means that each point of data represents one well and an interval of a well drilled The object of this research, the "K" reservoir, currently produces hydrocarbons using a secondary method (water injection). This is a very maturely explored and developed area, which allows significant improvements and upgrades of geological models in order to increase hydrocarbon recovery. This is especially reflected in the application of statistical and geomathematical analyses of numerous reservoir variables such as porosity and permeability (primary variables) and injected water volumes (secondary variable). Furthermore, most of the existing reservoirs maps (structural, porosity, etc.) should be redrawn in much more detail, using ever-advancing mapping algorithms.
Mathematical algorithms, such as variograms and Kriging, have been the main tools in analysing spatial dependence and the possible anisotropy of the selected variables in the CPBS. However, for a small input dataset, like the injected volume of water (three points of data) and permeability (18 points of data), simpler interpolation methods such as: Inverse Distance, Nearest and Natural Neighbourhood are more effective, e.g., 15 points of data were sufficient to apply the Inverse Distance and Nearest Neighbourhood methods, and even for reservoir characterisation in cases of very low numbers of datasets (e.g., [22]). However, the problem of visually expressed local values ("bull's-eye", "butterflies") permanently exists and can lead to misinterpretation. That is the reason for testing the modified Sheppard's methods in this study, as a replacement for using the IDW method, leading to similar cross-validation results and, visually, a lower number of non-acceptable geological shapes.
The interpolation of porosity, permeability and thickness of the Neogene sandstone reservoir "K" in northern Croatia, is valid for all similar reservoirs in the CPBS. This also means that the applied methods, interpretations and conclusions could be applied in the Middle Miocene (Badenian) as well as in the Upper Miocene (Pannonian/Pontian) sandstones in this area, even though their origin is partially different.

Basics of the Applied Interpolation Methods
Interpolation is an estimation of the values of variables in places where no values have been measured. When interpolating, one (primary) or multiple (primary + secondary) variables can be used. The interpolation methods applied in reservoir "K" are Inverse Distance Weighting (IDW) and the Modified Shepard's Method (MSM). Only primary variables (porosity, permeability, thickness) were used.
This was a bi-dimensional interpolation, standardised for hydrocarbon reservoirs. "Standardisation" means that each point of data represents one well and an interval of a well drilled through a selected reservoir. The selected value is an average across the reservoir interval, and Figure 3 also schematically depicts average depths, lithology and structural position, e.g., in [23].
The average values are calculated from logs or cores, which were not a part of this study, i.e., such average values were taken from available reports. 2D reservoir interpolation is the conventional procedure for reservoirs of low to medium thickness (in general, less than 50 m) and which are of mostly homogenous composition (the majority of sandstone reservoirs are often approximated as homogeneous, much like the one being analysed). Examples of different reservoir interpolations in northern Croatia have been vastly made in 2D space and have been published in numerous papers. Here, selections were made in geological depositional environments similar to the Sava Depression (as in the situation here) and the Drava Depression (just to the north), which are the two largest depressions in the CPBS. The choice to apply 2D interpolation (vs. the 3D possibility) had been made without regard to method selection (hand-made, software with statistical or geostatistical algorithms. All case reservoirs of lithostratigraphic units had been considered as vertically uniform sections, regardless of the selected variable. Such regional and local structural interpolation can be clearly observed on the map published for the Bjelovar Subdepression (hand-made interpolation given in [24]) or the Kloštar Structure in the Sava Depression [25]. In the same structure, even geostatistical sequential indicator simulation of porosity and thickness had been made with 2D interpolation [26]. Even back in the Bjelovar Subdepression, when the thickness had been remapped by Kriging, it was done in 2D space (e.g., [27,28]). Finally, 2D interpolation in the same region as selected in this analysis, had been especially favourable when the input dataset was small, regardless of whether the interpolation was done by statistical or geostatistical methods [3,19].
As previously mentioned, 3D interpolation is a useful approach when the mapped unit is heterogeneous, and can be precisely defined from the available data (like logs or cores). If such data are ambiguous, then it is much more appropriate to apply 2D interpolation and avoid any heterogeneity which cannot be accurately interpreted. Moreover, in a large number of cases, the mapped variable (thickness, porosity, permeability) does not depend on indirect variables such as topography or well length, because such variables are collected in a single unit with the same depositional history, which is a critical element. Even if the depths are mapped, the "problem" of topography is often easily solved using a datum plane (in the entire CPBS it is +100 m). All these facts outline why 2D mapping is a valid procedure for this analysis, and not 3D.

Inverse Distance Weighting (IDW)
The Inverse Distance Weighting method is a mathematically simple interpolation method, where the unknown value of a variable is estimated according to the known values inside the searching radius (radii). The radii of such an ellipsoid (see Figure 6) or the radius of a circle can include some data points, but generally, most software includes all points that are geologically (or in some other way) connected by default.
In any case, a searching area could be a radius (a searching circle) or radii (a searching ellipsoid), which covers both options available in 2D interpolation. We did not discuss a 3D option, because it would produce a spheroid, which was not used here, and software mostly approximates a spheroid with two perpendicular ellipsoids (thus software does not construct a real spheroid, but only three perpendicular axes).
Generally, each point is weighted by distance (see Figure 6). A single datum is inversely proportional, dependent on the distance between the unknown and the measured values of a variable. Stats 2020, 3 FOR PEER REVIEW 7 of 16 Figure 6. Points included within the IDW interpolation area (from [3]).

Modified Shepard's Method (MSM)
MSM interpolation is a modification of the IDW method, with the aim of reducing the expressive local values (outliers, extremes) that could cause "bull-eyeing" or "butterfly shape" effects. The method was developed by [31], sometimes named as Shepard's method, and modification of the method was carried out in the works of, for example, [32][33][34]. An estimation (e.g., [33,34]) using the MSM function is shown by Equation (2): where: F-MSM function; W-relative weights; Qk-bivariate quadratic function; x, y-data coordinates; n-number of data.
A bivariate quadratic function is a second-degree polynomial of the form ( , ) = 2 + 2 + + + + where x and y are variables, a,b,c,d,e variables, f is constant. Such a function describes a quadratic surface. Setting f(x,y)=0 describes the intersection of the surface with the plane z=0, which is a locus of points equivalent to a conic section.

Modified Shepard's Method (MSM)
MSM interpolation is a modification of the IDW method, with the aim of reducing the expressive local values (outliers, extremes) that could cause "bull-eyeing" or "butterfly shape" effects. The method was developed by [31], sometimes named as Shepard's method, and modification of the method was carried out in the works of, for example, [32][33][34]. An estimation (e.g., [33,34]) using the MSM function is shown by Equation (2): where: F-MSM function; W-relative weights; Q k -bivariate quadratic function; x, y-data coordinates; n-number of data.
A bivariate quadratic function is a second-degree polynomial of the form f (x, y) = ax 2 + by 2 + cxy + dx + ey + f where x and y are variables, a,b,c,d,e variables, f is constant. Such a function describes a quadratic surface. Setting f(x,y) = 0 describes the intersection of the surface with the plane z = 0, which is a locus of points equivalent to a conic section.
The MSM uses so-called relative weights determined (e.g., [33]) with Equation (3): Stats 2020, 3 75 where: W-relative weights; d k -Euclidean distance between points at locations (x, y) and (x k , y k ); R ω -radius of influence around node (x k , y k ). Figure 7 shows the searching radii of the MSM within the measured and estimated values (grid). Here it can be noted that each unknown point is surrounded with local searching radii, instead of global, like in the IDW method.
where: W-relative weights; dk-Euclidean distance between points at locations (x, y) and (xk, yk); Rω-radius of influence around node (xk, yk). Figure 7 shows the searching radii of the MSM within the measured and estimated values (grid). Here it can be noted that each unknown point is surrounded with local searching radii, instead of global, like in the IDW method.

Cross-Validation (CV)
The algorithm or the cross-validation is rather simple (see Equation 4), where some data (or a cluster of data) has been "deleted" and the estimation is performed from the remaining data. This procedure is often called "Mean Square Error" (MSE), due to the squaring data differences. Such a re-sampling method is characterized by dimensionless standard error. The famous pioneer work on cross-validation was published by [35] and one classical application of MSE in the CPBS has been published in [6].
In geological mapping, the MSE has been applied to get a value that shows a numerical error obtained from the same dataset interpolated by a different method (algorithm). In our case, if the porosity, permeability or thickness are mapped while simultaneously using IDW and the MSM, each map, although constructed from the same data, will have a different visual outcome (shape) and a different MSE value. In general, a lower MSE means better interpolation, but a numerical value cannot be the only criteria for quality evaluation. The calculation of MSE can be applied to all geological variables.

Cross-Validation (CV)
The algorithm or the cross-validation is rather simple (see Equation (4)), where some data (or a cluster of data) has been "deleted" and the estimation is performed from the remaining data. This procedure is often called "Mean Square Error" (MSE), due to the squaring data differences. Such a re-sampling method is characterized by dimensionless standard error. The famous pioneer work on cross-validation was published by [35] and one classical application of MSE in the CPBS has been published in [6].
In geological mapping, the MSE has been applied to get a value that shows a numerical error obtained from the same dataset interpolated by a different method (algorithm). In our case, if the porosity, permeability or thickness are mapped while simultaneously using IDW and the MSM, each map, although constructed from the same data, will have a different visual outcome (shape) and a different MSE value. In general, a lower MSE means better interpolation, but a numerical value cannot be the only criteria for quality evaluation. The calculation of MSE can be applied to all geological variables.
where: MSE-value of cross-validation, i.e., mean square error; measured-value measured d at location "I"; estimated-value estimated at location "I"; n-number of locations. The cross-validation applied in this analysis has been included as one of the software (Golden/Surfer) features. According to the given definition, all points were deleted one-by-one, i.e., sequentially. This means that in the "n-th" iteration, only one point was deleted, as the equation condition allows (see Equation (4)).

Interpolation Results in Reservoir "K" for Variables of Porosity, Thickness and Permeability
Here are three interpolated geological variables, measured in reservoir "K", namely: porosity, permeability and thickness. Tables 1 and 2 shows the part of their descriptive statistics, calculated from values measured from the laboratory measurement and/or interpretation of the logging curves (resistivity, density and neutron). According to Table 1, the analysed sets belong to a small geological input data set [8,36]. So, interpolation has been done with mathematically simpler methods. The porosity maps obtained by the IDW and MSM methods are shown in Figure 8. The value cross-validation (porosity) for IDW is 0.00119 and MSM is 0.00345. The other crucial parameters for interpolation presented in Figures 8-10 are: power exponent 2, grid size 50 × 50 cells, searching radii 1/2 of primary and secondary reservoir axes (the boxed/interpolated area in Figure 2). According to Table 1, the analysed sets belong to a small geological input data set [8,36]. So, interpolation has been done with mathematically simpler methods. The porosity maps obtained by the IDW and MSM methods are shown in Figure 8. The value cross-validation (porosity) for IDW is 0.00119 and MSM is 0.00345. The other crucial parameters for interpolation presented in Figures 8, 9 and 10 are: power exponent 2, grid size 50 × 50 cells, searching radii 1/2 of primary and secondary reservoir axes (the boxed/interpolated area in Figure 2). Permeability maps obtained by the IDW and MSM methods are shown in Figure 9. The accompanied cross-validation (permeability) for IDW is 480.8 and MSM 516.1. Thickness maps which have been interpolated by the IDW and MSM methods are shown in Figure 10. The value cross-validation (thickness) for IDW is 40.7 and MSM is 60.5. According to Table 1, the analysed sets belong to a small geological input data set [8,36]. So, interpolation has been done with mathematically simpler methods. The porosity maps obtained by the IDW and MSM methods are shown in Figure 8. The value cross-validation (porosity) for IDW is 0.00119 and MSM is 0.00345. The other crucial parameters for interpolation presented in Figures 8, 9 and 10 are: power exponent 2, grid size 50 × 50 cells, searching radii 1/2 of primary and secondary reservoir axes (the boxed/interpolated area in Figure 2). Permeability maps obtained by the IDW and MSM methods are shown in Figure 9. The accompanied cross-validation (permeability) for IDW is 480.8 and MSM 516.1. Thickness maps which have been interpolated by the IDW and MSM methods are shown in Figure 10. The value cross-validation (thickness) for IDW is 40.7 and MSM is 60.5. Permeability maps obtained by the IDW and MSM methods are shown in Figure 9. The accompanied cross-validation (permeability) for IDW is 480.8 and MSM 516.1.
Thickness maps which have been interpolated by the IDW and MSM methods are shown in Figure 10. The value cross-validation (thickness) for IDW is 40.7 and MSM is 60.5.
The maps obtained through the IDW and MSM methods are evaluated in a quick-look, according to the most observable feature of a highly expressed local value, i.e., bull-eye or butterfly shape effects. The MSM is clearly advantageous because it mitigates the mentioned effects, however a numerical review using the cross-validation values strongly favours the IDW method. Stats 2020, 3 FOR PEER REVIEW 11 of 16 The maps obtained through the IDW and MSM methods are evaluated in a quick-look, according to the most observable feature of a highly expressed local value, i.e., bull-eye or butterfly shape effects. The MSM is clearly advantageous because it mitigates the mentioned effects, however a numerical review using the cross-validation values strongly favours the IDW method.

Discussion
To consider the differences between the results obtained with the IDW and MSM methods, it is necessary to discuss the main differences in their equations. The IDW method does not use a weighting coefficient, i.e., each value is "weighted" by a simple (powered) inversely proportional distance from the measured point. The general form of the IDW method is Equation (5): where: u(x)-interpolated value in the location "x"; ui-i-th known value inside searching radius(es) centred in the location "x"; i-no. of samples (known/measured value); wi-weight for the "i-th" sample with known value; N-number of samples with known value.
So, the weighting coefficients for the location "x" could generally be defined as Equation (6): where: wi-weight for location "x"; d(x,xi)-is a given distance (metric operator) from the known point xi to the unknown point x, p-power factor (positive real number).
This is a basic IDW function as defined by Shepard ([31], where weights decrease as distances to the known values increase. On the other hand, higher values of power ("p") favour a larger influence of measured points closer to the interpolated point, and eventually ended up in a zonal interpolation (like a Voronoi diagram), where the inside zone is used as a constant value of a single known point centred into such a polygon and no isolines are constructed. On the contrary, for a 2D interpolation, if the power is 2 (like we used here and based on the previous interpolation done in the CPBS) or less,

Discussion
To consider the differences between the results obtained with the IDW and MSM methods, it is necessary to discuss the main differences in their equations. The IDW method does not use a weighting coefficient, i.e., each value is "weighted" by a simple (powered) inversely proportional distance from the measured point. The general form of the IDW method is Equation (5): where: u(x)-interpolated value in the location "x"; u i -i-th known value inside searching radius(es) centred in the location "x"; i-no. of samples (known/measured value); w i -weight for the "i-th" sample with known value; N-number of samples with known value. So, the weighting coefficients for the location "x" could generally be defined as Equation (6): where: w i -weight for location "x"; d(x,x i )-is a given distance (metric operator) from the known point x i to the unknown point x, p-power factor (positive real number). This is a basic IDW function as defined by Shepard ([31], where weights decrease as distances to the known values increase. On the other hand, higher values of power ("p") favour a larger influence of measured points closer to the interpolated point, and eventually ended up in a zonal interpolation (like a Voronoi diagram), where the inside zone is used as a constant value of a single known point centred into such a polygon and no isolines are constructed. On the contrary, for a 2D interpolation, if the power is 2 (like we used here and based on the previous interpolation done in the CPBS) or less, more distanced measured points will have a greater influence. In the "M" dimensions, this statement is valid for p ≤ M. In all cases, IDW minimizes the function that measures deviation between the unknown and known points using the condition Equation (7): On the contrary, the MSM is a more advanced technique than IDW, including a weighting coefficient (see Equation (8)), and calculating the interpolated value using only the nearest known points within the R-sphere (instead of using a full set of measured samples, which is the default setting of IDW in numerous algorithms). Weights are modified through Equation (8): where: Rd(x,x k )-is local searching radius(es), so called "R sphere", centred in point "x" and including "k" measured points; W k (x)-weighting coefficient for known point "k", at distance "x" from interpolated point. By default, IDW takes into account all measured points or work with a general searching radius (or radii for an ellipsoid). However, the MSM, by default, works with local searching (both can be changed in most software, but those are theoretical defaults). Therefore, the "R sphere" modification would need to reduce anomalous local values, i.e., outliers or extremes. The maps in Figures 7-9 show this to be true.
However, the problem of calculated cross-validation values remains. They were significantly lower for IDW than the MSM (see Table 3). The differences are not negligible. For porosity, the MSM cross-validation is 289% larger, permeability 7%, and thickness 49%. Table 3. Summary results of application of the IDW and MSM methods in reservoir "K". So, let's analyse the Equations (2) (the MSM function), (3) (weighting coefficient in the MSM searching radius) or (8) (the more general form of Equation (3)) one more time. It is not surprising that the MSM will mostly give larger cross-validation values, because it works by default with local search and, consequently, includes a smaller number of known values in each cross-validation iteration. If the "one point out" method is applied with a small number of remaining known values, the difference between the real and the estimated value in such a point will be mostly higher than if all the known measurements were used. So, generally, the MSM will tend to give a larger cross-validation value for the same dataset, in comparison with IDW.

Cross-Validation
Secondly, the MSM approximates the unknown value with a quadratic function, i.e., a polynomial of degree 2, for one or more variables in which the 2nd degree is the highest one. Any quadratic polynomial with two variables may be given as Equation (9): f (x, y) = ax 2 + by 2 + cxy + dx + ey + f where: x,y-variables; a, b, c, d, e, f -coefficients.
Consequently, the MSM will force gradual interpolation much more strongly for outliers (extremes), i.e., it will perform smoothing instead of bull-eyeing. Finally, let us look at the selection criteria of the interpolation algorithm for the small datasets in the same spatial area as selected for this work (Table 4, the Sava Depression, northern Croatia [3]): Table 4. Recommended interpolation methods for small input data set [3]. So, the discussed 2D interpolation, i.e., application of the IDW/MSM methods, would obviously be an appropriate choice for a quick evaluation of reservoirs. However, it must be pointed out that quick insights in such cases are mostly the first/early step. Later, in medium or large sized reservoirs, the development phase leads to additional data collection. Consequently, newly acquired and validated data could result in a larger number of measured points, laterally and vertically. Such datasets would sometimes be numerous enough for 3D subsurface modelling of different variables (which most often includes thickness, porosity and permeability).

Applicability of Interpolation Method
Nowadays, such 3D interpolation can be easily carried out by different software, especially if 2D models (like the ones presented) offered promising results. Advanced 3D models are currently used for many applications, especially in the mining of mineral ores (e.g., the model in [37] where the authors also applied the IDW interpolation in a complex grid/mesh based on the octahedron/quasi-pentahedron/pentahedral method). In any case, geostatistics offer numerous algorithms and one of them is Block Kriging, a classic 3D estimation method.

Conclusions
The authors in [3] divided small datasets into three classes regarding included points: (a) 1-5, (b) 6-10 and (c) 11-19 points. It is concluded that if the main selection criteria could be cross-validation, then IDW would be the best choice for interpolation. However, they also warned that this method sometimes created numerous "bull's-eyes" or "half-butterfly" features. So, we can summarize our results in conclusion and offer recommendations in the following statements and Table 5: 1.
The MSM could be recommended for subsurface geological mapping of Neogene deposits in northern Croatia; 2.
This method is valid for datasets smaller than 20 measured values, i.e., for the early exploration phase of hydrocarbon reservoirs or the later development phase when the number of measurements of a selected property is small, but quick insight in the spatial distribution of such variables is necessary; 3.
The selected variables could be porosity, permeability and thickness, measured in Neogene lithostratigraphic units by laboratory analysed well data (cores) or interpreted logs; 4.
The interpreters need to be aware of: (a) the IDW method resulted in the lowest cross-validation of all applied methods without a spatial model (variogram), i.e., namely the Nearest neighbourhood and the Natural neighbour, (b) the MSM interpolation eliminated all the unwanted geological features and did not result in a cross-validation 250 % higher than in the IDW, for the same variable.

5.
Based on the visual (geometric) criteria, the sets with less than five measured points are too small for interpolation (using the IDW or MSM methods), and they are only appropriate for zonal estimation. Table 5. Recommended choice between the IDW and MSM methods for a small input data set (less than 20 points of data.

Criteria
(1) If the Nearest Neighbourhood, the Natural Neighbour or any available zonal method is tested and the IDW has the smallest cross-validation. (2) If other non-zonal methods are used, but cross-validation is a single criterion for a quick-look.
(1) If IDW has a smaller cross-validation than the MSM, but the outliers are strongly and unnaturally emphasised on the map ("bull-eye", "butterfly" features). (2) If the MSM is chosen, an expert geologist, with knowledge about the regional geology of the mapped volume must evaluate the map(s).

5-20 Yes Yes
The Golden Software licensed program Surfer 15 was used for the presented interpolations. This software, among others, supports the IDW and MSM algorithms, the calculation of cross-validation and defining of the searching radius (radii).