Geoenergy potential of the Croatian part of Pannonian Basin: insights from the reconstruction of the pre-Neogene basement unconformity

Presented work focuses on the importance of unconformity that separates the Neogene infill from older Palaeozoic and Mesozoic rocks in the Croatian part of Pannonian Basin. Structure map of this horizon nearly represents the thickness map of the Neogene and Quaternary basin fill. Rock formations just below the unconformity are significantly weathered, which results in favourable petrophysical properties, making them interesting from the aspect of geoenergy potential. The pre-Neogene surface was constructed in 1:400,000 scale using publicly available subsurface maps of different scale and different level of detail. Harmonization and compilation of these maps enabled construction of a structured surface with near-vertical fault planes. Supplemental maps were constructed via basin modelling, showing the temperature distribution in the subsurface, potential source rock maturity near the mapped horizon, surface heat flow and geothermal gradient distribution. Constructed maps illustrate the importance of the mapped interval for regional planning of future geoenergy-related research.. ARTICLE HISTORY Received 29 October 2018 Revised 7 June 2019 Accepted 15 July 2019


Introduction
Geoenergy today is not anymore related to fossil fuels and exploration for oil and gas. It is more and more oriented to using the developed methods for investigating the deep subsurface rock formations in order to map other resourcesdeep geothermal, shallow geothermal, CO 2 geological storage and energy storage. Some of these resources, like geothermal, have already emerged and are being hastily developed, while others, CO 2 geological storage and energy storage are still investigated throughout the World, and not only in the economically developed countries. When evaluating the geoenergy potential of a certain area, it is of utmost importance to understand its subsurface geology and thicknesses of the most important prospective formations. This kind of insight can only be obtained by a construction of a geological model from a set of maps that define the subsurface volume.
Study area covers the Croatian part of the Pannonian Basin (CPB, Main Map) which spreads over the entire northern part of Croatia. Commonly, publicly available subsurface maps of the study area are made either in a large scale, showing only a small portion of the area, naturally in greater detail (e. g. Baketarić & Cvetković, 2015;Novak Zelenika, Cvetković, Malvić, Velić, & Sremac, 2013;Špelić, Malvić, Saraf, & Zalović, 2016) or cover the whole study area but lack in detail, illustrating only the major features (e.g. Saftić, Velić, Sztanó, Juhász, & Ivković, 2003;Velić, Weisser, Saftić, Vrbanac, & Ivković, 2002). To create a regional scale map with sufficient detail necessary for evaluation of the subsurface potential, a compilation of both groups of maps is needed.
In this work, the constructed map (Main Map) shows the depth of the pre-Neogene basement, which is overlain by a thick succession of dominantly Neogene sediments (up to 7000 m thick in the Drava Depression; according to Saftić et al., 2003 andVelić, 2007) that are in most parts of the CPB covered with locally thick Quaternary deposits (>400 m thick; Brkić, 2017;Hernitz, Kovačević, Velić, & Urli, 1981;Prelogović & Velić, 1992). Pre-Neogene basement surface represents one of the most important correlation horizons in CPB. In contrast to other horizons within the Neogene sedimentary succession, it can be recognized on most of the well and seismic data as it marks a strong unconformity (Saftić et al., 2003;Velić, 2007).
Since the mapped unconformity surface presents the depth of a base of an interval comprising the Neogene-Quaternary basin infill in which majority of the Croatian hydrocarbon reserves are found, it also represents an isopach map on a regional scale. It marks an interesting horizon underneath which hydrocarbon accumulations and geothermal reservoirs can occur in the fractured and weathered basement rocks. The sediments of the Neogene age basin fill, situated just above or near the unconformity, can contain source rocks whose maturity level depends on the depth, age and temperature distribution in the subsurface. Several thematic supplemental maps were constructed that are important from the aspect of geoenergy potential. These present the surface heat flow distribution (mW/m 2 ), averaged geothermal gradient (°C/100 m), temperature distribution along the pre-Neogene unconformity (°C) and maturity of the source rocks near the pre-Neogene unconformity based on vitrinite reflectance values (%Ro). The maps were constructed by means of basin modelling. The Main Map, together with supplemental maps enables the general inferences on geoenergy potential of the study area.
In the CPB, the pre-Neogene basement rocks incorporate carbonate, igneous and metamorphic rock units of Mesozoic and Palaeozoic age that are occasionally overprinted by Early Cretaceous metamorphism (Matoš, 2014 with references). These rocks are found on the surface, on inselbergs in the western and central part of continental Croatia (Figure 1), but also in nearby mountains in northern Bosnia, north-western Serbia and southern Hungary (Pamić, 1998).
Sedimentary rocks of Palaeozoic and Mesozoic age are characterized by numerous lithotypes, which are associated with changes in paleogeographical environments from Devonian to Cretaceous time (Jamičić, Vragović, & Matičec, 1989). Palaeozoic and Mesozoic sedimentary successions are proven hydrocarbon reservoirs within a predominantly Neogene petroleum system in the southern part of PB. Triassic dolomites are known reservoirs in the Mura Depression, whereas mixed Palaeozoic-Mesozoic tectonized carbonates and breccia-conglomerates are the most important reservoirs in the western part of Drava Depression (Tadej, 2011;Velić, 2007). Mesozoic dolomite massif is the principal source of clastic material in certain reservoirs in the eastern part of Drava Depression, whereas Palaeozoic to Mesozoic sandstones, limestones and schists are known petroleum reservoirs in the Slavonia-Srijem Depression (Tišljar, 1993;Velić, 2007). Also, the reservoir of the Velika Ciglena geothermal field is situated in fractured Triassic carbonates in the area between Sava and Drava Depression.

