Kriging with a Small Number of Data Points Supported by Jack-Knifing , a Case Study in the Sava Depression ( Northern Croatia )

The semivariogram and the ordinary kriging analyses of porosity data from the Sava Depression (Northern Croatia), are presented relative to the Croatian part of the Pannonian Basin system. The data are taken from hydrocarbon reservoirs of the Lower Pontian (Upper Miocene) age, which belong to the Kloštar Ivanić Formation. The original datasets had been jack-knifed with the purpose of re-sampling and calculating the more reliable semivariograms. The results showed that such improvements can assist in the interpolation of more reliable maps. Both sets, made by the original and re-sampled data, need to be compared using geological recognition of isoline’s shapes (such as “bull-eye” or “butterfly” effects) as well as cross-validation results. This comparison made it possible to select the most appropriate porosity interpolation.


Introduction
The hydrocarbon reservoirs in the Sava Depression (Northern Croatia) are in the secondary or tertiary recovery phase.Most are supported by water injection.Two fields from the abovementioned depression with still significant production of hydrocarbons from selected reservoirs are analyzed here.These are field "A", reservoir "L", and field "B", reservoir "K", hydrodynamic unit "K1".
The analysis in the continuation is based on many previous works published for the Sava Depression as well as the entire Croatian part of the Pannonian Basin System (CPBS).Statistical data and mapping of previous regional works are represented in continuation, especially in the Sava Depression.Hence, Malvić [1] presents the use of descriptive statistics (porosity vs. depth) and tests (t-, F-tests) and the Pearson correlation for the interpretation of the influence of increasing depth on the decrease of porosity in the Bjelovar Sub-depression (the Kloštar Ivanić Formation, Pepelana and Poljana Members).The author of [2] identified the kriging method as the most appropriate method for porosity interpolation in the Sava Depression, based on the mapping results done in the Sava Depression sandstones and cross-validation results.For the first time in Croatian geology, [3] applied jack-knifing in the variogram modelling of the reservoir variables for data from the Stari Gradac-Barcs Nygat Field (the Drava Depression).In the same field, [4] applied stochastically improved probability of success calculation for the first time.Reference [5] described the most often applied kriging techniques-simple, ordinary and indicator-and provided the optimal Lagrange value in the ordinary kriging method for the CPBS.Reference [6] gave sequential Gaussian simulation (SGS) results for the Kloštar Field in the Sava Depression.Reference [7] described the criteria for choosing between kriging, co-kriging and simulations using examples from the CPBS.Reference [8] analyzed lithofacies from logs, porosity maps, and made indicator transformation and maps (the Kloštar Field, Sava Depression).In the same field, [9] analyzed reservoir porosity, depth and thickness using 100 realizations of the SGS.Reference [10] applied geostatistics for the most detailed geological model of the Kloštar Field and typical depositional model of the Upper Miocene sandstone reservoirs in the Sava Depression.Reference [11] described the four main applications of geostatistics and recommended for the CPBS such interpolation for the datasets with 20 or more points.Reference [12] pointed out that ordinary kriging was the most often geostatistical technique applied for mapping of the sandstone reservoirs in the Sava Depression.Reference [13] applied the sequential indicator simulations for porosity mapping (of the CPBS), using reservoir thickness maps as the secondary source of information.Reference [14] applied a stochastically improved probability model of geological risk as a new approach to the clastic reservoirs in the CPBS.Reference [15] used ordinary and indicator kriging, as well as SGS, for CO 2 sequestration in oil reservoirs.Deterministic methods for geological risk have been adapted for sequestration purpose.Reference [16] mapped the Bjelovar Sub-depression using ordinary kriging and a digitalized set of older maps.Reference [17] analyzed data from the Šandrovac Field (the Drava Depression) comparing ordinary kriging, inverse distance weighting and nearest neighborhood.Reference [18] used universal kriging for the first time in the Bjelovar Sub-depression and proved it the best interpolation method in the CPBS if a regional trend is visible.Reference [19] described the application of stochastic simulation for the Lower Pontian reservoir in general (the CBPS).Reference [20] presented the problem of the selection of an appropriate mapping method for a dataset with a small number of data (i.e., nearest neighborhood vs. inverse distance weighting).
Regarding the geostatistical theory, there are a large number of available papers and books where the theory and/or its application had been presented in a clear manner.We decided to select several works where the review of the variogram calculations and kriging algorithms had been presented.Thus, [21] proposed a two-step procedure for variogram interpretations in hydrocarbon reservoir modelling, i.e., for modelling spatial variability of variables, such as lithofacies, porosity, etc.The emphasis was given to total variance of the analyzed variable(s), divided into variance regions.The same authors [22] stated the strong connection between the reliable variogram model and depositional model, i.e., the continuity of geological features.It is especially important in the case of geological anisotropies as a reflection of the depositional environment.
Kriging techniques are numerous, and all are developed as variations of simple kriging.For example, [23] presented the application of the sequential indicator simulation algorithm with the goal to obtain a set of equiprobable realizations.They used a conditional variant, respecting the data and simulating other cells.The intention was to show the heterogeneity in the Berea Sandstone reservoir.An interesting application of the secondary variable in kriging matrices has been presented [24] using the collocated co-kriging approach in defining hydrocarbon reservoirs, i.e., to mapping the top of stratum.They correlated wells with the seismic data, and compared kriging with external drift and co-kriging, where both use secondary information as drift or additional variable.Recently, [25] presented a method for connecting porosities from cores and logs for consequent modelling of permeability.Ordinary kriging and co-kriging was applied on the multi-facies, shale and sandstone datasets (South Rumaila Field, Southern Iraq).The permeability had been estimated using porosity data and algorithms of ordinary kriging.Then, the cross-variogram analysis was done to model core permeability (primary) using core porosity (secondary variable), and both are functions of logged values (neutron porosity, water saturation, and their depth).The results had been checked using cross-validation and visual observation of the predicted and observed values.The co-kriging algorithm obtained better results.In addition, [26] applied Bayesian kriging for spatial prediction of the permeability in the hydrocarbon reservoir.All previously applied "conventional" kriging techniques did not correctly estimate the uncertainty in the co-variance structure, so Bayesian kriging was selected with the purpose of taking prior knowledge of the data, such as expert opinions.The uncertainty is given in the form of a posterior probability, and unbiased estimation led to a more realistic view of the reservoir (the small unrealistic "regions" are eliminated).The model had been applied in the sandstone of the tidal depositional environment (the South Rumaila Field), even using the Bayesian approach for stochastic simulations and selection of the P10, P50 and P90 realizations.
Two famous resampling techniques are jack-knifing and bootstrapping.Both are previously used in different geological analysis, especially those of hydrocarbon reservoirs.Reference [27] presented multi-plausible (probable) variograms, where the uncertainty is addressed to each point (lag).They can be basically used for one or more, equally probable, kriging realizations.Multiple plausible variograms may be fitted knowing the uncertainty at each variogram point, 2γ (h).Multiple geostatistical realizations may then be constructed and subjected to process assessment in order to measure the impact of this uncertainty.Theoretically, if all random variables are independent, the Gaussian distribution could be applied, and although it is not the case in reservoir dataset(s), such multi-Gaussianity simplifies the calculation and is mostly acceptable.Authors made several local simulations of four-point locations (so evaluated the 4 th -order moments) and eventually calculated the confidence interval for each variogram point (25, 50, 75 and 95%).Later, [28] presented a generalized bootstrap method for the analysis of semivariogram (i.e., covariance) uncertainties.As the semivariogram is based on measurements, bootstrap could be applied to increase knowledge about the reliability of calculated spatial continuity.Authors used several approaches and datasets, such as LU decomposition and non-Gaussian random field.The bootstrap percentile confidence intervals were not centered around the empirical semivariograms and did not use distributional assumptions for their construction (vs.symmetrical intervals based on standard error, estimated from parametric equation or bootstrap).
All these works have been analyzed with the purpose of solving the problem of a relatively small number of data points, characteristic for almost each particular (and selected) structure in the Sava Depression.Hence, here we first describe two such typical structures in the depression, i.e., their lithostratigraphic parts where there are hydrocarbon reservoirs proven, where most of the available geological variables (petrophysical, saturation, granulometry, injected water) have been measured.This is followed by a relatively simple geostatistical analysis of the reservoir porosity, which is also assumed to have a normal distribution.Such a statement is given as axiomatic, due to an overly small number of hard data to test normality.Unfortunately, new data cannot be additionally collected (most due to expense), so the statistically available methods (jack-knifing) for increasing the number of calculations per each variogram lag and for the "similar datasets" has been tested.

