Remapping of depth of e-log markers between Neogene basement and Lower/Upper Pannonian border in the Bjelovar Subdepression

Abstract

Within the Croatian part of the Pannonian Basin System, there are several depressions and subdepressions. A geostatistical analysis was performed in the Bjelovar Subdepression, which is part of the larger Drava Depression. Data for the analysis were obtained from structural maps of the subdepression, previously made by hand interpolation. Square grid cells (2 × 2 km) were laid over the maps. The analysis of cell data included a sufficient number of data (>700 per map) and the semivariogram analyses after which experimental and theoretical semivariograms were made. Consequently a new map set was interpolated by the Ordinary Kriging technique including the pre-Neogene palaeorelief border (‘Tg’) and e-log markers (‘Rs7’ and ‘Rs5’). New maps were subsequently compared with the older, hand-interpolated maps. It was shown that the method of transferring non-computer into computer interpolated maps via/using an artificial, regular grid of point data is an effective way of ‘digitizing’ older paper-based maps, preserving almost all of the regional information.

Keywords:

• geostatistics
• Kriging
• structural maps
• Bjelovar Subdepression
• Croatia

1. Introduction

Geostatistics is a scientific and professional subdiscipline of spatial analysis. The use of geostatistics began with modern reservoir characterization in the 1980s. But the overly simplistic ‘picture’ showed that new tools for a more detailed study and evaluation of geological structure were needed. Geostatistics brings order to datasets of different origin and quality, thereby assisting in the construction of realistic 3D geological models of certain heterogeneity and the quantification of uncertainty by generating a number of possible models (e.g. Dubrule, 1998; Isaaks & Srivastava, 1989).

Geostatistics can be used to evaluate geological variables such as porosity, permeability, depth, and thickness. One of the input data is the location of the sampling point of geological exposure, and is defined by two coordinates (x, y) and a height or depth value (z). Analysis of geological data can be accessed using two methods – deterministic and stochastic.

The deterministic approach implies that all parameters and variables of the conceptual model are known, and within such a model, all measured variables are known at each point. Such models are very rare in practice. The stochastic approaches always try to preserve uncertainties included in the data, due to the number of points, variables, distribution, statistics, etc. (e.g. Malvić, 2008). In this way the validity of results can be better explained, but the interpretation is weighted with a range of solutions (not only one like in determinism). When using this method, the result is the mean value and the standard deviation as a measure of uncertainty. A variogram can show the structure of a random field, which may depend on the distance between two points and the difference between them (e.g. Andričević, Gotovac, & Ljubenkov, 2007).

2. Short review of the Bjelovar Subdepression geological settings

The Bjelovar Subdepression (2900 km²) represents a large and geomorphologically closed southwestern part of the Drava Depression (12,000 km²). The location in WGS84 projection is given on the Main Map, and the positions of numerical data points with analysed marker/border depths on Map 2. Characteristically, the subdepression is surrounded by the following mountains – Kalnik to the northwest, Bilogora to the northeast, Papuk, Ravna gora and Psunj to the east and southeast, and Moslavačka gora to the south. The subsurface rock's lithologies are defined within two groups. The first group belongs to the younger clastic rocks of the Neogene-Quaternary periods, while the second group consists of Palaeozoic and Mesozoic aged magmatites, metamorphites and carbonates. These two groups are lithologically very different which means that directly measured samples (such as cores) or indirectly measured values (such as e-logs or seismic) can be easily distinguished. The very detailed review of the subdepression's geological settings and palaeoenviroments is given in Malvić (2003).

The younger clastics have cyclic features, because they were deposited during three main sedimentation megacycles (Saftić, Velić, Szatano, Juhasz, & Ivković, 2003; Velić, 2007; Velić, Weisser, Saftić, Vrbanac, & Ivković, 2002). Each cycle corresponds to certain lithostratigraphic formations. They are a product of regional tectonic events, transtensional and transpressional, as well as the changes of detritus sources during the last approx. 16 Ma (Map 3). Formations are distinguished based on lithological criterium, but simultaneously, they (as well as several other members) are mostly bordered with regional marker beds (markers) which are recognized as lacustrine marls (e.g. Malvić & Cvetković, 2013). They are also easily recognizable on e-logs (especially resistivity) and are, interestingly, mostly synchronous on the regional scale, i.e. in the entire Croatian part of the Pannonian Basin System (CPBS). Two such e-log markers are used here. The older is called ‘Rs7’, marking the border between Sarmatian and Lower Pannonian (i.e. Mosti and Križevci Members, Map 4). The younger is ‘Rs5’, marking the border between Lower and Upper Pannonian (i.e. Moslavačka gora and Ivanić-Grad Formations, Map 5). The exception is the physical regional border, unconformity, called ‘Tg’ (or sometimes ‘Pt’), which represents the border between the oldest Neogene clastics and pre-Neogene basement (Map 3). It is a result of the long-lasting regional emersion during the entire Palaeogene and Lower Neogene in the entire CPBS.