Materials and methods
Construction of a structured subsurface map and supplemental maps showing the geoenergy potential requires several steps. These include preparation of input data, simple surface construction, construction of a structured surface and basin modelling. ArcMap 10.1 (input data preparation, legacy data digitalization), Schlumberger Petrel (simple and structured surface construction) and Schlumberger PetroMod (basin modelling) software were used in process of construction of the Main Map and supplemental maps.

Input data
Within this work, an extensive input dataset was compiled from different sources, then evaluated and in part interpreted in order to construct the Main Map and the supplemental maps that are related to the characterization of the geoenergy potential.

Input data for the Main Map
Input data for the pre-Neogene basement surface construction can be divided into the dataset for the construction of simple, non-faulted surface and fault geometry information for the construction of a structured surface. Data for simple surface construction was gathered from various published sources, which differ greatly in quality in respect of the input data available in time of construction of these previous maps. Regarding the data quality, four categories can be defined (Figure 2).
First category of the input data relates to the regions where the information about the surface geometry was extracted directly from geological models which were constructed based on the interpretation of the seismic data acquired by recent exploration. These represent the most accurate and reliable data for the construction of a regional surface as the data points are densely spaced and the interpretation is based on the most up-to-date data augmented with some vintage data (interpreted seismic sections acquired from 1960 to 2000). This data were extracted from Matoš (2014), Rukavina (2015) and Rukavina, Matoš, Tomljenović, and Saftić (2016).
For subsequent structured surface construction, fault geometry information was extracted from all four category data sources. As previously mentioned, sources differ greatly in resolution and quality, thus the extracted fault geometry was reduced to relatively simple polylines representing the intersection of the fault plane with the Pre-Neogene basement surface.

Input data for the supplemental maps
Main objective of the supplemental maps' construction was to enable outlining of the prospective areas along the pre-Neogene basement unconformity that should be further explored for geoenergy potential. For the construction of supplemental data, a basin modelling process has been applied (Hantschel & Kauerauf, 2009), which requires additional input for subsurface characterization. Crucial input data for such a task consist of stratigraphic and lithological properties of the pre-Neogene basement and Neogene-Quaternary infill, thermal settings of the exploration area (heat flow distribution) and calibration data for the basin modelling.

Workflow
After the preparation of the legacy input data in Arc-Map 10.1 software, modelling was performed as a three-step process. The simple surface construction, structured surface construction and preparation of the basin model were performed in Schlumberger Petrel, while basin model simulation and calibration were performed in Schlumberger PetroMod.

Simple surface construction
For the construction of the simple pre-Neogene surface, only the information about the surface geometry was regarded, without the incorporation of the fault trace information. Data were organized separately according to its origin/source and georeferenced to the same coordinate system, which is in this case WGS84 projected in UTM 33N zone. As the data were gathered from different sources, the spatial data overlapping and differences in the quality/resolution of the input data were emphasized. So, the data had to be harmonized. High-resolution data (first and second category data) were downscaled from modelled surface (Figure 3(a)) firstly to polyline databy extraction of the contour lines with 100 m increment ( Figure  3(b)) with node spacing of 500 m along the polyline and subsequently reduced to point data (Figure 3(c)).
In the marginal parts of the areas covered by different data sources, certain differences occurred, resulting from different input datasets used in original interpretations and the fact that marginal parts are usually characterized by the uncertainty of interpretation due to the lower data density. Since surface morphology became erratic when both datasets in the marginal/ overlapping parts of two neighbouring data sources were included (Figure 4), one input dataset had to be deleted (Figure 4) to achieve more consistent resulting surface. Relevancy of input data for describing the surface morphology in the overlapping area was defined with respect to the data category, under the presumption that the higher category data provides better quality input.
Area representing outcrops of pre-Neogene basement rocks were subtracted from publicly available DEM models and used as an additional input in the surface generation to avoid surface inconsistencies in the marginal part of the PB or near outcrops of pre-Neogene rocks. Area within the polygons that describe the pre-Neogene outcrops were later removed in the final phase of the map construction as they represent the topographic relief rather than the contact of pre-Neogene basement with younger sediments but were necessary for the map construction.