Basic Equations of the Applied (Geo-)Statistical Techniques/Methods
Jack-knifing is a resampling statistical technique, later upgraded into, e.g., bootstrapping.It is useful for statistics estimation, sequentially leaving out each value from datasets and calculating the statistical parameters of remaining data.For example, if the estimated parameter is the population mean (x), using jack-knifing it is possible to calculate the mean of each sub-dataset that includes all but the i-th measurement (x i ) using Equation (1): where x i : mean of sub-dataset, where i-th data is skipped, n: number of data points, j: data currently analyzed.
These n estimations of the mean also form the distribution that characterizes sample statistics with a mean equal to the average of x i .
Variogram (2γ) is the basic geostatistical tool for spatial analysis of a selected regionalized variable.This is defined as the difference between two values of a selected variable at distance h (Equation (2)): where 2γ(h): variogram calculated at distance h; z: value of variable in selected location; z(h): value of variable at distance h from selected location.
Once the experimental variogram is approximated with a theoretical function it can be calculated in the kriging equations only based on distances between two points (not from values).For the N data pairs at distance h (i.e., lag or class), the variogram is calculated from Equation (3): where N: number of data pairs compared at distance h; z: value of variable in selected location; z(h): value of variable at distance h from selected location.The simplest way to present linear estimation is given by Equation (4): where Z k : values of regionalized variable at selected location (estimated from hard data); Z i : known values at the location i (hard-data); λ i : weighting coefficient for location i.
The values of Z i are assumed to have a Gaussian distribution, especially if the number of samples is large (i.e., random variable with a large number of independent samples respects the central limit theorem).Weighting coefficients are calculated using a system of linear equations, such as the kriging matrices (e.g., [29][30][31]).The main advantages of kriging (vs.simpler mathematical methods) are a better consideration of the data anisotropy (by variogram analysis) and possible clustering.Moreover, the values of weighting coefficients in kriging depend on distance, not on values, thanks to the variogram model.This is why such a distance is sometimes known as the "statistical distance".The matrix equitation of kriging is given here with an example of simple kriging (Equation ( 5)): where γ(Z 1 − Z 2 ): semivariogram values between values Z 1 and Z 2 (can be read from variogram curve); γ(Z n − Z n ): 0 (semivariogram cannot be calculated between the same values); γ(x n − x): semivariogram on distance between known value x n and estimated value x. λ n : weighting coefficient for location n.
Simple kriging is a basic technique, but is not unbiased estimation.All other kriging techniques have added some constraint in the equations and are described as BLUE (best linear unbiased estimator) techniques.

