www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society Article history: Manuscript recieved: May 8, 2024 Revised manuscript accepted: May 28, 2024 Available online: June 21, 2024 Keywords: fracture network, carbonate aquifer, digital outcrop model, discrete fracture network modelling, northern Croatia 2024 | 77/2 | 99–125 | 12 Figs. | 1 Tab. | 4 Suppl. Figs. | 9 Suppl. Tabs. | 1. INTRODUCTION Groundwater represents a strategic resource since it provides approximately half of the global volume of water withdrawn for domestic use and approximately one-fourth of the water for irrigation (UNITED NATIONS, 2022). Valuable ground- water resources are hosted in fractured bedrock aquifers (e.g., SINGHAL & GUPTA, 2010). It has been estimated that 10% of the world population relies on freshwater from carbonate aquifers (STEVANOVIĆ, 2019). This importance increases in the Mediterranean region where the karst groundwater- dependent ecosystems are extensively used, facing environ- men tal threats (SIEGEL et al., 2023). Site-specific management plans considering both the principal flow paths and the hydra- ulic properties of the aquifer are crucial for a more efficient conservation of groundwater resources. In rocks having low primary porosity and permeability, the bulk hydraulic properties are controlled by the network of fractures crosscutting the rock mass (SINGHAL & GUPTA, Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Ivan Kosović1, Bojan Matoš2,*, Stefano Casiraghi3, Gabriele Benedetti3, Tihomir Frangen1, Kosta Urumović1, Ivica Pavičić2, Andrea Bistacchi3, Silvia Mittempergher3, Marco Pola1 and Staša Borović1 1 Croatian Geological Survey, Department of Hydrogeology and Engineering Geology, Sachsova 2, 10000 Zagreb, Croatia 2 University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering, Department of Geology and Geological Engineering, Pierottijeva 6, 10000 Zagreb, Croatia; (*corresponding author: bojan.matos@rgn.unizg.hr) 3 Università degli Studi di Milano Bicocca, Dipartimento di Scienze dell'Ambiente e della Terra, Piazza della Scienza, 4, 20126 Milan, Italy doi: 10.4154/gc.2024.11 Abstract Highly fractured Mesozoic carbonate rocks are the main reservoir of many geothermal re- sources in northern Croatia, being of environmental, cultural, and economic value for the local and regional communities. The Daruvar thermal springs (temperatures < 50°C) rep- resent the outflow area of an intermediate scale, tectonically controlled, hydrothermal sys- tem hosted in Triassic carbonate rocks. Several investigations have been conducted in the Daruvar area detailing the architecture of regional and local fracture networks and quanti- fying the hydrogeological parameters of the thermal aquifer. In this work, an integrated approach based on structural and hydrogeological investigations was employed to model the network of fractures in the reservoir and quantify its impact on the hydraulic properties. Structural investigations were conducted in the Batinjska Rijeka quarry, considered as an outcrop analogue of the thermal aquifer, employing both a classical field approach and the virtual quantitative analysis of a 3D digital outcrop model. Structural analysis of the digital outcrop model allowed identification of two sub-vertical systems of discontinuities, dipping to the NW and the WSW respectively, in accordance with the data collected through direct field measurements. The main geometric features of the discontinuity network and their statistical distributions were employed to construct discrete fracture network models at both the outcrop scale (approximately 100 m) and the aquifer scale in Daruvar (approximately 700 m). Calibration of the input parameters allowed modelling of porosity and permeability values that reproduce the field values assessed through pumping tests, well tests, and well logging. This work highlights the importance of integrating geological and hydrogeological investigations to obtain a more reliable reconstruction and quantification of the processes driving the fluid flow in fractured aquifers and affecting the spatial distribution of their hydraulic properties. 2010). For fractures here we consider every kind of dis con ti- nuity including bedding discontinuities, schistosity, joints, shear fractures, veins, stylolites and other dissolution features, deformation and compaction bands, etc. (e.g., BISTACCHI et al., 2020; STORTI et al., 2022). The permeability field in the rock mass depends on the geometric features of the fractures (i.e., orientation, intensity, aperture, length, etc.) and the connectivity of the fracture network (e.g., FAYBISHENKO et al., 2000; DE DREUZY et al., 2001, 2002). Fault zones affect the groundwater circulation at both regional and local scales, primarily influencing the development of springs in hard rock terrains (e.g., KEEGAN-TRELOAR et al., 2022). The well- established conceptual model of fault zones generally includes a low permeable fault core surrounded by a highly fractured and permeable damage zone (e.g., CAINE et al., 1996; FAULKNER et al., 2010), although different hydraulic behaviours are documented (i.e., high permeability fault core; SMITH et al., 2013). The structural architecture of fault zones and the impact of fracture networks on the permeability field mailto:bojan.matos@rgn.unizg.hr G eo lo gi a C ro at ic a 100 Geologia Croatica 77/2 have been extensively investigated using both structural and hydrogeological approaches. Still, the integration between these two disciplines is limited due to different methods and scales of investigation (e.g., BENSE et al., 2013). Outcrop geological and structural mapping supported by numerical modelling and pumping tests are the most common approaches to quantify the permeability of fractured aquifers in the field (SINGHAL & GUPTA, 2010; BENSE et al., 2013). Outcrop characterisation determines the main kinematic, chronological, and geometric features of the fracture network and the statistical distribution of fracture parameters. The classical field investigations can be supported by photogrammetric reconstruction of the outcrop. This technique can be used to investigate extensive outcrops, which may not even be accessible in field investigations, obtaining large datasets for a solid statistical analysis of the main discontinuity sets through the virtual structural analysis of the digital outcrop model (e.g., HODGETTS, 2013; BISTACCHI et al., 2015; MARTINELLI et al., 2020). These statistical distributions can be used to stochastically reconstruct the fracture network, employing the discrete fracture network (DFN) approach. Considering each fracture as a discrete object with peculiar hydraulic properties surrounded by an impermeable matrix, DFN models are profitably used to determine the permeability of fractured aquifers and rock masses (e.g., VOECKLER & ALLEN, 2012; LEI et al., 2017; MEDICI et al., 2020; CECCATO et al., 2021; SMERAGLIA et al., 2021; MAMMOLITI et al., 2023; MEDICI et al., 2023). However, calibrations and uncertainty analysis of parameters (particularly aperture) with permeability values from field investigations are generally lacking. Pumping tests consist of pumping groundwater from a well at different flow rates and measuring the drawdown in observation wells. The hydraulic properties (mostly transmissivity and hydraulic conductivity) of the aquifer are assessed by modelling the variations of the drawdown using analytical solutions (e.g., KRUSEMAN et al., 2000). Due to the heterogeneous structure of fractured aquifers, the identification of flow-bearing paths is crucial for assessing the most suitable modelling approach (e.g., CARRERA & MARTINEZ-LANDA, 2000; GUIHÉNEUF et al., 2021), as well as for understanding and managing the site. In contrast, hydrogeological investigations through pumping and well tests generally focus on assessing the hydrogeological properties of the aquifer and the efficiency of the tested well, but they generally lack insights into the fracture network architecture. Among groundwater resources in fractured aquifers, geothermal resources are particularly important since they are potential renewable sources of raw materials and energy (e.g., FINSTER et al., 2015; SZANYI et al., 2023). Geothermal resources generally have a lower yield and higher management costs than groundwater resources. Therefore, assessing the hydraulic characteristics of the reservoir and the impact of regional and local fracture systems is crucial for the quantification of the renewable component of the system and the proposal of management plans for long-term sustainable exploitation (e.g., RYBACH & MONGILLO, 2006). Northern Croatia is rich in thermal springs that are mostly used for balneotherapy and tourism (BOROVIĆ & MARKOVIĆ, 2015). Here, connected geothermal resources have great potential since they share the favourable thermal features of the Pannonian area (HORVÁTH et al., 2015). The thermal spring area in Daruvar city (Figs. 1 and 2) is one of the most investigated thermal regions in Croatia (e.g., BOROVIĆ, 2015; BOROVIĆ et al., 2019; KOSOVIĆ et al., 2023; URUMOVIĆ et al., 2023; KOSOVIĆ et al., 2024). Thermal springs in Daruvar (temperature up to 50°C; Fig. 3) have been documented since the Roman era, representing an economic and cultural value for the area. The Daruvar thermal field is the outflow area of an intermediate scale geothermal system hosted in a Mesozoic fractured carbonate reservoir. Previous studies investigated the hydraulic properties of the reservoir and the architecture of regional and local fracture systems (e.g., BOROVIĆ et al., 2019; KOSOVIĆ et al., 2024). In this work, an integrated approach combining structural data analysis and well testing is applied to detail the hydraulic properties of the Daruvar thermal aquifer. Structural investigations were conducted in the Batinjska Rijeka quarry (BRQ; Fig. 2) located 4 km NE of Daruvar. There, the Mesozoic carbonate rock complex is highly fractured representing an outcrop analogue of the aquifer. The discontinuities deforming the rock mass were investigated through both classical field investigations and the analysis of the digital outcrop model (DOM) of the quarry. A DFN model reproducing the observed network of discontinuities was constructed, and its capability of assessing the hydraulic properties of the aquifer was tested using data collected from several well tests and logs conducted in the Daruvar thermal field. 2. GEOLOGICAL AND HYDROGEOLOGICAL SETTINGS 2.1. Regional tectono-structural evolution The study area is located along the western margin of Mount Papuk, which is part of the Slavonian Mts. (Fig. 1). The Slavonian Mts. are the best exposures of the Tisza continental block, i.e., a lithospheric fragment of the European foreland that formed beside the Adria microplate during the Middle Jurassic (Fig. 1; SCHMID et al., 2008). The polyphasic tectono- metamorphic evolution of the study area started with the Variscan and continued with the Alpine-Dinarides-Carpathian orogeny. The continuous convergence between the Tisza-Dacia mega block, the Adria microplate, and the European plate through the late Mesozoic-Cenozoic caused the stacking of the Mecsek, Bihor, and Codru nappe systems (BALEN et al., 2006; SCHMID et al., 2008). The Slavonian Mts. are part of the Bihor nappe system that also crops out in southern Hungary. Their structural architecture was significantly influenced by the proximity of the Sava suture zone (Fig. 1), which resembled the Cretaceous-Palaeogene regional suture zone between the Tisza block toward the NE and the Adria microplate block toward the SW (SCHMID et al., 2020). Furthermore, regional fault zones such as the Periadriatic lineament and the Mid-Hungarian fault line (Fig. 1) played important roles as tectonic boundaries between the main regional blocks accommodating the clockwise rotation of the Tisza-Dacia mega block, the counterclockwise rotation of the Adria plate, and the eastward lateral extrusion, at a scale of several hundreds of km (Fig. 1; PRELOGOVIĆ et G eologia C roatica 101Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing al., 1998; TARI et al., 1999; CSONTOS & VÖRÖS, 2004; USTASZEWSKI et al., 2010) After the climax of Cretaceous-Palaeogene tectonism, the opening of the Pannonian Basin System (PBS) between the Carpathians and the Dinarides occurred (BALEN et al., 2006). The tectonic evolution of the PBS was characterised by polyphase extension, compression, and tectonic inversion. Formed by Early-Middle Miocene NE-SW oriented back-arc type lithospheric extension, the PBS was characterised by systems of NW-SE striking normal faults that formed rift and wrench-type troughs filled with large amounts of syn-rift deposits (HORVÁTH et al., 2006; USTASZEWSKI et al., 2010). These basins and correlative structural highs were subjected to localised tectonic inversion during the Middle Miocene, but deepening and rapid thermal subsidence along existing NE-SW and NW-SE striking faults prevailed until the Early Pliocene (TOMLJENOVIĆ & CSONTOS, 2001; MALVIĆ & CVETKOVIĆ, 2013). Changes in the stress field due to the counterclockwise rotation of the Adria microplate and the consumption of the subducted European plate yielded a transition from extension to N-S compression/transpression, which resulted in regional tectonic inversion and fault reactivation (BADA et al., 2007; USTASZEWSKI et al., 2014). The Pliocene-Quaternary tectonic inversion accommodated regional folding and prevalent reverse and strike-slip faulting (JAROSINSKI et al., 2011). It yielded the uplift of the basement structural highs (e.g., Slavonian Mts.), the tectonic inversion of pre-existing basement structures and the formation of strike-slip fault-related positive f lower structures (PRELOGOVIĆ et al., 1998; TOMLJENOVIĆ & CSONTOS, 2001). 2.2. Geological and structural setting of the study area The geological setting of the NW part of the Slavonian Mts. was investigated by JAMIČIĆ et al. (1989). Due to its complexity, the lithological units of the study area (Fig. 2) were Figure 1. Regional tectonic map of the Alps, Dinarides, Tisza-Dacia mega block, and the European Foreland framework (after SCHMID et al., 2008, 2020). The map shows the principal regional faults and the outline of the Neogene basins (black line) including the Pannonian basin. The principal towns and the neighbouring states are depicted (acronyms: AUT – Austria, BiH – Bosnia and Herzegovina, HRV – Croatia, HUN – Hungary, ITA – Italy, SLO – Slovenia, SRB – Serbia). G eo lo gi a C ro at ic a 102 Geologia Croatica 77/2 reorganised by KOSOVIĆ et al. (2024) as follows: i) pre- Permian crystalline rocks, ii) Permian sedimentary units, iii) Triassic carbonate 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 pre-Permian crystalline basement extends across a large surface area in western Papuk (Fig. 2), encompassing magmatic and metamorphic rock complexes that are partly transgressively covered by younger units (JAMIČIĆ et al., 1989; PAMIĆ et al., 2003). The onset of the Permian succession began with transgressive, well- layered conglomerates and quartz sandstones. Triassic and Jurassic formations constitute the Mesozoic carbonate rock complex that represents the main reservoir of the Daruvar thermal waters. The Triassic succession is characterised by Lower to Middle Triassic clastic deposits (e.g., sandstones, siltstones, and laminated shales) sedimented within a shallow- water environment (JAMIČIĆ et al., 1989; ŠIKIĆ, 1981). A gradual transition towards carbonates (e.g., dolomites, limestones, and crinoid limestones with chert) is observed in the Middle to Late Triassic, while the Late Triassic stratigraphic top is characterised by a pure dolomite-limestone succession (occasionally incorporating dolomite and limestone breccias and marly limestones; Fig. 2). In a sequence, Jurassic deposits are observed being represented by grey platy limestones with cherts. The Cretaceous collision between the Tisza block and the Adria and Europe plates caused emersion in the Slavonian Mts. Sedimentation was partly restored during the Middle Miocene E-W regional extension of the PBS. It continued through the Late Miocene-Pliocene and Quaternary (Fig. 2), characterised by a transition from coastal to freshwater depositional environments (SAFTIĆ et al., 2003). The youngest sediments are Quaternary unconsolidated sediments (from Figure 2. Simplified geological map of the Daruvar area and lithostratigraphic succession (modified from KOSOVIĆ et al., 2024). The black polygon rep- resents the extent of Fig. 4A. Topographic peaks are denoted by white triangles, while the main towns are shown as teal polygons. The NW-SE striking geological profile (bottom) depicts the general structural relationships along the western margin of Mount Papuk, in the vicinity of the Batinjska Rijeka quarry (BRQ). Fault acronyms: BoF – Borki fault; DF – Daruvar fault; DeF – Dežanovac fault; GF – Gradina fault; PF – Pakrac fault; TF – Toplica fault. G eologia C roatica 103Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing sandy gravels to sandy clays), deposited by alluvial and gravitational processes (JAMIČIĆ et al., 1989). The Cretaceous to Quaternary tectonic evolution of the study area is affected by the structural architecture of the Slavonian Mts. and the surrounding basins (KOSOVIĆ et al., 2024). The complex structural pattern of generally E-W striking thrust faults (i.e., Voćin fault, Gradina fault, Deža- novac fault, and Daruvar fault; Fig. 2) and multi-folded structures gently dipping towards both the NNW and SSE or ENE and WSW suggests a polyphase evolution, which began with regional Cretaceous-Palaeogene E-W contraction. This contraction formed an initial NW-SE striking fault and NNW and SSE dipping fold pattern, which rotated counterclockwise approximately 40° towards the NE during the Palaeogene (TOMLJENOVIĆ & CSONTOS, 2001; USTASZEWSKI et al., 2008). Counterclockwise regional rotation resulted in the final structural emplacement of the existing faults and folds, reoriented as E-W striking thrust faults and an ENE and WSW dipping fold system. Though the Neogene was characterized by regional E-W extension and the formation of predominantly normal faults, inherited faults/fold systems in the study area experienced only partial tectonic stretching and gravitational sliding. Contemporary, sediment deposition within pre- existing synclines occurred (Fig. 2; (JAMIČIĆ, 1995; TARI et al., 1999; TOMLJENOVIĆ & CSONTOS, 2001). The Pliocene- Quaternary regional N-S compression/transpression promoted continuous shortening and re-folding processes of the pre- existing structures coeval with the formation of dextral/ sinistral faults (e.g., Toplica fault; TF in Fig. 2), which locally displaced existing E-W striking thrust faults. 2.3. Hydrogeological setting of the Daruvar thermal system The lithostratigraphic sequence in the Daruvar thermal area (KOSOVIĆ et al., 2023) was assessed through boreholes and wells: i) Quaternary alluvial deposits, ii) Miocene or Pliocene marls, iii) Miocene bioclastic limestone, and iv) Triassic dolomites and limestones with an investigated thickness of 130 m. Core analyses and well logging evidenced that the Triassic carbonates are characterised by a uniform background fracture network and by localised corridors (total thickness of 24 m in the deepest well), showing a higher frequency of fractures. The thermal waters are hosted in fractured Triassic carbonate rocks, and subthermal waters can be found in the sandy layers of the alluvial cover and the Miocene biocalcarenites (BOROVIĆ et al., 2019). The aquifer transmissivity assessed through pumping and well tests ranges from 0.015 to 0.03 m2/s (BOROVIĆ et al., 2019; URUMOVIĆ et al., 2023), while the average porosity measured through neutron borehole logging is 7.9%. Hydrochemical and isotope analyses of the thermal waters (BOROVIĆ, 2015) evidenced: i) a temperature from 18.2 to 49.8°C, nearly-neutral pH, and electrical conductivity of 550–700 μS/cm, ii) a calcium-bicarbonate hydrochemical facies indicating flow through a carbonate aquifer, iii) the meteoric origin of the water as depicted by O and H stable isotope ratios, iv) a reservoir equilibrium temperature of 80°C calculated using SiO2 geothermometers, and v) a mean residence time between 11 and 15 ka assessed through 14C activity. These data, combined with the results of regional and local scale structural and geophysical investigations, were used to propose a conceptual model of the Daruvar hydrothermal system (BOROVIĆ et al., 2019; KOSOVIĆ et al., 2024). The recharge area of the system is located in the topographically higher Petrov vrh, approximately 8 km E of Daruvar (Fig. 2). Here, Mesozoic carbonate rocks extensively crop out and are deformed by N-S and E-W striking fractures that are potential paths for the infiltration of the meteoric waters (KOSOVIĆ et al., 2024). The general westward dip of the strata and E-W tensional open fractures favour the westward flow of the infiltrated waters in the Mesozoic reservoir. The waters are warmed to 80°C by the regional heat flow and convection processes in the aquifer. In the Daruvar area, they intercept an asymmetric anticline cogenetic to the Daruvar reverse fault (DF in Fig. 2). The fold is structurally reactivated by the regional polyphase tectonic evolution, increasing the fracturing of the bedrock. The thermal waters rise to shallow depths in the damage zone of the Daruvar fault and the fracture network in the hinge of the Daruvar anticline. The fold hinge could be affected by a local extensional regime increasing both the aperture of the fractures and the permeability field of the aquifer. The main outflow occurs within the interaction zone of local scale faults/fractures mapped through geophysical investigations (KOSOVIĆ et al., 2023). This structure further localises the flow of the thermal waters resulting in four thermal springs (Fig. 3) with temperatures between 38 and 50°C. Currently, the thermal waters are exploited from a well and two springs providing approximately 10 L/s. 3. METHODS 3.1. Structural-geological investigations and analysis Geological field investigations were carried out from 2021 to 2023. Structural and geological data were collected NE of Daruvar in the Batinjska Rijeka quarry (BRQ) and the surrounding area (Fig. 2). This area was selected as it provides an extensive outcrop of the Mesozoic carbonate rock complex and it can be considered as an outcrop analogue of the Daruvar thermal aquifer. Data were collected at 59 locations and were spatially georeferenced using GIS software (Fig. 4A). Structural measurements included dip direction and dip angle of bedding, fractures, and faults. For faults, lineation and sense of slip, i.e. striations and slickenside orientations with their sense of movement, were also recorded. The resulting cumulative dataset incorporated 109, 284, and 50 measurements of bedding, fractures, and faults, respectively. For the structural data analysis, the software Stereonet v.11 was used to plot bedding and fracture/joint data (ALL- MEN DINGER et al., 2011; CARDOZO & ALLMENDINGER, 2013). Determination of the most representative fracture sets was conducted by plotting the poles of fracture planes and constructing the contour plots of the poles’ distribution using the 1% area contour method (ALLMENDINGER et al., 2011; CARDOZO & ALLMENDINGER, 2013). Fault planes and shear planes were used to determine the palaeostress field based on the geometric properties of the fault and its kinematics G eo lo gi a C ro at ic a 104 Geologia Croatica 77/2 (DOBLAS, 1998). Based on kinematic criteria, the fault data were analysed by the software Win-Tensor v. 5.9.2 and TectonicsFP v. 1.7.9 (ORTNER et al., 2002; DELVAUX & SPERNER, 2003). Fault data were separated into compatible fault groups, and theoretical maximum (σ1), mean (σ2), and minimum stress axes (σ3) were calculated using the P–T axis method (TURNER, 1953; MARRETT & ALLMENDINGER, 1990). The palaeostress field was analysed by constructing the synthetic focal mechanisms for each fault group using the Right Dihedra Method (ANGELIER & MECHLER, 1977). 3.2. Digital structural analysis and DFN modelling The reconstruction of the network of discontinuities for quantifying the hydraulic parameters of the aquifer analogue was conducted through: i) imaging of the BRQ obtaining a 3D digital outcrop model (DOM) of the quarry, ii) digital structural analysis of the DOM reconstructing the main systems of discontinuities (i.e., both bedding and fractures), iii) statistical analysis of the main geometrical parameters of the systems, and iv) discrete fracture network (DFN) modelling of the discontinuity network for the hydraulic parameterisation. The imaging of BRQ was performed in April 2023 using DJI Matrice 300 RTK unmanned aerial vehicle (UAV). The UAV was equipped with a Zenmuse P1 camera featuring a resolution of 45 megapixels, a full-frame size sensor, and a 35 mm focal length lens. The location of the UAV was determined with centimetre accuracy during the whole flight using the RTK (real-time kinematics) method and employing a supplementary GNSS device on a fixed position. Due to the complex quarry geometry, manual flying was used to control the UAV while capturing images resulting in a camera-wall distance of 15 to 65 m. 702 images were acquired with an image overlap up to of 90%. They were processed into a 3D DOM of the quarry (Fig. 4B) by means of “Structure from motion” algorithms using the photogrammetry software Pix4Dmapper v. 4.7.5 (URL 2). The resulting DOM was in sub-centimetre resolution (0.29-0.56 px/cm) covering a surface of 222,135 m2. The software CloudCompare v. 2.13.0 (URL 3) was employed for the digital structural analysis of the DOM. Sections of the walls (QS2-QS6; Fig. 4B), where discontinuity sets were more visible, were selected for analysis. Discontinuity Figure 3. Map of the Daruvar thermal field area showing the thermal springs (blue dots) and the main wells and exploration boreholes (red circle). Figure 4. (A) Investigated area of the Batinjska Rijeka quarry (BRQ; image source: URL 1) with main faults (acronyms in Fig. 2). Observations (purple points) were conducted within the quarry area and its surroundings. (B) Digital outcrop model of the BRQ showing the sections (polygons) used for the digital structural analysis of the discontinuity systems. The section QS3A (red polygon) was chosen to detail the length of the main sets of discontinuities being the base for the DFN modelling. G eologia C roatica 105Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing mapping was done employing the manual sampling technique “Trace tool” of the plugin Compass in the CloudCompare software (THIELE et al., 2017). This tool was used for the manual selection of the DOM’s points that belong to a discontinuity surface and for the calculation of the 3D plane that best fits the points’ selection. The dip direction and dip of the discontinuity were calculated from the fitted plane. The software Stereonet v.11 (ALLMENDINGER et al., 2011; CARDOZO & ALLMENDINGER, 2013) was used to determine the most representative sets of discontinuities in each section of the quarry following the methodology described in the field structural investigations. The analysis of the DOM showed that the rock mass in level 3 of the quarry was highly fractured. Therefore, the QS3A (Fig. 4B) was selected to detail the geometric analysis of the discontinuity sets, since it is the longest section in level 3 and has the greatest amount of data. The DomStudioOrientation script (URL 4) for clustering analysis and Fisher concentration parameter calculations were used to separate and assess the statistics of the main sets of discontinuities. The code was run using the software Matlab v. 9.11.0 (URL 5). A 2D front-view orthophotograph of the QS3A outcrop was performed and imported into the software QGIS v. 3.34.5 (URL 6) to manually vectorise the 2D traces of the discontinuities. The traces were used to: i) calculate the areal fracture intensity P21 (length of fractures over area) for each discontinuity set, and ii) assess the statistical distribution of the heights of the sets in the outcrop using the FracAbility package (URL 7), which allows fitting a statistical distribution to length data affected by censoring bias due to the finite size of the outcrops by applying survival analysis (BENEDETTI et al., 2024). The data obtained from the structural analyses of the DOM and the QS3A section were used to build the DFN model of the investigated outcrop. The input parameters for each set of discontinuities were: mean dip, mean dip direction, Fisher K, mean length, standard deviation of length, length/height ratio, volumetric fracture density P32 (area of fractures over volume), and aperture. All parameters were determined from the previous analyses except for P32 and aperture, which were calibrated through the DFN modelling. P32 cannot be directly measured and needs to be estimated from P21. A linear relationship is expected between P21 and P32 for domains with similar properties (DERSHOWITZ & HERDA, 1992; MAULDON, 1994; WANG, 2005): P32 C21*P21= (1) Here, P32 was calibrated using the DFNworks Python package (HYMAN et al., 2015) with a procedure similar to that described in ANTONELLINI et al. (2014), KORNEVA et al. (2015) and ZAMBRANO et al. (2016). For each discontinuity set, it is necessary to target at least 3 different increasing P32 and run simulations for each target. The resulting DFN is sliced along the reference outcrop’s mean plane to calculate the corresponding P21 value. The overall mean of the P21 values is the simulated value for the chosen target P32. At the end of each calibration procedure, it is possible to calculate the theoretical linear relationship of the P32/P21 ratio for each fracture set in a given homogeneous volume. This calibrated linear coefficient C21 can be then used to obtain the P32 value corresponding to the P21 measured on the DOM. The calibration of the aperture and the calculation of the distributions of porosity (Φ) and permeability (k) were con- ducted employing the software Petroleum Experts Move v. 2019.1 (URL 8). This software has been profitably used to cal- culate the k field of fault zones in fractured carbonates (ANTO- NELLINI et al., 2014; KORNEVA et al., 2015; ZAMBRANO et al., 2016; GIUFFRIDA et al., 2019; ROMANO et al., 2020; SMERAGLIA et al., 2021; MAMMOLITI et al., 2023). Different approaches used by the software to assign the aperture to the discontinuity were tested. They consist of: i) a constant aperture, ii) an aperture proportional to the root length of the fracture, and iii) an aperture proportional to the length of the fracture. The proportionality constant in the second and third approaches is calculated by the software inputting the average aperture of the fracture set. The software discretises the modelling domain using a grid (i.e., geocellular volume), and then the Φ and the k tensor of the cell are computed. The Φ is obtained as the ratio between the total fracture volume in the cell and the cell volume, while the k tensor calculation is based on the geometric methodology described by ODA (1985). 3.3. Hydrogeological parameterisation The hydrogeological parameterisation of the thermal aquifer was performed by analysing the results of a pumping test conducted on 27 February 2022 in Daruvar. The water was pumped from the well Dar 1 at different flow rates (Q). The water level variation was measured in the wells Dar 1, D 1, and DS 1 and in the spring Antunovo vrelo (Fig. 3) using a HOBO Water level data logger with a range of 9 m. The resolution of the water level measurement was 1 cm, while the temporal interval between the measurements was 60 s in Dar 1 and 10 s in the other objects. Manual measurements of the water level were conducted to constrain the measured data. The flow rate was measured at least three times during every single step. The results of this test are partially reported in URUMOVIĆ et al. (2023) where part of the dataset measured in Dar 1 was used to develop a new methodology for interpreting step- drawdown well tests. In this study, a double approach was used: i) the dataset was interpreted as a pumping test (i.e., exploitation and observation well), and ii) the data collected in Dar 1 were reinterpreted using the classical approach for interpreting step-drawdown well tests. Since the Daruvar thermal aquifer is highly fractured, it can be considered as a quasi-homogenous porous medium at the investigated scale. The aquifer is confined because the water level is generally a few metres below the ground while the top of the aquifer is at least 5 m below the ground level. Therefore, the Theis-Hantush analytical solution (THEIS, 1935; HANTUSH, 1961a, b) for pumping tests in nonleaky confined aquifers was applied. The analysis of the collected data was performed using the software AQTESOLV v. 4.5 (DUFFIELD, 2007). It consisted of manually matching the Theis-type curve to the drawdown (Δ) plotted as a function of time on double logarithmic axes. Afterward, the software optimised the curve through a nonlinear weighted least- G eo lo gi a C ro at ic a 106 Geologia Croatica 77/2 squares parameter estimation, employing the Gauss-Newton linearisation method procedure. The estimated parameters were the transmissivity (T), the storativity (S), the hydraulic conductivity anisotropy ratio (Kz/Kr), and the saturated thickness (b). Besides the pumping test analysis, the data collected in Dar 1 were further investigated as a step-drawdown well test. Details on the interpretations of well tests can be found in (KRUSEMAN et al., 2000). The collected data were used to determine the theoretical drawdown vs. flow rate curve of the well, the well efficiency (W), and its specific capacity (SC = Q/Δ). The specific capacity can be in turn used to determine the T of the aquifer using experimental relationships: 2 1.07T (m /day) = 0.85*SC (2) 2 1.07T (m /s) = 2.39*SC (3) These equations were considered since they were developed, respectively, for a carbonate thermal aquifer in Italy (FABBRI, 1997) and dolomite aquifers in Slovenia (VERBOVŠEK, 2008). A similar methodological approach was profitably used to assess the transmissivity of a thermal aquifer in central Croatia (PAVIĆ et al., 2023). 4. RESULTS 4.1. Structural analysis of field data The Batinjska Rijeka quarry (BRQ) is part of the Petrov vrh structural domain (KOSOVIĆ et al., 2024). Structurally, it is located at the intersection of the E-W striking Gradina fault to the NE and Daruvar fault to the SW (Fig. 4A). These two reverse faults are laterally displaced by a NE-SW striking fault (Fig. 4A) by approximately 570 m. Lithologically, the quarry is predominantly composed of a Middle Triassic dolomite and limestone succession (Fig. 2). At the southern face, interlayers of thin layered shale and siltstone/sandstone occur. Strata orientation and fracture system Measured bedding orientations indicate subvertical bedding within the dolomites and limestone strata (Fig. 5A) that are dipping either NE/SW (29 data) or SE/NW (39 data) at the dip angle >60° (Supplementary Tab. 1). The subvertical bedding in the quarry indicate at least two generations of isoclinal folds (Fig. 6A) that are characterised by fold axes dipping either towards the SE (138°) or the SW (228°) at dip angles of 42° and 16°, respectively (Fig. 6A). If the fold generations are considered as second-order structures, they can be interpreted as related to a first-order folded structure that generally dips towards the N at an angle of >65°. The structural analysis encompassed 284 fractures (Fig. 6B). Three representative fracture systems were delineated (Supplementary Tab. 2 and Fig. 6B). The Fs1 fracture system (59 measurements) includes NNE-SSW striking subvertical fracture sets that dip to either the ESE or the WNW. The Fs2 fracture system counts 42 measured data. With a NNE-SSW strike similar to Fs1, this fracture system dips towards the ESE at an angle of 45°. The fracture system Fs3 (54 data), composed of two subsets, resembles two coupled subvertical fractures that dip either towards the NNE (22°) or the NNW (338°). Faults and shear fracture analysis Fault analysis focused on addressing the structural architecture and identification of the faults and shear fractures within the area. Due to heavily tectonised carbonate quarry faces and the active exploration of the quarry, 50 shear fractures/fault plane data were collected. Considering the stress field and kinematic criteria, three principal fault categories were identified and subdivided by fault groups and subsets (Fig. 6C and Supplementary Tab. 3). Reverse fault planes (20 data; Figs. 6C and 7 and Supplementary Tab. 3) were separated into two fault sets (RFa and RFb). RFa subset shows an average SE dip direction (dip Figure 5. (A) Tectonised Triassic dolomites and dolomitised limestones in the BRQ at location 45.603°N, 17.272°E. Strata are steeply dipping towards the SW (So: 233/70) with a normal fault (198/60) at the base of carbonate beds. (B) Coupled fracture systems in the BRQ at location 45.606°N, 17.273°E. Measured fracture systems dip towards the ESE (Fs2: 100/45), NW (Fs3: 338/73), W (Fs1b: 288/82). G eologia C roatica 107Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Figure 6. (A) Stereonets of the bedding orientations (small black circles) in the BRQ. Their interpretation highlights the occurrence of generally two gen- erations of folds (fold limbs are shown as small blue and red circles) characterised by different fold hinges (blue and red points). (B) Structural diagrams of the main fracture systems (black small circles) in the study area. The orientations of the representative fracture systems (FS) are reported in Supple- mentary Tab. 2. They are obtained from the distribution of the measured discontinuities represented here as poles (black points). (C) Structural diagrams for the interpreted fault systems (RF - reverse fault group; NF - normal fault group; SSF-1 and SSF-2 - strike-slip fault groups). In all plots, N indicates the number of measurements. G eo lo gi a C ro at ic a 108 Geologia Croatica 77/2 angle of 56°), whereas the RFb fault subset dips towards the WNW at an angle of 36°. Structural analysis and computation of the representative palaeostress field using the P–T axes method and derived synthetic structural focal mechanism show that the palaeostress compressional field was associated with a NW-SE trending P-axis (the T-axis is subvertical, steeply dipping towards the NNW), hence NW-SE compression. This compressional palaeostress field may have resulted in the formation of the fold generation dipping towards SE (138°; Fig. 6A) Normal fault planes (10 data; Supplementary Tab. 3 and Fig. 6C) showed orientations similar to the reverse faults. They were divided into two fault sets (NFa and NFb) characterised by NE‒SW striking fault planes dipping towards either the SE or the NW, at angles of 41° and 50°, respectively. Kinematic and stress analysis indicates that observed normal faults were formed due to a subvertical P-axis steeply dipping towards the NNE (P-axis orientation is 14/87) and subhorizontal T-axis trending NW‒SE. This stress field resulted in the NW–SE extension of the existing faulted structures. Besides the reverse and normal faults, 19 strike-slip fault planes were measured (Supplementary Tab. 3; Figs. 6C and 7B). Measured data were separated based on geometry and kinematic compatibility into two fault systems (SSF-1 and SSF-2) and four group subsets. The SSF-1 group included 10 measured dextral/sinistral faults (dip angle between 53° and 73°) striking either NE–SW or NW–SE. The SSF-2 group incorporated conjugate fault pairs with NNE–SSW and SE-NW strike (dip angle between 66° and 71°). Mapped strike- slip fault planes often showed structural reactivation with both dextral/sinistral motions. Similar orientations of two NW-SE striking fault subsets (SSF-1b and SSF-2b) indicate fault structural reactivation and reorientation of the principal stress axes σ1 and σ3 within a relatively short timeframe. Both strike-slip fault groups were formed in two slightly different palaeostress fields. The SSF-1 fault group indicates a palaeostress field associated with the NNE–SSW trending P-axis, whereas SSF-2 fault group indicates a palaeostress field associated with an ENE-WSW trending P-axis. This implies that both strike-slip fault groups were tectonically active within a general NE-SW transpressional stress field. 4.2. Digital structural analysis of the DOM The discontinuities digitalised in the selected sections of the DOM are shown in Fig. 8. The contour plots of the poles were used to identify the main sets of discontinuities in each section (Supplementary Tab. 4). Quarry section 2 (length: 145 m, height: 17 m, represen- ta tive dip direction and dip: 282/68) is located on the lowest level of the quarry. The DOM analysis reveals that three discontinuity systems crosscut the rock mass. The discontinuity system D1 is composed of two NE-SW striking subvertical discontinuities dipping at high angles (76° and 83°, respectively) toward the NW (D1a; 20 data) or ESE (D1b; 17 data). The D2 system is divided into D2a and D2b sets (23 and 15 data, respectively) and includes predominantly NNW-SSE striking discontinuities dipping towards the WSW or W with dip angles of 53° and 57°, respectively. The D3 system (21 data) is characterised by discontinuities that are steeply dipping at high angles (≥87°) towards the W. Quarry face 3 is divided into sections A (length: 97 m, height: 30 m, representative dip direction and dip: 259/70) and B (length: 63 m, height: 31 m, representative dip direction and dip: 297/68). Analysis of the DOM in section 3A shows two discontinuity systems D1 and D2. The D1 system is characterised by two discontinuity sets (D1a and D1b; 132 and 43 measurements, respectively) dipping towards the WNW with dip angles of 68° and 83°. The D2 system comprises two sets dipping towards the WSW but at different dip angles. D2a (140 data) is steeply dipping towards the WSW at 75°, while set D2b (74 data) is gently dipping with a dip angle of 41°. The structural fabric of section 3B is similar to section 3A being deformed by two discontinuity systems (D1 and D2). The D1 system is characterised by NE-SW striking discontinuity sets. Figure 7. (A) Fault plane (167/45) observable across most of the BRQ faces, at location 45.606°N; 17.273°E. Fault plane shows reverse kinematics with indication of tectonic transport to the NW, while in the same time this fault plane experienced structural reactivation and tectonic inversion characterised by NW-SE oriented extension. (B) Reverse fault planes with average orientation 142/30 (Fp). Besides NW tectonic transport, we observed subvertical bedding with average strata orientation 318/75. In the right section, a subvertical bedding plane (Fp2) was structurally reactivated as a dextral fault. G eologia C roatica 109Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Set D1a (124 data) dips towards the NNW (dip angle of 70°), whereas the D1b set (73 data) includes discontinuities that steeply dip towards the NW at a dip angle of 77°. The second discontinuity system (D2; 115 data) encompasses the NNW- SSE striking discontinuities that dip steeply towards the WSW at an angle of 72°. Quarry face 4 (length: 108 m, height: 30 m, representative dip direction and dip: 290/66) is similar to quarry face 3 being deformed by two discontinuity systems D1 and D2. The D1 system is characterised by two sets (D1a and D1b; 103 and 86 data, respectively) steeply dipping towards the NNW or NW. The D2 system (89 measurements) includes predominantly NNW-SSE striking discontinuities steeply dipping towards the WSW (dip angle of 72°). Quarry section 5 (length: 78 m, height: 18 m, representative dip direction and dip: 286/69) is deformed by two discontinuity systems (D1 and D2) that dip towards either the NW or W at angles of 71° and 72°, respectively. Two discontinuity systems D1 and D2 were delineated in the highest face of the quarry (section 6; length: 53 m, height: 11 m, representative dip direction and dip: 298/65). The D1 system includes the D1a and D1b sets (16 and 9 data, respectively), characterised by discontinuities steeply dipping towards the NNW with a dip angle of approximately 70°. The D2 system encompasses NE-SW striking discontinuities that are steeply dipping towards the W at an angle of 60° (D2a set; 5 measurements) and 81° (D2b; 8 data), respectively. Based on the structural analysis of the DOM, two subvertical systems (D1 and D2) of discontinuities can be distinguished. The D1 system generally dips toward the NW, while the D2 system shows a general NNW-SSE strike dipping towards the WSW. The section of the quarry QS3A (Fig. 4B) was selected for a detailed analysis of the discontinuity sets to derive the input parameters for the discrete fracture network (DFN) modelling. The high-resolution orthophotograph of the section was digitised in QGIS obtaining a shapefile of discontinuity traces, i.e. the intersections of the discontinuity planes with the quarry face (Fig. 9). Since the face dips 70° towards W, it intersects all the discontinuity systems identified in this sector of the quarry (Fig. 8). Sets D1a and D1b (subvertical in the orthophoto; red set in Fig. 9) differ by just 20° in strike and it was impossible to separate them based on their traces on the orthophotograph. They were merged into one set (named set V). The traces of systems D2a and D2b (at medium and low angle in the orthophoto; green set in Fig. 9) show a continuous range of plunges, without a clear cutoff angle between them. The medium to low angle traces were therefore considered as a merged system (named set H). The parameters calculated for the two merged sets visible in section QS3A are listed in Supplementary Tab. 5. The merged sets show a dispersed distribution of their spatial orientation, as evidenced by the relatively low values of the Fisher K (13.53 for set H and 18.83 for set V). Since section QS3A is subvertical, the length of the traces corresponds to fracture heights. The percentage of censored heights is 20% and 11.08% in sets V and H, respectively (Supplementary Tab. 5). A log-normal distribution best fits the distribution of the measured heights, accounting for censoring bias (Supplementary Fig. 1). The Figure 8. Structural diagrams of the main discontinuity systems in the different sections of the quarry (N: number of data). The orientations of the rep- resentative discontinuity systems are reported in Supplementary Tab. 4. They are obtained from the distribution of the measured discontinuities in the DOM of BRQ here represented as poles (black hollow points). G eo lo gi a C ro at ic a 110 Geologia Croatica 77/2 cumulative fracture lengths were used to calculate the areal fracture intensity P21 being 0.47 and 0.37 for set H and V, respectively. 4.3. DFN modelling of the aquifer outcrop analogue The results obtained from the DOM and QS3A virtual structural analyses were used to build the DFN model of the outcrop. The model size was 100×100×30 m (XYZ; Fig. 10A) corresponding to the length and height of the QS3A section. The input parameters used for the construction of DFN models in both the software DFNworks and Move are reported in Supplementary Tab. 6 including the P32 and aperture values obtained after the calibration procedure. The calibration of P32 was based on the iterative building of DFNs in the DFNworks software, pointing to a ratio C1 (Eq. 1) of 0.24 and 0.21 for set H and set V, respectively. These linear coefficients were used to calculate the P32 from the P21 obtained through the digitalisation of the fracture traces in the QS3A section resulting in P32 values of 1.98 and 1.18 m2/m3 for set H and set V, respectively. The aperture was calibrated to fit the porosity (Φ) distribution in the DFN model with the Φ field measured during well logging of Dar 1 well. The analysis of the well log showed a large variation in the Φ ranging from 0.03 to 42% (average = 7.9%). The occurrence of outliers was indicated by the fact that the 10th and 90th percentiles of the distribution are 0.9 and 17%, respectively, suggesting a narrower Φ field. A closer inspection of the caliper well log showed a variation in the diameter of the borehole from 35 to 45 cm (nominal diameter of the drill bit was 31.1 cm), with the high Φ occurring in the well sections with larger diameters. These data suggest that the wall of the well locally collapsed (most probably in the more fractured part of the reservoir), affecting the measurement of the Φ. The section of the well with both a smaller and more homogeneous diameter was selected to obtain a more realistic assessment of the Φ. These conditions should point to a more compact wall of the well, resulting in Φ values that should better reflect the Φ field in the aquifer. The selected section was located at a depth of 100-130 m showing a homogenous diameter between 35 and 36 cm. The Φ within this interval is between 0.03 and 9.1% with an average value of 2.7%. The 10th and 90th percentiles of the distribution are 0.2 and 6.3%, respectively, suggesting a more homogeneous distribution around the average than the whole Φ dataset. The Φ calibration was conducted testing values of the fracture aperture from 0.1 to 10 mm and combining them with the different approaches for assigning the aperture in the Move software. The modelling domain was discretised using a grid spacing of 0.5 m considering that the Φ in the Dar 1 well was calculated from the neutron log generally having a vertical resolution of 1 m and depth of investigation of 25 cm (LAI et al., 2024). The distributions of Φ and k obtained for the tested aperture values and employing the different methodological approaches are shown in Supplementary Figs. 2 and 3, respectively, while their statistical analyses are in Supplementary Tabs. 7 and 8. A general increase of both the Φ and k fields with the aperture was observed except for the minimum and 10th percentile of the distributions, which were always 0. This value is the outcome of choosing a small grid size since the software assigns null Φ/k to those cells that are not cut by any structure. Considering the same aperture, the highest values were obtained with the method “aperture proportional to length”, while the method “constant aperture” showed the lowest values. This discrepancy can be explained by the different methodological approaches for calculating the aperture. In the “aperture proportional to length” approach, long fractures have larger apertures than short fractures resulting in both higher Φ and k values. In the “constant aperture” approach, the same aperture is applied to all fractures diminishing the effect of long and potentially wider fractures on the hydraulic parameters. This effect is also visible by analysing the skewness and kurtosis values of the distributions. Skewness and kurtosis are higher than 0.9 and 3.8, respectively, suggesting right-skewed and leptokurtic distributions. These data suggest a clustering of the datasets toward low values with a long right tail toward high values. The highest skewness and kurtosis were observed in the “aperture proportional to length” approach enforcing the Figure 9. Orthophotograph of section QS3A digitised in QGIS to vectorise the traces of the discontinuity systems. Traces are divided into two disconti- nuity systems: subvertical (red; set V in Supplementary Tab. 5) and subhorizontal (green; set H in Supplementary Tab. 5). G eologia C roatica 111Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing occurrence of high Φ/k values associated with wider fractures. It should be noted that a similar asymmetric distribution was depicted in the experimental Φ of the Dar 1 well. The statistical analysis of the Φ in the DFN models was further detailed to calibrate the aperture. A linear relationship (Fig. 10B) between the aperture and all statistical indexes is observed, except for the minimum and 10th percentile values. The slopes of the average values for the different methodological approaches were calculated to obtain the aperture fitting with the average Φ measured in Dar 1 well at a depth of 100–130 m. The resulting apertures were 8.9, 5.1, and 4.1 mm for the constant, root length, and length approaches, respectively. The aperture values obtained with the second and third approaches seem more realistic and applicable for a fractured aquifer, while wider apertures could be more representative of karst aquifers. The model was run with an aperture of 4.1 mm and employing the “aperture proportional to the fracture length” approach. The resulting Φ ranged from 0 to 22.7%, but the statistical analysis showed that the Φ distribution (red dots in Fig. 10C) fits the experimental values at the depth of 100–130 m (red lines in Fig. 10C) with only a minor discrepancy observed for the 75th percentile. The resulting k value ranged from 0 to 4.4×106 D and an average value of 2.7×105 D. The statistical distributions of both Φ and k are right-skewed and leptokurtic (skewness of 1.1 and 2.4, respectively, and kurtosis of 4.3 and 11.2, respectively) indicating a highly asymmetric distributions. Figure 10. (A) DFN model obtained using the input parameters in Supplementary Tab. 5. The red and green discontinuities represent the vertical and horizontal discontinuity sets, respectively. (B) Mean (black lines) and maximum (red lines) porosity values obtained using different approaches for as- signing the aperture to fractures in the DFN model (dotted lines: constant aperture; dashed lines: aperture proportional to the root of the length; full lines: aperture proportional to the length). (C) Porosity distribution (black lines) obtained using the “aperture proportional to length” approach. It was considered the best approach for fitting the experimental porosity distribution (dashed red lines) of the Dar 1 well at a depth of 100-130 m. The results obtained from the best model with aperture of 4.1 mm are shown as red dots suggesting a good fit with the experimental data. G eo lo gi a C ro at ic a 112 Geologia Croatica 77/2 4.4. Interpretation of pumping and well tests The pumping test in the Daruvar thermal field lasted approxi- mately 12 h. It was conducted performing 3 steps (i.e., step- drawdown test) with flow rates of 12.1, 13.4, and 7.2 L/s, respectively, lasting 1 h each (Fig. 11A). The pumping was stopped for 1 h measuring the recovery of the water level. Afterwards, a flow rate of 12.3 L/s was set for approximately 7 h (i.e., constant test). The subsequent recovery of the water level was measured for 1.5 h reaching the initial water level. A maximum drawdown (Δ) of 0.59 m was observed in the Dar 1 well at the end of the second pumping step, while the highest Δ in the observation wells D 1 and DS 1 was measured during the constant test reaching 0.09 m (Fig. 11A). The water level in the Antunovo vrelo spring was slightly affected by the pumping showing a Δ up to 0.01 m. The Δ values measured in Dar 1, D1, and DS 1 were used for determining the hydrogeological parameters of the aquifer using both the Theis-Hantush analytical solution for pumping tests. The solution was calculated using the data from the whole pumping test (i.e., step-drawdown and constant tests; Fig. 11B) and only from the constant test (Fig. 11C). The obtained parameters are summarised in Supplementary Tab. 9. The results of the two methodological approaches are com parable. The transmissivity (T) is 3.82×10–2 m2/s and 3.46×10–2 m2/s, while the saturated thickness of the aquifer (b) is 61.5 and 68.4 m for the whole and constant tests, respectively. The saturated thickness value slightly corresponds to the length of the screened section of the Dar 1 well (53 m, and up to 59 m if a short part of casing between the two well screens is included). Furthermore, the data measured in Dar 1 were used to assess the drawdown vs. flow rate curve of the well, the well efficiency (W), and the specific capacity (SC). The obtained results are summarised in Supplementary Fig. 4. The calculated drawdown vs. flow rate curve is: 21 4.6 2186.6 Q Q∆ = + (4) W is between 33% and 48% and decreases by increasing the flow rate. Although W is low showing a poor efficiency, it is acceptable considering that the Dar 1 well is relatively old and it has practically never been used. The SC was calculated considering both the experimental and theoretical Δ showing minor differences. The resulting T values are 3.8×10–2 m2/s and 4.8×10–2 m2/s for the equations 2 and 3, respectively. Considering that these values are obtained using experimental relationships developed in geological and hydrogeological settings different from the Daruvar thermal field, the slight discrepancy with the T calculated from the pumping test approach is acceptable. 5. DISCUSSION 5.1. Integrated interpretation of field and virtual structural data Fractures are the main f low path in fractured aquifers (FAULKNER et al., 2010; BENSE et al., 2013). Reconstructing the geometric characteristics of the fracture network and the stress field deforming the rock mass can be useful in determining the preferential flow paths fostering the conceptual Figure 11. (A) Drawdown measured in the wells Dar 1, D 1, and DS 1 and the Antunovo vrelo spring during the pumping test conducted in 2022. The flow rate is shown as a step plot with a dashed line. (B-C) Drawdown data for Dar 1 (black circle), D 1 (red cross), and DS 1 (blue square) fitted with the theoretical Theis curves (black, red, and blue lines for Dar 1, D 1, and DS 1, respectively) for the whole pumping test (B) and the constant flow rate test (C). G eologia C roatica 113Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing model of the groundwater resource (SINGHAL & GUPTA, 2010). Integrated structural mapping based on field mea- surements and digital outcrop model (DOM) analysis was employed in this work to analyse the fracture network in the Batinjska Rijeka quarry (BRQ). The results of the field measurements indicate that the Mesozoic carbonate complex in the quarry is deformed by systems of folds and fractures being a representative outcrop analogue for the fractured carbonate reservoir of the Daruvar hydrothermal system. Folded structures in the BRQ are characterised by subvertical bedding (dipping either NE/SW or SE/NW; Fig. 6A) that have undergone extensive brittle deformation with the formation of three systems of fractures (Fs1 and Fs3 including fractures steeply dipping toward ESE or NNE, and Fs2 including low angle fractures dipping toward the ESE; Fig. 6B). The geometry, structural position, and subvertical orientation of these systems indicate a cogenetic relationship with the fold hinge zones that have an predominantly NE-SW and ESE- WNW strike. Structural observations of shear fractures/fault planes (Fig. 6C) corroborate the results of the bedding and fracture system analyses, suggesting that the present-day fracture pattern observed in the BRQ results from the multiphase regional tectonics (KOSOVIĆ et al., 2024). Cretaceous-Palaeogene E-W (or NE-SW) compression rotated by 40º counterclockwise is observed on NNE-striking reverse fault planes associated with a NW-SE oriented P-axis (Supplementary Tab. 3 and Fig. 6C). NE-striking normal faults indicate a NW-SE oriented extension and clockwise rotation of the existing NNE-striking reverse fault due to a tectonic inversion. The observed extension may be related to the collapse of the hanging wall structures of the Cretaceous- Palaeogene reverse faults due to either gravitational slip (e.g. TAVANI et al., 2012) or the initial extension of the PBS in the late Oligocene and Miocene (TOMLJENOVIĆ & CSONTOS, 2001). The N-S (NE-SW) compressional/transpressional phase from the Pliocene-Quaternary is observed on strike-slip fault planes (Supplementary Tab. 3), showing polyphase structural reactivation and repeated movements in both a dextral and sinistral sense. Considering the cross-cutting relationships between SSF-1 and SSF-2 fault systems and mapped reverse/normal faults, it indicates that the strike-slip fault planes are younger and probably associated with the most recent Pliocene-Quaternary tectonic phase in the study area (e.g., TOMLJENOVIĆ & CSONTOS, 2001; SCHMID et al., 2008). This NE-SW oriented stress field may also influence the bedding reorientation and the formation of a younger fold system that gently dips towards the SE (Fig. 6A). Results of the DOM virtual structural analysis depicted two dominant subvertical systems (D1 and D2; Supplementary Tab. 4) of discontinuities. The D1 system generally dips toward the NW, while the D2 system dips towards the WSW. These results correspond well with the field observations (Supple- mentary Tabs. 1 and 2). The D1 system has the same orientation as the bedding So3 and the fracture set Fs3b. Discontinuity system D2 has a similar orientation to the bedding So2 and fracture set Fs1b. Other sets of bedding and fractures observed by field measurements don’t clearly correspond with the discontinuity systems delineated in the DOM. Only the fracture set Fs1a (Supplementary Tabs. 1 and 2) could have an analogy with the discontinuity set D1b measured in quarry section 2 (Supplementary Tab. 4). However, this discontinuity set is depicted only in one section of the quarry and is obtained from a limited dataset, decreasing the reliability of its geometric characterisation in the DOM. The lack of correlation between field data and the DOM results could be associated with several factors. Fracture sampling in the DOM highly depends on the outcrop inclination and the fracture 3D spatial orientation (PRIEST, 1993; ZHANG & EINSTEIN, 1998; ZEEB et al., 2013). This difference significantly reduces the identification of potential fracture systems that are either parallel or steeply dipping in opposite directions to the quarry face during the DOM analysis. Furthermore, field investigations and virtual outcrop analysis have different resolution scales (TZIAVOU et al., 2018; MARTINELLI et al., 2020): disconti- nuities at the centimetre scale can be measured in the field, while they are impossible to digitalise in decametre scale DOMs. Bearing in mind the different resolutions of the methods, only metre scale discontinuities were measured during the field investigations in the BRQ. Still, a discrepancy between the results occurred highlighting the intrinsic difference between the two methodological approaches. However, the digital outcrop analysis of the BRQ covered a relatively large area favouring: i) reconstruction of the discontinuity network in remote or inaccessible portions of the outcrop, and ii) collection of a large dataset that was used for a detailed statistical analysis of the geometric features of the main discontinuity sets. 5.2. Hydrogeological parameterisation of the Daruvar thermal aquifer Besides their primary role in the circulation of groundwater, networks of fractures can influence the distribution of the bulk permeability (k) and porosity (Φ) in fractured aquifers (e.g., FAYBISHENKO et al., 2000; SINGHAL & GUPTA, 2010; BENSE et al., 2013). In this study, both classical hydrogeological investigations and discrete fracture network (DFN) modelling were conducted for the hydraulic parameterisation of the Daruvar carbonate thermal reservoir. Several pumping and well tests have been conducted in the Daruvar thermal field to assess the transmissivity (T) of the aquifer (BOROVIĆ, 2015; BOROVIĆ et al., 2019; URUMOVIĆ et al., 2023). In particular, BOROVIĆ (2015) critically revised and reanalysed the results from old pumping tests (1971, 1980, 1982, 2009), providing a solid estimation of T. Here, the results of a recent pumping test (2022) were analysed using several methodological approaches (THEIS, 1935; HANTUSH, 1961a, b; KRUSEMAN et al., 2000). The T values obtained in this work and from the literature range from 1.58×10-2 to 4.83×10-2 m2/s showing a good concordance among the different methods (Tab. 1). The average T was used to calculate the aquifer k value considering: i) the lower density and viscosity of the thermal waters due to their temperature and salinity, and ii) both the thicknesses of the dolomite aquifer and the fracture corridors in the Dar 1 well (d and f, respectively; Tab. 1) and the saturated thicknesses calculated from the pumping test analysis (b; Tab. 1). The resulting average values of k (Tab. 1) are 16.4, 32.7, and 88.7 D considering d, b, and f as the aquifer thickness, respectively. G eo lo gi a C ro at ic a 114 Geologia Croatica 77/2 These values show good consistency within the same order of magnitude and are within the range of k values for fractured and karst carbonate aquifers (FREEZE & CHERRY, 1979). The Φ of the aquifer was obtained through well logging of the Dar 1 well. Since the distribution of Φ values showed a large variability (Tab. 1), a section of the well with tight walls (depth of 100–130 m) was considered to obtain a narrower, more realistic, range of Φ. The average Φ for this section was 2.7% (Tab. 1) being within the range of Φ in carbonates (NG & SANTAMARINA, 2023). The distribution of Φ at a depth of 100–130 m was used for the calibration of the DFN model of the BRQ outcrop analogue. The best fit of the model was obtained using an aperture of 4.1 mm resulting in an average Φ of 2.7% and an average k of 2.7×105 D. Although the Φ fits with the experimental values, the k value is 3-4 orders of magnitude higher than the measured k (Tab. 1). The reasons for this discrepancy could be that: i) the Φ used for the calibration corresponds to the total Φ, while effective Φ is relevant for the fluid flow, and ii) the size and the discretisation of the DFN model do not correspond to the area investigated by the pumping test and the hydrogeological representative elementary volume. In order to overcome these issues, a larger DFN model was conducted. The model size was increased to 700×700×150 m (XYZ; Fig. 12A) being more representative of the aquifer volume investigated by the pumping test. It was considered that the radius of influence of the pumping test was between the DS 1 and D 1 observation points and the Antunovo vrelo spring, because the experimental data showed variations of the water level in the observation wells (Figs. 3 and 11A). Therefore, the horizontal size of the model was calculated to be representative of twice the average distance between the exploitation well Dar 1 and the observation points (Dar 1 - D 1: 265 m, Dar 1 – Antunovo vrelo: 450 m; Fig. 3). The vertical size was chosen to be representative of the maximum investigated thickness of the dolomite in the Daruvar area (d in Tab. 1). The model was conducted using the software Move and it was discretised using: i) a geocellular volume with a cell size as the DFN model of the outcrop analogue for generating the fracture network, and ii) a single cell (i.e., the whole volume) for the hydraulic parameterisation. This “double-sized” geocellular volume permits: i) construction of the fractures using the input parameters of the outcrop DFN model (Supplementary Tab. 6), and ii) calculation of the hydraulic parameters over a volume that is representative of the investigated section of the aquifer. Through this approach, the fracture geometry at the aquifer scale DFN is constrained using the geometric parameters obtained from the DOM virtual structural analysis, and the bulk hydraulic parameters are calculated obtaining an estimation at the scale of the pumping test. The model was run using the “aperture proportional to the fracture length” approach since it seemed more reliable for fractured aquifers. Different aperture values were tested obtaining a linear and power correlation with Φ and k (Figs. 12B-C). The aperture of 4.1 mm (i.e., the best fit in the outcrop analogue DFN model) results in a Φ and k of 3.7% and 3.9×105 D, respectively. These values are slightly higher than the results of the outcrop scale DFN model, probably due to the different geocellular volume used for the calculations. The linear relation- ship between the aperture and Φ was used to recalibrate the aperture value fitting the average Φ at the depth of 100–130 m (Tab. 1). The resulting aperture of 3 mm was run in the model obtaining a k value of 1.5×105 D, which is 3–4 orders of magni- tude higher than the experimental values. This result suggests once more that two different aperture values are needed for fitting the experimental Φ and k values. In particular, the hydra- ulic aperture of the fractures, which is smaller than the mecha- nical aperture (RENSHAW, 1995; POURASKARPARAST et al., 2024), should be used for reconstructing the k value in DFN models (e.g., SMERAGLIA et al., 2021). Furthermore, the effective Φ should be considered since fractured aquifers usually have dual porosity, with most of the storage volume being in the matrix and most of the flow being in the fracture network (WORTHINGTON et al., 2019). The Φ dataset used for the calibration of the DFN models was calculated from the neutron log and describes the total Φ of the Daruvar aquifer. Considering that effective Φ in carbonate rocks is at least one order of magnitude less than total Φ (WORTHINGTON et al., 2019), the fracture aperture for obtaining a Φ of 0.27% (i.e., one order of magnitude lower than the average Φ at the depth of 100-130 m; Tab. 1) was calculated and input to the model. The calculated fracture aperture was 0.3 mm resulting in a k value of 149.4 D, which is slightly higher than the experimental values. The effective Φ value was refined on the basis of the experimental Φ dataset, and the 10th percentile of the Φ distribution was tested (Tab. 1). The corresponding fracture aperture (0.22 mm) was input to the model resulting in Φ and k values of 0.2% and 60.5 D, respectively. The k value fits with the experimental values being in the range of the k calculated from f and slightly higher than the values obtained from d and b (Tab. 1). The results of the conducted numerical investigations highlight the importance of employing a “dual aperture” approach for fitting both experimental Φ and k values in DFN models. This approach is necessary since fractured aquifers are a dual continuum medium where the groundwater is Table 1. Hydrogeological parameterisation of the Daruvar thermal aquifer obtained from: i) the interpretation of pumping and well tests in this study, ii) the well logging of Dar 1 well, and iii) the literature (BOROVIĆ et al., 2019; URUMOVIĆ et al., 2023). Average value Minimum – Maximum 10th – 90th percentile Investigated thickness of dolomite d (m) 130 – – Saturated thickness b (m) 65 61.5 – 68.4 – Thickness of fracture corridors f (m) 24 Transmissivity T (m2/s) 3.5e-2 1.6e-2 – 4.8e-2 2.5e-2 – 4.3e-2 Storativity S (–) 1.6e-3 1.1e-3 – 2e-3 – Hydraulic conductivity anisotropy ratio Kz/Kr (–) 1.3e-2 6.8e-3 – 2e-2 – Permeability from d k d (D) 16.4 7.4 – 22.7 11.8 – 20.3 Permeability from b k b (D) 32.7 14.8 – 45.3 23.6 – 40.6 Permeability from f k f (D) 88.7 40 – 122.8 63.8 – 109.9 Porosity Φ (%) 7.9 0.03 – 42 0.9 – 17.4 Porosity at depth 100-130 m Φ* (%) 2.7 0.03 – 9.1 0.2 – 6.3 G eologia C roatica 115Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing mostly stored in the rock matrix, affecting the total porosity of the aquifer, while the fluid flow occurs in the fractures, being influenced by their hydraulic aperture and the effective porosity of the aquifer. This difference should also be considered for the analytical interpretation of pumping tests. Here, a generally well accepted analytical method reproducing the aquifer as an equivalent porous medium was used (SINGHAL & GUPTA, 2010), but it provides only a rough estimation of the k value. The estimation of k in dual porosity media requires more sophisticated field investigations (i.e., flowmeter and quantitative visual logging, tracer tests) and analytical methods for modelling the pumping test results (i.e., MOENCH, 1984). 6. CONCLUSIONS This paper describes the results of an integrated approach for reconstructing the network of fractures and quantifying the hydraulic parameters in fractured aquifers. Structural and hydrogeological investigations were conducted in the Daruvar thermal area and its hinterland detailing the structural architecture and the hydraulic parameterisation of the carbonate thermal aquifer. The results obtained through field and virtual structural analyses were used as input parameters to construct discrete fracture network models at different scales to estimate the aquifer properties. The modelling results were compared with the results of hydrogeological field investigations obtaining a good correlation. These results permit the quantification of the impact of the fracture network on the permeability field of the Daruvar thermal aquifer as observed in fractured aquifers. The results obtained in the investigation of the Daruvar thermal area demonstrate that an integrated approach is invaluable in investigating and modelling fractured media. The utilisation of such an integrated approach permits a Figure 12. (A) Aquifer scale DFN model obtained using the input parameters in Supplementary Tab. 5. The red and green discontinuities (1% of the generated discontinuities is shown) represent the vertical and horizontal discontinuity set, respectively. (B-C) Porosity and permeability values (black dots) obtained from the aquifer scale DFN using different fracture apertures. The results of the best model (aperture of 0.22 mm) are shown as red dots pointing to a good fit with the experimental data (dashed lines). G eo lo gi a C ro at ic a 116 Geologia Croatica 77/2 quantitative structural and hydrogeological reconstruction of the aquifer and the network of fractures favouring the fluid flow accounting for the advantages of both techniques. Field analysis can be used to determine the geological relationship among structures and the stress field developing the dis conti- nuities. Remote sensing photogrammetry allows the construc- tion of a 3D digital model of the outcrop favouring the collection of large datasets for a detailed quantitative description and statistical analysis of geological objects' geometry and their spatial relationships (HODGETTS, 2013; BISTACCHI et al., 2015; MARTINELLI et al., 2020). Field research can add a conceptual link to the results of the data analysis in the DOM obtaining a more comprehensive geological interpretation. Furthermore, the obtained results highlight the importance of correlating the results from structural and hydrogeological investigations. Structural data can be used to construct discrete fracture network models of the fractured aquifer that need to be calibrated using the data from hydrogeological investigations. Depending on the available experimental dataset of hydrogeological properties, the model calibration needs to account for the dual porosity of fractured aquifers and the difference between the mechanical and hydraulic aperture of the fractures. ACKNOWLEDGMENT The authors would like to thank the editorial team for the support during the preparation and the revision of the manuscript, the two anonymous reviewers for their useful comments, the Daruvarske toplice–Special Hospital for Medical Rehabilitation and Pivovara Daruvar d.o.o. for logistic help on-site and sharing of existing materials and Matej MANJARIĆ and INTER-PROMET d.o.o. for allowing the access to the Batinjska Rijeka quarry and the logistic support during the investigations. This research was funded by the HyTheC project of the Croatian Science Foundation (HRZZ), grant number UIP-2019-04-1218. REFERENCES ALLMENDINGER, R.W., CARDOZO, N. & FISHER, D.M. (2011): Struc- tural geology algorithms: Vectors and tensors.– Cambridge University Press, 289 p. doi: 10.1017/CBO9780511920202 ANGELIER, J. & MECHLER, P. (1977): Sur une methode graphique de re- cherche des contraintes principales egalement utilisables en tectonique et en seismologie: la methode des diedres droits.– Bulletin de la Société Géologique de France, S7-XIX(6), 1309–1318. doi: 10.2113/gssgfbull.S7- XIX.6.1309 ANTONELLINI, M., CILONA, A., TONDI, E., ZAMBRANO, M. & AGO- STA, F. (2014): Fluid flow numerical experiments of faulted porous car- bonates, Northwest Sicily (Italy).– Marine and Petroleum Geology, 55, 186–201. doi: 10.1016/j.marpetgeo.2013.12.003. BADA, G., HORVÁTH, F., DÖVÉNYI, P., SZAFIÁN, P., WINDHOFFER, G. & CLOETINGH, S. (2007): Present-day stress field and tectonic in- version in the Pannonian basin.– Global and Planetary Change, 58/1–4, 165–180. doi: 10.1016/j.gloplacha.2007.01.007 BALEN, D., HORVÁTH, P., TOMLJENOVIĆ, B., FINGER, F., HUMER, B., PAMIĆ, J. & ÁRKAI, P. (2006): A record of pre-Variscan Barrovian re- gional metamorphism in the eastern part of the Slavonian Mountains (NE Croatia).– Mineralogy and Petrology, 87/1–2, 143–162. doi: 10.1007/ s00710-006-0120-1 BENSE, V. F., GLEESON, T., LOVELESS, S. E., BOUR, O. & SCIBEK, J. (2013): Fault zone hydrogeology.– Earth-Science Reviews, 127, 171–192. doi: 10.1016/j.earscirev.2013.09.008 BENEDETTI, G., CASIRAGHI, S., BISTACCHI, A., & BERTACCHI, D. (2024): FracAbility: A python toolbox for survival analysis in fractured rock systems.– In: EGU General Assembly, Vienna, Austria, 14–19 Apr 2024. EGU24-22156. doi: 10.5194/egusphere-egu24-22156, 2024. BISTACCHI, A., BALSAMO, F., STORTI, F., MOZAFARI, M., SWENNEN, R., SOLUM, J., TUECKMANTEL, C. & TABERNER, C. (2015): Pho- togrammetric digital outcrop reconstruction, visualization with textured surfaces, and three-dimensional structural analysis and modeling: In- novative methodologies applied to fault-related dolomitization (Vajont Limestone, Southern Alps, Italy).– Geosphere, 11/6, 2031–2048. doi: 10.1130/GES01005.1 BISTACCHI, A., MITTEMPERGHER, S., MARTINELLI, M. & STORTI, F. (2020): On a new robust workflow for the statistical and spatial analysis of fracture data collected with scanlines (or the importance of stationar- ity).– Solid Earth. doi: 10.5194/se-2020-83 BOROVIĆ, S. (2015): Integrirani hidrogeološko - hidrogeokemijski model Daruvarskog geotermalnog vodonosnika [Integrated hydrogeological- hydrogeochemical model of Daruvar geothermal aquifer – in Croatian].– Unpubl. PhD Thesis, Faculty of Mining, Geology and Petroleum Engi- neering, University of Zagreb, 228 p. BOROVIĆ, S. & MARKOVIĆ, I. (2015): Utilization and tourism valorisation of geothermal waters in Croatia.– Renewable and Sustainable Energy Reviews, 44, 52–63. doi: 10.1016/j.rser.2014.12.022 BOROVIĆ, S., POLA, M., BAČANI, A. & URUMOVIĆ, K. (2019): Constrain- ing the recharge area of a hydrothermal system in fractured carbonates by numerical modelling.– Geothermics, 82, 128–149. doi: 10.1016/j.geo- thermics.2019.05.017 CAINE, J. S., EVANS, J. P. & FORSTER, C. B. (1996): Fault zone architec- ture and permeability structure.– Geology, 24/11, 1025–1028. doi: 10.1130/0091-7613(1996)024<1025:FZAAPS>2.3.CO CARDOZO, N. & ALLMENDINGER, R. W. (2013): Spherical projections with OSXStereonet.– Computers and Geosciences, 51, 193–205. doi: 10.1016/j.cageo.2012.07.021. CARRERA, J. & MARTINEZ-LANDA, L. (2000): Mixed discrete-continu- um models: A summary of experiences in test interpretation and model prediction.– In: FAYBISHENKO, B., WITHERSPOON, P.A. & BEN- SON, S.M. (eds.): Dynamics of Fluid in Fractured Rock, 251–265. doi: 10.1029/GM122p0251 CSONTOS, L. & VÖRÖS, A. (2004): Mesozoic plate tectonic reconstruction of the Carpathian region.– Palaeogeography, Palaeoclimatology, Palaeo- ecology, 210/1, 1–56. doi: 10.1016/j.palaeo.2004.02.033 CECCATO, A., VIOLA, G., ANTONELLINI, M., TARTAGLIA, G. & RY- AN, E.J. (2021): Constraints upon fault zone properties by combined structural analysis of virtual outcrop models and discrete fracture net- work modelling.– Journal of Structural Geology, 152, 104444. doi: 10.1016/j.jsg.2021.104444 DELVAUX, D. & SPERNER, B. (2003): New aspects of tectonic stress inversion with reference to the TENSOR program.– Geological Society, London, Special Publications, 212/1, 75–100. doi: 10.1144/GSL. SP.2003.212.01.06 DERSHOWITZ, W.S. & HERDA, H.H. (1992): Interpretation of fracture spac- ing and intensity.– All Days, ARMA-92-0757. DOBLAS, M. (1998): Slickenside kinematic indicators.– Tectonophysics, 295/1–2, 187–197. doi: 10.1016/S0040-1951(98)00120-6 DE DREUZY, J.‐R., DAVY, P. & BOUR, O. (2002): Hydraulic properties of two‐dimensional random fracture networks following power law distri- butions of length and aperture.– Water Resources Research, 38/12, 1276. doi: 10.1029/2001WR001009 DE DREUZY, J., DAVY, P. & BOUR, O. (2001): Hydraulic properties of two‐ dimensional random fracture networks following a power law length distribution: 2. Permeability of networks based on lognormal distribu- tion of apertures.– Water Resources Research, 37/8, 2079–2095. doi: 10.1029/2001WR900010 G eologia C roatica 117Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing DUFFIELD, G.M. (2007): AQTESOLV for Windows Version 4.5 User's Guide.– HydroSOLVE, Inc., Reston, VA. FABBRI, P. (1997): Transmissivity in the Geothermal Euganean Basin: A Geostatistical Analysis.– Ground Water, 355, 881–887. doi: 10.1111/j.1745-6584.1997.tb00156.x FAULKNER, D.R., JACKSON, C.A.L., LUNN, R.J., SCHLISCHE, R.W., SHIPTON, Z.K., WIBBERLEY, C.A.J. & WITHJACK, M.O. (2010): A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones.– Journal of Structural Geology, 32/11, 1557–1575. doi: 10.1016/j.jsg.2010.06.009 FAYBISHENKO, B., WITHERSPOON, P.A. & BENSON, S.M. (2000): Dynamics of Fluids in Fractured Rock. Geophysical Monograph Series. Vol. 122. – American Geophysical Union, Washington, D. C. doi:10.1029/ GM122 FREEZE, R.A. & CHERRY, J.A. (1979): Groundwater.– Englewood Cliffs, NJ: Prentice-Hall, 604 p. FINSTER, M., CLARK, C., SCHROEDER, J. & MARTINO, L. (2015): Geo- thermal produced fluids: Characteristics, treatment technologies, and management options.– Renewable and Sustainable Energy Reviews, 50, 952–966. doi: 10.1016/j.rser.2015.05.059 GIUFFRIDA, A., LA BRUNA, V., CASTELLUCCIO, P., PANZA, E., RUS- TICHELLI, A., TONDI, E., GIORGIONI, M. & AGOSTA, F. (2019): Fracture simulation parameters of fractured reservoirs: Analogy with outcropping carbonates of the Inner Apulian Platform, southern Italy.– Journal of Structural Geology, 123, 18–41. doi: 10.1016/j.jsg.2019.02.007 GUIHÉNEUF, N., DAUSSE, A., DE DREUZY, J.R. & PARKER, B.L. (2021): Flow-bearing structures of fractured rocks: Insights from hydraulic property scalings revealed by a pumping test.– Journal of Hydrology, 598. doi: 10.1016/j.jhydrol.2020.125715 HANTUSH, M.S. (1961a): Aquifer Tests on Partially Penetrating Wells.– Journal of the Hydraulics Division, 87/5, 171–195. doi: 10.1061/ JYCEAJ.0000639 HANTUSH, M. S. (1961b): Drawdown Around a Partially Penetrating Well.– Journal of the Hydraulics Division, 87/4, 83–98. doi: 10.1061/JYCE- AJ.0000633 HODGETTS, D. (2013): Laser scanning and digital outcrop geology in the petroleum industry: A review.– Marine and Petroleum Geology, 46, 335–354. doi: 10.1016/j.marpetgeo.2013.02.014 HORVÁTH, F., MUSITZ, B., BALÁZS, A., VÉGH, A., UHRIN, A., NÁDOR, A., KOROKNAI, B., PAP, N., TÓTH, T. & WÓRUM, G. (2015): Evolu- tion of the Pannonian basin and its geothermal resources.– Geothermics, 53, 328–352. doi: 10.1016/j.geothermics.2014.07.009 HORVÁTH, F., BADA, G., SZAFIÁN, P., TARI, G., ÁDÁM, A. & CLOE- TINGH, S. (2006): Formation and deformation of the Pannonian Basin: Constraints from observational data.– Geological Society Memoir, 32, 191–206. doi: 10.1144/GSL.MEM.2006.032.01.11 HYMAN, J.D., KARRA, S., MAKEDONSKA, N., GABLE, C.W., PAIN- TER, S.L. & VISWANATHAN, H.S. (2015): dfnWorks: A discrete fra- cture network framework for modeling subsurface flow and transport.– Computers & Geosciences, 84, 10–19. doi: 10.1016/j.cageo.2015.08.001 JAMIČIĆ, D. (1995): The role of sinistral strike-slip faults in the formation of the structural fabric of the Slavonian Mts. (eastern Croatia).– Geologia Croatica, 48/2, 155–160. JAMIČIĆ, D., VRAGOVIĆ, M. & MATIČEC, D. (1989): Osnovna geološka karta SFRJ 1:100 000. Tumač za list Daruvar (Basic geological map of SFRJ 1:100 000. Explanatory notes for sheet Daruvar). – Geol. zavod, Zagreb, Sav. geol. zavod Beograd, Beograd, 55 p. JAROSINSKI, M., BEEKMAN, F., MATENCO, L. & CLOETINGH, S. (2011): Mechanics of basin inversion: Finite element modelling of the Pannonian Basin System.– Tectonophysics, 502/1–2, 121–145. doi: 10.1016/j.tecto.2009.09.015 NG, K. & SANTAMARINA, J. C. (2023): Mechanical and hydraulic proper- ties of carbonate rock: The critical role of porosity.– Journal of Rock Mechanics and Geotechnical Engineering, 15/4, 814–825. doi: 10.1016/j. jrmge.2022.07.017 KEEGAN-TRELOAR, R., IRVINE, D.J., SOLÓRZANO-RIVAS, S.C., WERNER, A.D., BANKS, E.W. & CURRELL, M.J. (2022): Fault-con- trolled springs: A review.– Earth-Science Reviews, 230, 104058. doi: 10.1016/j.earscirev.2022.104058 KORNEVA, I., CILONA, A., TONDI, E., AGOSTA, F. & GIORGIONI, M. (2015): Characterisation of the permeability anisotropy of Cretaceous platform carbonates by using 3D fracture modeling: the case study of Agri Valley fault zones (southern Italy).– Italian Journal of Geosciences, 134/3, 396–408. doi: 10.3301/IJG.2014.26 KOSOVIĆ, I., BRIŠKI, M., PAVIĆ, M., PADOVAN, B., PAVIČIĆ, I., MATOŠ, B., POLA, M. & BOROVIĆ, S. (2023): Reconstruction of Fault Archi- tecture in the Natural Thermal Spring Area of Daruvar Hydrothermal System Using Surface Geophysical Investigations (Croatia).– Sustaina- bility, 15/16, 12134. doi: 10.3390/su151612134 KOSOVIĆ, I., MATOŠ, B., PAVIČIĆ, I., POLA, M., MILEUSNIĆ, M., PA- VIĆ, M. & BOROVIĆ, S. (2024): Geological modeling of a tectonically controlled hydrothermal system in the southwestern part of the Pannonian basin (Croatia).– Frontiers in Earth Science. 12:1401935. doi: 10.3389/feart.2024.1401935 KRUSEMAN, G.P., DE RIDDER, N.A. & VERWEIJ, J.M. (2000): Analysis and evaluation of pumping test data.Vol. 47.– International institute for land reclamation and improvement, Wageningen, 378 p. LAI, J., SU, Y., XIAO, L., ZHAO, F., BAI, T., LI, Y., LI, H., HUANG, Y., WANG, G. & QIN, Z. (2024): Application of geophysical well logs in solving geologic issues: Past, present and future prospect.– Geoscience Frontiers, 15/3, 101779. doi: 10.1016/j.gsf.2024.101779 LEI, Q., LATHAM, J.-P. & TSANG, C.-F. (2017): The use of discrete fracture networks for modelling coupled geomechanical and hydrological behav- iour of fractured rocks.– Computers and Geotechnics, 85, 151–176. doi: 10.1016/j.compgeo.2016.12.024 MAMMOLITI, E., PEPI, A., FRONZI, D., MORELLI, S., VOLATILI, T., TAZIOLI, A. & FRANCIONI, M. (2023): 3D Discrete Fracture Network Modelling from UAV Imagery Coupled with Tracer Tests to Assess Frac- ture Conductivity in an Unstable Rock Slope: Implications for Rockfall Phenomena.– Remote Sensing, 15/5, 1222. doi: 10.3390/rs15051222 MALVIĆ, T. & CVETKOVIĆ, M. (2013): Lithostratigraphic units in the Dra- va Depression (Croatian and Hungarian parts) – a correlation.– Nafta, 64/1, 27–33. MARRETT, R. & ALLMENDINGER, R.W. (1990): Kinematic analysis of fault-slip data.– Journal of Structural Geology, 12/8, 973–986. doi: 10.1016/0191-8141(90)90093-E MARTINELLI, M., BISTACCHI, A., MITTEMPERGHER, S., BONNEAU, F., BALSAMO, F., CAUMON, G. & MEDA, M. (2020): Damage zone characterization combining scan-line and scan-area analysis on a km- scale Digital Outcrop Model: The Qala Fault (Gozo).– Journal of Struc- tural Geology, 140, 104144. doi: 10.1016/j.jsg.2020.104144 MAULDON, M. (1994): Intersection probabilities of impersistent joints.– In- ternational Journal of Rock Mechanics and Mining Sciences & Geome- chanics Abstracts, 31/2, 107–115. doi: 10.1016/0148-9062(94)92800-2 MEDICI, G., LING, F. & SHANG, J. (2023): Review of discrete fracture net- work characterization for geothermal energy extraction.– Frontiers in Earth Science, 11, 1–17. doi: 10.3389/feart.2023.1328397 MEDICI, G., SMERAGLIA, L., TORABI, A. & BOTTER, C. (2020): Review of Modeling Approaches to Groundwater Flow in Deformed Carbonate Aquifers.– Groundwater, 593, 334–351. doi: 10.1111/gwat.13069 MOENCH, A.F. (1984): Double-porosity models for a fissured groundwater reservoir with fracture skin.– Water Resources Research, 207, 831–846. doi: 10.1029/WR020i007p00831 ODA, M. (1985): Permeability tensor for discontinuous rock masses.– Géo- technique, 35/4, 483–495. doi: 10.1680/geot.1985.35.4.483 ORTNER, H., REITER, F. & ACS, P. (2002): Easy handling of tectonic data: the programs TectonicVB for Mac and TectonicsFP for WindowsTM.– Computers & Geosciences, 28/10, 1193–1200. doi: 10.1016/S0098- 3004(02)00038-9 PAMIĆ, J., RADONIĆ, G. & PAVIĆ, G. (2003): Geološki vodič kroz Park prirode Papuk (Geological guide through the Papuk nature park).– Park prirode Papuk, Požega, 66 p. G eo lo gi a C ro at ic a 118 Geologia Croatica 77/2 PAVIĆ, M., KOSOVIĆ, I., POLA, M., URUMOVIĆ, K., BRIŠKI, M. & BOROVIĆ, S. (2023): Multidisciplinary Research of Thermal Springs Area in Topusko (Croatia).– Sustainability, 15/6, 5498. doi: 10.3390/ su15065498 POURASKARPARAST, Z., AGHAEI, H., COLOMBERA, L., MASOERO, E. & GHAEDI, M. (2024): Fracture aperture: A review on fundamental concepts, estimation methods, applications, and research gaps.– Marine and Petroleum Geology, 164, 106818. doi: 10.1016/j.marpet- geo.2024.106818 PRELOGOVIĆ, E., SAFTIĆ, B., KUK, V., VELIĆ, J., DRAGAŠ, M. & LUČIĆ, D. (1998): Tectonic activity in the Croatian part of the Pannon- ian basin.– Tectonophysics, 297/1–4, 283–293. doi: 10.1016/S0040- 1951(98)00173-5 PRIEST, S.D. (1993): Discontinuity Analysis for Rock Engineering. Vol. 30. – Springer Netherlands, Dordrecht, 323–324 p. doi:10.1007/978-94-011- 1498-1 RENSHAW, C.E. (1995): On the relationship between mechanical and hydraulic apertures in rough‐walled fractures.– Journal of Geophysical Research: Solid Earth, 100/B12, 24629–24636. doi: 10.1029/95JB02159 ROMANO, V., BIGI, S., CARNEVALE, F., DE’HAVEN HYMAN, J., KAR- RA, S., VALOCCHI, A.J., TARTARELLO, M.C. & BATTAGLIA, M. (2020): Hydraulic characterization of a fault zone from fracture distribu- tion.– Journal of Structural Geology, 135, 104036. doi: 10.1016/j. jsg.2020.104036 RYBACH, L. & MONGILLO, M. (2006): Geothermal Sustainability-A Review with Identified Research Needs.– GRC Transactions, 30, 1083–1090. SAFTIĆ, B., VELIĆ, J., SZTANO, O., JUHASZ, G. & IVKOVIĆ, Ž. (2003): Tertiary Subsurface Facies, Source Rocks and Hydrocarbon Reservoirs in the SW Part of the Pannonian Basin (Northern Croatia and South- Western Hungary).– Geologia Croatica, 56/1, 101–122. doi: 10.4154/232 SCHMID, S.M., FÜGENSCHUH, B., KOUNOV, A., MAŢENCO, L., NIEV- ERGELT, P., OBERHÄNSLI, R., PLEUGER, J., SCHEFER, S., SCHUS- TER, R., TOMLJENOVIĆ, B., USTASZEWSKI, K. & VAN HINSBER- GEN, D.J.J. (2020): Tectonic units of the Alpine collision zone between Eastern Alps and western Turkey.– Gondwana Research, 78, 308–374. doi: 10.1016/j.gr.2019.07.005 SCHMID, S.M., BERNOULLI, D., FÜGENSCHUH, B., MATENCO, L., SCHEFER, S., SCHUSTER, R., TISCHLER, M. & USTASZEWSKI, K. (2008): The Alpine-Carpathian-Dinaridic orogenic system: correla- tion and evolution of tectonic units.– Swiss Journal of Geosciences, 101/1, 139–183. doi: 10.1007/s00015-008-1247-3 SIEGEL, L., GOLDSCHEIDER, N., PETITTA, M., XANKE, J., ANDREO, B., et al. (2023): Distribution, threats and protection of selected karst groundwater-dependent ecosystems in the Mediterranean region.– Hy- drogeology Journal, 31/8, 2231–2249. doi: 10.1007/s10040-023-02711-9 SINGHAL, B.B.S. & GUPTA, R.P. (2010): Applied Hydrogeology of Fractured Rocks. – Springer Netherlands, Dordrecht. doi: 10.1007/978- 90-481-8799-7 SMERAGLIA, L., MERCURI, M., TAVANI, S., PIGNALOSA, A., KETTER- MANN, M., BILLI, A. & CARMINATI, E. (2021): 3D Discrete Fracture Network (DFN) models of damage zone fluid corridors within a reser- voir-scale normal fault in carbonates: Multiscale approach using field data and UAV imagery.– Marine and Petroleum Geology, 126, 104902. doi: 10.1016/j.marpetgeo.2021.104902 SMITH, S.A.F., BISTACCHI, A., MITCHELL, T.M., MITTEMPERGHER, S. & DI TORO, G. (2013). The structure of an exhumed intraplate seis- mogenic fault in crystalline basement.– Tectonophysics, 599, 29–44. doi: 10.1016/j.tecto.2013.03.031 STEVANOVIĆ, Z. (2019): Karst waters in potable water supply: a global scale overview.– Environmental Earth Sciences, 78/23, 662. doi: 10.1007/ s12665-019-8670-9 STORTI, F., BISTACCHI, A., BORSANI, A., BALSAMO, F., FETTER, M., OGATA, K. (2022). Spatial and spacing distribution of joints at (over-) saturation in the turbidite sandstones of the Marnoso-Arenacea Fm. (Northern Apennines, Italy).– Journal of Structural Geology, 156, 104551. doi: 10.1016/j.jsg.2022.104551 SZANYI, J., RYBACH, L. & ABDULHAQ, H.A. (2023): Geothermal Energy and Its Potential for Critical Metal Extraction – A Review.– Energies, 16/20, 1–28. doi: 10.3390/en16207168 ŠIKIĆ, K. (1981): Facijesi mezozoika Papuckog gorja [Facies of the Mesozoic of Mount Papuk – in Croatian]. – Unpubl. PhD Thesis, Faculty of Sci- ence, University of Zagreb. TARI, G., DÖVÉNYI, P., DUNKL, I., HORVÁTH, F., LENKEY, L., STEFA- NESCU, M., SZAFIÁN, P. & TÓTH, T. (1999): Lithospheric structure of the Pannonian basin derived from seismic, gravity and geothermal data.– Geological Society, London, Special Publications, 156/1, 215–250. doi: 10.1144/GSL.SP.1999.156.01.12 TAVANI, S., STORTI, F., BAUSÀ, J. & MUÑOZ, J. A. (2012): Late thrusting extensional collapse at the mountain front of the northern Apennines (Italy).– Tectonics, 31/4. doi: 10.1029/2011TC003059 THEIS, C.V. (1935): The relation between the lowering of the Piezometric surface and the rate and duration of discharge of a well using ground‐ water storage.– Eos, Transactions American Geophysical Union, 16/2, 519–524. doi: 10.1029/TR016i002p00519 THIELE, S.T., GROSE, L., SAMSU, A., MICKLETHWAITE, S., VOLLG- GER, S.A. & CRUDEN, A.R. (2017): Rapid, semi-automatic fracture and contact mapping for point clouds, images and geophysical data.– So- lid Earth, 8/6, 1241–1253. doi: 10.5194/se-8-1241-2017 TOMLJENOVIĆ, B. & CSONTOS, L. (2001): Neogene-quaternary structures in the border zone between Alps, Dinarides and Pannonian Basin (Hr- vatsko zgorje and Karlovac basins, Croatia).– International Journal of Earth Sciences, 90/3, 560–578. doi: 10.1007/s005310000176 TURNER, F.J. (1953): Nature and dynamic interpretation of deformation la- mellae in calcite of three marbles.– American Journal of Science, 251/4, 276–298. doi: 10.2475/ajs.251.4.276 TZIAVOU, O., PYTHAROULI, S. & SOUTER, J. (2018): Unmanned Aerial Vehicle (UAV) based mapping in engineering geological surveys: Con- siderations for optimum results.– Eng. Geol., 232, 12–21. doi: 10.1016/j. enggeo.2017.11.004 UNITED NATIONS (2022): The United Nations World Water Development Report 2022: Groundwater: Making the invisible visible. UNESCO, Pa- ris, 225 p. URUMOVIĆ, K., TERZIĆ, J., KOPIĆ, J. & KOSOVIĆ, I. (2023): Identifica- tion of Aquifer and Pumped Well Parameters Using the Data Hidden in Non-Linear Losses.– Sustainability (Switzerland), 15/14, 11170. doi: 10.3390/su151411170 USTASZEWSKI, K., HERAK, M., TOMLJENOVIĆ, B., HERAK, D. & MA- TEJ, S. (2014): Neotectonics of the Dinarides-Pannonian Basin transiti- on and possible earthquake sources in the Banja Luka epicentral area.– Journal of Geodynamics, 82, 52–68. doi: 10.1016/j.jog.2014.04.006 USTASZEWSKI, K., SCHMID, S. M., FÜGENSCHUH, B., TISCHLER, M., KISSLING, E. & SPAKMAN, W. (2008): A map-view restoration of the Alpine-Carpathian-Dinaridic system for the Early Miocene.– Swiss Jo- urnal of Geosciences, 101/SUPPL. 1. doi: 10.1007/s00015-008-1288-7 USTASZEWSKI, K., KOUNOV, A., SCHMID, S. M., SCHALTEGGER, U., KRENN, E., FRANK, W. & FÜGENSCHUH, B. (2010): Evolution of the Adria-Europe plate boundary in the northern Dinarides: From con- tinent-continent collision to back-arc extension.– Tectonics, 29/6, TC6017. doi: 10.1029/2010TC002668 VERBOVŠEK, T. (2008): Estimation of Transmissivity and Hydraulic Con- ductivity from Specific Capacity and Specific Capacity Index in Dolo- mite Aquifers.– Journal of Hydrologic Engineering, 13/9, 817–823. doi: 10.1061/(ASCE)1084-0699(2008)13:9(817) VOECKLER, H. & ALLEN, D. M. (2012): Estimating regional-scale fractu- red bedrock hydraulic conductivity using discrete fracture network (DFN) modeling.– Hydrogeology Journal, 20/6, 1081–1100. doi: 10.1007/ s10040-012-0858-y WANG, X. (2005): Stereological interpretation of rock fracture traces on bo- rehole walls and other cylindrical surfaces.– Virginia Techp. WORTHINGTON, S.R.H., FOLEY, A.E. & SOLEY, R.W.N. (2019): Transient characteristics of effective porosity and specific yield in bedrock https://doi.org/10.1016/j.enggeo.2017.11.004 https://doi.org/10.1016/j.enggeo.2017.11.004 G eologia C roatica 119Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing aquifers.– Journal of Hydrology, 578, 124129. doi: 10.1016/j.jhy- drol.2019.124129 ZAMBRANO, M., TONDI, E., KORNEVA, I., PANZA, E., AGOSTA, F., JANISECK, J.M. & GIORGIONI, M. (2016): Fracture properties analysis and discrete fracture network modelling of faulted tight limestones, Murge Plateau, Italy.– Italian Journal of Geosciences, 135/1, 55–67. doi: 10.3301/IJG.2014.42 ZEEB, C., GOMEZ-RIVAS, E., BONS, P. D. & BLUM, P. (2013): Evaluation of sampling methods for fracture network characterization using out- crops.– AAPG Bulletin, 97/9, 1545–1566. doi: 10.1306/02131312042 ZHANG, L. & EINSTEIN, H. H. (1998): Estimating the Mean Trace Length of Rock Discontinuities.– Rock Mechanics and Rock Engineering, 31/4, 217–235. doi: 10.1007/s006030050022 URL 1: https://earth.google.com/web/@45.53570123,17.29653113,232.316029 05a,6291.71769557d,35y,-12.560923h,2.44982084t,0r/data=OgMKATA (visited April 2024) URL 2: https://www.pix4d.com/product/pix4dmapper-photogrammetry-sof- tware/ (visited April 2024) URL 3: https://www.cloudcompare.org/ (visited April 2024) URL 4: https://doi.org/10.5281/zenodo.7890077 (visited April 2024) URL 5: https://www.mathworks.com/ (visited April 2024) URL 6: http://qgis.org (visited April 2024) URL 7: https://doi.org/10.5281/zenodo.11032168 (visited April 2024) URL 8: https://www.petex.com/products/move-suite/move/ (visited April 2024) https://earth.google.com/web/@45.53570123,17.29653113,232.31602905a,6291.71769557d,35y,-12.560923h,2.44982084t,0r/data=OgMKATA https://earth.google.com/web/@45.53570123,17.29653113,232.31602905a,6291.71769557d,35y,-12.560923h,2.44982084t,0r/data=OgMKATA https://www.pix4d.com/product/pix4dmapper-photogrammetry-software/ https://www.pix4d.com/product/pix4dmapper-photogrammetry-software/ https://www.cloudcompare.org/ https://doi.org/10.5281/zenodo.7890077 https://www.mathworks.com/ http://qgis.org https://doi.org/ https://zenodo.org/doi/10.5281/zenodo.11032168 https://www.petex.com/products/move-suite/move/ G eo lo gi a C ro at ic a 120 Geologia Croatica 77/2 Supplementary Table 1. Average orientations of the measured carbonate beds with computed fold axes. Strata orientation Data Dip direction (°) Dip angle (°) Strike Fold axis orientation Dip direction (°) Dip angle (°) So1 29 57 80 147-327 138 42 So2 226 88 136-316 So3 39 310 66 40-220 228 16 So4 138 88 48-228 Supplementary Table 2. Average orientation of the three fracture sets in BRQ. Fracture system Data Dip direction (°) Dip angle (°) Strike Fs1a 31 113 81 23-203 Fs1b 28 288 82 18-198 Fs2 42 100 45 10-190 Fs3a 32 22 85 112-292 Fs3b 22 338 73 68-248 Supplementary Table 3. Mean geometrical properties of the observed fault planes with calculated kinematic indicators and parameters. Fault planes (Fault types: R – reverse; N – normal; SSF – strike-slip were grouped following their geometrical properties. Orientations of the P- and T-axes are based on constructed synthetic structural beach-ball diagrams (Fig. 6). Fault group Subset Data Dip direction (°) Dip angle (°) Pitch (°) Striation P-axis T-axis Trend (°) Plunge (°) Trend (°) Plunge (°) Trend (°) Plunge (°) RF RFa 20 128 56 72 115 43 128 12 336 76 RFb 289 36 76 261 66 NF NFa 11 153 41 81 126 65 14 87 325 11 NFb 319 50 79 211 68 SSF-1 SSF-1a 10 150 52 24 78 20 13 17 275 23 SSF-1b 45 73 17 308 16 SSF-2 SSF-2a 9 113 66 18 67 11 252 1 345 30 SSF-2b 211 71 19 124 17 Supplementary Table 4. Average orientations of the main discontinuity sets in the quarry sections (QS2-QS6) and amount of data constituting the sets. Quarry section Set Data Dip direction (°) Dip angle (°) Strike QS2 D1a 20 315 76 45-225 D1b 17 116 83 26-206 D2a 23 248 53 158-338 D2b 15 277 57 137-317 D3 21 272 87 2-182 QS 3A D1a 132 304 68 34-214 D1b 43 285 83 15-195 D2a 140 241 75 151-331 D2b 74 240 41 150-330 QS 3B D1a 124 329 70 59-239 D1b 73 303 77 33-213 D2 115 258 72 168-348 QS 4 D1a 103 328 70 58-238 D1b 86 304 69 34-214 D2 89 254 72 164-344 QS 5 D1 182 308 71 38-218 D2 115 267 72 177-357 QS 6 D1a 16 331 70 61-241 D1b 9 313 72 43-223 D2a 5 270 60 0-180 D2b 8 265 81 175-355 G eologia C roatica 121Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Supplementary Table 5. Results of statistical analysis for the horizontal (H) and vertical (V) discontinuity network. Discontinuity set - H Fracture network - V Total number of fractures 352 260 % censored 11.08 20 Data Fit Data Fit Mean 3.5306 3.8836 3.6673 5.0660 Median 2.0381 2.3812 1.8787 2.4461 Mode 0.2029 0.8952 0.3225 0.5703 B5 0.5387 0.4680 0.4375 0.3361 B95 12.4790 12.1160 13.1371 17.8025 Std 4.5250 5.0036 4.5258 9.1879 Var 20.4755 25.0358 20.4825 84.4173 Supplementary Figure 1. Statistical distribution of the heights of the horizontal (H) and vertical (V) fracture network (Quarry Section - 3A). Statistical analysis of the main geometrical parameter (height) was performed by using FracAbility. Supplementary Table 6. Input parameters used for the construction of the DFN model of the QS3A section. The measurement units of the param- eters are reported in brackets. The obtained model is shown in Fig. 10. Parameter Set - H Set - V Dip direction 240.86 295.86 Dip 64.92 74.79 Fisher K 13.53 18.83 Height – mean (m) 3.88 5.06 Height – STD (m) 5.00 9.18 Length/Height ratio 0.25 0.25 Sum of fractures heights (m) 1244.96 961.19 P21 (m/m2) 0.47 0.37 P32 (m2/m3) 1.98 1.18 Aperture (mm) 4.15 4.15 G eo lo gi a C ro at ic a 122 Geologia Croatica 77/2 Supplementary Figure 2. Distributions of porosity (Φ) obtained for the tested aperture values and employing the different methodological approaches. G eologia C roatica 123Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Supplementary Figure 3. Distributions of permeability (k) obtained for the tested aperture values and employing the different methodological approaches. G eo lo gi a C ro at ic a 124 Geologia Croatica 77/2 Supplementary Table 7. Statistical distributions of the porosity (Φ) in the DFN models conducted considering different methodological approaches for assigning the fracture aperture and using different aperture values. Aperture (mm) Porosity – Φ (%) Mean Median Standard deviation Min Max Skewness Kurtosis 10th perc. 25th perc. 50th perc. 75th perc. 90th perc. Co ns ta nt 0.1 0.03 0.03 0.02 0.00 0.24 0.91 3.85 0.00 0.01 0.03 0.05 0.06 0.2 0.06 0.05 0.05 0.00 0.47 0.91 3.85 0.00 0.02 0.05 0.09 0.13 0.5 0.15 0.13 0.12 0.00 1.18 0.91 3.85 0.00 0.05 0.13 0.23 0.32 1 0.30 0.26 0.25 0.00 2.37 0.91 3.85 0.00 0.11 0.26 0.45 0.64 2 0.61 0.51 0.49 0.00 4.73 0.91 3.85 0.00 0.22 0.51 0.90 1.29 5 1.52 1.28 1.23 0.00 11.83 0.91 3.85 0.00 0.55 1.28 2.25 3.21 10 3.04 2.56 2.46 0.00 23.65 0.91 3.85 0.00 1.10 2.56 4.50 6.43 Ro ot s qu ar e le ng th 0.1 0.05 0.05 0.04 0.00 0.43 0.97 4.00 0.00 0.02 0.05 0.08 0.11 0.2 0.11 0.09 0.09 0.00 0.86 0.97 4.00 0.00 0.04 0.09 0.16 0.23 0.5 0.27 0.23 0.22 0.00 2.14 0.97 4.00 0.00 0.09 0.23 0.40 0.57 1 0.53 0.46 0.44 0.00 4.28 0.97 4.00 0.00 0.18 0.46 0.79 1.14 2 1.07 0.91 0.88 0.00 8.57 0.97 4.00 0.00 0.36 0.91 1.59 2.28 5 2.67 2.29 2.21 0.00 21.42 0.97 4.00 0.00 0.91 2.29 3.97 5.70 10 5.34 4.57 4.42 0.00 42.84 0.97 4.00 0.00 1.81 4.57 7.94 11.41 Le ng th 0.1 0.07 0.05 0.06 0.00 0.53 1.09 4.34 0.00 0.02 0.05 0.10 0.14 0.2 0.13 0.11 0.11 0.00 1.07 1.09 4.34 0.00 0.04 0.11 0.20 0.29 0.5 0.33 0.27 0.28 0.00 2.67 1.09 4.34 0.00 0.10 0.27 0.49 0.72 1 0.65 0.55 0.57 0.00 5.35 1.09 4.34 0.00 0.19 0.55 0.98 1.44 2 1.31 1.09 1.14 0.00 10.69 1.09 4.34 0.00 0.39 1.09 1.95 2.88 5 3.27 2.73 2.85 0.00 26.73 1.09 4.34 0.00 0.97 2.73 4.88 7.19 10 6.54 5.47 5.70 0.00 53.46 1.09 4.34 0.00 1.93 5.47 9.75 14.38 Supplementary Table 8. Statistical distribution of the permeability (k) in the DFN models conducted considering different methodological approaches for assigning the fracture aperture and using different aperture values. Aperture (mm) Permeability – k (D) Mean Median Standard deviation Min Max Skewness Kurtosis 10th perc. 25th perc. 50th perc. 75th perc. 90th perc. Co ns ta nt 0.1 2.5e-1 2.2e-1 2.e-1 0.00 1.9 0.85 3.69 0.00 9.3e-2 2.2e-1 3.8e-1 5.3e-1 0.2 2 1.7 1.6 0.00 1.5e1 0.85 3.69 0.00 7.4e-1 1.7 3 4.3 0.5 3.2e1 2.7e1 2.5e1 0.00 2.3e2 0.85 3.69 0.00 1.2e1 2.7e1 4.7e1 6.6e1 1 2.5e2 2.2e2 2.e2 0.00 1.9e3 0.85 3.69 0.00 9.3e1 2.2e2 3.8e2 5.3e2 2 2.e3 1.7e3 1.6e3 0.00 1.5e4 0.85 3.69 0.00 7.4e2 1.7e3 3.e3 4.3e3 5 3.2e4 2.7e4 2.5e4 0.00 2.3e5 0.85 3.69 0.00 1.2e4 2.7e4 4.7e4 6.6e4 10 2.5e5 2.2e5 2.e5 0.00 1.9e6 0.85 3.69 0.00 9.3e4 2.2e5 3.8e5 5.3e5 Ro ot s qu ar e le ng th 0.1 1.6 1.2 1.5 0.00 1.4e1 1.28 4.99 0.00 3.8e-1 1.2 2.3 3.6 0.2 1.2e1 9.8e 1.2e1 0.00 1.1e2 1.28 4.99 0.00 3.1 9.8 1.9e1 2.9e1 0.5 2.e2 1.5e2 1.8e2 0.00 1.7e3 1.28 4.99 0.00 4.8e1 1.5e2 2.9e2 4.5e2 1 1.6e3 1.2e3 1.5e3 0.00 1.4e4 1.28 4.99 0.00 3.8e2 1.2e3 2.3e3 3.6e3 2 1.2e4 9.8e3 1.2e4 0.00 1.1e5 1.28 4.99 0.00 3.1e3 9.8e3 1.9e4 2.9e4 5 2.e5 1.5e5 1.8e5 0.00 1.7e6 1.28 4.99 0.00 4.8e4 1.5e5 2.9e5 4.5e5 10 1.6e6 1.2e6 1.5e6 0.00 1.4e7 1.28 4.99 0.00 3.8e5 1.2e6 2.3e6 3.6e6 Le ng th 0.1 3.8 2.2 4.8 0.00 6.1e1 2.40 11.22 0.00 4.7e-1 2.2 5.4 9.7 0.2 3.1e1 1.7e1 3.8e1 0.00 4.9e2 2.40 11.22 0.00 3.8 1.7e1 4.3e1 7.8e1 0.5 4.8e2 2.7e2 6.e2 0.00 7.7e3 2.40 11.22 0.00 5.9e1 2.7e2 6.7e2 1.2e3 1 3.8e3 2.2e3 4.8e3 0.00 6.1e4 2.40 11.22 0.00 4.7e2 2.2e3 5.4e3 9.7e3 2 3.1e4 1.7e4 3.8e4 0.00 4.9e5 2.40 11.22 0.00 3.8e3 1.7e4 4.3e4 7.8e4 5 4.8e5 2.7e5 6.e5 0.00 7.7e6 2.40 11.22 0.00 5.9e4 2.7e5 6.7e5 1.2e6 10 3.8e6 2.2e6 4.8e6 0.00 6.1e7 2.40 11.22 0.00 4.7e5 2.2e6 5.4e6 9.7e6 G eologia C roatica 125Kosović et al.: Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing Supplementary Figure 4. Interpretation of the step-drawdown test results conducted in the Dar 1 well. Supplementary Table 9. Hydrogeological parameters obtained from the interpretation of the pumping test using Dar 1 as exploitation well and D 1 and DS 1 as observation wells. The estimated parameters are: transmis- sivity (T), storativity (S), hydraulic conductivity anisotropy ratio (Kz/Kr), and saturated thickness (b). Parameter Whole test Constant test T (m2/s) 3.82e-2 3.46e-2 S (-) 2.05e-3 1.14e-3 Kz/Kr (-) 6.81e-3 1.98e-2 b (m) 61.53 68.42