Geological modeling of a tectonically controlled hydrothermal system in the southwestern part of the Pannonian basin (Croatia)

Kosović, Ivan; Matoš, Bojan; Pavičić, Ivica; Pola, Marco; Mileusnić, Morena; Pavić, Mirja; Borović, Staša

doi:10.3389/feart.2024.1401935

ORIGINAL RESEARCH article

Front. Earth Sci., 07 June 2024

Sec. Structural Geology and Tectonics

Volume 12 - 2024 | https://doi.org/10.3389/feart.2024.1401935

Geological modeling of a tectonically controlled hydrothermal system in the southwestern part of the Pannonian basin (Croatia)

Ivan Kosović¹

Mirja Pavić¹

Staša Borović¹

¹Croatian Geological Survey, Zagreb, Croatia
²Faculty of Mining, Geology and Petroleum Engineering, Department of Geology and Geological Engineering, University of Zagreb, Zagreb, Croatia
³State Geodetic Administration, Zagreb, Croatia

Geothermal energy is an important resource in the green economy transition. For the preservation of a geothermal resource it is crucial to assess its renewability and the sustainability of the exploitation. These aspects are influenced by the interaction among the physical, chemical, geological, and hydrogeological processes. The reconstruction of the geological assemblage allows the detailing of the geometries of the reservoir and fracture systems that influence the fluid flow and the water/rock interaction. The control of regional/local scale fault and fold systems on the development of the Daruvar hydrothermal system (DHS), located in Croatian part of the Pannonian basin, is detailed in this work. Field investigations were conducted to collect structural data on strata orientation and fault/fracture systems. The dataset was integrated with geological and geophysical data to develop composite geological profiles and a 3D geological model. Results display a pattern of generally N-S and E-W striking folds and cogenetic fracture systems with orientations parallel to the fold axes. The geological reconstruction was integrated with geophysical, hydrogeological, and geochemical data to propose a conceptual model of the DHS. The DHS is a topographically driven system hosted in a Mesozoic carbonate reservoir where E-W striking fracture systems are regional flow paths that enable infiltration of meteoric water to 1 km depth and its reheating in its reservoir area. In Daruvar, an anticline and fault/fracture systems accommodate the uplift of reservoir to shallow depths, promoting the bedrock fracturing and increase of the permeability field. These conditions favor the localized upwelling of thermal water resulting in four thermal springs (38°C and 50°C) in Daruvar city area. This work highlights the importance of employing a multidisciplinary approach to detail the complex interaction among the processes driving the geothermal resource.

1 Introduction

The sustainable management of natural resources is one of the most important challenges in the 21st century (ECE, 2021). Natural resources developed by geological processes are very important since they support many industrial activities. Geothermal resources have a pivotal role in the current global economy since they are potential renewable sources of both raw materials and energy (e.g., Finster et al., 2015; Szanyi et al., 2023). The sustainable utilization of geothermal resources is crucial since their development and renewability depend on a delicate balance between physical and chemical processes (e.g., Rybach and Mongillo, 2006; Axelsson, 2010; Rman, 2014; Shortall et al., 2015; Fabbri et al., 2017).

The most relevant processes affecting the characteristics of a geothermal resource are the conduction of heat and the convection of fluids (Moeck, 2014). The magnitude of these processes mostly depends on the geological and hydrogeological settings of the geothermal system associated with the resource (Kühn and Gessner, 2009; Bundschuh and César Suárez A., 2010; Pasquale et al., 2014). The conductive component of the heat transfer is mostly related to the regional subsurface geological setting that affects: i) the distribution of lithologies and their thermal properties, and ii) the occurrence of deep structures favoring an increased heat flow from the deeper part of the Earth’s crust. Convection encompasses the transfer of mass and heat occurring by bulk fluid motion and depends on the subsurface permeability field (Ataie-Ashtiani et al., 2018). In bedrock aquifers, the original permeability can be enhanced by the fracturing with fault damage zones being preferential paths for the fluid flow (e.g., Faulkner et al., 2010; Bense et al., 2013). In particular, thermal springs are generally associated with systems of faults that enable or enhance the outflow of thermal waters (e.g., Curewitz and Karson, 1997; Nelson et al., 2009; Pola et al., 2014; Keegan-Treloar et al., 2022; Pavić et al., 2023). Therefore, a detailed reconstruction of the geological, tectonic, and hydrogeological settings is crucial to determine the processes favoring the circulation of thermal waters and influencing the renewability of the geothermal resource and its exploitation sustainability (Magri et al., 2010; Faulds et al., 2013; Scheck-Wenderoth et al., 2014; Brehme et al., 2016; Pola et al., 2020; Torresan et al., 2021).

The reconstruction of the geo-tectonic settings of a certain area is generally conducted by integration of field investigations and geophysical data, at regional and local scales. While field investigations enable the reconstruction of the surficial geometry of geological formations and fractures and the assessment of the kinematics of the principal faults, geophysical data can add crucial subsurface information to the geological reconstruction. Geological and geophysical datasets can be integrated to construct a 3D geological model of the subsurface (Pavičić et al., 2018; Olierook, 2020; Jia et al., 2021; Panzera et al., 2022). 3D modeling, as a backbone in geological and hydrogeological applications, provides a useful tool for the interpretation and visualization of geological structures, especially when their geometrical complexity cannot be fully represented through 2D sections (Caumon et al., 2009; Pavičić et al., 2018; Wellmann and Caumon, 2018; Pan et al., 2020; Lyu et al., 2021). Assessing the tectonic setting and the stress regimes can improve the geological reconstruction. It permits understanding the kinematics of the local fracture networks that influence the permeability field in the reservoir (Wang et al., 2014; Santilano et al., 2016; Xie et al., 2017; Li Y. et al., 2018; Price et al., 2018; Pan et al., 2020). Furthermore, it is crucial to quantify the hydrogeological properties of the thermal aquifer that can be used to forecast the quantity of exploitable waters. These reconstructions are fundamental for a detailed hydrogeological modeling, which explains the processes driving the formation of a thermal resource (Moeck et al., 2010; Calcagno et al., 2014; Fulignati et al., 2014; Mroczek et al., 2016; Montanari et al., 2017; Torresan et al., 2020).

According to Borović and Marković (2015), Northern Croatia is rich in geothermal resources sharing the favorable thermal properties of the Pannonian basin area (Horváth et al., 2015). Most of the thermal springs in Croatia are used for balneotherapy and tourism, however, they also have a great potential for additional utilization (e.g., district heating, industrial processes). One of the most investigated thermal regions in Croatia is the Daruvar city area (Figure 1), with thermal springs documented since the Roman era. Systematic geological, hydrogeological, and geophysical investigations in the Daruvar area have been conducted since the 1970s (e.g., Babić et al., 1971; Mraz, 1983; Larva and Mraz, 2008; Borović, 2015; Borović et al., 2019; Kosović et al., 2023; Urumović et al., 2023). The thermal spring area is the outflow area of an intermediate scale, tectonically-controlled, topographically driven, geothermal system hosted in a Mesozoic carbonate rock complex. The thermal waters are exploited from a thermal well and two springs that provide approximately 10 l/s.

Figure 1

Figure 1. Regional map (modified after Schmid et al., 2008; 2020) showing: i) the main tectonic units, ii) the regional fault structures, and iii) the border line of the Pannonian basin (white line). The study area is located in the Tisza-Dacia Mega Unit at the contact with the Sava Suture Zone and is a part of the Pannonian basin. The red rectangle marks the extent of the structural map in Figure 2. A schematic geographic map is reported in the upper left corner showing Daruvar, Zagreb (the capital city of Croatia), Croatia, and the neighboring states (acronyms: AT–Austria, BA–Bosnia and Herzegovina, HR–Croatia, HU–Hungary, IT–Italy, Sl–Slovenia, RS–Serbia).

Despite several multidisciplinary investigations in Daruvar, the detailed quantitative geological and structural reconstructions of the recharge area and flow-through parts of the Daruvar hydrothermal system (DHS) are still unreliable and outdated. This study aims to present results of field investigations conducted in the hinterland area of Daruvar to detail its geological setting, structural framework, and tectonic evolution. Geological and structural data were combined with available seismic reflection data and stratigraphic logs of deep wells to perform a 3D geological model. Results were furthermore used to refine the hydrogeological conceptual model of the DHS focusing on the correlation of the geological and tectonic models with the preferential flow paths in the system.

2 Geological and hydrogeological settings