Geographical Locations of the Analyzed Hydrocarbon Fields (Northern Croatia)
For the fields "A" (reservoir "L") and "B" (reservoir "K"/HD "K1") two variables were collected (porosity, injected water) and mapped with several methods (nearest neighborhood, inverse distance weighting, kriging).Both fields are located in the Sava Depression, about 90 km south-east from the capital Zagreb (Figure 1) and close to the highway A3 (5 km) and regional railroad (4 km).

Geological History of the Western Part of the Sava Depression
The CPBS is located at the very south-western margin of the PBS.Its marginal location, regional position of the pre-Neogene basement and later tectonics caused structural elongation northwestsoutheast.The entire CPBS is divided into four depressions-Sava, Drava, Mura and Slavonia-Srijem (Figures 1 and 3).Regional geological units of the 2nd order (depressions in the CPBS) and location of the fields "A" and "B" (boxed) in the Sava Depression (from [32]).
The smaller location maps of analyzed fields are presented at Figure 2 (the fields "A" and "B" with blue and red borders, respectively).Both fields are on hilly terrain."A" is an elongated structure with strike NW-SE and altitude 113-216 m, crossed by numerous landslides, and partially forested.Similarly, the altitude of field "B" is 120-231 m.The hills are steep and locally forested.

Geological History of the Western Part of the Sava Depression
The CPBS is located at the very south-western margin of the PBS.Its marginal location, regional position of the pre-Neogene basement and later tectonics caused structural elongation northwestsoutheast.The entire CPBS is divided into four depressions-Sava, Drava, Mura and Slavonia-Srijem (Figures 1 and 3).