3. Geostatistical input data and calculated Semivariograms

The data for geostatistical analysis were obtained from hand-made (hand-interpolated) structural maps of the Bjelovar Subdepression (Malvić, 2011). Depth values for 6 e-log markers and 1 e-log border are observed. Namely, the values (from the oldest) of ‘Tg’, ‘Rs7’, and ‘Rs5’ are given on the subdepression palaeostructural maps. These maps represent the periods of 1st transtension and 1st transpression in the CPBS, i.e. approximately from the Lower Badenian to the end of Lower Pannonian. The input data for the geostatistical analysis represents the new values for the particular depths related to analysed markers or the border. In fact, each palaeostructural map was re-sampled using the regular 2 × 2 km grid (constructed in AutoCAD 2013). The grids were laid over existing maps and the following values were read for each point (the centre of each grid square): (1) X and Y coordinates (Map 2) and (2) value from hand-interpolated map. These datasets (3 of them) were used for the re-interpolation of analysed palaeostructural maps, this time, however, from the set of 700 re-sampled points per map. The first step in re-mapping was the calculation of experimental semivariograms and approximation with the specific theoretical model (Figures 1^{Figure 2}–3). In geological practice, three models are mostly used – spherical, exponential, and Gaussian (e.g. Malvić, 2008). In this case, the Gaussian (Figure 1) and exponential (Figures 2 and 3) are applied. Both are a logical reflection of mapped lithologies. In the pre-Neogene basement depth changes relatively fast over small distances, due to frequent and different sized buried hills along the palaeosea bottom of the Central Paratethys in this area (e.g. Malvić, 2012). However, the depth of Neogene clastics gradually uplifts from central to marginal subdepression zones. Such gradual transition was well reflected in the exponential model. The semivariograms were calculated in Variowin 2.21 21 (Pannatier, 1996). The maps were interpolated in Surfer 2011 using Ordinary Kriging.

Figure 1. (a – up) the experimental semivariogram obtained from structural palaeorelief maps e-log border ‘Tg’ of the main axis (120°); (b – down) approximated semivariogram for the main axis of the Gaussian model.

Figure 2. (a – up) the experimental semivariogram obtained from structural palaeorelief maps e-log marker ‘Rs7’ of the main axis (120°); (b – down) approximated semivariogram for the main axis of the exponential model.

Figure 3. (a – up) the experimental semivariogram obtained from structural maps palaeorelief e-log markers ‘Rs5’ to the main axis (120°); (b – down) approximated semivariogram for the main axis of the exponential model.

The semivariogram showed the variability of analysed variables, which depended on the distance between measured points (it was a regular grid, which was favourable) and the difference value between them (which depended on previous maps). As always, experimental semivariograms included all imperfections that are subjected to measurements and variogram modelling itself (as example, see Andričević et al., 2007). However, a significant help to the variogram modelling was the database of the characteristics variograms calculated in different lithofacies for the CPBS (Malvić, Velić, & Cvetković, 2010). Variogram analysis was accompanied with the calculation of variogram surface, which had been used for the definition of the variogram model axis. As expected, they corresponded highly with the regional structural axis of the subdepression. All semivariograms (Figures 1^{Figure 2}–3) had been calculated with the purpose to reach values of lag distance and ranges, both for the major and minor axis. The following values are used in all of them:

1. On the main axis: lag (h) = 2000 m, the number of lags = 20, direction = 120°,

2. On the minor axis: lag (h) = 2500 m, the number of lags = 15, direction = 210°.

3.1. Semivariograms for e-log border ‘Tg’

‘Tg’ reflects the shape of pre-Neogene palaeorelief, i.e. the beginning of the regional transgression. The theoretical model was Gaussian (Figure 1). The range was set at 30 km, i.e. at the highest distance where autocorrelation could be observed. Experimental results are: nugget = 0; range = 30,000 m; sill = 600,000. The smallest square error of theoretical fitting was 4.0778 × 10⁻³.

3.2. Semivariograms for e-log marker ‘Rs7’

The theoretical model was exponential (Figure 2). The range was 24 km. Reduction of the range can be explained by the reduction of depositional environment, transition from marine into large lacustrine, where continuity of strata was also decreased. Experimental results are: nugget = 0; range = 23,918 m, sill = 500,000. The smallest square error of theoretical fitting was 2.8125 × 10⁻².

3.3. Semivariograms for e-log marker ‘Rs5’

The theoretical model was exponential (Figure 3). The input map showed a large unconformity at the eastern part. It reduced the number of re-sampled data and, consequently, the range. The reduction of the environment continued, revealing that was the final phase of changing from marine to brackish lake (e.g. Pavelić, 2001). All events coincided with the calculated semivariogram. Experimental results are: nugget = 0, range = 15,000 m, sill = 550,000. The smallest square error of theoretical fitting was 1.0769 × 10⁻¹.

4. Interpretation of kriged maps