2.1 Regional tectonics evolution

The DHS area is located in the western part of the Slavonian mountains (Mount Papuk), which are one of the best exposures of the Tisza-Dacia Mega-Unit, a lithospheric fragment formed between the European and Adria plates during the Middle Jurassic (Figure 1; e.g.; Balen et al., 2006; Schmid et al., 2008 with references). This area experienced a complex tectono-metamorphic evolution that started with the Variscan and continued through the Alpine-Dinarides-Carpathian orogeny. The Variscan events were characterized by structural stacking of the Mecsek, Bihor, and Codru nappe systems (Figure 1), with various metamorphic overprints of the preexisting crystalline basement, intruded by granitoids and migmatites (Balen et al., 2006 with references). According to Schmid et al. (2008, 2020), Bihor nappe system exposures can be found in southern Hungary and the Slavonian mountains. Balen et al. (2006) suggested that metamorphic and plutonic basement rocks of the Slavonian mountains were locally covered by a Paleozoic-Mesozoic carbonate/clastic succession deposited at the passive margin of the Adria plate. Neogene-Quaternary deposits further concealed the older rock sequence after the intensive Cretaceous-Paleogene tectonism (Balen et al., 2006 with references). As a part of the Alpine-Dinarides-Carpathian orogenic system, the study area is positioned in the vicinity of the Sava Zone, i.e., the Cretaceous-Paleogene suture zone between the Tisza-Dacia Mega-Unit in the NE and the Adria plate in the SW. The Sava Zone is composed of complex assemblage of ophiolitic, magmatic-metamorphic, and deep-water sedimentary rocks formed in the Neotethys and the Sava Ocean (Schmid et al., 2008 with references). These rock complexes are often displaced by regional fault zones (i.e., Mid-Hungarian fault line), which are characterized by a polyphase tectonic evolution (Figure 1; Figure 2; Schmid et al., 2008).

Figure 2

Figure 2. Simplified structural map of the regional tectonic framework with regional fault systems within northern Croatia. The yellow rectangle corresponds to the extent of the geological map in Figure 3. Fault systems are compiled after: Jamičić, 1995; Prelogović et al., 1998; Tomljenović and Csontos, 2001; Schmid et al., 2008; Šolaja, 2010; Matoš et al., 2016; Wacha et al., 2018; Baize et al., 2022; Herak and Herak, 2023.

The Neogene-Quaternary tectonic evolution of the Pannonian Basin System (PBS), which is characterized by repeated extension, compression, and tectonic inversion, further affected the structural assemblage of the study area. In particular, it was conditioned by the Adria-Europe collision, the eastward lateral-extrusion of the continental blocks between these plates, and the clockwise rotation of the Tisza-Dacia Mega-Unit (Prelogović et al., 1998; Tari et al., 1999; Csontos and Vörös, 2004; Schmid et al., 2008).

PBS was formed by Early to Late Miocene (c. 26–11.5 Ma) NNE-SSW directed “back-arc type” lithospheric extension along the NNW-striking listric normal faults (Figure 2; Fodor et al., 2005; Horváth et al., 2006; Schmid et al., 2008; Brückl et al., 2010). Rift and wrench-related troughs were filled with large amounts of syn-rift deposits (Tari and Pamić, 1998; Horváth and Tari, 1999; Steininger and Wessely, 2000; Ustaszewski et al., 2010). In the Croatian part of the PBS, deposition commenced along the listric faults forming the basins and subbasins (Pavelić et al., 2001; Ćorić et al., 2009). Though local structural tectonic inversion occurred at the end of the Middle Miocene (c. 13.0–11.6 Ma), Croatian part of PBS was characterized by continuous deepening and rapid thermal subsidence along the existing faults until the Late Miocene - Early Pliocene (c. 11.5–5.3 Ma; Csontos et al., 1992; Horváth and Tari, 1999; Tomljenović and Csontos, 2001; Fodor et al., 2005; Malvić and Velić, 2011). A significant change in the stress field, with N-S trending compressional and/or transpressional P-axes, occurred during the Pliocene. Translation and counterclockwise rotation of the Adria plate in combination with consumption of the subducted European plate lithosphere led to regional tectonic inversion (Horváth and Tari, 1999; Grenerczy et al., 2005; Dolton, 2006; Jarosiński et al., 2006; Bada et al., 2007; Jarosinski et al., 2011; Ustaszewski et al., 2014). The Pliocene-Quaternary tectonic inversion accommodated large-scale lithospheric folding, vertical/horizontal motions along the existing regional faults, and horizontal/vertical displacement along the co-genetic reverse and strike-slip faults (e.g., Periadriatic fault, Mid-Hungarian fault line, Sava and Drava fault zones; Figure 2; Horváth and Cloetingh, 1996; Prelogović et al., 1998; Horváth and Tari, 1999; Fodor et al., 2005; Dolton, 2006; Bada et al., 2007; Jarosinski et al., 2011). In the Croatian part of the PBS, the tectonic uplift yielded final uplift of pre-Neogene basement highs (e.g., Slavonian mountains), which caused tectonic overprint of basement structures, block rotations, and formation of the positive flower structures with kilometer-scale folds along the reactivated and newly formed strike-slip faults (Figure 2; Jamičić, 1995; Prelogović et al., 1998; Tomljenović and Csontos, 2001; Balen et al., 2006).

2.2 Geological setting

The geological setting of the Daruvar area was extensively investigated by Jamičić et al. (1989). Since the scope of this work is the geological reconstruction of the DHS for a detailed hydrogeological conceptual modeling, the original geological map (Supplementary Figure S1) was simplified considering the hydrogeological properties of the lithological units together with their age (Figure 3). The units (Supplementary Figure S2) were reorganized as follows: i) pre-Permian crystalline rocks, ii) Permian sedimentary units, iii) Triassic carbonate rock complex, iv) Jurassic limestones, v) Miocene sedimentary and magmatic rock complex, vi) Pliocene clastic sediments, vii) Pleistocene unconsolidated sediments, and viii) Holocene alluvial and colluvial unconsolidated sediments.

Figure 3

Figure 3. Simplified geological map of the Daruvar area (modified from Jamičić et al., 1989; Šolaja, 2010). Fault acronyms: DF, Daruvar fault; DeF, Dežanovac fault; GF, Gradina fault; PF, Pakrac fault; TF, Toplica fault; VF, Voćin fault. Topographic peaks are denoted by white triangles. White polygons indicate larger settlements and towns.

The oldest rocks in the study area are the pre-Permian crystalline rocks. They cover the largest area in western Papuk and are composed of migmatites, granitoids, pegmatites, gneisses, and chlorite schists. Granitoids are the most common lithology and they are S-type granites concordant with migmatite bodies (Jamičić et al., 1989; Pamić et al., 2003). The crystalline rocks are generally in transgressive contact with Permian or Miocene units. Permian rocks are composed of well-layered conglomerates and quartz sandstones. The conglomerates contain clasts with a variable lithological composition depending on the underlying basement. Locally, these sediments can show a low-grade metamorphism. The thickness of the Permian unit is approximately 400 m. The Triassic sedimentary rock complex was continuously sedimented over the Permian deposits. The formation of clastic deposits prevailed during the Lower Triassic, while shallow-water carbonate sedimentation with occasional clastic sediment deposition occurred in the Middle and Upper Triassic (Šikić, 1981; Jamičić et al., 1989). The result is a sedimentary complex composed of: i) Lower Triassic sandstones, siltstones, and laminated shales, ii) Middle Triassic dolomites, limestones, and crinoid limestones with chert, and iii) Upper Triassic dolomites and limestones. The total thickness of the Triassic unit is approximately 500 m. During the Jurassic, western Papuk as the contact zone of the Adria plate and the Tisza block was characterized by deep-sea basin sedimentation. Jurassic deposits are preserved exclusively in western Papuk (Jamičić, 2009). They are represented by platy limestones with cherts (Jamičić et al., 1989), with thickness approximately 100 m. From the Lower Cretaceous, western Papuk experienced tectonic uplift and emersion. It was characterized by coastal environments with frequent sea-level oscillations and alternations of marine, brackish, and freshwater sedimentation. The sedimentation was partly restored in the Middle Miocene due to the PBS E-W extension and regional transgression. Sea level oscillations and alternations in the Neogene-Quaternary deposition environments resulted in marine, brackish, and freshwater sedimentation in the structural lows of the previously formed structures. Miocene sediments are mostly composed of conglomerates, sandstones, marls, marly and bioclastic limestones, and loose clayey-sandy sediments. Furthermore, pyroclastic and effusive rocks can be found, with andesites occurring in the western part of the study area. The thickness of the Miocene succession is in the range of 600–650 m. The clastic sedimentation continued in the Pliocene with sandstones and marls in different proportions, as well as sands and gravels, with a total thickness between 700 and 900 m. The youngest unconsolidated Quaternary sediments can often be found as a “transgressive” cover of the older units. Quaternary sediments consist of sandy gravels, quartz sands, silty sands, and sandy clays (Jamičić et al., 1989). Pleistocene deposits are characterized by alluvial and loess-like deposits, while alluvial and slope sediments occur during the Holocene. The Pleistocene and Holocene units are up to 25 and 5 m thick, respectively.