Geological History of the Western Part of the Sava Depression
The CPBS is located at the very south-western margin of the PBS.Its marginal location, regional position of the pre-Neogene basement and later tectonics caused structural elongation northwest-southeast.The entire CPBS is divided into four depressions-Sava, Drava, Mura and Slavonia-Srijem (Figures 1 and 3).There are several characteristic lithofacies for the Sava Depression (as well as for the CPBS).Carbonates and coarse-grained clastics dominated in the Middle Miocene (Badenian-Sarmatian), i.e., had been deposited in the marine environments.The Upper Miocene (Lower Pannonian-Upper Pontian) is characterized by two clastic lithotypes-sandstones and marls.Both are deposited in the lacustrine environments [37].The analyzed hydrocarbon reservoirs, in the Sava Depression, belong to that depositional phase (Figure 5), when huge volumes of sand and silt were transported by turbidites [33,36,37].During "calm" periods, different carbonate-rich mud was deposited.Such alterations of clastic is typical for "high" and "low" energy environments and characteristic of the Upper Miocene for the entire CPBS [32].There are several characteristic lithofacies for the Sava Depression (as well as for the CPBS).Carbonates and coarse-grained clastics dominated in the Middle Miocene (Badenian-Sarmatian), i.e., had been deposited in the marine environments.The Upper Miocene (Lower Pannonian-Upper Pontian) is characterized by two clastic lithotypes-sandstones and marls.Both are deposited in the lacustrine environments [37].The analyzed hydrocarbon reservoirs, in the Sava Depression, belong to that depositional phase (Figure 5), when huge volumes of sand and silt were transported by turbidites [33,36,37].During "calm" periods, different carbonate-rich mud was deposited.Such alterations of clastic is typical for "high" and "low" energy environments and characteristic of the Upper Miocene for the entire CPBS [32].There are several characteristic lithofacies for the Sava Depression (as well as for the CPBS).Carbonates and coarse-grained clastics dominated in the Middle Miocene (Badenian-Sarmatian), i.e., had been deposited in the marine environments.The Upper Miocene (Lower Pannonian-Upper Pontian) is characterized by two clastic lithotypes-sandstones and marls.Both are deposited in the lacustrine environments [37].The analyzed hydrocarbon reservoirs, in the Sava Depression, belong to that depositional phase (Figure 5), when huge volumes of sand and silt were transported by turbidites [33,36,37].During "calm" periods, different carbonate-rich mud was deposited.
Such alterations of clastic is typical for "high" and "low" energy environments and characteristic of the Upper Miocene for the entire CPBS [32].
About 95% of Upper Miocene sandstone reservoirs in the CPBS are of turbidite origin.They include several lithofacies.The central part of the structure is "pure" medium-grained sandstones, which are gradually changed into silty or clayey sandstones and siltstones on structural margins.Such heterogeneity determined secondary migration paths, recovery and production regime.
Both analyzed fields ("A" and "B") are located in the western part of the Sava Depression (Figure 5).Present-day production is supported by water injection in the Lower Pontian reservoirs (Figure 6), lithostratigraphically belonging to the Kloštar Ivanić Formation, i.e., the Pepelana Member.The sandstone and marl sequences are determined from spontaneous potential and resistivity logs (Figures 7 and 8).The marls are also isolator rocks for Lower Pontian reservoirs.
Figure 5. Schematic view of regional lacustrine environments at the border of the Sava and Drava Depressions re-formatted from [36].The huge amount of pelitic (marlstone) and turbiditic (sandstone, siltstone) clastics had been deposited [33,37].About 95% of Upper Miocene sandstone reservoirs in the CPBS are of turbidite origin.They include several lithofacies.The central part of the structure is "pure" medium-grained sandstones, which are gradually changed into silty or clayey sandstones and siltstones on structural margins.Such heterogeneity determined secondary migration paths, recovery and production regime.
Both analyzed fields ("A" and "B") are located in the western part of the Sava Depression (Figure 5).Present-day production is supported by water injection in the Lower Pontian reservoirs (Figure 6), lithostratigraphically belonging to the Kloštar Ivanić Formation, i.e., the Pepelana Member.The sandstone and marl sequences are determined from spontaneous potential and resistivity logs (Figures 7 and 8).The marls are also isolator rocks for Lower Pontian reservoirs.

Results
Exploration drilling in the analyzed fields dates to the 1960's.Consequently, several types of reservoir data are available.

Reservoir Type
According to classification [38] the reservoirs in the field "A" are those of bedding type, capped by isolators, and margined with lithological and/or tectonic screens.The quartz dominated sandstones are enriched by mica.The reservoir sandstones are rhythmically interbedded with thin marls or sandy marls.The reservoir "L" is a single hydrodynamic unit.

Granulometry, Porosity, Saturation and Permeability of Sandstones
Granulometry has been calculated from cores and CaCO 3 using calcimetry analysis (Table 1).Based on 619 data from cores taken in 10 wells, the porosity varied from 22.3 to 24.7%.The porosity estimated from logs showed slightly different values: 16.7-20.1% in gas saturated and 15.6-23.9% in oil saturated parts.The saturation is also calculated dually-in 42 wells from cores, and in 36 from logs.The accepted values are 59.8-77.0%for S g and 57.3-71.6 for S o .The permeability has been calculated from the data taken in 7 wells: 8-27 × 10 −3 µm 2 .