The Map 3 is palaeorelief (e-log border ‘Tg’) of the Neogene basement. The orange colour represents the shallowest parts of the basement depth while the blue represents the deepest basement up to 3600 m. The map shows that the elongation of the basin is northwest-southeast. In the central part there are several remnants of smaller depositional centres which generally follow the strike of the basin. To the south, ‘Tg’ gradually becomes shallower. Faults from original maps were not transferred onto the new grid, because resolution was not high enough. However, the main fault zone could be interpreted from the stratoisohypse shape and areal density and is transversal with the main subdepression axis.

The Pannonian basement is represented by the e-log marker ‘Rs7’ shown in Map 4. The deepest parts are located in the central part and are oriented NE-SW, with the shallower ones on the north side. This is explained by a fault whose displacement is almost 1000 m. The assumed fault zone from Map 3 is also visible here. Two depositional centres are visible. The smaller is in the central part and the larger one extends from the central towards southern part and is separated by anticline along Psunj Mt. The deepest part is up to 3200 m, and the shallowest is exposed at the surface along the marginal mountains.

The Map 5 reflects structures along e-log marker ‘Rs5’. The essential difference from previous maps is the regional unconformity that covers the eastern part. It reflects the reduction of the sedimentary basin that caused a considerably smaller range of the variogram. The regional synclines in the central part are heritage structures from Map 4.

5. Discussion and conclusion

In practice, geostatistical methods are used mostly for mapping different point data of one or more spatial variables, where normal distribution is assumed. Point data, regardless of the variables, were considered as hard data. But that's not entirely true. Uncertainties in data occur because of technical disadvantages in logging equipment or the drilling process. Such technical disadvantages decrease accuracy and affect knowing the right position. However, all these problems do not generally influence data sources in an analysed area. Generally, any dataset larger than 30 data points is statistically representative. This number of data was overreached several times in the present mapping. These input data are a key factor that influences the results obtained by using the Kriging method. In this case, more than 700 depth values were used as input data for the map, with the exception of Map 5.

The first uncertainties appeared in the determination of semivariogram axes. In this case, variogram surface maps compared with palaeostructural maps confirmed that main and secondary axes in those two cases were the same. Moreover, the second uncertainty was part of the range and nugget determination, followed by theoretical models. The largest range of 30 km is calculated for the basement of the Neogene sediments. The smallest range of 15 km was calculated for the e-log marker ‘Rs5’ which was interpreted as the first stage of a reduced, brackish lake located in the subdepression. Such a gradual reduction in the range coincides with the progressive transtension and local uplifting, reducing the ‘continuous’ properties of rocks, thereby decreasing values of autocorrelation in rock or structural variables.

The third uncertainty was the Kriging algorithm itself. All variogram values were transferred in/during the interpolation, logically including possible ‘errors’ or subjective determinations. Next was the selection of the Kriging techniques. The Ordinary Kriging has been proven in several papers as the ‘best’ interpolation approach for subsurface geological variables in the CPBS (e.g. Balić, Velić, & Malvić, 2008; Husanović & Malvić, 2014).

The proposed approach for the final evaluation is considered valuable for practical mapping. The first was a visual comparison of the previous, hand-interpolated maps with the new set of Kriging maps. It was concluded that, on the kilometre scale, structures observed in the original set are also preserved on the new maps. The only differences are the fault zones that are not interpreted, but could be easily recognized regionally from the stratoisohypses density per area unit. The second was the interpretation of mapped structure in accordance with the regional model of geological evolution. Maps 3–5 clearly show locally developed depositional centres and changes in structural forms, which can be followed in previous works. For the sake of consistency, and consequently resolution, computer interpolated maps are ‘non-biased’ due to interpolation on small scale where two or more isoline shape variations are partially possible. Regardless of how many times the same input is provided, a deterministic, computer algorithm will always produce the same output. This is not a condition or a rule for interpolation done by a human being. Both approaches have advantages and disadvantages, but if the goal is to ‘digitize’ ‘old’ maps with a new point set, the presented approach can be favourable. However, geostatistical maps are not perfect and still include uncertainties connected with variogram models and sampling set. Therefore, they need to be expertly evaluated for each map.

Software

The semivariograms were calculated and approximated with open-source (freeware) Variowin 2.2, created by Y. Pannatier (1996). The interpolation was done in licensed software Surfer 8 and AutoCAD 2013 (installed at University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering, Department of Geology and Geological Engineering).

Acknowledgement

This analysis was performed during the ‘Development of geomathematical methods for analysis of Neogene depositional environments on the Croatian part of the Pannonian Basin System’ project (led by T. Malvić), funded from University of Zagreb as part of the ‘Support to researching II in the 2013’ programme, as well as a thematic work in research project titled ‘Stratigraphical and geomathematical research of petroleum geological systems in Croatia’ (project no. 195-1951293-0237, led by J. Velić), financed by the Ministry of Science, Education and Sports of the Republic of Croatia.