Structured surface construction
After the completion of the simple surface, the next step included the construction of the fault planes from fault traces in order to generate displacement of the surface along the fault planes. As the presentation of pre-Neogene structured surface was the principal aim of the research and mapping is of regional scale, the fault planes geometry was considerably simplified. This means that the fault planes were modelled as subvertical planes from polylines representing the fault intersection with the pre-Neogene surface. Plane geometry was extrapolated vertically so the fault could cover the vertical extent of the modelled surface for computation of the displacement within the structural modelling process. In some instances, the physical aspects of fault intersections/displacements from legacy data could not be honoured due to the subsurface volumes that needed to be displaced. This was only the case in non-model derived data, i.e. 2nd to 4th category data. Once the faults are modelled without volume errors, the simple surface that was constructed in the previous step could now be re-constructed honouring displacements in the near-fault areas, thus resulting in a valid structured surface ( Figure 5).

Basin modelling
As geothermal gradient can differ in vertical scale due to the differences in conductivity of the rocks (Eppelbaum, Kutasov, & Pilchin, 2014), simple estimation of temperature at depth using averaged geothermal gradient values can yield significant errors. In order to determine the subsurface temperature values, 3D basin modelling procedure has been used. 3D basin modelling was also necessary to provide the information about the possible maturity, presented as possible vitrinite reflectance (%Ro) values, of the source rocks situated near the pre-Neogene unconformity. For this process, the subsurface rock formations were grouped into four layers representing four distinct stratigraphic and lithological intervals. The first layer represents the pre-Neogene basement rocks constituted of magmatic and metamorphic rocks (dominantly) or Palaeozoic and Mesozoic sedimentary rocks (subordinate). Neogene infill was subdivided in three layers which represent sediments belonging to Lower/Middle Miocene, Upper Miocene and Pliocene, Pleistocene and Holocene. Thickness of each layer was approximated based on the assumed share of each layer in the thickness of the Neogene-Quaternary sequence. The percentage itself was determined from spatial analysis of thickness distribution of each interval from Saftić et al. (2003). After the input of all parameters, including average lithological composition per layer, heat flow (HF) values, PWD values and sediment-water interference temperature (SWIT, Wygrala, 1989), the model was calibrated based on temperatures obtained from drill-stem-tests (DST), bottom hole static temperatures (BHST) and %Ro values. Three values the temperature on pre-Neogene surface, vitrinite reflectance and surface heat flow corrected after calibration process were extracted from the basin model, while the geothermal gradient was calculated. All of them are generalized and presented as supplemental maps below the Main Map.