Structural Settings
The isopach map of reservoir "L" is given in Figure 9.The selected cross-sections are presented in Figures 10 and 11.The thickness of reservoirs in the Kloštar Ivanić Formation (Figure 10) are mostly constant in the central part (reservoir "N" about 150 m, "B" 20 m, "K" and "JL" 80 m, and "L" 45 m).However, the analyzed reservoir "L" thinned from 70 m (L-161) to only 20 m (L-153).Similar is seen at the section SW-NE (Figure 11), where the thickness of the reservoir "L" from 20 m (L-63, L-57, L-27) gradually increases to 39 m (L-131alfa).

Variogram Analysis and Ordinary Kriging of Porosity
The porosity data was available in 25 wells (Figure 12).The previously interpolated porosity maps were too simplified, mostly including one average value for a single tectonic block.Moreover, kriging has been proven as the most appropriate method for interpolation of different variables in the Upper Miocene sandstone reservoirs of the CPBS (e.g., [2,8,10,18]).This is why, the for selected reservoir, the experimental semivariogram for available porosities had been calculated (Figure 13), and kriging was applied.

Variogram Analysis and Ordinary Kriging of Porosity
The porosity data was available in 25 wells (Figure 12).The previously interpolated porosity maps were too simplified, mostly including one average value for a single tectonic block.Moreover, kriging has been proven as the most appropriate method for interpolation of different variables in the Upper Miocene sandstone reservoirs of the CPBS (e.g., [2,8,10,18]).This is why, the for selected reservoir, the experimental semivariogram for available porosities had been calculated (Figure 13), and kriging was applied.

Variogram Analysis and Ordinary Kriging of Porosity
The porosity data was available in 25 wells (Figure 12).The previously interpolated porosity maps were too simplified, mostly including one average value for a single tectonic block.Moreover, kriging has been proven as the most appropriate method for interpolation of different variables in the Upper Miocene sandstone reservoirs of the CPBS (e.g., [2,8,10,18]).This is why, the for selected reservoir, the experimental semivariogram for available porosities had been calculated (Figure 13), and kriging was applied.Unfortunately, the semivariogram is characterized by large oscillations and a low number of data pairs per class.Thus, the approximation was possible only with the simplest, linear, model.The ordinary kriging map is shown in Figure 14.
The cross-validation (where porosity is expressed as parts of units, not as percentages) for the kriged map is 0.000676.However, the variogram model is obviously characterized with large uncertainties (estimation of nugget, number of data pairs).This is why it has been decided to apply the jack-knifing statistical method with the purpose of artificially increasing the number of data points.Jack-knifing had been a previously applied or described in the CPBS [3,39].Basically, this method has a similar algorithm as cross-validation, where data or groups of data have been "deleted" and the estimation is done from the rest.Such a re-sampling method is characterized by standard error (e.g., [40][41][42]) and can be used as one of the improvements in the basic kriging techniques or data intended to interpolate using kriging (e.g., digital terrain model (DTM) data, as in [43]).
In this study, jackknifing had been applied for each variogram class (step), which was calculated several times.The final value was the approximation of each calculation set for the same class (lag).The lag (of class of data) has been determined experimentally, i.e., in several iterations.The selected lag tolerance was half the lag value, so there was no overlapping of data in two classes.The value of the lag did not change significantly in the (a) shape of oscillation or (b) number of data pairs per class, which indicated that the number of original data points was too small for an advanced and more reliable spatial (variogram) modelling.
Consequently, the number of data pairs included in the calculation of class was much larger.The "jack-knifed" semivariogram for the same reservoir porosity is given in Figure 15.Oscillations are largely decreased in several points, and the experimental semivariogram is approximated with the exponential theoretical model, significantly decreasing the range.Unfortunately, the semivariogram is characterized by large oscillations and a low number of data pairs per class.Thus, the approximation was possible only with the simplest, linear, model.The ordinary kriging map is shown in Figure 14.The cross-validation (where porosity is expressed as parts of units, not as percentages) for the kriged map is 0.000676.However, the variogram model is obviously characterized with large uncertainties (estimation of nugget, number of data pairs).This is why it has been decided to apply the jack-knifing statistical method with the purpose of artificially increasing the number of data points.Jack-knifing had been a previously applied or described in the CPBS [3,39].Basically, this method has a similar algorithm as cross-validation, where data or groups of data have been "deleted" and the estimation is done from the rest.Such a re-sampling method is characterized by standard error (e.g., [40][41][42]) and can be used as one of the improvements in the basic kriging techniques or data intended to interpolate using kriging (e.g., digital terrain model (DTM) data, as in [43]).
In this study, jackknifing had been applied for each variogram class (step), which was calculated several times.The final value was the approximation of each calculation set for the same class (lag).The lag (of class of data) has been determined experimentally, i.e., in several iterations.The selected lag tolerance was half the lag value, so there was no overlapping of data in two classes.The value of the lag did not change significantly in the (a) shape of oscillation or (b) number of data pairs per class, which indicated that the number of original data points was too small for an advanced and more reliable spatial (variogram) modelling.
Consequently, the number of data pairs included in the calculation of class was much larger.The "jack-knifed" semivariogram for the same reservoir porosity is given in Figure 15.Oscillations are largely decreased in several points, and the experimental semivariogram is approximated with the exponential theoretical model, significantly decreasing the range.Using the new semivariogram, the new porosity map was interpolated (Figure 16).Due to a much smaller range, the new porosity map is characterized with numerous "bull-eye" shapes.
Using the new semivariogram, the new porosity map was interpolated (Figure 16).Due to a much smaller range, the new porosity map is characterized with numerous "bull-eye" shapes.3.2.The Field "B", Reservoir "K", Hydrodynamic Unit "K1"

