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).

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).

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.

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).

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.

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.

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).

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.

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).

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.

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).

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.

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.