Hydrogeological parameterisation of the Daruvar thermal aquifer: integration of fracture network analysis and well testing

Highly fractured Mesozoic carbonate rocks are the main reservoir of many geothermal resources in northern Croatia, being of environmental, cultural, and economic value for the local and regional communities. The Daruvar thermal springs (temperatures < 50°C) represent 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 quantifying 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.


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 groundwater 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 groundwaterdependent ecosystems are extensively used, facing environmen tal threats (SIEGEL et al., 2023).Site-specific management plans considering both the principal flow paths and the hydraulic 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, 2010).For fractures here we consider every kind of dis con tinuity 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 wellestablished 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 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.

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

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(after SCHMID et al., , 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).
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, welllayered 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 shallowwater 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 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žanovac 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 preexisting 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 preexisting 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.

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 m 2 /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 SiO 2 geothermometers, and v) a mean residence time between 11 and 15 ka assessed through 14 C 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.

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

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 m 2 .
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  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): 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 conducted employing the software Petroleum Experts Move v. 2019.1 (URL 8).This software has been profitably used to calculate 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).

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 stepdrawdown 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-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: T (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).

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 tectonized 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

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, representa 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.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 cumulative fracture lengths were used to calculate the areal fracture intensity P21 being 0.47 and 0.37 for set H and V, respectively.

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 m 2 /m 3 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 10 th and 90 th 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 10 th and 90 th 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 10 th 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 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 10 th 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 75 th percentile.The resulting k value ranged from 0 to 4.4×10 6 D and an average value of 2.7×10 5 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.

Interpretation of pumping and well tests
The pumping test in the Daruvar thermal field lasted approximately 12 h.It was conducted performing 3 steps (i.e., stepdrawdown 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 m 2 /s and 3.46×10 -2 m 2 /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: 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 m 2 /s and 4.8×10 -2 m 2 /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.

Integrated interpretation of field and virtual structural data
Fractures are the main flow 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 model of the groundwater resource (SINGHAL & GUPTA, 2010).Integrated structural mapping based on field measurements 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 (Supplementary 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): discontinuities 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.

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 m 2 /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.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×10 5 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×10 5 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 relationship 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×10 5 D, which is 3-4 orders of magnitude 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 hydraulic aperture of the fractures, which is smaller than the mechanical 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 10 th 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 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).

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 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 continuities.Remote sensing photogrammetry allows the construction 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.

Figure 2 .
Figure 2. Simplified geological map of the Daruvar area and lithostratigraphic succession (modified from KOSOVIĆ et al., 2024).The black polygon represents 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.

Figure 3 .
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 .
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.

Figure 6 .
Figure 6.(A) Stereonets of the bedding orientations (small black circles) in the BRQ.Their interpretation highlights the occurrence of generally two generations 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 (F S ) are reported in Supplementary 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.
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 strikeslip 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.

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

Figure 8 .
Figure 8. Structural diagrams of the main discontinuity systems in the different sections of the quarry (N: number of data).The orientations of the representative 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).

Figure 9 .
Figure 9. Orthophotograph of section Q3A digitised in QGIS to vectorise the traces of the discontinuity systems.Traces are divided into two discontinuity systems: subvertical (red; set V in Supplementary Tab. 5) and subhorizontal (green; set H in Supplementary Tab. 5).

Figure 10 .
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 assigning 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.

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

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

Table 1 .
(BOROVIĆ et al., 2019;URUMOVIĆ et al., 2023)var 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).

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.th perc.25 th perc.50 th perc.75 th perc.90 th perc.

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.th perc.25 th perc.50 th perc.75 th perc.90 th perc.