Reservoir Type
According to [37], the reservoirs in the field "B" are of bedding type, capped by isolators, and margined with lithological and/or tectonic screens.In that brachianticline, with a northwestsoutheast strike, sandstone reservoirs are enriched by mica minerals and quartz and occasionally interbedded with marls and (shallower) sandy marls.The "K" reservoir is divided into several tectonic blocks, where each is a mostly separated hydrodynamic (HD) unit.This is the largest such unit analyzed here (the HD "K1").

Figure 16.
Porosity map for reservoir "L" interpolated by ordinary kriging with jack-knifed semivariogram.

Reservoir Type
According to [37], the reservoirs in the field "B" are of bedding type, capped by isolators, and margined with lithological and/or tectonic screens.
In that brachianticline, with a northwest-southeast strike, sandstone reservoirs are enriched by mica minerals and quartz and occasionally interbedded with marls and (shallower) sandy marls.The "K" reservoir is divided into several tectonic blocks, where each is a mostly separated hydrodynamic (HD) unit.This is the largest such unit analyzed here (the HD "K1").

Structural Settings
The HD "K1" unit is the largest tectonic block and HD in the reservoir "K".The part where the reservoir pressure is supported by water injection is given in Figure 17 The HD "K1" unit is the largest tectonic block and HD in the reservoir "K".The part where the reservoir pressure is supported by water injection is given in Figure 17.This reservoir has been crossed with two sections (Figure 17).The first has strike approximately NW-SE (Figure 18), and the second NE-SW (Figure 19).The section in Figure 18 shows that the reservoirs have similar thickness in both wells ("N" about 120 m, "K" 100 m, "Z" 40 m).Some marls pinch out, which could define separated hydrodynamic units.The next section (Figure 19) clearly shows partial structural inversion in the Kloštar Ivanić Formation, i.e., transition from monocline (deeper) to fold (shallower).This reservoir has been crossed with two sections (Figure 17).The first has strike approximately NW-SE (Figure 18), and the second NE-SW (Figure 19).The section in Figure 18 shows that the reservoirs have similar thickness in both wells ("N" about 120 m, "K" 100 m, "Z" 40 m).Some marls pinch out, which could define separated hydrodynamic units.The next section (Figure 19) clearly shows partial structural inversion in the Kloštar Ivanić Formation, i.e., transition from monocline (deeper) to fold (shallower).

Variogram Analysis and Ordinary Kriging of Porosity
Porosity data were available in 19 wells (Figure 20).The porosity data were analyzed using a semivariogram, and the experimental one is shown in Figure 21.