The youngest Plio-Quaternary tectonic phase significantly affected the structural assemblage of the study area. This tectonic activity is associated with the proximal compressional/transpressional stresses accommodated along the Drava and Sava fault zones (Figure 2) which led to both the formation of folds, new faults, and the structural reactivation with local inversion of inherited structures (Figure 3; Jamičić, 1995). The principal mapped faults in the study area are mostly E-W, NE-SW, and NW-SE striking faults with cogenetic N-S, NW-SE, and NE-SW striking fold axes (Jamičić et al., 1989; Jamičić, 1995; Šolaja, 2010).

2.3 Hydrogeological setting

Four thermal springs with temperatures between 38°C and 50°C occur in Daruvar. Furthermore, approximately 100 shallow boreholes and a few deep wells have been drilled since the 1970s. The lithostratigraphic sequence was detailed through two wells deeper than 100 m (Kosović et al., 2023). It consists of: i) Quaternary alluvial deposits with a thickness of up to 20 m, ii) Miocene or Pliocene marls up to 30 m thick, iii) Miocene bioclastic limestone with a thickness of 10 m, and iv) Triassic dolomites and limestones with thickness of 130 m. Where the thermal springs occur, the stratigraphic logs show that the Triassic carbonates are in direct contact with the Quaternary cover. The Triassic carbonates are moderately to highly fractured. This formation represents the primary thermal aquifer, while secondary, colder, thermal aquifers are found in the sandy layers of the alluvial cover and the Miocene biocalcarenite (Borović et al., 2019). The transmissivity of the Triassic carbonate aquifer was assessed through pumping and well tests ranging from 0.015 to 0.03 m²/s (Borović et al., 2019; Urumović et al., 2023). The main physico-chemical characteristics of thermal waters in springs and wells (Borović, 2015) are: i) temperature from 18.2°C to 49.8°C, ii) nearly neutral pH with values between 6.7 and 7.8, iii) EC between 550 and 700 μS/cm, and iv) calcium-bicarbonate hydrochemical facies. O and H stable isotope ratios suggest a meteoric origin of the Daruvar waters, while ¹⁴C activity points to a mean residence time between 11 and 15 ka (Borović, 2015).

These data were used to propose an initial conceptual model of the DHS (Borović et al., 2019). The recharge area of the system is located in the topographically highest part of the eastern hinterland of Daruvar (i.e., western Papuk and Petrov vrh area) where the Triassic carbonates are uplifted by a regional fault system (Figure 3). The deep infiltration of the meteoric waters is favored by the dip of the layers, the karstification, and the fracturing of the rock mass. The low permeable Permian and pre-Permian units at the base act as a barrier for further downward circulation. Locally, shallow cold and deep warm waters mix developing hypothermal and subthermal springs in the Daruvar hinterland. In the Daruvar area, existing regional thrust and strike-slip faults (Figure 3) and their polyphase tectonic history enhanced the deformation of the aquifer and its permeability field. Furthermore, faults juxtapose the aquifer with low permeable Miocene-Pliocene deposits forming a barrier for the fluid flow toward W.

3 Materials and methods

3.1 Structural-geological investigation and analysis

Structural-geological field investigation was carried out from 2021 to 2024. Data were collected at 305 field points, digitized, and spatially georeferenced using GIS software. The structural investigation included measurements of the geometrical properties of strata bedding and fractures/joints (i.e., dip direction and dip angle). Furthermore, outcrop-scale structural data on fault geometrical properties, fault slip, and kinematics were also collected. The resulting dataset comprised 134, 659, and 91 measurements of bedding, fractures, and faults, respectively. Where local scale folds were visible, the fold axis was determined. Following the principal structural units, the Daruvar structural framework (Figure 4) was divided into three structural domains separated by regional E-W or NE-SW striking reverse faults: i) Dujanova kosa (DK), ii) Petrov vrh (PV), and iii) Sirač (SI). Structural observations in these domains were subdivided into eastern, central, and western sections obtaining a detailed analysis of the structural style within each block composing the structural fabric of the Daruvar hinterland.

Figure 4

Figure 4. Distribution of the field measurement points (white dots)¸ and traces of constructed composite geological profiles (red line) in the Daruvar area. Measured data, coupled with seismic reflection profiles (light yellow lines) and stratigraphic logs of deep wells, were used to construct the 3D model (blue rectangle). The structural domains (Dujanova kosa, Petrov vrh, and Sirač; orange polygons), the main geological structures, and towns (light grey polygons) are also shown. The acronyms of the main faults are reported in Figure 3.

In this study, data were plotted by the Stereonet v.11 software (Allmendinger et al., 2011; Cardozo and Allmendinger, 2013). A representative bedding was calculated for every measurement point. The results were plotted to graphically determine the most common orientations within each domain. Poles of fracture planes were used to construct contour plots of the poles distribution using the 1% area contour method (Allmendinger et al., 2011; Cardozo and Allmendinger, 2013). The plots were used to determine the most representative sets of fractures in each domain. For the determination of fault kinematics, we used data of dip direction and dip angle of fault planes and kinematic indicators, i.e., orientation of slickensides defined by azimuth and plunge, and, their sense of movement (Doblas, 1998). Based on kinematic criteria, the shear fracture data were analyzed by Win-Tensor v. 5.9.2 software (Delvaux and Sperner, 2003). The obtained data was separated into compatible fault groups and processed by TectonicsFP v. 1.7.9 software (Ortner et al., 2002). Theoretical maximum (σ1), mean (σ2) and minimum stress axes (σ3) were calculated using the P–T axis method (Turner, 1953; Marrett and Allmendinger, 1990). For the analyzed fault groups, synthetic focal mechanisms were calculated using the Right Dihedra Method (Angelier and Mechler, 1977).

3.2 Composite geological profiles and 3D geological modeling

Composite geological profiles are key components in the interpretation of the 2D/3D subsurface relations and structures. The composite geological profiles were constructed using: i) field data collected within this study, ii) existing geological data (i.e., geological maps, explanatory notes, published geological data), iii) seismic reflection profiles collected for hydrocarbon exploration, and iv) stratigraphic logs of deep wells drilled during the hydrocarbon exploration campaign and geophysical well-logging data. Here, five NW-SE composite geological profiles were constructed (DHS-1 to DHS-5; Figure 4), representing surface/subsurface 2D models of the Lonjsko-Ilovska depression and the western part of Mount Papuk.