Discussion and conclusion
Constructed map shows a structured surface of the pre-Neogene unconformity representing the depth to the base of the interval in which most of northern Croatia's geoenergy potential is located. The Map depicts complex fault systems and associated map-scale structures (i.e. depressions and structural highs) within the CPB that were formed during the Neogene-Quaternary tectonic evolution of the study area. Constructed surface indicates inherited pre-Neogene paleorelief that is the result of Late Cretaceous-Paleogene compressional/ transpressional tectonic features (Csontos & Nagymarosy, 1998;Ratschbacher et al., 1991;Ustaszewski et al., 2010 with references). Majority of the NW-SE and NE-SW striking faults mapped both parallel and perpendicular to CPB depressions are normal faults with listric geometries that have accommodated an E-W extension, i.e. opening and deepening of the CPB depressions during Early to Middle Miocene (Matoš, 2014;Tomljenović & Csontos, 2001 with references). In the same time, NW-SE and NE-SW elongated depressions with depth ≥3000 m (e.g. Drava Depression with max. depth of 7000 m; see Main Map) are mainly bordered by pre-Neogene structural highs, i.e. inselbergs (e.g. Mt. Medvednica, Mt. Moslavačka gora, Mt. Kalnik, etc.) and their subsurface prolongations, which affected sedimentation rates and lithofacies distributions within the CPB.
Fault traces represent important features on the surface as they can be fairways for hydrocarbon and hydrothermal fluid migration from deeper formations or can be boundaries of hydrocarbon traps primarily influenced by near-fault petrophysical properties of rocks. Additionally, they can also be very valuable in determining the possibility of CO 2 storage as larger faults can be related to an active fault area and possible paths of leakage of injected CO 2 .
Supplemental maps are a result of basin modelling and do not represent only the legacy data values as structured surface on the Main map. The three 'modelled maps' presented confirm the regional settings of the studied PB area as having elevated heat flow values (Lenkey et al., 2002), higher than average 65 mW/m 2 estimated for continents (Jaupart, Labrosse, Lucazeau, & Mareschal, 2015;Pollack, Hurter, & Johnson, 1993), resulting in higher geothermal gradient than Earth's average. Although, it is worth noting that, surface heat flow values determined in the calibration process of the 3D basin modelling in certain places showed significantly lower results than those presented in Lenkey et al. (2002) and Békési et al. (2017). Also, estimated values of geothermal gradient were in certain areas lower than those presented in Jelić, Kevrić, and Krasić (1995) and Kurevija et al. (2014).
The first direct consequence of the elevated heat flow is the source rock maturity reached at relatively shallow interval (early oil window onset at only −1600 m) making it a proliferous hydrocarbon exploration area. Modelled Ro values and hydrocarbon maturity windows are in accordance with the up to date discoveries in the CPB. Favourable areas for directing future petroleum explorations can easily be outlined based on the modelled %Ro values presented in the Map.
Supplemental map of averaged geothermal gradient along with the Main map can be used to geographically focus the exploration of geothermal resources. It can be used to delineate areas prospective for exploration of hydrothermal systems as well as areas that are potentially interesting for exploration of enhanced geothermal systems (hot-dry-rocks). As it was previously mentioned, the presented map of averaged geothermal gradient shows certain differences in comparison with previously published maps (Jelić et al., 1995;Kurevija et al., 2014), mainly related to generally lower averaged geothermal gradient values in the Sava Depression. It should be noted that the more prominent difference of averaged geothermal gradient values in Sava and Drava Depression visible in this supplemental map is in accordance with the difference of Earth's crust thickness between the Sava and Drava Depression (after Horváth et al., 2006).
The importance of the Main Map from the aspect of exploration aiming to CO 2 geological storage capacity estimation is manifested mainly in the definition of depth of pre-Neogene basement unconformity as well as identification of deep discontinuities that represent potential leakage paths for CO 2 from geological storage objects. So far, regional exploration of CO 2 storage capacity estimates has been mainly focused on Upper Miocene sandstones as potential regional storage formations, while storage in fractured basement reservoirs has not been assessed. Furthermore, estimates were based on interpretation of well data to enable correlation and mapping of storage unit thickness (e.g. Kolenković, Saftić, & Perešin, 2013) and this approach partly disregarded structural analysis. In that respect, Main map can be used to direct further exploration to areas where pre-Neogene basement unconformity is situated shallower than 2500 m and where faults are not densely situated.
Presented Map is a most detailed publicly available structured surface of the Base Neogene unconformity covering the CPB. Constructed maps are important for planning of regional research of the geoenergy potential, including conventional and unconventional hydrocarbon resources, geothermal energy, CO 2 geological storage capacity and of the subsurface energy storage potential as well. Implications of developing this research are two-fold. Firstly, northern Croatia is a mature petroleum province with many oil and gas fields discovered and partly already depleted. This means that either the quest for the by-passed oil in the existing exploitation blocks is inevitable (which is ongoing) or one should also consider if some regional elements of petroleum systems have been omitted or underexplored, which would render additional HC resources, both conventional and unconventional. Secondly, the already developed subsurface exploration and results offer the opportunity for using this as an impetus to more quickly develop emerging resources that are becoming increasingly important in the light of energy transition to renewables or needs for decarbonization of both the energy and industrial sector. Regarding the renewables, deep geothermal potential is a non-intermittent source and as such is highly desirable in the energy mix which means that if Croatia will be able to speed up the development of this resource it is bound to have a significant impact on the national economy.

Software
For the preparation of the data for the map construction, ArcGIS ArcMap 10.1 and Schlumberger Petrel 2017 software were used. Construction of the Main Map and maps with supplemental data were performed using Schlumberger Petrel 2017 and PetroMod 2017 software. Final map layout was drafted in CorelDraw X4 software.