1. INTRODUCTION Geothermal systems, which represent a renewable resource for energy and raw material production, vary considerably world- wide due to the different mechanisms governing their forma- tion. They are classified according to their principal properties and characteristics, including their geological, hydrogeologi- cal, geochemical, and thermal aspects (MOECK, 2014). A sub- set of geothermal systems is referred to as hydrothermal when the heat transfer mechanism involves circulating water, whether as liquid or vapour (OJHA et al., 2021; KHODAYAR & BJÖRNSSON, 2024). In the study of hydrothermal systems, it is necessary to determine the origin of the fluid and the area of recharge, the heat transfer mechanism, the direction of fluid flow and the depth to which it descends, the geometry of the aquifer and its hydrogeological and thermal properties, and the conditions favouring the outflow of the thermal water. Research of hydrothermal systems usually includes the application of an integrated multidisciplinary approach. The geological framework and tectonic evolution of the area influ- enced by the circulation of thermal fluids are typically recon- A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system Mirja Pavić1, Marco Pola1, Bojan Matoš2,*, Katarina Mišić2, Ivan Kosović1, Ivica Pavičić2 and Staša Borović1 1 Croatian Geological Survey, Sachsova 2, 10000 Zagreb, Croatia 2 University of Zagreb, Faculty of Mining, Geology and Petroleum Engineering, Pierottijeva ulica 6, 10000 Zagreb, Croatia; (*corresponding author: bojan.matos@rgn.unizg.hr) doi: 10.4154/gc.2024.14 Abstract A comprehensive understanding of hydrothermal systems is often obtained through the in- tegration of conceptual and numerical modelling. This integrated approach provides a struc- tured framework for the reconstruction and quantification of fluid dynamics in the reservoir, thereby facilitating informed decision-making for sustainable utilisation and environmental protection of the hydrothermal system. In this study, an updated conceptual model of the Topusko hydrothermal system (THS), central Croatia, is proposed based on structural, geo- chemical, and hydrogeological analyses. The stratigraphic sequence and the structural framework of the THS were defined based on geological maps and field investigations. As depicted by hydrochemical and isotope analyses, the thermal waters in the Topusko system (temperatures < 65 °C) are of meteoric origin and circulate in a carbonate aquifer. The THS receives diffuse recharge approximately 13 km S of Topusko, where Triassic carbonates crop out. Gravity-driven regional groundwater circulation is favoured by regional thrusts that tectonically uplifted Palaeozoic rocks of low permeability. These structures confine the flu- id flow in the permeable, fractured and karstified Triassic carbonates, favouring the north- ward circulation of the water. A regional anticline lifts the aquifer closer to the surface in Topusko. Open fractures in the anticline hinge zone increase the fracturing and permeabil- ity field of the aquifer, promoting the rapid upwelling of thermal water resulting in the Topus- ko thermal springs. Numerical simulations of fluid flow and heat transport corroborate the proposed conceptual model. In particular, a thermal anomaly was modelled in the Topusko subsurface with temperature values of 31.3 °C and 59.5 °C at the surface and at the base of the thermal aquifer, respectively, approaching the field observations. These findings show that the circulation of Topusko thermal water is influenced by regional and local geological structures suggesting that the enhanced permeability field in the discharge area enables the formation of the natural thermal springs. structed by combining regional and local field investigations and geophysical data (e.g., MUFFLER & CATALDI, 1978; FLÓVENZ et al., 2012; KOSOVIĆ et al., 2023, 2024). They provide insights into the surface geometry of geological for- mations and fracture networks and the kinematics of fault sys- tems that are consequently used for subsurface geological re- constructions, supported by 2D or 3D geological modelling. Hydrogeochemical research is the requisite for understanding and managing geothermal aquifers. It involves continuous monitoring of the thermal water to evaluate baseline levels, track their changes, and assess the impact of water abstraction, which is crucial for resource protection and legislative compli- ance (HOUNSLOW, 1995; MARINI, 2000; MAZOR, 2004; HEASLER, 2009). Analysis of groundwater chemistry and isotopic content aids in characterising hydrothermal systems, identifying water sources, estimating reservoir temperatures, and assessing potential mixing. Hydrogeological research helps in understanding the overall groundwater regime and distribution by providing data (i.e., hydraulic parameters of the aquifer and surrounding rocks, water flow velocities) on sub- 2024 | 77/3 | 291–310 | 9 Figs. | 1 Tab. | www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society Article history: Manuscript received: May 16, 2024 Revised manuscript accepted: June 11, 2024 Available online: August 28, 2024 Keywords: thermal spring, conceptual modelling, numerical modelling, recharge area, fractured carbonates, SW Pannonian basin, central Croatia mailto:bojan.matos@rgn.unizg.hr G eo lo gi a C ro at ic a 292 Geologia Croatica 77/3 surface conditions influencing the water circulation (FETTER 2001; GOLDSCHEIDER et al.2010; SZANYI & KOVÁCS, 2010; LEI & ZHU, 2013; RMAN, 2014; FABBRI et al., 2017). Quantifying and monitoring the hydrogeological parameters of the thermal aquifer and water is necessary for predicting the exploitable water volumes with an acceptable drawdown and identifying detrimental effects on the system. Addition- ally, thermal parametrisation of the geological units involved in the thermal fluid flow helps define the changes in the tem- perature field and fluid distribution across the system (FUCHS & BALLING, 2016; XIONG et al., 2020). The investigations mentioned above aid in the construc- tion of a conceptual model of the studied hydrothermal system. Developing the conceptual model of a hydrothermal system consolidates the existing understanding by integrating multi- disciplinary and multiscale datasets. Conceptual models de- scribe the main processes governing both fluid flow and heat transport, which influence the volume of the hydrothermal re- source and its geochemical and thermal characteristics. The aforementioned geological reconstructions serve as the foun- dation for a hydrogeological conceptual model of the hydro- thermal system, elucidating the mechanisms governing hydro- thermal resource formation (MOECK et al., 2014; CALCAGNO et al., 2014, MROCZEK et al., 2016). The physical reliability of the conceptual model can be constrained by developing variable-density fluid flow and heat transport numerical simu- lations of the system. Numerical models can be used for test- ing and quantifying the importance of different processes in the development of the geothermal resource and its physico- chemical characteristics (e.g., MÁDL-SZŐNYI & TÓTH, 2015; HAVRIL et al., 2016; MONTANARI et al., 2017 BO RO- VIĆ et al., 2019; POLA et al., 2020; TORRESAN et al., 2022). Furthermore, they can be used to reconstruct the historical and current state of the system and to forecast future impacts (AN- DERSON et al., 2015). Thermal springs in Croatia are generally part of interme- diate-scale hydrothermal systems, including recharge areas in the nearby mountainous hinterlands and geothermal aquifers mainly hosted in Mesozoic carbonate rocks (GOLDSCHEI- DER et al., 2010; BOROVIĆ et al., 2016). Their occurrence is favoured by the regional thermal characteristics in central and northern Croatia that are part of the Pannonian Basin System (PBS). The tectonic setting of the PBS is characterised by the thinned lithosphere, which enables an above-average heat flow from the asthenosphere (HORVÁTH et al., 2015). In the PBS, three levels of the regional flow of thermal water have been identified: i) gravity flows in the Neogene-Quaternary clastic rocks and sediments of the basin fill (the shallowest), ii) grav- ity flows in pre-Neogene confined carbonate aquifers below them, and iii) flow caused by overpressure in the deepest Mes- ozoic aquifers (HORVÁTH et al., 2015; VASS et al., 2018). The Mesozoic carbonate rocks, representing the deepest geother- mal aquifers, usually crop out either as inselbergs or along the margins of the basin. This could imply the occurrence of greater recharge in the marginal parts of the PBS, where the aquifer is shallower and covered by thinner Neogene deposits, favouring the development of local to intermediate scale hy- drothermal systems (STEVANOVIĆ, 2015; HAVRIL et al., 2016). The artesian thermal springs of Topusko have been re- nowned since Roman times, ranking as the second warmest in Croatia (BOROVIĆ et al., 2016; ŠIMUNIĆ, 2008). Thermal water with temperatures of up to 65°C has been used since the 1980s for health and recreational purposes and district heat- ing. Despite this fact, the geological features driving the de- velopment of the Topusko hydrothermal system (THS) and regulating the regional groundwater flow direction remained uncertain. In the initial stages of THS research, the potential recharge area was determined by defining outcrops of perme- able rocks in topographically prominent areas. The lack of previous systematic and detailed structural-geological and hy- drogeological investigations hindered the reconstruction of the regional geological evolution and understanding of how the THS functions. Though available publications and unpub- lished reports suggest contradictory hypotheses, the most commonly used conceptual model is from the early 2000s (ŠIMUNIĆ, 2008). In this study, we propose a novel concep- tual model of the THS, detailing the recharge area and main circulation paths using a collection of structural, geochemical, and hydrogeological field data. The second objective of this research involves conducting 2D numerical modelling to ana- lyse fluid flow and heat transport within the THS. Here, the 2D numerical modelling served as a physical validation for the proposed conceptual model, supporting it with the quantifica- tion of the main processes governing the development of the Topusko geothermal resource. 2. MATERIALS AND METHODS 2.1. Tectonic setting The THS formed in the pre-Neogene basement units of the Internal Dinarides (Fig. 1). The Internal Dinarides, as an integral part of the Adria Microplate, convey Adria’s eastern passive margin that was involved in a complex tectonic collision between the Adria Microplate and European foreland during the Cretaceous-Paleogene period (SCHMID et al., 2020 with references). This tectonic contraction (the recent convergence rate between the Adria indenter and Europe is ≤ 4.17 mm/yr according to D’AGOSTINO et al. 2008) resulted in the formation of a Dinaridic orogen system, tectonic suture zone (i.e., Sava Suture Zone), and 400 km eastward extrusion of the ALCAPA block (i.e., Eastern Alps, West Carpathians and Transdanubian ranges; TARI et al., 1999; CSONTOS & VÖRÖS, 2004). Besides the formation of an orogen-parallel thrust fault system, tectonic contraction accommodated the formation of regional dextral/sinistral faults (e.g., the Split- Karlovac fault, Periadriatic fault that extends into the Mid- Hungarian fault zone further to the E; Fig. 1), which enabled CCW/CW rotation and partial tectonic exhumation of the nearby tectonic blocks of the Adria Microplate and the Tisza- Dacia Mega-Unit (e.g. TOMLJENOVIĆ, 2002; TOMLJE NO- VIĆ et al., 2008; USTASZEWSKI et al., 2010; SCHMID et al., 2020). At the same time, as the study area is positioned in the immediate vicinity of the transient zone between the Adria G eologia C roatica 293Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system Microplate (W) and the Tisza Mega-Unit (E), the inherited fault system is characterised by polyphase tectonic evolution and persistent structural reactivation of the faults (Fig. 2; SCHMID et al., 2008). As a result, the THS bedrock units (Fig. 3) resemble a complex lithostratigraphic mosaic of Palaeozoic- Triassic clastic and carbonate units, that are often seen in tectonic contact with younger Jurassic-Cretaceous ophiolitic mélange units (SCHMID et al., 2008), and its Paleogene- Neogene cover. The Neogene-Quaternary tectonic evolution of the study area, on the other hand, was further affected by back-arc type formation of the PBS (ROYDEN & HORVÁTH, 1988; HORVÁTH et al., 2006; CLOETINGH at al., 2006). Formation of the PBS was characterised by repeated Early-Middle Miocene E-W oriented lithospheric extension (c. 26-11.5 Ma), along the NNW-striking normal listric faults (i.e., the Sava fault; Fig. 2), which was followed by its Late Miocene- -Pliocene-Quaternary tectonic inversion due to N – S com- pression (e.g., PRELOGOVIĆ et al., 1998; TARI et al., 1999; TOMLJENOVIĆ & CSONTOS, 2001; CLOETINGH et al. 2006; SCHMID et al., 2008; BRÜCKL et al., 2010). In the Croatian part of the PBS, a Neogene-Quaternary sediment succession (Fig. 3) is associated with the dominant NNW-striking Sava, Karlovac and Glina basins and subbasins, which were tectonically inverted and highly deformed during the Late Miocene-Pliocene-Quaternary periods (PAVELIĆ, 2001; PAVELIĆ et al., 2003; TOMLJENOVIĆ & CSONTOS, 2001). Tectonic deformation of the Neogene-Quaternary Figure 1. Regional tectonostratigraphic units that surround the study area of the THS (red dashed polygon corresponds to the extent of Fig. 3). Map shows the main regional fault systems that accommodated the tectonic collision of the Adria Microplate and European Foreland during Cretaceous- Paleogene time. The study area is located at the eastern margin of the Adria Microplate, within the Internal Dinarides, close to the Sava Suture Zone (modified after SCHMID et al., 2008; 2020). G eo lo gi a C ro at ic a 294 Geologia Croatica 77/3 succession in the study area is especially pronounced along the contact with the Palaeozoic-Triassic anticlinal core of the Petrova gora Mt. (HORVÁTH & TARI 1999; TOMLJENOVIĆ & CSONTOS, 2001). Though tectonic uplift of the Petrova gora Mt. had probably already started during the Cretaceous-Paleogene contraction (similar to the other PBS pre-Neogene basement highs e.g., Trgovska gora, Slavonian Mts.), the final uplift commenced during the Late Miocene-Pliocene-Quaternary compression/transpression phase that resulted in tectonic exhumation, tectonic overprint of its Palaeozoic-Mesozoic structures, and formation of kilometre-scale folds along the reactivated and newly formed faults in the area (Fig. 2; PRELOGOVIĆ et al., 1998; TOMLJENOVIĆ & CSONTOS, 2001). Ongoing, local NNE-SSW regional compression in the study area is driven by residual Adria indentation shortening locally at the scale of 1-2 mm/yr and is accommodated along the inherited faults with slip rates below 0.1 mm/yr (GRENERCZY et al., 2005; KASTELIC & CARAFA, 2012; USTASZEWSKI et al., 2014). 2.2. Geological setting A composite geological map covering the study area (Fig. 3) was constructed using basic geological maps of the former Yugoslavia at a scale of 1:100.000, sheets Karlovac (BENČEK et al., 2014), Sisak (PIKIJA, 1987), Slunj (KOROLIJA et al., 1980), and Bosanski Novi (ŠIKIĆ, 1990), as well as the 1:500.000 scale geological map (FEDERAL GEOLOGICAL SURVEY, 1970). A description of the subsurface geological composition was compiled, including a composite geological column that outlines the lithostratigraphic and chrono strati- graphic sequence of deposits in the area. The synthesis of existing data was undertaken using GIS and graphical editing tools. The geological setting of the study area predominantly comprises Late Palaeozoic, Triassic, and Plio-Quaternary to Quaternary deposits, which cover older Variscan bedrock and structures. In the SE part, Carboniferous deposits (C) are the oldest exposed rocks (Fig. 3). They are predominantly com- posed of clastic and subordinately carbonate deposits, includ- Figure 2. The structural map shows a simplified tectonic framework of the fault systems at the SW margin of the PBS. Fault abbreviations: SF – Sava fault; PF- Pokupsko fault; GF – Glina fault. Fault systems are compiled after KOROLIJA et al. (1980), VELIĆ & SOKAČ (1982), BUKOVAC et al. (1984), PIKIJA (1987), ŠIKIĆ (1990), PRELOGOVIĆ et al. (1998), TOMLJENOVIĆ & CSONTOS (2001), BENČEK et al. (2014), and HERAK & HERAK (2023). G eologia C roatica 295Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system ing shales, siltites, sandstones (greywacke, quartz-greywacke), and dolomitised and ankeritised limestones. Together with the Devonian deposits (D), characterised by thick lenses of clay, limonitised and schisty limestones within interbeds of shales, siltstone and sandstones form the pre-Permian low-grade meta morphic basement (ŠIKIĆ, 1990). Continuation of the post-Carboniferous deposition se- quence is characterised by a clastic sequence of younger Pala- eozoic age and is found on the southern slopes of Petrova gora Mt. (Fig. 3). Permian deposits (P) are represented by schists, quartz-greywacke sandstones, shales, and fine-grained con- glomerates (KOROLIJA et al., 1980). This turbidite-like com- plex was developed under conditions of rapid and constant in- filling within the existing basement basinal structures, reaching thicknesses of 500 m. At the same time, Permo-Tri- assic deposits (P, T) are consistent and composed of fine- grained brick-red sandstones and sandy-clay shales (KO- ROLIJA et al., 1981). The quartz-rich greywackes of the Upper Permian, observed along the southern slopes of the Petrova gora Mt., are a marker that represents transitional strata from the Upper Permian to the Lower Triassic (Fig. 3), which were deposited in shallower marine environments due to orogenic uplift. Lower Triassic (T1) deposition of mica sandstones, siltites, and shales with a gradual transition to carbonate marls and eventually limestones and dolomites continuously followed the Upper Palaeozoic normal superposition (KOROLIJA, 1981). These deposits can be observed on the eastern slopes of Petrova gora Mt., concordantly overlapping older Permian de- posits (Fig. 3). During the Middle Triassic (T2), evidence of consistent limestone and dolomite deposition can be observed. At the local scale, within the carbonate facies, tuffs, fine- grained sandstones, and interlayers of sheet limestone with chert alongside the shales were observed. The diminishing presence of clastic and pyroclastic components within the up- permost Middle Triassic succession suggests a reduction in tectonic activity and volcanism. During the Late Triassic pe- riod (T3), stable marine conditions enabled further massive Figure 3. Geological map of the study area of the THS according to the basic geological maps of the SFRY at a scale of 1:100.000, sheets Karlovac (BENČEK et al., 2014), Sisak (PIKIJA, 1987), Slunj (KOROLIJA et al., 1980), Bosanski Novi (ŠIKIĆ, 1988) and geological maps of SFR Yugoslavia, scale 1:500.000 (FEDERAL GEOLOGICAL SURVEY, 1970). Acronyms of the lithostratigraphic units: Q – Quaternary; Pl, Pl,Q – Pliocene, Plio – Quaternary; M - Miocene; Pc,E-Ol – Paleo- cene and Eocene-Oligocene; K, Pc – Cretaceous-Paleogene K – Cretaceous; J - Jurassic; T - Triassic; P,T – Permian - Triassic; P – Permian; C – Carboniferous; D – Devonian. Thermal springs’ locations in Bosnia and Herzegovina after HRVATOVIĆ (2005). The extent of the study area is shown in Figure 1. G eo lo gi a C ro at ic a 296 Geologia Croatica 77/3 carbonate deposition in the area. The Upper Triassic dolomites (W of Hrvatsko Žarište, Fig. 3) are usually in tectonic contact with the Lower – Upper Cretaceous limestones. According to SCHMID et al. (2008), Late Jurassic re- gional movements prompted the intraoceanic subduction of the Neotethys oceanic realm, which coexisted with the frag- mentation of the Adria Microplate. This led to the Adria Mi- croplate terrain differentiation, which resulted in variable sed- imentation patterns and continuous Jurassic-Cretaceous carbonate sedimentation in the area of the Adriatic Carbonate Platform (VLAHOVIĆ et al., 2005), while clastic-carbonate- volcanic sedimentation prevailed along its passive margins. The Jurassic rock complex in the study area (J2,3), part of this passive margin, primarily comprises a magmatic-sedimentary ophiolitic complex within the Central Dinaridic Ophiolitic Zone (SCHMID et al., 2008). This approximately 800 m thick complex consists of low-metamorphosed sedimentary rocks, (i.e., sandstones, shales, cherts, and some siltites, marly shales, and fine-grained limestones) and various magmatic rocks (i.e., basalts, gabbros, diabases; KOROLIJA et al., 1981; ŠIKIĆ, 1990). The Cretaceous-Paleogene transgressive sequence (Fig. 3) indicates passive margin sedimentation (KOROLIJA et al., 1981). Cretaceous deposits (K1) include limestones, greenish- gray marls, silicified sandstones, and carbonate breccias, fol- lowed by a flysch-turbidite like succession (K1,2) (VLAHOVIĆ et al., 2005). Paleogene deposits are mainly flysch-like, with carbonate marls, marly limestones, and thin layers of shales, marls, and fine-grained sandstones (K, Pc). Palaeocene and Eocene-Oligocene flysch-like deposits (Pc, E-Ol) crop out to a small extent in the eastern part of the study area. The Neo- gene-Quaternary sediment succession, linked to the tectonic evolution of the Croatian part of the PBS, was deposited in half-graben structures formed during the Early-Middle Mio- cene (PRELOGOVIĆ et al., 1998; TOMLJENOVIĆ & CSON- TOS, 2001; SAFTIĆ et al., 2003). Extensional tectonics led to basins filled with marine, lacustrine, and freshwater sediments (PAVELIĆ et al., 2003). This succession, about 1.1 km thick, includes conglomerates, sandstones, gravels, clays, marls, and limestones (M; Pl, Q; Fig. 3) (ŠIKIĆ, 1990). The youngest Pli- ocene and Quaternary deposits (Pl; Q), approximately 200 m thick, consist of conglomerates, sandstones, siltstones, sands, gravels, clays, and occasional interlayers of clay and coal, forming a final terrigenous/alluvial cover. 2.3. Hydrogeological setting In the Topusko area, there are three natural artesian thermal springs with a total capacity of approx. 25 l/s (BAĆ & HERAK, 1962; BAHUN & RALJEVIĆ, 1969) and temperatures ranging from 46 °C to 53 °C. Three exploitation wells were drilled to depths of up to 250 m in the immediate vicinity of the natural springs during the 1980s. They are currently used for heating, recreational, and medicinal purposes, with a total yield of 200 l/s and pressure of around 1.4 bar (PAVIĆ et al., 2023). A decrease in pressure from 2.18 bar (1978) to 1.52 bar (1982) was observed and attributed to overexploitation (ŠEGOTIĆ & ŠMIT, 2007). The Topusko thermal water shows a slightly acidic character, with an average pH of 6.5 – 6.8. The electrical conductivity ranges from 582 μS/cm to 680 μS/cm, and the total dissolved solids are approx. 500 mg/l indicating a medium to low mineralised fresh water. Hydrochemical analyses show a Ca-HCO3 hydrochemical facies (PAVIĆ et al., 2023, 2024) indicating water flow in carbonate rocks. Stable water isotopes δ2H and δ18O suggest a meteoric origin for thermal water and recharge during colder climatic conditions The uniform δ18O values indicate deep circulation and a large, thick aquifer, which homogenises seasonal precipitation variations over longer residence times. The estimated mean residence time of approximately 8.5 – 9.5 kyr is based on 14C content in DIC (PAVIĆ et al., 2024). Tritium activity in the thermal water is generally below the detection limit, but trace levels were detected after extensive abstraction for district heating during winter months, suggesting the potential infiltration of modern precipitation through the semi- permeable Neogene cover sequence of the aquifer. The most plausible equilibrium aquifer temperature was estimated at 90 °C using quartz geothermometers (PAVIĆ et al., 2023, 2024). Hydrodynamic measurements in the Topusko discharge area estimated an aquifer transmissivity of approx. 2 x 10-2 m2/s (PAVIĆ et al., 2023). According to ŠIMUNIĆ (2008), a set of faults forming a block in the shape of a three-sided prism enabled the uplifting of the aquifer. The result of an electrical resistivity survey identified the damaged zones of these faults in the spring area (PAVIĆ et al., 2023). From a hydrogeological point of view, fractured and karstified, highly permeable, Middle and Upper Triassic carbonates are the main geothermal aquifer of the THS. Lower Triassic and Permian-Triassic carbonate and clastic rocks represent the semi-confining units at the base of the reservoir, while the Palaeozoic complex is the impervious basement at the bottom of the stratigraphic sequence due to its low-grade metamorphism. Furthermore, low permeability Jurassic ophiolitic complex and Neogene-Quaternary deposits cover the geothermal aquifer. However, a limited number of deep boreholes hinders a more comprehensive understanding of local and regional scale hydrogeological conditions. Different conceptual models of the THS have been previously proposed: 1) In the early 20th century, it was considered that the thermal water at the Topusko springs was heated by a magmatic body based on the existence of basalt outcrops in the vicinity (e.g., Lasinja GORJANOVIĆ-KRAMBERGER, 1905, 1917). To support the volcanic origin of the thermal water in Topusko, GORJANOVIĆ-KRAMBERGER (1917) postulated that mountains S of Topusko could not produce such hydrostatic pressure that would raise water from a depth of approx. 1.5 km, according to the thermal water temperature at the surface. However, these basalts are of Mesozoic age (MAJER, 1978; 1993) and Quaternary magmatic bodies behaving as a recent heat source for hydrothermal systems are absent within the PBS. Therefore, this conceptual model was abandoned. 2) BAĆ & HERAK (1962) interpreted the Palaeozoic rocks as behaving as an impermeable rock complex that regulates regional water circulation and defines the watershed of the THS together with Lower Triassic and Neogene formations. According to the same authors, permeable Middle G eologia C roatica 297Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system and Upper Triassic carbonate deposits facilitate deep groundwater circulation, especially in tectonically active areas. The thermal aquifer receives recharge in the Glina river headwaters, where Triassic carbonates crop out. The proposed primary regional flow direction is from the S to N, from the Glina area towards the Topusko depression. 3) ŠIMUNIĆ (2008) proposed that the THS recharge area is located W of the Petrova gora Mt., where Triassic carbonates crop out. According to the proposed model, the water circulates from W to E below the Petrova gora nappe, heats up by the geothermal gradient, and discharges in Topusko due to the highly permeable fault damaged zones. 2.4. Structural-geological research Structural-geological fieldwork was conducted through field campaigns in 2021 and 2022. Field investigations involved structural and lithological data acquisition at 162 locations covering an area of approx. 2,000 km2 (Figs. 2 and 3). During the fieldwork, observation points were archived and processed in the Avenza PDF Maps application (URL 1) and MS Excel field database, while database creation and geospatial positioning of the collected geological data were performed using the software ArcMap 10.1 (URL 2). The structural investigations detailed the stratigraphic relationships among the main lithological units and the analysis of the tectonic and structural settings in the study area. They encompassed geometric and structural measurements of strata bedding, fractures/joint systems, and fault/shear planes. In particular, measurements of fault/shear planes included the identification of fault kinematic indicators and determination of the tectonic relationships between the observed structures. Collected structural data were used in the construction of three NE-SW striking regional structural-geological profiles (THS-1 to THS-3; Figs. 2–4). These profiles represent surface/sub- surface 2D models of the distribution of geological units and structures in the area of Petrova gora Mt. (SW) and the Glina depression (NE) and they were used here for recon struc tion of the geological and structural settings of the THS. 2.5. Conceptual and numerical modelling In this study, we initially tested and compared the validity of the conceptual model proposed by ŠIMUNIĆ (2008) with the results of structural investigations and the constructed geological profiles. Furthermore, these data combined with the results of previous investigations in the Topusko area were used to update the conceptual model of the THS, which served as a base for a 2D numerical model of heat flow and fluid transport. The conceptual model was developed including: lithological field data, geological maps, hydrodynamic and electrical resistivity survey data acquired in the discharge area of Topusko, chemical and isotopic compositions of both the thermal water and the precipitation in the assumed recharge area. The proposed conceptual model was tested by performing numerical simulations of f luid f low and heat transport accounting for the variation of the thermal water density due to the temperature distribution in the subsurface (BUNDSCHUH & CÉSAR SUÁREZ, 2010; DIERSCH, 2014). The numerical model solves the governing equations of fluid flow and heat transport, obtaining the distribution of the primary variables (i.e., hydraulic head and temperature) within the modelling domain and over time. In this research, we used the FEFLOW 7.4 software (DIERSCH, 2014), a Finite Element- based simulator specifically designed to address fluid flow and transport phenomena within porous and fractured media. The fundamental equations governing fluid flow and heat transport in these media are formulated based on conservation principles for fluid mass, momentum, and thermal energy (DIERSCH, 2014). The equations are coupled by setting a functional form behaviour for fluid density and viscosity depending on the respective primary variable of the problem. For the specific problem at hand, fluid viscosity was kept constant. In contrast, the density of the fluid varies linearly with temperature, following r = r0 (1 – aDT) with a being the volumetric thermal expansion coefficient of the fluid (POLA et al., 2020). 3. RESULTS 3.1. Structural-geological profiles Structural – geological profiles THS-1 to THS-3 (Fig. 4) represent an interpretation of the surface/subsurface relationships in the study area down to an investigation depth of approx. 6 km. Constructed geological profiles include three structural domains that convey a system of thrust faults separating the External Dinaridic lithological units in the SW parts of the profiles from the Internal Dinaridic lithological units. The Internal Dinaridic units can be further subdivided into two subdomains. The SE and central part of the profiles included the Petrova gora area, while the Glina depression characterises the NE profile domains (Fig. 4). The profiles (Figs. 2 and 4) pinpoint complex faulted and folded structures that are predominantly affected by low angle thrust faults (dip angle ≤30°; faults R1, R2, and R3 in Fig. 4) and reverse faults in their immediate hanging wall (faults R5 and R2 in Fig. 4). Besides these faults, higher angle normal and tectonically inverted normal faults are observed (dip angle ≥45°; faults R4 and R6 in Fig. 4). Constructed profiles show that the THS subsurface is composed of folded structures, with the largest Topusko anticline formed in the hanging wall of the R3 thrust fault and R5 blind fault (profile THS-3 in Fig. 4). The Petrova gora structure also resembles a larger-scale anticlinal structure that is tectonically uplifted in the hanging walls of the R1 and R2 thrust faults (profile THS-1 in Fig. 4). The Petrova gora structure, composed of a Palaeozoic clastic sequence, is tectonically uplifted by approx. 3.5 km in relation to its R1 footwall units, which resemble an undeformed, 6 km thick, Palaeozoic-Mesozoic succession. Towards the NE, the Petrova gora structure is additionally vertically displaced by tectonically inverted normal faults (faults R3 and R4; profile THS-1 in Fig. 4). Here, vertical displacement along the R3 and R4 normal faults accommodates NE-SW oriented extension, with approx. 1.5 – 2 km vertical subsidence. Further to the NE, the R4 fault delineates a folded system that incorporates the area of the Topusko anticline. With an approx. 4 km thick Palaeozoic-Mesozoic succession, this folded system (profile THS-3 in Fig. 4) is mainly composed of the Carboniferous-Permian clastic succession that is covered G eo lo gi a C ro at ic a 298 Geologia Croatica 77/3 by Triassic clastic-carbonate rocks. The Topusko anticline hinge zone is partly covered by the Jurassic ophiolitic complex and its transgressive Cretaceous clastic-carbonate succession (see Figs. 3 and 4, profile THS-3). In the area of the Topusko anticline, both the Cretaceous and Jurassic units are significantly reduced due to extensive uplift along the R5 blind Figure 4. Structural-geological profiles in the THS area. The NE-SW striking profiles are perpendicular to the geological structures in the research area. The profiles show the Topusko anticline formed in the hanging wall of the thrust fault R3 that behaves as a low-angle detachment surface. The anticline is composed of a thick Palaeozoic-Mesozoic sequence, which is partly eroded in the immediate area of Topusko town. Towards the NE, it is covered by the thick Paleogene-Neogene sedimentary succession of the Glina depression. Horizontal and vertical scale ratio is 1:1. G eologia C roatica 299Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system fault, which resulted in hanging wall tectonic erosion. The final NE domains of the THS geological profiles incorporate the gently dipping NE limb of the Topusko anticline that is faulted by the NE dipping R6 normal fault. The R6 normal fault (vertical displacement here is approx. 500 m; profile THS-3 in Fig. 4) is part of the extensional structure of the Glina depression, which is one of the extensional basins formed during the Neogene extension of the Croatian part of the PBS. As a result, the NE segments of the THS profiles are covered by an additional approx. 2 km thick clastic-carbonate succession deposited during the Paleogene and Neogene within the Glina depression. Reconstructed THS subsurface relationships along the profiles indicate that the observed faults and cogenetic folded structures have a history of polyphase deformation. The same asymmetric anticlinal structures (e.g., Petrova gora) in the hanging wall of the low angle detachment faults (e.g., R3 and R4 faults, Fig. 4) show indications of tectonic transport towards the SW, localised uplift and erosion (e.g., Topusko anticline), as well as NE – SW oriented extension. This implies that the interpreted faults accommodated both initial Cretaceous- Paleogene NE – SW compression and the following Neogene NE – SW extension. 3.2. Conceptual and numerical modelling of groundwater flow and heat transport in the THS 3.2.1. Novel conceptual model of the THS The new conceptual model of the THS is established based on constructed geological profiles (Fig. 4) and previous hydrogeological and geophysical research (PAVIĆ et al., 2023, 2024). These data lead to development of an interpolated schematic geological profile, which refines the understanding of the THS (Figs. 5 & 6). Analysis of the geological profiles Figure. 5 Trace of the interpolated schematic geological profile of the THS (black line). Sinkholes (white circles) in the SW part of the profile indicate higher karstification of the carbonate formations, representing the potential recharge area of the THS. G eo lo gi a C ro at ic a 300 Geologia Croatica 77/3 here, challenges the previously proposed hypothesis by ŠIMUNIĆ (2008) about the recharge area of the THS being W of the Petrova gora nappe. The Triassic carbonates cropping out W of the Petrova gora Mt. are hydrogeologically isolated from the discharge area in Topusko by the impermeable Petrova gora structure, which is composed of a Palaeozoic clastic sequence with low-grade metamorphism. Therefore, a regional groundwater flow from W to E becomes highly unlikely. Instead, in our conceptual model, attention is directed towards a new hypothesis: the recharge occurs S from Topusko, where Triassic carbonates are partly exposed at the surface. The surface manifestation of thermal water in the Topusko spring area can be understood through the concept of gravity- driven groundwater flow, determined by: i) Middle-Upper Triassic carbonate complex rocks, which are the main thermal aquifer of the THS, ii) Palaeozoic-Lower Triassic and the Jurassic-Neogene formations representing the semi-confining layers at the bottom and top of the reservoir, respectively, iii) a large area characterised by intense karstification occuring approx. 13 km S of Topusko within the Triassic carbonates (Fig. 5), facilitating the infiltration and deep circulation of the meteoric water with a long residence time, and iv) the fault- thrusted regional tectonic setting and the existence of substantial hydraulic boundaries (e.g., the R3 and R4 thrust faults and near-surface fault zones), which influence the direction of the groundwater flow. The Topusko anticline within the hanging wall of the R3 thrust fault facilitates the uplift of the aquifer, bringing it closer to the surface in the thermal springs area. Cogenetic faults and fracture networks in the hinge zone of the Topusko anticline, probably affected by an extensional regime, increase the fracturing of the bedrock and the permeability field in the aquifer enabling the thermal water outflow. The occurrence of fault damage zones in the Topusko subsurface was determined by the electrical resistivity surveys (PAVIĆ et al., 2023). Hydrogeochemical research corroborates the meteoric origin of the Topusko thermal water and the interaction with the carbonate aquifer. A residence time of approx. 9 kyr is evidenced by the 14C analysis of the thermal water (PAVIĆ et al., 2024). Additionally, the consistent major ion composition of the thermal water over a two-year period suggests a large and stable system. Based on all findings, it is suggested that the novel THS conceptual model favours a regional groundwater flow direction from S to N. 3.2.2. Setup of the numerical model The physical validity of the proposed conceptual model of the THS was tested by conducting 2D numerical modelling of both the regional fluid flow and heat transport. To simulate gravity-driven groundwater flow in the THS, we adopted a shortened version of the schematic geological profile (Fig. 6), excluding the SW portion of the R3 fault. The recharge area of the THS is located NE of R3, where numerous sinkholes occur (Fig. 5). Therefore, the SW portion of the section is not strictly connected to the hydrothermal system, and was not reproduced in the numerical model to diminish the computa- tional effort for solving the numerical simulations. The profile was digitised in GIS software and used to generate a su- Figure 6. Interpolated geological profile of the THS. The profile is perpendicular to the Topusko anticline formed in the hanging wall of the thrust fault R3 that behaves as a low angle detachment surface. The same fault, combined with the R4 and R5 reverse faults, lifts the reservoir up to shallow depths in the immediate area of Topusko town. Horizontal and vertical scale is 1:1. G eologia C roatica 301Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system permesh in the FEFLOW software outlining the geometry of the geological units within the modelling domain. The domain was discretised through a triangular mesh employing the Tri- angle triangulation code (DIERSCH, 2014), which permits honouring the complex geometrical relationships of the units. The mesh was refined in the hinge zone of the Topusko anti- cline, in the aquifer unit, and along the main faults. Such re- finement allows a better discretisation in the parts of the model where numerical instability is expected. The obtained mesh (Fig. 7) comprised 19,241 triangular elements and 9,831 nodes. Geological units with similar hydrogeological and thermal properties were grouped, obtaining seven hydrostratigraphic units: i) Palaeozoic low-metamorphic sedimentary complex (B), ii) Permian and Permian-Triassic clastic units (P), iii) Early Triassic sedimentary rocks (T1), iv) Middle and Late Tri- assic carbonates (T2,3), v) Jurassic magmatic – sedimentary ophiolitic complex and Jurassic and Cretaceous carbonates (J), vi) Neogene siliciclastic units (N), and vii) Quaternary alluvial cover (Q). In particular, the T2,3 hydrostratigraphic unit repre- sents the main thermal aquifer of the THS. The equivalent po- rous medium approach was employed for the hydrogeological and thermal parametrisations of the hydrostratigraphic units (DIERSCH, 2014; ANDERSON et al., 2015). This approach considers the unit as a porous medium with homogeneous and isotropic properties being suitable for the numerical modelling of regional groundwater flow systems in heavily fractured and karstified carbonate rocks (TEUTSCH & SAUTER, 1991; SCANLON et al., 2003; GHASEMIZADEH et al., 2012; MÁDL-SZŐNYI & TÓTH, 2015). In a regional numerical model, the scale of these discontinuities is smaller than the representative elementary volume (i.e., the mesh size), result- ing in a constant and homogeneous value of the considered parameter. For assigning the hydraulic and thermal properties (Table 1), we considered the main lithology of the hydrostrati- graphic units and relied on data from the literature (CERMAK & RYBACH, 1982; DOMENICO & SCHWARTZ, 1997; FET- TER, 2001; STYLIANOU et al., 2016; BOROVIĆ et al., 2018; STOBER & BUCHER, 2021; LALOUI & ROTTA LORIA, 2020), parametrisation of similar units in hydrothermal sys- tems within the PBS (RMAN & TÓTH, 2011; MÁDL- SZŐNYI & TÓTH, 2015; HAVRIL et al., 2016), and datasets obtained from field investigations in Topusko (PAVIĆ et al., 2023). The hydraulic conductivity (K) and the porosity (φ) were assigned accounting for the role of the units in the THS. Excluding the loose sediments of the Quaternary cover (Q unit; Fig. 3), the highest K and φ values were assigned to the T2,3 hydrostratigraphic unit representing the thermal aquifer. Conversely, the impervious basement B was reproduced using the lowest K and φ values. In the hinge zone of the Topusko anticline, K and φ were increased simulating the impact of the local scale fracturing that favours the upwelling of the thermal water. Open fractures can increase the permeability and po- rosity fields enhancing the fluid flow in the aquifer (i.e., FAULKNER et al., 2010; WORTHINGTON et al., 2019; POURASKARPARAST et al., 2024). The regional values of K and φ were increased by: i) two orders of magnitude and 10%, respectively, for the T2,3, T1, and P units, and ii) one or- der of magnitude and 5%, respectively, for the B unit consid- ering the decrease of fracture apertures with depth. In particu- lar, the K value of the T2,3 unit in the hinge zone of the Topusko anticline was calculated from a collection of transmissivity values obtained from well tests conducted in Topusko (PAVIĆ et al., 2023). Furthermore, the thermal conductivities (λ) of the N and T2,3 units were derived from RMAN & TÓTH (2011), MÁDL-SZŐNYI & TÓTH (2015), and HAVRIL et al. (2016), who investigated regional groundwater flow and heat transport patterns in the PBS. Boundary conditions (BCs) for fluid flow and heat trans- port were applied at the border of the modelling domain fol- lowing the conceptual model of the THS. The fluid flow BCs (Fig. 7A) included: i) a 2nd kind (Neumann) BC at the SW part of the top domain, and ii) a 3rd kind (Cauchy) BC at the top of the domain in the Topusko area. The Neumann BC produces an inflow in the modelling domain and was employed to sim- ulate the recharge of the THS. The condition spanned for a length of 2 km, which corresponds to the extension of the area with a higher density of sinkholes along the modelling domain. An inflow of 90 mm/yr was imposed considering an effective infiltration of 10 % (as being used for the simulation of similar hydrothermal systems in a carbonate reservoir in Croatia; BOROVIĆ et al., 2019) and an average annual precipitation of 900 mm (DHMZ, 2021; MARTINSEN et al., 2022). The Cauchy BC simulates a fluid flux through the model boundary depending on: i) the difference between the hydraulic head imposed at the boundary and the value calculated within the domain, and ii) a transfer-rate coefficient (DIERSCH, 2014). Since it is generally used to simulate the variable outflow in springs (ANDERSON et al., 2015), it was employed in the THS numerical model to reproduce the occurrence of thermal springs in the Topusko area. The imposed hydraulic head was set as the average ground elevation (122 masl), while the trans- fer rate was calculated as the K of the Q unit (Table 1) divided by the thickness of the alluvial cover in Topusko obtained from the stratigraphic well logs. The heat transport BCs (Fig. 7B) included: i) a 1st kind (Dirichlet) BC at the SW and NE verti- cal boundaries of the modelling domain and at the top except for the Topusko area, ii) a 2nd kind (Neumann) BC at the bot- tom, and iii) a 3rd kind (Cauchy) BC at the top of the domain in the Topusko area. The Dirichlet BC imposes a constant value of the primary variable. In the THS model, it was used to reproduce both the ground temperature at the surface and the increasing temperature with depth due to the regional ge- othermal gradient at the vertical boundaries. The temperature imposed at the top was set to 10 °C, which corresponds to the mean average annual air temperature, while the values at the lateral boundary were obtained from the initial distribution of temperature. The Neumann BC was used to simulate the re- gional inflow of heat from the deeper parts of the crust. A value of 100 mW/m2 was imposed being consistent with the average heat flow in the Croatian part of the PBS (LENKEY et al., 2002; HORVÁTH et al., 2015). Similarly to the fluid flow, the Cauchy BC reproduces a heat flux through the boundary depending on: i) the difference between the reference and the calculated temperature values, and ii) a transfer-rate coefficient (DIERSCH, 2014). The reference temperature was set to 10 °C following the value imposed for the Dirichlet BC, while the transfer rate was calculated as the λ of the Q unit (Table 1) di- vided by its thickness. G eo lo gi a C ro at ic a 302 Geologia Croatica 77/3 The THS numerical modelling was conducted performing transient state numerical simulations that account for the var- iations of the hydraulic head and temperature distributions over time. Therefore, it is important to set appropriate initial conditions for these variables achieving a better convergence of the numerical solution at the initial stages of the simulation. The initial value of the hydraulic head was set as 122 m, cor- responding to the hydraulic head imposed at the fluid flow Cauchy BC. The initial distribution of the temperature was obtained through a steady-state simulation using: i) the de- scribed set of hydrogeological and thermal parameters for the hydrostratigraphic units (Table 1), ii) a 1st kind (Dirichlet) BC with a constant temperature value of 10 °C at the top of the domain, and iii) a 2nd kind (Neumann) BC at the bottom with imposed inflow of 100 mW/m2. The resulting temperature dis- tribution reproduces a regional geothermal gradient of 35.8 °C/km, being consistent with available regional values in the study area (30 – 40 °C/km; MACENIĆ et al., 2020). 3.2.3. Simulation results The best simulation of the THS was conducted employing the described set of hydrogeological and thermal parameters for the hydrostratigraphic units (Table 1) and the imposed boundary and initial conditions (Fig. 7). The transient simulation was run for 100 kyr obtaining a quasi-stationary distribution of the primary variables. This condition reproduces the long-term natural state of the hydrothermal system (i.e., GARG et al., 2007; KAISER et al., 2013; HAVRIL et al., 2016; Figure 7. Fluid flow (A) and heat transport (B) boundary conditions (BCs) in the THS numerical model. Blue circles – 1st kind (Dirichlet) BC; Red crosses – 2nd kind (Neumann) BC; Green circles with crosses – 3rd kind (Cauchy) BC. Horizontal and vertical scale ratio is 1:1. Table 1. Hydrogeological and thermal parameters assigned to the hydrostratigraphic units in the THS numerical model. Unit Age Lithology Thick K Ss j l rc [m] [m s-1] [m-1] [%] [W m-1 K-1] [MJ m-3 K-1] Q Quaternary Gravel, sand, silt 90 1*10-4 1*10-4 20 1.6 1.3 N Neogene Sand, clay, silt, marls < 2,000 2*10-7 1*10-4 5 1.8 1.6 J Jurassic Magmatic-sedimentary ophiolitic complex 750 – 1,600 2*10-10 3.3*10-7 5 2 2.2 T2,3 Middle-Upper Triassic Limestone, dolostone 400 – 1,200 2*10-6 3.3*10-7 10 2.5 2.52 T1 Lower Triassic Sandstone, marl, limestone, dolomite 600 2*10-9 3.3*10-7 5 2 1.6 P Permian-Triassic Siliciclastic rocks 650 2*10-11 3.3*10-7 5 3 2.1 B Palaeozoic Schists shales, sandstones, conglomerates with limestone and dolomite lenses 2,000 – 3,000 2*10-12 3.3*10-7 2 3.5 2.1 Note: K: Hydraulic conductivity (DOMENICO & SCHWARTZ, 1997; FETTER, 2001; PAVIĆ et al., 2023); φ: Porosity (DOMENICO & SCHWARTZ, 1997); λ: Thermal conductivity (CERMAK & RYBACH, 1982, RMAN & TÓTH (2011), MÁDL-SZŐNYI & TÓTH (2015), STYLIANOU et al., 2016, BOROVIĆ et al., 2018; STOBER & BUCHER, 2021); Ss: specific storativity (DOMENICO & SCHWARTZ, 1997); ρc: volumetric heat capacity (LALOUI & ROTTA LORIA, 2020) G eologia C roatica 303Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system BOROVIĆ et al., 2019; POLA et al., 2020; TORRESAN et al., 2022), which is the goal of the THS numerical modelling. The temperature variations over time at different depths in the subsurface of Topusko are represented in Fig. 8A. The simulated temperature gradually increases from the initial distribution reaching maximum values between 18 and 30 kyr depending on the considered depth. The temperature at the surface peaks to 42.6 °C at 18 kyr (initial value of 10 °C; depth = 0 km in Figs. 8A and 8B), while it reaches the maximum value of 77.4 °C at the same simulation time at base of the aquifer (initial value of 55.8 °C; depth = 1 km in Figs. 8A and 8B). The temperature differential progressively decreases with depth, being up to 2 °C in the deeper part of the modelling domain (depth = 3.5 km). After this increasing phase, the simulated temperatures progressively decrease. The final modelled temperatures at the surface and base of the aquifer are 31.3 °C and 59.5 °C (Figs. 8A and 8B), respectively, de pict- ing a drop of 11.3 °C and 17.9 °C, respectively. The drop from the peak value decreases with the depth as well and is 3 °C in the deeper part of the modelling domain. Despite these variations, a practically constant temperature distribution is observed between 80 and 100 kyr with a maximum variation of 1.4 °C at the aquifer base. This result suggests that a quasi- stationary state of the solution was achieved. The regional and local scale spatial distributions of temperature at the end of the simulation period (Figs. 9A and 9B, respectively) were further detailed. The temperature distribution within the aquifer at both regional and local scales exhibits significant spatial variability. In particular, a decrease from the initial values up to 14 °C is observed in the recharge area of the system where the infiltration of cold waters (10 °C as imposed by the Dirichlet BC for heat transport) occurs. Similar behaviour is observed in the recharge area of regional, topographically-driven groundwater systems (i.e., DOMENICO & PALCIAUSKAS, 1973; ANDERSON, 2005; AN et al., 2015). However, the temperature within the aquifer (pink polygon in Fig. 9) is up to 25 °C (depth of 0.7 km) as expected in shallow carbonate aquifers. The temperature in the aquifer progressively decreases eastward as depicted by the 20 °C isotherm located at a depth of approx. 1.6 km in the central part of the modelling domain. This decrease is connected to the deep circulation of the infiltrated waters in the aquifer driven by both the topographic gradient and the deepening of the strata to the SW of Topusko (Fig. 6). The temperature increases in the hinge zone of the Topusko anticline where the upwelling of the thermal water favours the rise of the isotherms toward the surface (Fig. 9B). In particular, the isotherm of 35 °C is located approx. at the top of aquifer (depth of 0.1 km), while the isotherm 50 °C approaching to the water temperature of the Topusko thermal springs is located at a depth of 0.4 km. The modelled temperature is: i) up to 44 °C in the central part of the thermal anomaly at the depth investigated by the thermal wells in Topusko (0.2 km), ii) between 56 and 74 °C at the base of the aquifer, and iii) between 80 and 100 °C at the base of the T1 hydrostratigraphic unit, which represents the semi-confining unit at the bottom of the aquifer (depth of approx. 1.8 km; Fig. 9). The circulation of the fluid is generally directed from the recharge area of the model to the outflow area in Topusko. The flow mostly occurs in the T2,3 hydrostratigraphic unit where an average Darcy velocity of 0.18 m/yr (5.02×10-4 m/d) is modelled. It drops to 0.01 m/yr (2.9×10-5 m/d) in the T1 and P hydrostratigraphic units at the bottom of the aquifer, and to 3.91×10-7 m/yr (1.07×10-9 m/d) in the impervious basement (B hydrostratigraphic unit). The direction and velocity of the flow in the aquifer show different distributions from the recharge to the outf low area depending on the K value of the hydrostratigraphic unit and the thickness of the aquifer. In the recharge area, the flow is almost horizontal with an average Darcy velocity of 0.27 m/yr (7.3×10-4 m/d). The velocity in the recharge area shows a high variability spanning over almost three orders of magnitude (minimum and maximum of 7.06×10-4 m/yr to 1.77 m/yr, respectively). In particular, the highest values are observed NE of the R3 fault (Fig. 6), where the thickness of the aquifer decreases (Fig. 9A) due to the anticline between the R3 and R4 faults. In the flow-through part of the system, the flow (average velocity of 0.17 m/yr corresponding to 4.54×10-4 m/d) is generally downward following the deepening of the aquifer to the SW of Topusko and resulting in the local drop of the isotherms (as depicted by Figure 8. Modelled temperature variations over time (A) and final temperature distribution (B) in the Topusko area. G eo lo gi a C ro at ic a 304 Geologia Croatica 77/3 the 20 °C isotherm; Fig. 9A). The flow is generally upward in the southern limb of the Topusko anticline, but with similar magnitudes of the horizontal and vertical components of the Darcy velocity (average value of 0.15 m/yr corresponding to 4.16×10-4 m/d). In the Topusko area, the flow is upward and with a prevalent vertical component of the Darcy velocity. The average velocity value is of 0.13 m/yr (3.5×10-4 m/d), being slightly lower than the values observed in other parts of the aquifer, but the velocity field is more evenly distributed with minimum and maximum values of 1.21×10-3 and 0.28 m/yr, respectively. The Darcy velocity increases up to 0.28 m/yr in the upper part of the aquifer (i.e., the depth investigated by the thermal wells in Topusko) pointing to the quick rise of the thermal waters in the final part of the circulation path. The outflow modelled by the 3rd kind (Cauchy) BC is 0.18 m3/yr. 4. DISCUSSION The Topusko hydrothermal system area is located at the tec- tonic boundary between the External and Internal Dinarides in the vicinity of the Sava suture zone (SCHMID et al., 2020). The Palaeozoic-Mesozoic terrains in the study area experi- enced a complex polyphase tectonic evolution through the Cre- taceous-Paleogene and the Neogene-Quaternary, which re- sulted in tectonic overprints of both regional shortening and extension. Observable regional shortening is associated with several NW-striking low-angle thrust faults and reverse faults (Fig. 4) that were probably formed due to the Cretaceous-Pale- ogene Adria Microplate and European Foreland collision (SCHMID et al., 2020). This regional NE – SW oriented short- ening (TOMLJENOVIĆ & CSONTOS, 2001; SCHMID et al., 2008; USTASZEWSKI et al., 2008) resulted in the tectonic uplift of the Petrova gora Mt. and the formation of a folded system that also incorporates the Topusko anticline (Fig. 6). The polyphase tectonic history is further evidenced by struc- tural reactivations of these faults and their tectonic inversions, which is observable in the tectonically inverted normal faults R3 and R4 (Fig. 4) indicating a NE – SW extension in the THS area. According to TOMLJENOVIĆ & CSONTOS (2001), NE – SW stretching of the inherited structures within the Croatian part of the PBS occurred during the Early-Middle Miocene due to back-arc type extension. Tectonic rejuvenation of the THS area through N – S regional contraction is observable at Figure 9. Regional temperature distribution in the modelling domain (A) and in the hinge zone of the Topusko anticline (B). A local temperature anomaly in the Topusko subsurface occurs reaching values of up to 74 °C at the base of the THS carbonate aquifer (pink polygon). The horizontal and vertical scale ratio is 1:1. G eologia C roatica 305Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system the outcrop scale, with evidence of re-folding processes and mapped reverse faults within the Pliocene-Quaternary clastic sequence (e.g., USTASZEWSKI et al., 2008; HERAK et al., 2009; HERAK & HERAK, 2023). The interpretation of the geological profiles constructed in this study suggests that the most used THS conceptual model (ŠIMUNIĆ, 2008) was based on an out-dated geologi- cal reconstruction. Furthermore, Šimunić’s model was not constructed using a series of stuctural-geoloical profiles and a detailed reconstruction of subsurface relationships, but only spatial distribution of the lithological units. These shortcom- ings hindered the detailed understanding of complex subsur- face relationships in the Topusko area. The new structural- geological interpretation was combined with geochemical and hydrogeological data (PAVIĆ et al., 2023, 2024) to propose a novel conceptual model of the THS. This model suggests that the THS is an intermediate-scale, gravity-driven, tectonically – controlled, hydrothermal system hosted in a Middle-Upper Triassic carbonate rock complex. The thermal water is of me- teoric origin infiltrating approx. 13 km to the S of Topusko, where an area with intense karstification depicted by a high density of sinkholes occurs (Fig. 5). These sinkholes play a significant role in facilitating the diffuse recharge of the THS aquifer. Our interpretation in general coincides with BAĆ & HERAK's (1962) hypothesis, which identifies Middle and Up- per Triassic carbonate deposits as the aquifer units facilitating deep groundwater circulation. They proposed that the thermal aquifer receives recharge from the Glina river headwaters, where Triassic carbonates are exposed, with a primary re- gional flow direction from S to N, aligning with our findings regarding the recharge area S of Topusko. Regional folded structures favour both the deep circulation of the infiltrated water, warming due to increased heat f low of the PBS (HORVÁTH et al., 2015), and the upwelling of the thermal wa- ter in the Topusko area. Fault and fracture systems in the hinge zone of the Topusko anticline promote the fracturing of the bedrock, the increase of the permeability field in the aquifer, and the quick rise of the thermal water from the deeper part of the reservoir. The extensional regime in the anticline hinge zone, combined with cogenetic fault and fracture systems and their permeable damage zones (PAVIĆ et al., 2023), creates a synergistic effect. This synergy enhances the overall perme- ability field in the discharge area, leading to the formation of abundant natural thermal springs. The physical validity of the THS model was constrained conducting variable-density numerical simulations of fluid flow and heat transport. The numerical model was populated using a set of hydrogeological and thermal parameters from both the literature and field datasets. Boundary conditions were im- posed at the borders of the modelling domain following the conceptual model of the THS. The best simulation was run for 100 kyr developing a local increase of the temperature distri- bution in the subsurface of Topusko (Fig. 9). In particular, the modelled temperatures were 31.3 and 44 °C at the surface and at the depth investigated by wells in Topusko, respectively. De- spite these values are approx. 20 °C lower than the observed temperatures at the same depths, they can be considered as ac- ceptable due to the simplifications imposed during the con- struction of the numerical model. The thermal aquifer was re- produced using an equivalent porous medium approach, which assigns homogeneous values of the hydrogeological properties, while localised highly permeable damage zones channelling the outflow of the thermal water were not considered. The im- plementation of such structures in the numerical model would increase the fluid flow promoting the local increase of the tem- perature. In addition, a temperature between 80 and 100 °C was modelled in the semi-confining unit at the bottom of the aqui- fer being comparable with the reservoir temperature of 90 °C calculated with geothermometers. As discussed for the surficial temperature, the effect of localised highly permeable fault zones was not included in the conducted numerical model. These areas could favour the raising of the isotherms in the deeper part of the reservoir resulting in the higher temperature values calculated with geothermometers. The simulated temperature distribution at the local scale reproduced the field values and was almost constant during the final 20 kyr of the simulation time (Fig. 8A). This result suggested that a quasi-stationary condition reproducing the natural state of the THS was achieved for a time span that is almost twice the mean residence time of the Topusko thermal water. The simulation results indicate a predominant ground- water flow pattern towards the Topusko discharge area. They also highlight the significant influence of the hydraulic con- ductivity field on groundwater flow patterns within the system. The modelled Darcy velocities ranged from 0.13 m/yr (To- pusko) to 0.27 m/yr (recharge area) in the aquifer and de- creased within the hydrostratigraphic units below it, directly reflecting the contrasting hydraulic conductivity values of dif- ferent hydrostratigraphic units. This emphasises the impor- tance of incorporating accurate K data into groundwater flow models for reliable predictions. An outflow of 0.18 m3/yr was modelled in the Topusko area, reproducing the thermal springs. Historical data on the flow rate of the springs before the drill- ing of the wells are not available and a more detailed numeri- cal model would be needed to reproduce the effect of exploita- tion of the thermal water. The localised thermal anomaly in the anticline hinge zone corroborates the structural-causative processes enhancing fluid flow and heat transfer in the To- pusko subsurface. Therefore, the results of the numerical mod- elling corroborate the validity of the conceptual model of the THS. In order to constrain the processes favouring the devel- opment of the modelled temperature distribution, the Rayleigh number (Ra) was calculated following the procedure described in TURCOTTE & SCHUBERT (1982). Ra is a dimensionless parameter indicating the threshold for the onset of free con- vection occurring when it exceeds the critical value of approx- imately 40 (TURCOTTE & SCHUBERT, 1982; NIELD & BEJAN, 1999; PESTOV, 2000; MÁDL-SZŐNYI & TÓTH, 2015). Considering the thickness of the reservoir (approx. 1 km) and its hydrogeological and thermal parametrisations, the conditions for the development of convection cells in the THS predominantly exist in the more permeable hinge zones of the Topusko anticline (Ra = 5,500 – 8,600 at temperature range simulated in the aquifer). Conduction is the predominant heat transfer mechanism in the rest of the THS aquifer (Ra < 60). G eo lo gi a C ro at ic a 306 Geologia Croatica 77/3 The used modelling approach proved to be a valuable initial step for understanding regional hydrogeological systems despite its simplifications. The lack of extensive boreholes or geophysical data limits the understanding of the regional geological and hydrogeological settings. Complex modelling with such limited data can lead to numerous assumptions, increasing the uncertainty of the simulations. Therefore, a simple numerical approach with a detailed hydrogeological and thermal parameterisation or robust geological reconstruction is initially preferable to reduce uncertainty. In the absence of a local dataset of hydrogeological and thermal properties, the parameterisation of hydrostratigraphic units in THS was conducted considering relevant regional studies and literature collections. Incorporating well-documented information from the literature enhanced the credibility and robustness of the model, facilitating a more comprehensive understanding of the hydrothermal system's dynamics. As fluid flow and heat transport are three-dimensional processes, the use of a 2D profile in the conducted THS representation is likely the most important simplification. This approach does not account for the proper local and regional characterisation of faults and related subsidiary structures, which are critical for fluid flow in fractured rocks. In the wider area of the THS, a dome-like anticline structure affects the geological and structural architecture in the subsurface, as illustrated in the presented profiles (Fig. 4, profiles THS-1-3). Describing this complex geological setting with the 2D model is inherently problematic, as it fails to capture the total spatial variability and structural-geological intricacies. This limitation adds to the discrepancies observed between modelled and actual temperatures and fluid flows, emphasising the need for more sophisticated modelling approaches capable of accounting for the geological complexity of the THS. The construction of a 3D model, supported by a more extensive dataset of hydrogeological, structural-geological, and thermal parameters, could improve the understanding of the thermal system dynamics. Considering the 3D geological architecture will significantly enhance the accuracy and reliability of future numerical simulations. Future investigations will focus on conducting a comprehensive 3D geological reconstruction of the subsurface based on detailed structural-geological and geophysical investigations that are necessary to characterise and justify the proposed local and regional structural- geological setting of the THS. This consequently includes the reconstruction of regional and local fault meshes for a more sophisticated numerical model. 5. CONCLUSION The presented study utilised structural-geological field investigations and numerical modelling to develop a novel conceptual model of the THS. By constructing new geological profiles and integrating them with hydrogeological and geochemical data, we were able to provide new insights into the THS dynamics proposing a novel conceptual model. The recharge area of THS lies S of Topusko, where Triassic carbonates crop out. Regional and local scale fold and fault systems favour the deep circulation of the thermal water and the upwelling in Topusko forming a valuable geothermal resource. Numerical modelling corroborates the proposed con cep- tual model of the THS. A local temperature anomaly was simulated in the Topusko subsurface reproducing the temperature field observations. The occurrence of such a localised thermal anomaly in the hinge zone of the Topusko anticline suggests an enhanced fluid flow and heat transport pointing to the impact of terrestrial heat f low and local geological structures on the system. In particular, the synergetic effect of the anticline hinge, which was determined by structural investigations, and cogenetic fault/fracture systems, which were identified by surface geophysical research, results in an enhanced permeability field in the discharge area, thereby enabling the formation of abundant natural thermal springs. This study both advances the understanding of the THS and highlights the importance of integrating field observations with numerical modelling to physically validate and refine conceptual models. This integrated approach permits a more detailed reconstruction of the hydrogeological and thermal processes driving the development of geothermal resources. These aspects are crucial for developing site-specific exploitation plans of the resource, favouring its long-term utilisation. ACKNOWLEDGEMENT This research was funded by the Croatian Science Foundation (HRZZ), grant number UIP-2019-04-1218. REFERENCES ANDERSON, M.P. (2005): Heat as a Ground Water Tracer.– Groundwater, 43/6, 951–968. doi: 10.1111/j.1745-6584.2005.00052.x ANDERSON, M., WOESSNER, W. & HUNT, R. (2015): Applied Ground Water Modeling: Simulation of Flow and Advective Transport.– Else- vier, London, 533 p. AN, R., JIANG, X.-W., WANG, J.-Z., WAN, L., WANG, X.-S. & LI, H. (2015): A theoretical analysis of basin-scale groundwater temperature distribu- tion.– Hydrogeol. J., 23/2, 397–404. doi: 10.1007/s10040-014-1197-y BAĆ, J. & HERAK, M. (1962): Prijedlog određivanja užih i širih zaštitnih zona termomineralnih izvora u Hrvatskoj [Recommendation for deter- mination of wider and narrow protection zones for thermo-mineral springs in Croatia – in Croatian].– Unpublished report, Institute for ge- ological research, Zagreb, 147 p. BAHUN, S. & RALJEVIĆ, B. (1969): Mineralna, termalna i ljekovita vrela [Mineral and Thermal Springs – in Croatian].– Unpublished report, In- stitute for geological research, Zagreb. BÉKÉSI, E., LENKEY, L., LIMBERGER, J., PORKOLÁB, K., BALÁZS, A., BONTÉ, D., VRIJLANDT, M., HORVÁTH, F., CLOETINGH, S. & VAN WEES, J.D. (2018): Subsurface temperature model of the Hungari- an part of the Pannonian Basin.– Global and Planetary Change, 171, 48–64. doi: 10.1016/j.gloplacha.2017.09.020 BENČEK, Đ., BUKOVAC, J., MAGAŠ, N. & ŠIMUNIĆ, A. (2014): Osnovna geološka karta Republike Hrvatske 1:100.000, List Karlovac L33-92 [Ba- sic Geological Map of the Republic of Croatia 1:100.000, Karlovac sheet – in Croatian].– Croatian Geological Survey, Zagreb. BENČEK, Đ., BUKOVAC, J., MAGAŠ, N. & ŠIMUNIĆ, A. (2014): Osnovna geološka karta Republike Hrvatske 1:100.000, Tumač za list Karlovac L33-92 [Basic Geological Map of the Republic of Croatia 1:100.000, Ex- planatory notes for Karlovac sheet – in Croatian].– Croatian Geological Survey, Zagreb, 85 p. BOROVIĆ, S., MARKOVIĆ, T., LARVA, O., BRKIĆ, Ž. & MRAZ, V. (2016): Mineral and Thermal Waters in the Croatian Part of the Pannonian Ba- G eologia C roatica 307Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system sin.– In: PAPIĆ, P. (ed.): Mineral and Thermal Waters of Southeastern Europe. Springer, Cham. BOROVIĆ, S., URUMOVIĆ, K., TERZIĆ, J. & PAVIČIĆ, I. (2018): Examin- ing thermal conductivities of shallow subsurface materials for ground source heat pump utilization in the Pannonian part of Croatia.– The Mining-Geological-Petroleum Engineering Bulletin, 33/5, 27–35. doi: 10.17794/rgn.2018.5.3 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 BRÜCKL, E., BEHM, M., DECKER, K., GRAD, M., GUTERCH, A., KELLER, G.R. & THYBO, H. (2010): Crustal structure and active tec- tonics in the Eastern Alps.– Tectonics, 29/2. https://doi.org/10.1029/ 2009TC002491 BUNDSCHUH, J. & CÉSAR SUÁREZ, A.M. (2010): Introduction to the Nu- merical Modeling of Groundwater and Geothermal Systems.– CRC Press, London. 522 p. doi: 10.1201/b10499 BUKOVAC, J., ŠUŠNJAR, M., POLJAK, M. & ČAKALO, M. (1984): Os- novna geološka karta SFRJ 1:100.000, List Črnomelj L33–91 [Basic Ge- ological Map of the Republic of Croatia, 1:100.000, Črnomelj sheet – in Croatian].– Federal Geological Survey, Belgrade. CALCAGNO, P., BAUJARD, C., GUILLOU-FROTTIER, L., DAGALLIER, A. & GENTER, A. (2014): Estimation of the deep geothermal potential within the Tertiary Limagne basin (French Massif Central): An integrat- ed 3D geological and thermal approach.– Geothermics, 51, 496–508. doi: 10.1016/j.geothermics.2014.02.002 CERMAK, V., & RYBACH, L. (1982): Thermal properties: Thermal conduc- tivity and specific heat of minerals and rocks.– In: ANGENEISTER, G. (ed.): Landolt-Börnstein: Zahlenwerte and Funktionen aus Naturwissen- schaften und Technik, Neue Serie, Physikalische Eigenschaften der Gesteine. Springer Verlag, Berlin, Heidelberg and New York, V/1a, 305–343. CLOETINGH, S., CORNU, T., ZIEGLER, P.A., BEEKMAN, F. & ENTEC Working Group (2006): Neotectonics and intraplate continental topo- graphy of the northern Alpine Foreland.– Earth. Sci. Rev., 74, 127–196. doi: 10.1016/j.earscirev.2005.06.001 CRITTENDEN, J.C., TRUSSELL, R.R., HAND, D.W., HOWE, K.J. & TC- HOBANOGLOUS, G. (2012): MWH’s Water Treatment.– John Wiley & Sons, Inc. doi: 10.1002/9781118131473 CSONTOS, L. & VÖRÖS, A. (2004): Mesozoic plate tectonic reconstruction of the Carpathian region.– Palaeogeogr. Palaeoclimatol. Palaeoecol., 210/1, 1–56. doi: 10.1016/j.palaeo.2004.02.033 D’AGOSTINO, N., AVALLONE, A., CHELONI, D., D’ANASTASIO, E., MANTENUTO, S. & SELVAGGI, G. (2008): Active tectonics of the Adriatic region from GPS and earthquake slip vectors.– Journal of Geo physical Research: Solid Earth, 113/B12. https://doi.org/10.1029/ 2008JB005860 DHMZ. (2021): Monthly values and extremes.– Državni Hidrometeorološki Zavod [Croatian Meteorological and Hydrological Service – in Croa- tian]. (retrieved January 18, 2021) from https://meteo.hr/klima. php?section=klima_podaci¶m=k2_1 DIERSCH, H.J.G. (2014): FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media.– Springer Berlin, Heidelberg. doi: 10.1007/978‐3‐642‐38739‐5 DI NAPOLI, R., MARTORANA, R., ORSI, G., AIUPPA, A., CAMARDA, M., DE GREGORIO, S., GAGLIANO CANDELA, E., LUZIO, D., MESSINA, N., PECORAINO, G. et al. (2011): The Structure of a Hy- drothermal System from an Integrated Geochemical, Geophysical, and Geological Approach: The Ischia Island Case Study.– Geochem. Geo- phys. Geosyst., 12/7. https://doi.org/10.1029/2010GC003476 DOMENICO, P.A. & SCHWARTZ, F.W. (1997): Physical and Chemical Hy- drogeology (2nd edition).– Wiley, New York. 506 p. DOMENICO, P. A. & PALCIAUSKAS, V. V. (1973): Theoretical Analysis of Forced Convective Heat Transfer in Regional Ground-Water Flow.– GSA Bulletin, 84/12, 3803–3814. doi: 10.1130/0016-7606(1973)84<3803:TAO FCH>2.0.CO;2 EVANS, M.A. & FISCHER, M.P. (2012): On the distribution of fluids in folds: A review of controlling factors and processes.– J. Struct. Geol., 44, 2–24. doi: 10.1016/j.jsg.2012.08.003 FABBRI, P., POLA, M., PICCININI, L., ZAMPIERI, D., ROGHEL, A. & DALLA LIBERA, N. (2017): Monitoring, utilization and sustainable de- velopment of a low-temperature geothermal resource: A case study of the Euganean Geothermal Field (NE, Italy).– Geothermics, 70, 281–294. doi: 10.1016/j.geothermics.2017.07.002. 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 f low properties of fault zones.– J. Struct. Geol., 32/11, 1557–1575. doi: 10.1016/j.jsg.2010.06.009 FEDERAL GEOLOGICAL SURVEY. (1970): Geological Map of SFRY 1:500.000.– Federal Geological Survey, Beograd. FETTER, C.W. (2001): Applied Hydrogeology (4th edition).– Prentice Hall, Upper Saddle River, New Jersey. FLÓVENZ, Ó.G., HERSIR, G.P., SÆMUNDSSON, K., ÁRMANNSSON, H. & FRIDRIKSSON, T. (2012): Geothermal Energy Exploration Tech- niques.– In: SAYIGH, A. (ed.): Comprehensive Renewable Energy, 51–95. doi: 10.1016/B978-0-08-087872-0.00705-8 FUCHS, S. & BALLING, N. (2016): Improving the temperature predictions of subsurface thermal models by using high-quality input data. Part 1: Uncertainty analysis of the thermal-conductivity parameterization.– Geo thermics, 64, 42–54. doi: 10.1016/j.geothermics.2016.04.010 GARG, S.K., PRITCHETT, J.W., WANNAMAKER, P.E. & COMBS, J. (2007): Characterization of geothermal reservoirs with electrical sur- veys: Beowawe geothermal field.– Geothermics, 36/6, 487–517. doi: 10.1016/j.geothermics.2007.07.005 GHASEMIZADEH, R., HELLWEGER, F., BUTSCHER, C. et al. (2012): Review: Groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer system of Puerto Rico.– Hydrogeol J., 20, 1441–1461. doi: 10.1007/s10040-012- 0897-4 GOLDSCHEIDER, N., MÁDL-SZŐNYI, J., ERŐSS, A. & SCHILL, E. (2010): Review: Thermal water resources in carbonate rock aquifers.– Hydrogeol. J., 18/6, 1303–1318. doi: 10.1007/s10040-010-0611-3 GORJANOVIĆ-KRAMBERGER, D. (1905): Geologijske i hidrografijske pri- like oko Topuskoga s obzirom na topuske terme [Geological and hydro- graphic conditions in the area of Topusko with respect to thermal spring – in Croatian].– Works of the Yugoslav Academy of Sciences and Arts – Mathematics and Science Division, 37, 2–34. GORJANOVIĆ-KRAMBERGER, D. (1917): Geologijske i hidrografijske pri- like oko Topuskog s osobitim obzirom na Topuske terme [Geological and hydrographic conditions in the area of Topusko with particular interest in thermal springs – in Croatian].– Reports of discussions in Mathe matics and Science Division, Yugoslav Academy of Sciences and Arts, Zagreb. GRENERCZY, G., SELLA, G.F., STEIN, S. & KENYERES, A., (2005): Tec- tonic implications of the GPS velocity field in the northern Adriatic re- gion.– Geophys. Res. Lett., 32/16. doi: 10.1029/2005GL022947 HAVRIL, T., MOLSON, J.W. & MÁDL-SZŐNYI, J. (2016): Evolution of flu- id flow and heat distribution over geological time scales at the margin of unconfined and confined carbonate sequences - A numerical investiga- tion based on the Buda Thermal Karst analogue.– Mar. and Pet. Geol., 78, 738–749. doi: 10.1016/j.marpetgeo.2016.10.001 HEASLER, H.P., JAWOROWSKI, C. & FOLEY, D. (2009): Geothermal sys- tems and monitoring hydrothermal features.– In: YOUNG, R. & NOR- BY, L. (eds.): Geological Monitoring. Geological Society of America, Boulder. doi:/10.1130/2009.monitoring(05) HERAK, D., HERAK, M. & TOMLJENOVIĆ, B. (2009): Seismicity and earthquake focal mechanisms in North-Western Croatia.– Tectonophys- ics, 465/1–4, 212–220. doi: 10.1016/j.tecto.2008.12.005 https://doi.org/10.1016/j.marpetgeo.2016.10.001 https://doi.org/10.1130/2009.monitoring(05) G eo lo gi a C ro at ic a 308 Geologia Croatica 77/3 HERAK, M. & HERAK, D. (2023): Properties of the Petrinja (Croatia) earth- quake sequence of 2020–2021 – Results of seismological research for the first six months of activity.– Tectonophysics, 858, 229885. doi: 10.1016/j. tecto.2023.229885 HORVÁTH, F. & TARI, G. (1999): The IBS Pannonian basin project: A re- view of the main results and their bearings on hydrocarbon exploration.– Geol. Soc. London Spec. Publ., 156, 195–213. doi: 10.1144/GSL. SP.1999.156.01.11 HORVÁTH, F., BADA, G., SZAFIÁN, P., TARI, G., ÁDÁM, A. & CLOET- INGH, 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 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 HOUNSLOW, A.W. (1995): Water Quality Data: Analysis and Interpretation (1st ed.).– CRC Press. doi: 10.1201/9780203734117 HRVATOVIĆ, H. (2005): Geological Guidebook through Bosnia and Herze- govina.– Geological Survey Sarajevo-Geology Department,163 p. KAISER, B.O., CACACE, M., & SCHECK-WENDEROTH, M. (2013): 3D coupled fluid and heat transport simulations of the Northeast German Basin and their sensitivity to the spatial discretization: different sensiti- vities for different mechanisms of heat transport.– Environ. Earth Sci., 70/8, 3643–3659. doi: 10.1007/s12665-013-2249-7 KASTELIC, V. & CARAFA, M.M.C. (2012): Fault slip rates for the active External Dinarides thrust‐and‐fold belt.– Tectonics, 31/3. doi: 10.1029/2011TC003022 KHODAYAR, M. & BJÖRNSSON, S. (2024): Conventional Geothermal Sys- tems and Unconventional Geothermal Developments: An Overview.– Open J. Geol., 14/2, 196–246. doi: 10.4236/ojg.2024.142012 KOROLIJA, B., ŽIVALJEVIĆC, T. & ŠIMUNIĆ, A. (1980): Osnovna geološka karta SFRJ 1:100.000, List Slunj L 33–104 [Basic Geological Map of SFRY 1:100.000, Slunj sheet – in Croatian].– Institute for Geo- logical Research, Zagreb, Geological Survey, Sarajevo; Federal Geolog- ical Survey, Belgrade. KOROLIJA, B., ŽIVALJEVIĆ, T. & ŠIMUNIĆ, A. (1981): Osnovna geološka karta SFRJ 1:100.000, Tumač za list Slunj L33–104 [Basic Geological Map of SFRY 1:100.000, Explanatory notes for Slunj sheet – in Croa- tian].– Institute for geological research, Zagreb; Geological survey, Sa- rajevo; Federal geological survey: Belgrade. KOSOVIĆ, I., BRIŠKI, M., PAVIĆ, M., PADOVAN, B., PAVIČIĆ, I., MA- TOŠ, B., POLA M. & BOROVIĆ, S. (2023): Reconstruction of Fault Ar- chitecture 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 Panno- nian basin (Croatia).– Front. Earth Sci., 12. doi: 10.3389/feart.2024.1401935 LALOUI, L. & ROTTA LORIA, A. F. (2020): Heat and mass transfers in the context of energy geostructures.– In: Analysis and Design of Energy Geostructures. Academic Press., 69–135. doi: 10.1016/B978-0-12- 816223-1.00003-5 LEI, H. & ZHU, J. (2013): Numerical modeling of exploitation and reinjection of the Guantao geothermal reservoir in Tanggu District, Tianjin, China.– Geothermics, 48, 60–68. doi: 10.1016/j.geothermics.2013.03.008 LENKEY, L., DÖVÉNYI, P., HORVÁTH, F. & CLOETINGH, S.A.P.L. (2002): Geothermics of the Pannonian basin and its bearing on the neo- tectonics.– In: CLOETINGH, S.A.P.L, HORVÁTH, F., BADA, G. & LANKREIJER, A.C. (eds.): EGU Stephan Mueller Special Publication Series, 3, 29–40. doi: 10.5194/smsps-3-29-2002 MACENIĆ, M., KUREVIJA, T. & MEDVED, I. (2020): Novel geothermal gradient map of the Croatian part of the Pannonian basin system based on data interpretation from 154 deep exploration wells.– Renew. Sustain. Energy Rev., 132. doi: 10.1016/J.RSER.2020.110069 MÁDL-SZŐNYI, J. & TÓTH, Á. (2015): Basin-scale conceptual groundwater flow model for an unconfined and confined thick carbonate region.– Hy- drogeol. J., 23/7, 1359–1380. doi: 10.1007/s10040-015-1274-x MAJER, V. (1978): Stijene “dijabaz-spilit-keratofirske asocijacije” u području Abez-Lasinja u Pokuplju i Baniji (Hrvatska, Jugoslavija) [Rocks of dia- base-spilite-keratophyre association in the area of Abez-Lasinja in Pokuplje and Banije (Croatia, Yugoslavia) – in Croatian].– Acta Geol., 9/4, 42, 137–158. MAJER, V. (1993): Ofiolitni kompleks Banije s Pokupljem u Hrvatskoj i Pa- stirevo u Bosni [Ophiolite complexes of Banija and Pokuplje in Croatia and Pastirevo in Bosnia – in Croatian].– Acta Geol., 23/2, 39–84. MANDAL, A., BASANTARAY, A.K., CHANDROTH, A. & MISHRA, U. (2019): Integrated Geophysical Investigation to Map Shallow Surface Alteration/Fracture Zones of Atri and Tarabalo Hot Springs, Odisha, In- dia.– Geothermics, 77, 24–33. doi: 10.1016/j.geothermics.2018.08.007 MARINI, L. (2000): Geochemical techniques for the exploration and exploi- tation of geothermal energy.– Dipartamento per lo Studio del Territorio e delle sue Risorse, Universita degli Studi di Genova, Italy. MARTINSEN, G., BESSIERE, H., CABALLERO, Y., KOCH, J., COLLA- DOS-LARA, A. J., MANSOUR, M., SALLASMAA, O., PULIDO-VE- LAZQUEZ, D., WILLIAMS, N.H., ZAADNOORDIJK, W.J. & STI- SEN, S. (2022). Developing a pan-European high-resolution groundwater recharge map – Combining satellite data and national sur- vey data using machine learning.– Science of the Total Environment, 822 p. doi: 10.1016/j.scitotenv.2022.153464 MAZOR, E. (2004): Chemical and Isotopic Groundwater Hydrology, 3rd ed.– Marcel Dekker, New York, 13–179 p. MOECK, I.S. (2014): Catalog of geothermal play types based on geologic con- trols.– Renewable and Sustainable Energy Reviews, 37, 867–882. doi: 10.1016/j.rser.2014.05.032 MONTANARI, D., MINISSALE, A., DOVERI, M., GOLA, G., TRUMPY, E., SANTILANO, A. & MANZELLA, A. (2017): Geothermal resources within carbonate reservoirs in western Sicily (Italy): A review.– Earth- Science Reviews, 169, 180–201. doi: 10.1016/j.earscirev.2017.04.016 MROCZEK, E.K., MILICICH, S.D., BIXLEY, P.F., SEPULVEDA, F., BER- TRAND, E.A., SOENGKONO, S. et al. (2016): Ohaaki geothermal sys- tem: Refinement of a conceptual reservoir model.– Geothermics, 59, 311–324. doi: 10.1016/j.geothermics.2015.09.002 MUFFLER, P. & CATALDI, R. (1978): Methods for regional assessment of geothermal resources.– Geothermics, 7/2–4, 53–89. doi: 10.1016/0375- 6505(78)90002-0 NIELD, D.A. & BEJAN, A. (1999): Convection in Porous Media (2nd Edi- tion).– Springer, New York. doi:10.1007/978-1-4757-3033-3 OJHA, L., KARUNATILLAKE, S., KARIMI, S. & BUFFO, J. (2021): Ama- gmatic hydrothermal systems on Mars from radiogenic heat.– Nat. Com- mun., 12/1. doi: 10.1038/s41467-021-21762-8 PAVELIĆ, D. (2001): Tectonostratigraphic model for the North Croatian and North Bosnian sector of the Miocene Pannonian Basin System.– Basin Res., 13, 359–376. doi: 10.1046/j.0950-091x.2001.00155.x PAVELIĆ, D., KOVAČIĆ, M., MIKNIĆ, M., AVANIĆ, R., VRSALJKO, D., BAKRAČ, K., TIŠLJAR, J., GALOVIĆ, I. & BORTEK, Ž. (2003): The Evolution of the Miocene Environments in the Slavonian Mts. Area (Northern Croatia).– In: VLAHOVIĆ, I. & TIŠLJAR, J. (eds.): 22nd IAS Meeting of Sedimentology – Opatija 2003, Field Trip Guidebook. Cro- atian Geological Survey, Zagreb. POLA, M., CACACE, M., FABBRI, P., PICCININI, L., ZAMPIERI, D. & TORRESAN, F. (2020): Fault Control on a Thermal Anomaly: Con- ceptual and Numerical Modeling of a Low-Temperature Geothermal System in the Southern Alps Foreland Basin (NE Italy).– J. Geophys. Res. Solid Earth, 125/5, e2019JBO17394, doi: 10.1029/2019JB017394 PAVIĆ, M., KOSOVIĆ, I., POLA, M., URUMOVIĆ, K., BRIŠKI, M., & BO- ROVIĆ, S. (2023): Multidisciplinary Research of Thermal Springs Area in Topusko (Croatia).– Sustainability, 15/6, 5498. doi: 10.3390/ su15065498 PAVIĆ, M., BRIŠKI, M., POLA, M., & BOROVIĆ, S. (2024): Hydrogeoche- mical and environmental isotope study of Topusko thermal waters, Cro- https://doi.org/10.4236/ojg.2024.142012 https://doi.org/10.1038/s41467-021-21762-8 G eologia C roatica 309Pavić et al.: A conceptual and numerical model of fluid flow and heat transport in the Topusko hydrothermal system atia.– Environ. Geochem. Health, 46/4, 133. doi: 10.1007/s10653-024- 01904-9 POURASKARPARAST, Z., AGHAEI, H., COLOMBERA, L., MASOERO, E. & GHAEDI, M. (2024): Fracture aperture: A review on fundamental concepts, estimation methods, applications, and research gaps.– Mar. Pet. Geol., 164, 106818. doi: 10.1016/j.marpetgeo.2024.106818 PRELOGOVIĆ, E., SAFTIĆ, B., KUK, V., VELIĆ, J., DRAGAŠ, M. & LU- ČIĆ, D. (1998): Tectonic activity in the Croatian part of the Pannonian basin.– Tectonophysics, 297/1–4, 283–293. doi: 10.1016/S0040- 1951(98)00173-5 RMAN, N. & TÓTH, G. (2011): Hydrogeological conceptual model.– Geolo- gical Survey of Slovenia, Ljubljana; Geological Institute of Hungary, Budapest [in Hungarian], 25 p. RMAN, N. (2014): Analysis of long-term thermal water abstraction and its impact on low-temperature intergranular geothermal aquifers in the Mura-Zala basin, NE Slovenia.– Geothermics, 51, 214–227. doi: 10.1016/j.geothermics.2014.01.011 RMAN, N., BĂLAN, L.L., BOBOVEČKI, I., GÁL, N., JOLOVIĆ, B., LA- PANJE, A. et al. (2020): Geothermal sources and utilization practice in six countries along the southern part of the Pannonian basin.– Environ. Earth Sci., 79, 1–12. doi: 10.1007/s12665-019-8746-6 ROYDEN, L.H., & HORVÁTH, F. (1988): The Pannonian Basin: A Study in Basin Evolution.– American Association of Petroleum Geologists. doi: 10.1306/M45474 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-We- stern Hungary). Geologia Croatica, 56/1, 101–122. doi:10.4154/232 SCANLON, B., & DUTTON, A. (2000): Groundwater Recharge in Texas.– The University of Texas at Austin, and Marios Sophocleous, Kansas Ge- ological Survey, Lawrence, KS. SCANLON, B.R., MACE, R.E., BARRETT, M.E. & SMITH, B. (2003): Can we simulate regional groundwater flow in a karst system using equiva- lent porous media models? Case study, Barton Springs Edwards aquifer, USA.– Journal of Hydrology, 2761–4, 137–158. doi: 10.1016/S0022- 1694(03)00064-7 SCHMID, S.M., BERNOULLI, D., FÜGENSCHUH, B., MATENCO, L., SC- HEFER, S., SCHUSTER, R., TISCHLER, M. & USTASZEWSKI, K. (2008): The Alpine-Carpathian-Dinaridic orogenic system: Correlation and evolution of tectonic units.– Swiss J. Geosci., 101, 139–183. doi: 10.1007/s00015-008-1247-3 SCHMID, S. M., FÜGENSCHUH, B., KOUNOV, A., MAŢENCO, L., NIE- VERGELT, P., OBERHÄNSLI, R., PLEUGER, J., SCHEFER, S., SC- HUSTER, R., TOMLJENOVIĆ, B., USTASZEWSKI, K. & VAN HIN- SBERGEN, 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 STEVANOVIĆ, Z., DULIĆ, I. & DUNČIĆ, M. (2015): Some experiences in tapping deep thermal waters of the Triassic karstic aquifer in the Panno- nian Basin of Serbia.– Central European Geology, 581–2, 50–61. doi: 10.1556/24.58.2015.1-2.3 STYLIANOU, I., TASSOU, S., CHRISTODOULIDES, P., PANAYIDES, I. & FLORIDES, G. (2016): Measurement and analysis of thermal proper- ties of rocks for the compilation of geothermal maps of Cyprus.– Rene- wable Energy, 88, 418–429. doi:10.1016/j.renene.2015.10.058" SZANYI, J. & KOVÁCS, B. (2010): Utilization of geothermal systems in So- uth-East Hungary.– Geothermics, 39/4, 357–364. doi: 10.1016/j.geother- mics.2010.09.004 ŠEGOTIĆ, B. & ŠMIT, I. (2007): Studija optimirane energetske učinkovitosti korištenja geotermalnih voda [Study of Optimized Energy Efficiency of Geothermal Water Use – in Croatian].– Unpublished report, Termoinže- njering-projektiranje, Zagreb. ŠIKIĆ, K. (1990): Osnovna geološka karta Republike Hrvatske 1:100.000, list Bosanski Novi L 33–105 [Basic Geological Map of the Republic of Cro- atia 1:100.000, Bosanski Novi sheet – in Croatian].– Croatian Geological Survey, 2014. ŠIKIĆ, K. (1990): Osnovna geološka karta Republike Hrvatske 1:100.000, Tu- mač za list Bosanski Novi L 33-70 [Basic Geological Map of the Repu- blic of Croatia 1:100.000, Explanatory notes for Bosanski Novi sheet – in Croatian].– Croatian Geological Survey Zagreb, 2014. ŠIMUNIĆ, A. (2008): Topusko.– In: ŠIMUNIĆ, A. & HEĆIMOVIĆ, I. (eds.): Geotermalne i mineralne vode Republike Hrvatske [Mineral and Ther- mal Waters of the Republic of Croatia – in Croatian]. Croatian Geologi- cal Survey, 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.– Geol. Soc. London Spec. Publ., 156, 215–250. doi: 10.1144/GSL. SP.1999.156.01.12 TEUTSCH, G. & SAUTER, M. (1991): Groundwater modeling in karst terra- nes: Scale effects, data acquisition and field validation.– In: Third conference on hydrogeology, ecology, monitoring, and management of ground water in Karst Terranes.– National Ground Water Association, Dublin, Ohio,17–35. TOMLJENOVIĆ, B. & CSONTOS, L. (2001): Neogene–Quaternary structu- res in the border zone between Alps, Dinarides and Pannonian Basin (Hrvatsko zagorje and Karlovac Basins, Croatia).– Int. J. Earth Sci., 90/3, 560–578. doi: 10.1007/s005310000176 TOMLJENOVIĆ, B. (2002): Strukturne značajke Medvednice i Samoborskog gorja [Structural features of Medvednica and Samobor hills – in Croa- tian].– Unpubl. PhD Thesis, Faculty of Mining, Geology and Petroleum Engineering, University of Zagreb, Zagreb. TOMLJENOVIĆ, B., CSONTOS, L., MÁRTON, E. & MÁRTON, P. (2008): Tectonic evolution of the northwestern Internal Dinarides as constrained by structures and rotation of Medvednica Mountains, North Croatia.– Geol. Soc. London Spec. Publ., 298/1, 145–167. doi: 10.1144/SP298.8 TORRESAN, F., PICCININI, L., CACACE, M., POLA, M., ZAMPIERI, D. & FABBRI, P. (2022): Numerical modeling as a tool for evaluating the renewability of geothermal resources: the case study of the Euganean Geothermal System (NE Italy).– Environ. Geochem. Health, 44/7, 2135–2162. doi: 10.1007/s10653-021-01028-4 TÓTH, J. (2009): Gravitational systems of groundwater flow: theory, evalua- tion, utilization.– Cambridge University Press, Cambridge. doi: 10.1017/ CBO9780511576546 TURCOTTE, D.L. & SCHUBERT, G. (1982): Geodynamics: Applications of continuum mechanics to geological problems.– Wiley, New York, 464 p. 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 J. Geosci., 101/Suppl 1, 273-294. 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. doi: 10.1029/2010TC002668 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.– J. Geodyn., 82, 52–68. doi: 10.1016/j.jog.2014.04.006 VASS, I., TÓTH, T.M., SZANYI, J. & KOVÁCS, B. (2018): Hybrid numerical modelling of fluid and heat transport between the overpressured and gravitational flow systems of the Pannonian Basin.– Geothermics, 72, 268–276. doi: 10.1016/j.geothermics.2017.11.013 VELIĆ, I. & SOKAČ, B. (1982): Novi nalazl naslaga donjeg i srednjeg trijasa u zapadnom Kordunu (središnja Hrvatska) [New discoveries of the Lower and Middle Triassic in the western part of Kordun area (central Croatia) – in Croatian].– Geološki vjesnik, 35, 47–57. VLAHOVIĆ, I., TIŠLJAR, J., VELIĆ, I. & MATIČEC, D. (2005): Evolution of the Adriatic Carbonate Platform: Palaeogeography, main events and depositional dynamics.– Palaeogeogr. Palaeoclimatol. Palaeoecol., 220/3–4, 333–360. doi: 10.1016/j.palaeo.2005.01.011 https://doi.org/10.1556/24.58.2015.1-2.3 https://doi.org/10.1556/24.58.2015.1-2.3 https://doi.org/10.1016/j.geothermics.2017.11.013 G eo lo gi a C ro at ic a 310 Geologia Croatica 77/3 WORTHINGTON, S.R.H., FOLEY, A.E., & SOLEY, R.W.N. (2019): Transi- ent characteristics of effective porosity and specific yield in bedrock aqu- ifers.– J. Hydrol., 578, 124129. doi: 10.1016/j.jhydrol.2019.124129 XIONG, J., LIN, H., DING, H., PEI, H., RONG, C. & LIAO, W. (2020): Inve- stigation on thermal property parameters characteristics of rocks and its influence factors.– Natural Gas Industry B, 7/3, 298–308. doi: 10.1016/j. ngib.2020.04.001 URL 1: https://www.avenza.com/avenza-maps/. Accessed on 29 April, 2024. URL 2: https://www.esri.com/news/arcnews/spring12articles/introducing-arc- gis-101.html. Accessed on 29 April, 2024. https://www.esri.com/news/arcnews/spring12articles/introducing-arcgis-101.html https://www.esri.com/news/arcnews/spring12articles/introducing-arcgis-101.html