Constructed composite geological profiles combined with field dataset were further integrated to develop the 3D geological model of the DHS area (blue rectangle in Figure 4). Petroleum Experts Move 2019.1 (https://www.petex.com/products/move-suite/move/) software package was used to build the subsurface model. The workflow for the model construction is shown in Figure 5. Since the geological setting of the study area is very complex (i.e., fault with variable architecture, regional folds, disconformities, angular unconformities, nonconformities), 15 additional geological sections mostly perpendicular and longitudinal to the geological structures were constructed. Furthermore, 20 smaller auxiliary (temporary) sections were made to obtain a more detailed reconstruction at the local scale. Fault surfaces were obtained by extending fault traces from geological profiles and maps into the subsurface using detailed measurements from the fieldwork and fitting the geological sections. Horizon surfaces were constructed by the ordinary kriging interpolation algorithm with “Use Meshed Alpha Shape as Input Points” option. This option allows a better interpolation of the complex surfaces based on sparse and irregular data. A construction mesh with vertices exactly at the XY location of irregularly spaced data points was performed (Petroleum Expert, 2019). The algorithm then used a convex hull to create a surface. Grid geometry provided an option (i.e., Honour Points) to create an interpolated surface with Z elevations at the XY locations of existing data points and is geomathematically predicted at other vertex locations of the construction mesh (Petroleum Expert, 2019). Surface Sampling controls the triangle size of the mesh, and it was fitted based on the density of input data.

Figure 5

Figure 5. Schematic flow chart used in the construction of the 3D geological model of the DHS area.

4 Results

4.1 Analysis of bedding and fracture system

4.1.1 Dujanova kosa (DK)

The structural domain of Dujanova kosa (DK) encompasses the northern segment of the study area (Figure 4). It is positioned between two E-W striking low-angle thrust faults (Voćin and Gradina faults; VF and GF, respectively, in Figure 3; Figure 4). Structural measurements were conducted on 68 geological stations (Figure 4) with collected data about 48 strata orientations and 158 fracture planes.

Field observations of the basement section evidenced that the granitoids are heavily weathered. They are structurally missing geometrical properties, i.e., foliations/bedding orientations, but they are deformed by two systems of fractures (Table 1; Figure 8). The system Fs₁ is characterized by steeply dipping (dip angle 80°) NNW-SSE striking fractures, whereas Fs₂ is composed of discontinuities that dip towards SSE at a dip angle of 52°.

Table 1

Table 1. Average orientations of the main fracture systems in the structural domains of the study area and number of data constituting the sets.

Figure 8

Figure 8. Structural diagrams of the main fracture systems (black circles) in the structural domains of the study area. The orientations of the representative fractures (F_S) in the systems are reported in Table 1. They are obtained from the distribution of the measured discontinuities represented here as poles (black points). The distribution is shown as a contour plot using a rainbow color ramp with values depending on the number of available data (N).

The eastern section of the DK encompasses the Permian clastic and the Triassic carbonate succession. It is mostly characterized by folded structures with fold axes gently dipping towards SW or NW (dip angle of 12° and 5°, respectively; Figure 6). In the first fold set beds are gently dipping SE and WNW at average angles of respectively 48° and 37°, or towards ENE and WSW at 28° and 52° for the second fold set. These observations indicate a polyphase tectonic evolution of the Permian-Triassic succession that enabled the formation of at least two generations of folds. Four fracture systems were further evidenced (Table 1; Figure 8). The Fs₁ fracture system is characterized by discontinuities steeply dipping towards N with a dip angle of 71°. The system Fs₂ includes predominantly N-S striking fractures that are gently dipping towards E at 36°. Fracture system Fs₃ shows a strike parallel to Fs₂, but these discontinuities are dipping towards W at the angle of 58°. The Fs₄ system is characterized by subvertical NE-SW striking fractures that are dipping at the angle of 76° towards SE.

Figure 6

Figure 6. Stereonets of the bedding orientations (thin black circles) in the different domains of the study area. Their interpretation highlights the polyphase tectonic evolution of the study area with processes of refolding and occurrence of two generations of folds (thick blue and red circles) characterized by different fold axes (blue and red points). N indicates the number of measurements available in each structural domain.

The structural measurements in the central section of the DK were mainly conducted in Upper Triassic and Jurassic carbonates (Figure 3; Figure 4). Beds are dipping either towards NE or W at the average angle of 43° and 64°, respectively, resulting in folded structures with fold axes gently dipping towards NNW at the angle of 15° (Figure 6). Fracture data suggest the existence of three fracture systems (Table 1; Figure 8). The system Fs₁ is characterized by fractures steeply dipping towards WSW with an average dip angle of 63°. The Fs₂ set is composed of NE-SW striking discontinuities that are steeply dipping towards SE at the angle of 73°. Fs₃ has a strike similar to Fs₁, but these discontinuities are dipping towards NE at 61°. All these fracture systems are characterized by subvertical geometry implying their structural position with respect to the fold hinge zones.

The western section of the DK includes field investigations that were conducted in Triassic carbonates and its Neogene sedimentary cover. Bedding measurements show similar structural pattern to the central section of the DK. NW-SE striking beds dipping towards NE and SW (dip angle of 81° and 63°, respectively) compose isoclinal folded structures with axes gently dipping towards NW, at the angle of 19° (Figure 6). Measured fracture discontinuities suggest the presence of two systems (Table 1; Figure 8). The Fs₁ fracture system is characterized by E-W striking discontinuities dipping towards N with a dip angle of 53°. The Fs₂ fracture system includes NNE-SSW striking discontinuities dipping towards ESE at the angle of 59°.

4.1.2 Petrov vrh (PV)

The PV structural domain includes the central segment of the study area (Figure 4). It is bounded by the E-W striking Gradina fault (GF) and the NW-SE striking Pakrac fault (PF). Here, the rock complex has been intensively folded by a series of kilometer-scale displaced anticlines and synclines. Structural measurements were conducted on 184 geological stations (Figure 4) with collected data about 57 strata orientations and 332 fracture planes.

As observed in DK, the pre-Permian granitoids of the basement are structurally missing recognizable foliation/bedding orientations. However, three systems of fractures (Table 1; Figure 8) were observed. Fracture system Fs₁ is characterized by discontinuities that dip towards NE at an average dip angle of 45°. The system Fs₂ is composed of discontinuities that dip towards W at a dip angle of 68°. The group Fs₃ is a WNW-ESE striking system of fractures that are steeply dipping (dip angle 87°) towards SSW.

The eastern section of the PV, which includes the Permian clastic and Triassic carbonate successions, is mostly characterized by folded structures with fold axes gently dipping towards S at an angle of 4° (Figure 6). Field measurements in these folded clastic and carbonate layers indicate that beds are dipping either E or W at the average angle of 48° and 31°, respectively. Three fracture systems were determined (Table 1; Figure 8). The first fracture system (Fs₁) is characterized by discontinuities steeply dipping towards N (dip angle of 88°). The Fs₂ system includes predominantly NNE-SSW striking discontinuities that are dipping towards ESE at the angle of 50°. The fracture system Fs₃ has a NNW-SSE strike and these discontinuities are dipping towards ENE at 61°.

The structural measurements in the central section of the PV, similarly to the central section of DK, were mainly conducted in fractured Upper Triassic-Jurassic carbonates (Figure 7A). Folded structures are characterized by fold axes gently dipping either towards SSE or SW (dip angle of 21° and 25°, respectively; Figure 6). Beds of the first group of folds are steeply dipping towards E and WSW at the average angle of 66° and 80°. For the second fold series, beds are gently dipping towards W and S at the angle of 33° and 32°, respectively. These observations indicate a polyphase tectonic evolution in the domain that enabled the formation of at least two generations of folds. At the same time, two main fracture systems can be delineated (Table 1; Figure 8). The system Fs₁ is characterized with fractures steeply dipping towards ESE at a dip angle of 82°. The second fracture system Fs₂ is characterized by WNW-ESE striking discontinuities that are steeply dipping towards SSW at the angle of 87°.

Figure 7

Figure 7. (A) Heavily fractured Triassic dolomites intercalated with dolomitized limestones and shales at the quarry of Batinjska rĳeka steeply dipping toward WNW (central part of PV structural domain, NE of Daruvar; 45°36′20.76″N, 17°16′19.77″E). (B) Permian layered sandstones (eastern section of the SI domain; 45°33′14.52″N, 17°21′38.74″E) with bedding orientation of 237/32 (dip direction/dip angle) deformed by two sets of N-S and E-W striking subvertical fractures system (average dip direction/dip angle of 123/64 and 354/70, respectively).

The western section of the PV encompasses the Triassic carbonates and their transgressive Neogene sedimentary cover (Figure 3; Figure 4). Similarly, to the central section, bedding measurements indicate a polyphase tectonics, with at least two generations of folds. Folded structures are characterized by fold axes gently dipping towards SSW or ENE (dip angle of 23° and 14°, respectively; Figure 6). Beds in the first fold group are dipping either ESE or W at the angle of 60° and 65°, respectively. For the second generation of folds, they are dipping towards SE and NNW at the angle of 26° and 53°, respectively. The analyses of the measured fracture discontinuities (Table 1; Figure 8) suggested the presence of three fracture systems. The Fs₁ fracture system includes NNE-SSW striking discontinuities that are dipping towards ESE with a dip angle of 52°. The Fs₂ system is dipping steeply towards N (average dip angle of 84°), while fracture system Fs₃ has a strike parallel to the strike of Fs₁, but these discontinuities are steeply dipping towards WNW at the angle of 82°.

4.1.3 Sirač (SI)

The SI structural domain encompasses the southern part of the study area and is mostly composed of Mesozoic carbonate succession (Figure 3). Its northern boundary follows the NW-SE striking low-angle Pakrac thrust (PF; Figure 3). Field investigations at the 58 geological stations resulted in measurements of 29 strata orientations and 243 fracture planes.

The basement section of the SI refers to a small area located in the farthest NE part. The structural fabric of the pre-Permian granitoids is similar to the DK and PV domains being deformed by two fracture systems (Table 1; Figure 8). Fs₁ is an NNW-SSE striking fracture system dipping towards WSW (dip angle of 44°), whereas the second fracture system, Fs₂, includes discontinuities that gently dip towards NNE at a dip angle of 13°.

The eastern section of the SI encompasses the Permian sedimentary units (Figure 7B) and the Triassic carbonate succession that are intensively fractured and folded. Folded structures are characterized by fold axes gently dipping towards S or W (dip angle of 10° and 4°, respectively; Figure 6). Beds of the first fold generation are dipping either E or W at the angle of 69° and 32°, or towards N or S (dip angle of 62°and 53°, respectively) for the second group of folds. As observed in other domains, it suggests a polyphase tectonic evolution that enabled at least two deformation events, i.e., fold generations. Measured fracture discontinuities pinpoint the existence of at least three fracture systems (Table 1; Figure 8). The first fracture system Fs₁ is characterized by discontinuities steeply dipping towards S with a dip angle of 72°. The second set Fs₂ includes predominantly NNW-SSE striking discontinuities that are dipping towards WSW (dip angle of 49°). Fs₃ fracture system is characterized by ENE-WSW striking discontinuities that are dipping at the angle of 30° towards SSE.

The structural measurements in the central section of the SI were mainly conducted in Upper Triassic dolomites (Figure 3). NNW-SSE striking beds (Figure 6) point to folded structures with fold axes gently dipping towards SSE at the angle of 2°. Beds in this structures are dipping towards either NE or SW at the angle of 40°, and 60°, respectively. Another set of measurements indicates a homocline where beds are gently dipping towards NNW at the angle of 33°. Three fracture systems were delineated (Table 1; Figure 8). The first fracture system Fs₁ is characterized by fractures steeply dipping towards ENE with a dip angle of 84°. The second set Fs₂ encompasses NE-SW striking discontinuities that are steeply dipping towards SE at the angle of 61°. Fs₃ group has a WNW-ESE strike and dips towards SSW at the angle of 60°. These fracture systems are characterized by subvertical geometry that implies their structural position with respect to the fold hinge zones.

The western section of the SI included field investigations that were mainly conducted in Triassic carbonates (Figure 3). NE-SW striking beds compose folded structures with axes gently dipping towards SW at the angle of 2° (Figure 6). Measured fracture discontinuities (Table 1; Figure 8) suggested the presence of three fracture systems. The system Fs₁ is characterized by ESE-WNW striking discontinuities that are dipping towards SSW with a dip angle of 68°. The Fs₂ fracture system includes NE-SW striking discontinuities dipping towards SE at the angle of 65°. Fs₃ is characterized by subvertical N-S striking discontinuities that are steeply dipping at the angle of 77° towards E.

4.2 Fault system and shear fracture analysis

Structural investigations focused also on the identification of the principal faults/fault zones that built the structural assemblage of the study area. Approximately 90 shear fractures/fault plane data at 47 geological stations were collected. Considering the stress field and kinematic criteria, three principal fault categories were delineated and further subdivided into compatible fault groups and subsets (Table 2; Figure 9).

Table 2

Table 2. Mean geometrical properties of the observed fault planes with calculated kinematic indicators and parameters. Fault planes were grouped following their geometrical and kinematic properties (Figure 9). Fault types: RF–reverse faults; NF–normal faults; SSF–strike-slip faults. Orientations of the P- and T-axes are based on constructed synthetic structural beach-ball diagrams.

Figure 9

Figure 9. Structural diagrams for the interpreted fault systems in DHS area. RF1 and RF2 represent reverse fault groups. NF1, NF2, and NF3 represent normal fault groups, whereas SSF1(a,b,c) and SSF2 (a,b,c) represent strike-slip fault groups. The red points, white rectangles, and blue triangles indicate σ1, σ2, and σ3 stress axes, respectively.

4.2.1 Reverse faults

Reverse fault planes (17 measurements) were separated into RF/1 and RF/2 fault groups (Table 2; Figure 9). Measured dominantly in the Mesozoic carbonate succession, the RF/1 group is characterized by two fault subsets: RF1/a shows an average ENE dip direction (dip angle of 35°, Figure 9), whereas RF1/b steeply dips towards NW at a dip angle of 74°. The RF/2 reverse fault group is composed of two subsets striking both N-S and WNW-ESE. The RF2/a subset is characterized by an average E dipping direction (dip angle of 67°), whereas the RF2/b subset includes planes dipping towards SW (dip angle of 44°). Structural analysis of the representative paleostress field mechanisms indicate that the paleostress compressional field for the RF1 fault group is associated with a P-axis generally trending NW‒SE (Table 2; Figure 9). The computations for the RF2 fault group show a compressional paleostress field that is associated with a P-axis trending NE‒SW (Table 2; Figure 9). These compressional paleostress fields resulted in the formation of the cogenetic fault-related structures (Figure 10A) that exhibit tectonic transport dominantly to NW/SE or NE/SW respectively.

Figure 10

Figure 10. (A) Conjugate reverse and normal fault pairs with cogenetic fold axes. Reverse faults and asymmetric folds indicate tectonic transport top to the NW, whereas normal faults show post-folding extensional relaxation. Structures are observed within Permian sandstone (E of Daruvar in the PV structural domain; 45°35′17.24″N, 17°20′17.25″E). (B) Normal fault (F-2/37) measured within Triassic dolomites. Besides striations and slickensides fault kinematics are characterized by systematic tensional fractures filled with syntaxial minerals (E from Daruvar, in the PV structural domain; 45°35′31.55″N, 17°17′37.28″E). (C) Normal fault (F-293/88) measured within Jurassic limestone. Striations and slickensides indicate normal displacement with sinistral oblique movement (E of Daruvar in the PV structural domain; 45°35′34.87″N, 17°17′57.78″E). (D) Subvertical fault plane (F-340/88) measured within Permian sandstone. Striations and slickensides indicate sinistral/dextral movements (E of Daruvar in the PV structural domain; 45°35′16.61″N, 17°20′20.96″E).

4.2.2 Normal faults

Besides reverse fault planes, 31 normal fault planes were measured (Figure 10B,C). Normal faults were separated into NF1, NF2, and NF3 groups (Table 2; Figure 9). Group NF1 is characterized by NE‒SW striking fault planes with fault subsets that are dipping towards either NNW or SE (dip angle of 46° and 43°, respectively). Kinematic analysis shows that these normal fault planes were formed within a paleostress field characterized by a subvertical P-axis steeply dipping towards the SSW (P-axis orientation is 206/86; Table 2) and subhorizontal T-axis trending NW‒SE resulting NW–SE extension. The group NF2 encompasses ESE‒WNW striking fault planes with fault subsets that are dipping towards N and S (dip angle of 64° and 28°, respectively). Kinematic and paleostress field analyses show that measured normal fault planes were formed within the paleostress field characterized by a subvertical P-axis steeply dipping towards the SSW and subhorizontal T-axis trending N‒S (Table 2) that resulted in NNE–SSW directed extension. The third normal fault group NF3 encompasses N‒S striking fault planes with fault subsets that were dipping towards ESE and WNW (dip angle of 44° and 47°, respectively; Table 2; Figure 9). Paleostress field analyses for this group show that fault planes were formed within the ESE–WNW directed extension that was influenced by the subvertical P-axis steeply dipping towards N (P-axis orientation is 11/86; Table 2), whereas the subhorizontal T-axis is generally trending E‒W.

4.2.3 Strike-slip faults

Strike-slip fault planes (43 measurements) were separated according to their geometric properties and kinematic compatibility into two principal fault groups (SSF1 and SSF2; Table 2; Figure 9). The SSF1 group includes dextral/sinistral faults (dip angle between 79° and 89°) that are dominantly striking either NE–SW or NW–SE, while the SSF2 group resembles steeply dipping conjugate fault pairs with ESE–WNW and NNE–SSW strike (dip angle between 61° and 89°). Mapped strike-slip fault planes were often observed with structural reactivation features that encompassed slickenside overgrowths. Both dextral and sinistral movement indicators were visible (Figure 10D). This reactivation indicates the interchange and re-orientation of the principal stress axes σ1 and σ3 within the same stress field. Kinematic analysis also points to faults’ structural reactivation indicating that the mapped fault planes were formed within two slightly different paleostress fields (Table 2): i) paleostress field associated with the NE–SW trending P-axis (T-axis trending NW–SE), and ii) paleostress field associated with general N-S trending P-axis that bends towards NNW or NNE (T-axis trending NE-SW or NW-SE; Figure 9). In addition, field observations of cross-cutting relationships between the mapped reverse/normal faults and strike-slip faults indicate that the strike-slip fault planes usually cut across and offset the reverse/normal faults. Structural reactivation indicators further evidence that some reverse/normal fault planes were also reactivated as dextral/sinistral faults.