Variogram Analysis and Ordinary Kriging of Porosity
Porosity data were available in 19 wells (Figure 20).The porosity data were analyzed using a semivariogram, and the experimental one is shown in Figure 21.Due to large uncertainties and a small number of data pairs, the semivariogram is approximated with a simple linear model.Such a model did not allow the estimation of nugget and certain range.It was applied in ordinary kriging and resulting map is shown in Figure 22.The cross-validation is 0.0013197.Again, the new semivariogram was calculated with "re-sampled", i.e., jack-knifed, dataset.The "one-leave-out" technique was used here, i.e., the jack-knife estimation was done for each sub-dataset where the i-th data was omitted (i = 1, n; n is total number of data points).Practically, this means that each variogram value per lag was calculated n times, and averaged.The direct consequence is that the number of data pairs per lag (class) was drastically increased, e.g., for the distance 80 m it was increased from 2 to 32 data pairs, and secondly, the lag distance could be decreased.
The new experimental model is given in Figure 23.This could be approximated with the Gaussian theoretical model and the resulting kriging map is shown in Figure 24.The semivariogram sill is lower and range larger (based on goodness of fit) than in the previous variogram model based on the original dataset (19 points).The new kriged porosity map is shown in Figure 24 and the accompanying cross-validation value is 0.0009704.Again, the new semivariogram was calculated with "re-sampled", i.e., jack-knifed, dataset.The "one-leave-out" technique was used here, i.e., the jack-knife estimation was done for each sub-dataset where the i-th data was omitted (i = 1, n; n is total number of data points).Practically, this means that each variogram value per lag was calculated n times, and averaged.The direct consequence is that the number of data pairs per lag (class) was drastically increased, e.g., for the distance 80 m it was increased from 2 to 32 data pairs, and secondly, the lag distance could be decreased.
The new experimental model is given in Figure 23.This could be approximated with the Gaussian theoretical model and the resulting kriging map is shown in Figure 24.The semivariogram sill is lower and range larger (based on goodness of fit) than in the previous variogram model based on the original dataset (19 points).The new kriged porosity map is shown in Figure 24 and the accompanying cross-validation value is 0.0009704.Porosity map for reservoir/HD "K1" interpolated by ordinary kriging with jack-knifed semivariogram.

Discussion and Conclusions
The presented analysis is the first of such a kind in the Sava Depression (Northern Croatia).It represents the continuation of previous geostatistical analyses conducted in that depression and the entire CPBS.The permanent problem of small datasets is emphasized as crucial for the subsurface mapping of different geological variables.In numerous analyzes, new hard data can hardly be obtained, mostly due to the costs of drilling (logging or seismic).Consequently, any statistical method that could increase the reliability of existing datasets and their analytical results is more than welcome.Jack-knifing is one such method.Previously, it was successfully applied in the (adjacent) Drava Depression, and here, for the first time, in the Sava Depression.The method could be highly appropriate for datasets of 15-30 points, where the basic, descriptive statistics are more or less representative (variance and mean), and the Gaussian distribution can be assumed.
In both analyzed reservoirs, the original semivariogram results were highly uncertain, with large oscillations, a small number of data pairs per class and unknown nugget.Consequently, the linear model was the only acceptable theoretical model to use.This was why a much larger number of re-sampled data points was calculated using jack-knifing.As a result, the new semivariograms could be approximated more reliably and with mathematically advanced theoretical models (exponential and Gaussian).
The hypothesis that new kriged maps interpolated from "jack-knifed" semivariograms are more reliable was eventually tested (a) visually (maps without the "bull-eye" or "butterfly" effects are better) and (b) numerically, using cross-validation.The, so called "bull-eyes" or "butterflies" have the same origin, and both can indicate contradiction between the data and the spatial continuity model.Such contradiction is the result of the "computing algorithm", which allows us to have two contour lines of the same value in the same point.However, the probability for such an event is almost infinitesimally low.This is why two such effects (both are similar, but the bull-eye is mostly circular and the butterfly ellipsoidal) indicate errors in the applied algorithm, the majority resulting from small spatial continuity (such as variogram range).

Figure 1 .
Figure 1.Regional geological units of the 2nd order (depressions in the CPBS) and location of the fields "A" and "B" (boxed) in the Sava Depression (from[32]).

Figure 2 .
Figure 2. Satellite maps of the fields "A" and "B" (taken from Google).The figure represents the area boxed in Figure 1, and locations of analyzed fields.

Figure 1 .
Figure 1.Regional geological units of the 2nd order (depressions in the CPBS) and location of the fields "A" and "B" (boxed) in the Sava Depression (from[32]).

Figure 2 .
Figure 2. Satellite maps of the fields "A" and "B" (taken from Google).The figure represents the area boxed in Figure 1, and locations of analyzed fields.

Figure 2 .
Figure 2. Satellite maps of the fields "A" and "B" (taken from Google).The figure represents the area boxed in Figure 1, and locations of analyzed fields.

Figure 3 .
Figure 3. Regional geological units of the 1st order in the Pannonian Basin System (PBS) (modified from[33,34]).The Sava Depression is located at very southwestern margin (blue circle).

Figure 4 .
Figure 4. Timescale of the main tectonic and depositional events during Neogene and Quaternary in the CPBS [36].Analyzed sandstone reservoirs belong to Early Pontian stage, i.e., to late 2nd transtensional phase.

Figure 3 . 25 Figure 3 .
Figure 3. Regional geological units of the 1st order in the Pannonian Basin System (PBS) (modified from [33,34]).The Sava Depression is located at very southwestern margin (blue circle).The CPBS had been covered by shallow marine environments of the Paratethys during the Badenian and Sarmatian (16.4-11.5 Ma).It was reduced in the Pannonian Lake during the Lower Pannonian (11.5-9.3Ma), and disintegrated during the Upper Pannonian (9.3-7.1 Ma) and Pontian (7.1-5.7 Ma) into several smaller lakes.The next, relatively short, phase of a single Slavonian Lake [35] happened during the Pliocene.From the Upper Romanian the entire CPBS has been covered with continental facies, which also lasts today (Holocene).The main stages of geological evolution in the CPBS are shown on Figure 4.

Figure 4 .
Figure 4. Timescale of the main tectonic and depositional events during Neogene and Quaternary in the CPBS [36].Analyzed sandstone reservoirs belong to Early Pontian stage, i.e., to late 2nd transtensional phase.

Figure 4 .
Figure 4. Timescale of the main tectonic and depositional events during Neogene and Quaternary in the CPBS [36].Analyzed sandstone reservoirs belong to Early Pontian stage, i.e., to late 2nd transtensional phase.

Figure 6 .
Figure 6.Typical geological column for the fields "A" and "B".The formal and informal lithostratigraphic units are listed.

Figure 6 .
Figure 6.Typical geological column for the fields "A" and "B".The formal and informal lithostratigraphic units are listed.

Figure 7 .
Figure 7.The composite log-lithological column in the field "A" (from archive of INA Plc. and redrawn).Figure 7. The composite log-lithological column in the field "A" (from archive of INA Plc. and re-drawn).

Figure 7 .
Figure 7.The composite log-lithological column in the field "A" (from archive of INA Plc. and redrawn).Figure 7. The composite log-lithological column in the field "A" (from archive of INA Plc. and re-drawn).

Figure 8 .
Figure 8.The composite log-lithological column in the field "A" (from archive of INA Plc. and redrawn).

Figure 8 .
Figure 8.The composite log-lithological column in the field "A" (from archive of INA Plc. and re-drawn).

Figure 9 .
Figure 9. Isopach map of part of the reservoir "L", where reservoir pressure is supported with water injection (from archive of INA Plc. and re-drawn).

Figure 9 .
Figure 9. Isopach map of part of the reservoir "L", where reservoir pressure is supported with water injection (from archive of INA Plc. and re-drawn).

Geosciences 2019, 9 , 25 Figure 9 .
Figure 9. Isopach map of part of the reservoir "L", where reservoir pressure is supported with water injection (from archive of INA Plc. and re-drawn).

Figure 12 .
Figure 12.Porosity data available for reservoir "L" (from archive of INA Plc. and re-drawn).

Figure 12 .
Figure 12.Porosity data available for reservoir "L" (from archive of INA Plc. and re-drawn).Figure 12. Porosity data available for reservoir "L" (from archive of INA Plc. and re-drawn).

Figure 12 .
Figure 12.Porosity data available for reservoir "L" (from archive of INA Plc. and re-drawn).Figure 12. Porosity data available for reservoir "L" (from archive of INA Plc. and re-drawn).

Figure 14 .
Figure 14.Porosity map for reservoir "L" interpolated by ordinary kriging.Figure 14.Porosity map for reservoir "L" interpolated by ordinary kriging.

Figure 14 .
Figure 14.Porosity map for reservoir "L" interpolated by ordinary kriging.Figure 14.Porosity map for reservoir "L" interpolated by ordinary kriging.

Figure 17 .
Figure 17.Isopach map of the HD unit "K1", where reservoir pressure is supported with water injection (from archive of INA Plc.).

Figure 17 .
Figure 17.Isopach map of the HD unit "K1", where reservoir pressure is supported with water injection (from archive of INA Plc.).

Figure 21 .
Figure 21.Experimental semivariogram with approximation, reservoir/HD "K1", porosity.Due to large uncertainties and a small number of data pairs, the semivariogram is approximated with a simple linear model.Such a model did not allow the estimation of nugget and certain range.It was applied in ordinary kriging and resulting map is shown in Figure22.The cross-validation is 0.0013197.

Table 1 .
Granulometry of reservoir "L" (M d -mean grain diameter, S o -coefficient of sorting, S k -coefficient of asymmetry) (from archive of INA Plc.).

Table 2 .
Granulometry of reservoir "K" (M d -mean grain diameter, S o -coefficient of sorting, S k -coefficient of asymmetry) (from archive of INA Plc.). .