4.3 Composite geological profiles and subsurface geological model of DHS

Composite geological profiles DHS-1 to DHS-5 (Figure 4; Figure 11; Supplementary Figure S3) represent the proposed geological relationships in the surface and subsurface of the study area reaching approximately an investigation depth of 2.5 km. The model (Figure 12) extends on an area of approximately 1,000 km² (36 and 28 km in E-W and N-S directions, respectively) including the Lonja-Ilova subdepression and the western margin of Mount Papuk. It consists of eight fault plane surfaces representing the principal regional faults and 42 horizon surface segments representing the base of the six main mapped units. In particular, the Quaternary deposits were modeled together with the Pliocene unit due to their limited thickness, while the Miocene andesites were merged with the coeval sedimentary unit.

Figure 11

Figure 11. Composite NW-SE striking geological profile DHS-3 (trace in Figure 4) investigating the central part of the study area. It shows the Daruvar anticline deforming the subsurface in the vicinity of Daruvar. The main regional faults in the profile are: Dežanovac fault (DeF), Končanica fault (KF), Daruvar fault (DF), Toplica fault (TF), Pakrac fault (PF), Velika fault (VeF), Glavica fault (GF). A detailed lithological description of the units is in Figure 3 and Supplementary Figure S2. Horizontal and vertical scale is 1:1.

Figure 12

Figure 12. (A) Set of composite geological profiles and geological data used to construct the 3D geological model. (B) Finalized 3D model. It conveyed the stratigraphic surfaces and the regional faults that delineate the structural framework of the wider DHS area.

Structurally, the constructed profiles could be subdivided into two domains. Here, the DHS-3 geological profile is used as an example (see Supplementary Figure S3 for other profiles) to describe the tectonic styles, geological units, and structures since it is located in the central part of the study area (Figure 11).

The NW domain of the geological profile covers the area of the Lonja-Ilova subdepression. This subdepression, a part of the Bjelovar depression, is filled with Neogene-Quaternary sediments that are up to 1.5 km thick, while the pre-Neogene basement is composed of Permian-Triassic sedimentary units and pre-Permian crystalline rocks (Malvić and Velić, 2011). The Permian-Triassic sedimentary complex is generally following the pre-Permian basement paleorelief. Going from NW towards SE, the Permian-Triassic complex is shallowing reaching the surface in the vicinity of Daruvar (Figure 3). Several NE-striking normal faults (e.g., Munija one and Munija 2; DHS-5 in Supplementary Figure S3) pinpoint the Neogene extension in the PBS and the opening of accommodation space. Furthermore, differential thicknesses of the Neogene deposits in the fault’s hangingwall/footwall are observed in a few locations (e.g., Dežanovac fault in DHS-4; Supplementary Figure S3) suggesting structural reactivation and tectonic inversion. This implies that some of the interpreted faults are polyphase structures, accommodating extension through the Neogene, and tectonic inversion during the Pliocene-Quaternary. In the central part of the study area, the pre-Neogene complex of the Lonja-Ilova subdepression crops out forming the western slopes of Mount Papuk (Figure 3; Figure 11; Figure 12). Here, the contact between the Neogene-Quaternary sediments and the pre-Neogene rock complex is mainly transgressive, but dozens of mapped tectonic contacts indicate NW-SE contraction with cogenetic reverse faults. Low angle reverse faults are usually either: i) blind faults forming cogenetic asymmetric anticlines (e.g., Dežanovac and Daruvar faults; DeF and DF, respectively, in Figure 11; Figure 12), or ii) thrust faults with ramps and flats (e.g., Pakrac fault; PF in Figure 11) forming fault-bend fold systems in the immediate hanging wall (Fossen, 2016; Nabavi and Fossen, 2021). The cogenetic asymmetric folds generally show gently inclined NW limbs, while SE limbs are steeper and shorter. This peculiar geometry suggests structure tectonic transport towards N-NW. Daruvar anticline is an example of a gentle asymmetric anticline (Figure 11 and DHS-2 and DHS-4 in Supplementary Figure S3) associated with the Daruvar fault (DF). It generally resembles a remobilized pre-Permian structural high that was transgressively covered by Permian-Triassic sediments and faulted afterward. In this context, the subvertical Toplica fault (TF in Figure 11) could be interpreted as a tensional fracture system developed in the hinge zone of the Daruvar anticline that was probably later reactivated as a dextral strike-slip fault zone (Kosović et al., 2023). Furthermore, NE-striking subvertical backthrust faults associated to the regional reverse faults were observed (e.g., Barica and Borki faults; Supplementary Figure S3). Mapped regional reverse faults are characterized by average dip angles of 45° and 55° in their steepest segments, while flat fault segments are characterized by dip angles ≤20°. Their maximum relative displacements are in a range between 0.5–1 km (Supplementary Figure S3). The Daruvar fault shows a maximum displacement of approximately 0.3 km (profile DHS-1 in Supplementary Figure S3).

The SE parts of the constructed profiles reflect the structural architecture of the western margin of the Slavonian mountains. Here, the tectonic uplift of the crystalline basement resulted in overall exposure of pre-Permian basement due to the erosion of the Permian-Mesozoic cover, while at the local scale, we could find patches of transgressively deposited Neogene sediments (Figure 11). Structural architecture of this area is cogenetic with two principal reverse low-angle faults, i.e., NNE-striking Voćin and Gradina faults (Figure 12; Supplementary Figure S3) which accommodated regional N-S compression. Relative displacements along the Voćin and Gradina reverse faults are in a range of a few hundred meters (Supplementary Figure S3).

5 Discussion

5.1 Tectonic emplacement and structural evolution of the study area

The Daruvar area, located in the immediate vicinity of the collision zone (i.e., Sava Suture Zone) between the Internal Dinarides and the Slavonian mountains resembles a complex litho-tectonic terrain that experienced a polyphase tectonic evolution (e.g., Jamičić, 1995; Balen et al., 2006). The different tectonic phases affected the structural relations in the subsurface, often showing tectonic overprints, without straightforward indication of distinctive deformation phases and associated geological features. Tectonic embayment of the study area started with the Variscan and afterward the Alpine-Dinarides-Carpathian orogeny which due to Cretaceous-Paleogene collision conveyed formation of the Slavonian mountains, a integral part of the regional nappe system (Schmid et al., 2008). The Cretaceous-Paleogene regional E-W (or NE-SW) compression promoted the formation of regional NW-SE striking reverse fault systems, and en-échelon folds characterized by NW-SE oriented fold axes. These structures were afterwards rotated counterclockwise of approximately 40° towards NE during the Paleogene (Tomljenović and Csontos, 2001; Ustaszewski et al., 2008). The post Cretaceous-Paleogene tectonic evolution of Daruvar area was further complicated by the Neogene-Quaternary evolution of the Pannonian Basin System (PBS). Inherited structures were affected by the Neogene E-W extension in the PBS, which locally caused additional rotation, structural reactivation, and tectonic inversion (Prelogović et al., 1998; Tari et al., 1999; Tomljenović and Csontos, 2001; Csontos and Vörös, 2004; Schmid et al., 2008; Ustaszewski et al., 2008). As a part of PBS, the Daruvar area during Pliocene-Quaternary was finally affected by the tectonic inversion of the existing structures due to regional N-S compression/transpression stresses (Ustaszewski et al., 2008; Schmid et al., 2020). Rejuvenated regional contraction enhanced the N-S shortening and continuous folding/re-folding processes of the existing structures synchronously with the lateral displacement processes along the strike-slip faults (Jamičić, 1995; Ustaszewski et al., 2008).

Structural data analysis of strata orientation, fracture systems, and fault systems presented in this study undoubtedly supports and confirms the complexity of the tectonic evolution in the Daruvar area. Analyses of strata orientations evidenced a polyphase folding in the area. Observed folds with fold axes that are gently dipping towards N or NW and S or SE (red folds in Figure 6) correspond to the Cretaceous-Paleogene E-W (NE-SW) contraction phase that culminated with the counterclockwise rotation of structures during the Paleogene (e.g., Tomljenović and Csontos, 2001; Ustaszewski et al., 2008). At the same time, E-W (or NE-SW) striking folded structures (blue folds in Figure 6) developed by the Pliocene-Quaternary N-S contraction were also observed in this study (Jamičić, 1995; Tomljenović and Csontos, 2001).

Cogenetically with the formation of these folded and refolded systems, the Paleozoic-Mesozoic sedimentary complex experienced extensive brittle deformation and the formation of fracture systems with preferred orientations (Figure 8). Measured fractures resembled N-S (locally NNE-SW, NNW-SSE) and subordinately E-W striking subvertical tensional fractures that were subparallel with the observed fold hinge zones (Figure 8). Locally, especially in the Mesozoic carbonate complex, N-S striking fracture systems show shear reactivation features (e.g., slickenside overgrowths) characterized by dextral/sinistral motions. This reactivation is connected to the Pliocene-Quaternary N-S oriented P-axis (Herak et al., 2009), and in general widens the damage zone of the N-S striking folded structures increasing the fracturing of the bedrock. On the other hand, E-W striking discontinuities are less frequent and generally without indications of structural reactivation suggesting that E-W striking folded structures are less affected by ongoing tectonic deformation due to their structural position in respect to the low strain rates (<1–2 mm/y; Grenerczy et al., 2005) of the N-S oriented P-axis (Herak et al., 2009).

Field observations of shear fractures/fault planes with associated kinematics and cross-cutting relationships support the results of bedding and fracture system analyses. Correlative to the Cretaceous-Paleogene (E-W contraction and counterclockwise structural rotation) and the Pliocene-Quaternary (N-S contraction) tectonics, formed fold systems and analyzed reverse fault group subsets suggested paleostress field with P-axes generally trending either NW‒SE (RF1) or NE‒SW (RF2) (Table 2; Figure 9). Synthetic structural focal mechanisms and field data suggested that orientations of average fault subsets are in correspondence with fault systems that accommodated the Cretaceous-Paleogene E-W contraction in the study area (today NW-SE orientated P-axis that rotated 40° counterclockwise) as well as Pliocene-Quaternary N-S contraction (i.e., Voćin fault, Gradina fault, Dežanovac fault and Daruvar fault; see Figure 3). At the same time, the RF2 reverse fault planes and the computed structural focal mechanisms suggested also a NE‒SW oriented contraction. Our structural results for this fault subset are slightly different in relation to the regional N-S contraction suggesting that local structural differences and inherited structures may contribute to a local reorientation of the P-axis. Focal mechanisms of normal faults generally suggest NW–SE, NNE–SSW, and ESE–WNW extension with subvertical P-axes (Table 2; Figure 9). This extension may be related to the final stages of the thrusting and folding in the area during the Late Cretaceous-Paleogene as a result of gravitational sliding of existing structures (e.g., Tavani et al., 2012) or by structural hanging wall collapse due to Late Oligocene-Miocene extension of the PBS. The most preserved slickensides, striations, and other kinematic indicators were observed along the strike-slip fault planes (SSF1 and SSF2; Table 2; Figure 9). Paleostress analysis and structural focal mechanisms indicate generally N-S oriented transpression which locally deflects to NW-SE and NE-SW (Figure 9). With their subvertical geometry mapped strike-slip faults may be associated with the Pliocene-Quaternary structural reactivation of faults originated during the Neogene extension of the PBS (Jamičić, 1995) or with the partly-tectonic reactivation of inherited fracture systems deforming the fold system hinge zones.

5.2 Hydrogeological conceptual model of the Daruvar hydrothermal system

The structural setting of the Daruvar hinterland and its tectonic evolution are integrated with available geochemical, hydrogeological, and geophysical data (Borović, 2015; Borović et al., 2019; Kosović et al., 2023; Urumović et al., 2023) on the Daruvar thermal field and the thermal waters accounting for: i) the correlation between fault and thermal systems (e.g., Curewitz and Karson, 1997; Faulds et al., 2013; Moeck, 2014), and ii) the increase of the permeability field due to fracture zones acting as preferential flow paths (e.g., Faulkner et al., 2010; Bense et al., 2013).

Geochemical and isotope data of the Daruvar thermal waters (Borović, 2015) show that: i) the waters are relatively young (10–15 ka) and originate from the precipitation in the nearby mountainous hinterland, which provides the hydraulic potential for the fluid flow, and ii) they circulate and most likely infiltrate into the Mesozoic carbonate complex constituting the geological assemblage of the western Papuk due to their calcium-bicarbonate hydrochemical facies. As a matter of fact, Mesozoic carbonates are an important reservoir for thermal waters in the PBS (Horváth et al., 2015; Rman et al., 2020) hosting approximately 20 geothermal installations just in its Croatian part (mostly thermal spas and balneological therapy centers; Borović and Marković, 2015).

The potential DHS recharge area could span over the outcrops of Mesozoic carbonates in the structural domains of Dujanova kosa, Petrov vrh, and Sirač (Figure 3). Since these areas are comprised in different tectonic blocks separated by reverse faults, they can be considered separate hydrogeological compartments. The Sirač structural domain is in the hanging wall of the Pakrac reverse fault (Figure 11). The Mesozoic units are here mostly in contact with Miocene formations in the footwall reducing the northward continuity of the Mesozoic carbonate aquifer (Figure 3). Furthermore, subthermal springs occur in the vicinity of Sirač and an artesian aquifer in a quarry nearby is reported (personal communication). Due to their geographical positions, they likely drain the recharge in the Sirač domain diminishing the potential recharge to the Daruvar thermal area. The Dujanova kosa structural domain has two outcrops of Mesozoic carbonates that could act as recharge areas of the DHS. The carbonates in the southeastern part of the domain (i.e., south of the Dujanova kosa and Crni vrh peaks; Figure 3) are included in a syncline structure bounded by the low permeable pre-Permian basement both to the N and the S (section DHS-5 in Supplementary Figure S3). Therefore, they can be considered hydrogeologically isolated. The carbonates in the western part of the block are at the footwall of the Daruvar/Gradina reverse fault being deeper than the same formations in both the Petrov vrh domain and the Daruvar area (sections DHS-5 and DHS-3, respectively, in Supplementary Figure S3). Therefore, they could eventually drain the nearby blocks and it is unlikely that they could host waters with pressure high enough to flow into shallower aquifers. This fact is enforced by the low elevations of the Mesozoic outcrops in the study area, thus providing a similar hydraulic potential in all blocks.

The Mesozoic carbonate complex in the Petrov vrh structural domain is the most probable recharge area for the DHS due to its geographical and geological settings. The Mesozoic complex extends for approximately 15 km² eastward of Daruvar being approximately 4–8 km from the spring area. The elevation of the recharge area in Petrov vrh is generally the highest with an average elevation of 400 m a.s.l. And up to 613 m a.s.l. (average elevation of 397 and 349 m a.s.l. In Dujanova kosa and Sirač, respectively; maximum elevation of 570 and 552 m a.s.l. In Dujanova kosa and Sirač, respectively). The geological framework of Petrov vrh is characterized by a regional N-S striking overturned syncline (Figure 3) formed in the hanging wall of the Daruvar/Gradina fault (sections DHS-5 and DHS-3 in Supplementary Figure S3). The structural analysis conducted in the PV domain (Figure 8) revealed two sets of steeply dipping fractures that generally strike E-W (NE-SW) and N-S (NW-SE). Both sets are deformed by the last Pliocene-Quaternary N-S compression (Herak et al., 2009) and they can be considered tectonically active. The constant deformation of fractures is crucial since they maintains their aperture preventing the sealing by precipitation of minerals and preserving the permeability field of the bedrock. Due to its favorable structural position in respect to recent N-S oriented stress field, the N-S discontinuities could have a wider and highly deformed damage zone connected to their reactivation. The E-W discontinuities are less frequent, but they may result from a local transtensional regime in the hinge zone of the currently deformed E-W oriented folds. These conditions are favorable for the high permeability of the bedrock resulting in a high effective infiltration in the recharge zone.

The meteoric waters infiltrate due to the intense fracturing of the bedrock and flow westward favored by: i) the general westward dipping of the strata, ii) the pre-Permian and Permian units acting as the aquitard below the Mesozoic reservoir, and iii) the E-W tensional open fractures. The set of beds dipping W is the most frequent in both the eastern and central sectors of the PV domain (Figure 6). The pre-Permian crystalline rocks and the Permian sedimentary units have moderate to low permeability due to their lithologies (e.g., Domenico and Schwartz, 1998). Conversely, dissolution processes in the Mesozoic carbonates could add to the fracturing enhancing their permeability field (e.g., Goldscheider et al., 2010). This contrast prevents a deep infiltration of the meteoric waters that are more prone to flow in the carbonate reservoir. In addition, the dipping of the aquitard units toward W channels the fluid flow favoring the westward circulation of the infiltrated waters. Finally, the E-W fractures could act as preferential flow paths. The transtensional regime favors the opening of the fractures and increases the permeability that depends on the square of the fracture hydraulic aperture (e.g., Domenico and Schwartz, 1998).

The infiltrated waters are warmed by the heat flow of 80–100 mW/m² in this part of the PBS resulting in a local geothermal gradient of 35°C–40°C/km (Horváth et al., 2015; Borović et al., 2019). The reservoir equilibrium temperature of the thermal waters calculated using SiO₂ geothermometers is approximately 80°C (Borović, 2015). Considering an infiltration temperature of approximately 10°C and a purely conductive heat flow, the waters should reach a depth of approximately 2 km to approach the reservoir equilibrium temperature. However, the fracturing of the bedrock favors the occurrence of convective processes that increase the circulation and the temperature in the aquifer. As a matter of fact, a gradient of approximately 70°C/km is measured in the thermal wells of Daruvar (Borović et al., 2019) corroborating the impact of local convection on the temperature distribution. At the scale of the recharge and flow-through area of the DHS, it is reasonable to expect a gradient slightly higher than the regional value, which could permit to approach the reservoir equilibrium temperature at the maximum aquifer depth of the aquifer being 800–900 m below the ground level in the central part of Petrov vrh. The high geothermal gradient in the Daruvar area could add to the water temperature in the flow-through part of the system reaching the reservoir equilibrium temperature.

The thermal waters flow in the Mesozoic reservoir and reach the Daruvar city area located in the immediate vicinity of the Daruvar anticline (Figure 11; Figure 13). This NNE-SSW striking structure is an asymmetric fold (tectonic transport top to the W) cogenetic to the SE dipping Daruvar reverse fault. The polyphase tectonic evolution affecting the study area suggests structural reactivation of the Daruvar anticline favoring the continuous fracturing of the bedrock. A localized extensional regime is expected in the topmost section of the fold hinge zone increasing the fracture aperture and the permeability field (e.g., Frehner, 2011; Li N. et al., 2018). Field observations (Figure 13B) and core samples from the thermal wells showed that the Mesozoic reservoir is moderately fractured with localized zones of intense fracturing. The thermal waters rise to shallow depths in the damage zone of the Daruvar fault and cogenetic fractures that are deforming the hinge of the Daruvar anticline. In this context, the subvertical NE-striking dextral Toplica fault, which is a structurally reactivated and tectonically inverted tensional fracture system, could act as a preferential flow path for the quick rise of the thermal waters with a minor loss of temperature from the deeper part of the reservoir. The Toplica fault could be the fault F5 in the local scale model of the Daruvar spring area obtained through shallow geophysical investigations (Kosović et al., 2023). This fault borders Daruvar thermal field westwardly, and together with an E-W striking fault to the S (F1 in Kosović et al., 2023) accommodates the uplift of the Mesozoic thermal reservoir to shallower depths. The geophysical investigations highlighted that the thermal springs occur within the interaction zone of local scale faults/fractures. Interaction zones are preferential locations for thermal springs because they favor the kinematic transfer between faults increasing the rock fracturing and the permeability field (e.g., Curewitz and Karson, 1997; Faulds et al., 2013). This structure further localizes the flow of the Daruvar thermal waters resulting in four thermal springs with temperature between 38°C and 50°C.

Figure 13

Figure 13. (A) Schematic conceptual model of the structural assemblage in the subsurface of the Daruvar area (modified from Frehner, 2011; Li Y. et al., 2018) affecting the upwelling of theDaruvar thermal waters. The NNE-SSW striking Daruvar anticline structurally forms in the hangingwall of the Daruvar reverse fault (see Figure 3 for location). In the hinge zone, cogenetic tensional fractures were structurally reactivated during the Pliocene-Quaternary and tectonically inverted as strike-slip fault zone (i.e., Toplica fault). The localized extensional regime in the topmost part of the hinge zone and its polyphase deformation increase the fracturing of the bedrock and the permeability field favouring the outflow of the thermal waters in the Daruvar spring area. (B) Heavily tectonized Triassic dolomites in the Daruvar thermal spring area (45°35′41.74″N, 17°13′38.43″E).

6 Conclusion

This research focused on the reconstruction of the geological framework and the tectonic evolution of western Papuk to detail the impact of regional and local scale fold and fault/fracture systems on the development of the Daruvar hydrothermal system (DHS) and its geothermal resource. The reconstruction was conducted by integrating surficial field investigations and available surface and subsurface geological and geophysical data, both at regional and local scales. The structural data analysis evidences a complex pattern of folds deformed by faults and fracture systems in a manner compatible with the polyphase tectonic evolution of the Slavonian mountains and the SW part of the Pannonian Basin System. The construction of 2D composite geological profiles that were integrated into a 3D geological model favored the visualization of the geological assemblage detailing the structural relations among the different blocks in the study area. The geological and tectonic reconstructions were integrated with geochemical data on the Daruvar waters and local scale geological and geophysical data on the thermal field to propose a conceptual model of the DHS. The conceptual model highlighted the importance of regional and local structures (i.e., folds, faults, networks of fractures) that are causative factor for the regional to local flow of the Daruvar thermal waters.

The geological structure of the thermal system is characterized by a regional NNE-SSW striking asymmetric fold (tectonic transport top to the W) formed in the hangingwall of SE dipping Daruvar reverse fault, a western prolongation of the E-W striking Gradina thrust fault. This structure favors both the extensive outcropping of the Mesozoic carbonate rock complex (i.e., the reservoir of the system) in the recharge area and the fracturing of the bedrock. Deformed by the current Pliocene-Quaternary N-S compressional regime, cogenetic steeply dipping E-W and N-S striking fracture systems are potential paths for the infiltration, flow, and rise of the thermal waters. The E-W discontinuities could represent the principal regional flow paths, while N-S fractures (as the Toplica fault imaged in Daruvar through local scale geophysical investigations) could enable the local quick rise of the thermal water from the deeper part of the reservoir due to their structural position in the anticline hinge zone and the local extensional regime increasing the fractures aperture and the permeability field of the bedrock.

Insights gained through the application of these research methodologies can be utilized as an example for the 3D subsurface reconstruction of areas affected by the deep circulation of waters resulting in the development of a geothermal resource. Such kind of multidisciplinary reconstruction could foster the estimation of the potential of a geothermal resource aiding the assessment of the reservoir volumes and the development of hydrogeological numerical modeling of fluid flow and heat transport.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

IK: Conceptualization, Data curation, Formal Analysis, Investigation, Writing–original draft, Writing–review and editing. BM: Conceptualization, Formal Analysis, Investigation, Methodology, Writing–original draft, Writing–review and editing. IP: Investigation, Software, Validation, Visualization, Writing–review and editing. MP: Conceptualization, Investigation, Visualization, Writing–original draft, Writing–review and editing, Data curation. MM: Data curation, Formal Analysis, Writing–review and editing, Investigation. MiP: Investigation, Project administration, Writing–review and editing. SB: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing–review and editing, Resources.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The research was funded by the HyTheC project of the Croatian Science Foundation (HRZZ), grant number UIP-2019-04-1218, and the publication fee was paid by the GeoTwinn project (H2020-WIDESPREAD-05-2017-Twinning project), grant number 809943.

Acknowledgments

The authors would like to thank the Daruvarske toplice–Special Hospital for Medical Rehabilitation for logistic help on-site and sharing of existing materials as well as to the Croatian Hydrocarbon Agency for granting access to the legacy hydrocarbon research data (seismic reflection profiles and deep exploration wells).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2024.1401935/full#supplementary-material

References

Allmendinger, R. W., Cardozo, N., and Fisher, D. M. (2011). Structural geology algorithms: vectors and tensors. doi:10.1017/CBO9780511920202