Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression

The Bjelovar sub-depression is situated in the north-western part of Croatia, on the southwestern margin of the Pannonian Basin System where favourable geothermal conditions exist. Thermal waters are used for recreation, balneotherapy, space heating, directly as sanitary water and elec- tricity production. Geophysical, geological and borehole data were used to determine the types of geothermal reservoirs. In addition, several campaigns were conducted to sample geothermal wa- ters from the Daruvar spa, Velika-1 and Krečaves locations for isotope (δ 18 Ο and δ 2 H) and physi-co-chemical (EC, T, pH, DO, Na + , K + , Mg 2+ , NH 4+ , Ca 2+ , SO 42- , Cl - , Br - , F - , SiO 2 and H 2 S) analyses to determine their hydrochemical characteristics. Two major types of geothermal reservoirs were determined: (i) ‘basement – BM’ reservoir, (ii) ‘basin fill – BF’ reservoir. The BM reservoir consists of Palaeozoic and Mesozoic sediments composed of: (i) fractured/karstified carbonate sediments and/or (ii) fractured/fissured crystalline/metamorphic rocks. The BF reservoirs are ‘Lower Panno nian’ and ‘Upper Pannonian’ sediments composed of coarse and fine-grained sand, sandstones and marls. The BM geothermal aquifers are the most important ones in the study area. The stable isotope δ 2 H and δ 18 O indicate that the monitored thermal waters have a meteoric origin, but without recent replenishment. Monitored waters belong to mixed hydrochemical types, from Na-ClHCO 3 to Na-HCO 3 types in the deep basin thermal waters and a CaMg-HCO 3 type in the thermal waters from shallower parts. The study area has great geothermal potential. The estimated total available thermal power from Križevci, Velika-1 and the Daruvar spa is 70.5 MWt, but only 28 % of this ther mal power is used. Since the predominant activity in the study area is agriculture, the geothermal resources available could lead to modern agricultural development and consequently contribute to increasing the standard of living of the local population. However, additional geophysical, geo- logical, hydrogeological and hydrochemical investigations at a number of new potential locations are required to estimate the total available geothermal resources.


INTRODUCTION
The North-western part of Croatia represents the southwestern margin of the Pannonian Basin System (PBS), and it is characterized by a high average geothermal gradient (49°C/km) and surface heat flow (76 mW/m 2 ) (BOŠNJAK et al., 1998). Favourable geothermal conditions in this area are a consequence of the very shallow Moho (Mohorovičić seismic discontinuity), at between 25-28 km depth (ALJINOVIĆ & BLAŠKOVIĆ, 1984; dicts a 39% Renewable Energy Source (RES) share in electricity production and 19.6% in heating and cooling from geothermal sources. In addition, geothermal waters in Croatia contribute to economic growth through health and tourist resorts. Two different laws regulate the usage of thermal waters, based on the method of utilization. If water is used for bottling, swimming and balneology, then the Water Act is applicable (OFFICIAL 63/11,130/11,56/13,14/14). Alternatively, if the water is used for heating and electricity production, the Exploration and Exploitation of Hydrocarbons Act has been applicable since 2018 (OFFICIAL GAZETTE 52/18,52/19). Previously, the Mining Act was applicable (OFFICIAL GAZETTE 56/13,14/14). The sustainable use and management of geothermal resources is essential to ensure their long-term availability.
The Bjelovar sub-depression is located in the centre of northwestern Croatia. It is a part of the Drava depression and extends towards the Sava depression, bordered by mountains on all sides. Geothermal springs occur at the east -southeast edge of the subdepression. Due to several hydrocarbon exploration campaigns in the 1970's, 80's and 90's, numerous geothermal reservoirs were discovered.
The aim of this paper is to hydrochemically characterize the geothermal aquifers that are utilized in the Bjelovar sub-depression and to give an overview on potential areas suitable for their direct or indirect use.

Study area
The Bjelovar sub-depression is bordered by the Kalnik Mt. to the northwest, Bilogora Mt. to the northeast, the Papuk, Ravna gora and Psunj Mts. to the east and southeast, and Moslavačka Mt. to the south (Fig. 1). The local climate is categorized as part of the Cfb group of the Köppen-Geiger classification system with an average annual air temperature of 10.9°C and average annual precipitation around 800 mm (ZANINOVIĆ et al., 2008). The study area comprises territories of three counties: the whole territory of Bjelovarsko-bilogorska, around 50% of Koprivničkokriževačka, and 10% of the Zagrebačka County. Moreover, these counties support tourism, agriculture, small and medium enterprises, and one of the many important factors which ensures future economic growth is a sustainable source of energy for consumption (heat and electricity).

Data collection
Intensive geological data collection of the Bjelovar sub-depression was carried out as part of the Danube Region Leading Geothermal Energy project (DARLINGe, 2017(DARLINGe, -2019. In addition, chemical analyses of geothermal waters undertaken after the drilling of boreholes were used to improve the interpretation of the hydrochemical characteristics of geothermal waters. These data were provided by the Ministry of Economy and Energy (MINGO) in 2017.
The geothermal waters from the Daruvar, Krečaves and Ciglena locations were sampled between March 2018 to December 2019 in the scope of the DARLINGe and MAMEBio (A multiphasic approach for deciphering the microbial ecology and biotechnological potential of hot springs in Croatia) projects for chemical and isotopic analyses. Four sampling campaigns were carried out in Daruvar and Krečaves and two in Ciglena. In Križevci, it was technically impossible to take samples. Therefore, hydrochemical interpretation was based on available chemical analyses done during the pumping test in 1990. In addition, one chemical analysis which was performed during DTS tests after the drilling was used for Ciglena. Prior to sampling, the electrical conductivity (EC), pH, temperature (T) and dissolved oxygen content (O 2 ) were determined using portable WTW electrodes. The total alkalinity was measured by titration with 1.6 N H 2 SO 4 , using phenolphthalein and bromocresol green-methyl red as indicators. Dissolved cations (Na + , K + , Mg 2+ , NH 4 + , Ca 2+ ) and anions (SO 4 2-, Cl -, Br -, F -) were measured using Ion Chromatography on Dionex ICS-6000 DC, at the Hydrochemical laboratory of the Croatian Geological Survey. Silica concentrations (SiO 2 ) and H 2 S were measured spectrophotometrically at the same laboratory using HACH DR3900. The water quality analyses was checked for a cation-anion balance by measuring the relative deviation from charge balance (Δ meq = 100 x (Σ meq+ − Σ meq− )/(Σ meq+ + Σ meq− ) < ±5%) (MANDEL & SHIFTAN, 1980;DOMENICO & SCHWARTZ, 1990). The δ 18 Ο and δ 2 H were determined using Picarro L2130i (Santa Clara, USA) in the Hydrochemical Laboratory of the Croatian Geological Survey. The instrument uses Cavity Ring-Down Spectroscopy (CRDS) technology (http:// www.picarro.com/technology/cavity_ring_down_spectroscopy). All measurements were checked with Picarro's standards (Depleted −29.6 ± 0.2 δ 18 Ο; −235 ± 1.8 δ 2 H; Mid −20.6 ± 0.2 δ 18 Ο; −159 ± 1.3 δ 2 H; Zero 0.3 ± 0.2 δ 18 Ο; 1.8 ± 0.9 δ 2 H), which were checked periodically against the International Atomic Energy Agency (IAEA) standards: Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation 2 (SLAP2). Measurement precision was ± 0.3 ‰ for δ 18 Ο and ± 1 ‰ for δ 2 H.

GEOLOGICAL AND HYDROGEOLOGICAL SETTINGS OF GEOTHERMAL AQUIFERS
Two major types of geothermal reservoirs were determined using the methodology developed in the project DARLINGe (RO-TÁR-SZALKAI et al., 2018): (i) 'basement -BM' reservoirs, (ii) 'basin fill -BF' reservoirs.

The basement reservoirs
The BM reservoirs are Palaeozoic and Mesozoic sediments the formation of which was associated with the Alpine-Carpathian orogen, and represent the base of Tertiary sediments. Aquifers formed in this reservoir are the most important type in the study area. Two types of geological strata represent reservoir rocks: (i) fractured/karstified carbonate sediments (dolomite; dolomite breccia -mainly Triassic in age and limestones of Triassic/Cretaceous age) and/or (ii) fractured/fissured crystalline/metamorphic rocks (gabbro, schist, etc.). Triassic dolomite is the most favourable type of geothermal reservoir in Croatia due to its high permeability at some locations. Miocene Badenian sandstone and carbonate overlying Triassic deposits are often associated with these aquifers. These aquifers are mainly low enthalpy systems (temperatures < 150°C), except for the eastern part of the Bjelovar sub-depression (dolomite aquifers) which are high enthalpy (temperatures > 150°C). Water temperatures range from 40 to 170°C. At the northwest, west and southeast margins of Bjelovar subdepression, the BM reservoirs occur as outcrops. Towards the sub-depression interior, the depth to the top of the BM varies from 100 to over 2000 m bsl. The depth mainly depends on the reservoir position in the sub-depression, because they are arranged into knappes along thrust sheets, dissected by strike-slip, normal and transcurrent faults, associated with the multiphase tectonic development of the basin (PRELOGOVIĆ et al., 1998;SAFTIĆ et al., 2003;MALVIĆ & VELIĆ, 2010;CVETKOVIĆ et al., 2019) (Figs. 2-5). At present, the geothermal aquifer is utilized at two locations, and exploration is ongoing at a third location.
The first location is situated in the middle of the study area, between the villages of Patkovac and Ciglena near the town of Bjelovar, where the INA oil company drilled four boreholes during the 1990's (Ptk-1, VC-2, VC-1, VC-1A) (Fig. 1). All boreholes reached the carbonate rocks (dolomite, dolomite breccia) of Mesozoic age, where thermal water with a temperature of 175°C was observed at the bottom of VC-1 after the drilling. Depth to the top of carbonate sediments (dolomite, dolomite breccia) varies from 2390 to 2314 m bsl, depending on borehole positions. The thickness of the carbonate sediments ranges from 500 to 1300 m. The dolomite and dolomite breccia are highly porous (~ 36%). During drill stem tests (DST), average values for vertical and horizontal permeability between 13.3 mD and 350 mD were determined (INA, 1990). Due to the high permeability, the total yield of boreholes VC-1 and VC-1A is 114 L/s. The pressure in the aquifer is more than 300 bar.
The second location is situated at the east -southeast edge of the sub-depression in Daruvar town, where outcrops of dolomite exist. The upwelling of thermal waters occurs via a N-S strike transcurrent fault in the Toplica River valley (Figs. 2 and Figure 2. Schematic map of the depth to the top of the BM reservoir, based on gravity data. 5). The largest natural springs are the Ivanovo vrelo, Antunovo vrelo, and Marijina vrela with temperatures ranging from 40 to 48°C. Since the springs are situated in the river valley, alluvial deposits overlie the reservoir. Such hydrogeological conditions enable the mixing of 'cold' alluvial water with accumulated geothermal water in the underlying aquifer, which is further enhanced during periods of increased geothermal water abstraction at wells and springs. Since the spring zone is very dispersive, it is very difficult to determine its total discharge. Measurements performed in 1962 by BAĆ & HERAK determined that the total spring's discharge is 7.5 L/s, then 10.5 L/s in 1983 by MRAZ, and the latest measurements in 2008 by LARVA et al. are 12 L/s. Besides the natural springs, two boreholes exist, Dar-1 drilled in 2009, and D-1 drilled a few years later for ensuring additional capacity. Drilling reports of both boreholes are indicate that the thickness of the thermal aquifers (dolomite, limestone and dolomite breccia) is over 200 m. According to pumping tests and geophysical well-logging measurements (electrical logging, neutron natural radioactivity), which were conducted in the borehole Dar-1 in 2009, the porosity in the area where aquifer has low permeability (it is not highly tectonised) is 10%, but in the very permeable area is 40% (CROSCO, 2009). The horizontal hydraulic permeability varies between 7.58 x 10 -5 and 1.91 x 10 -5 m/s (CROSCO, 2009). The maximum pumping yield is 40 L/s, but at this yield, lowering of the water level in the aquifer affects the discharge of the springs. The optimum pumping yield of the borehole Dar-1 is 12.2 L/s and it is established as a steady state flow without endangerment of the spring's discharges (LARVA et al., 2009). The optimum pumping yield of borehole D-1 is 2 L/s (LARVA et al., 2008).
The third location is in Križevci town, where borehole Kža-1 was drilled in 1986. The borehole ended in Palaeozoic rocks (slate, sandstone, slaty breccia-conglomerates) at a depth of 1496 m. In 1990, a 2-month pumping test was conducted and the water temperature during the test was 70°C, with an optimal pumping yield of 2 L/s (ŠARIN, 1990). In 2008, a second phase of testing was conducted and the water temperature was 83.9°C (MALJKOVIĆ & GUDMUNDSSON, 2017). Unfortunately, the original testing report was unavailable, so the duration of pumping and pumping yield are not known. Nevertheless, the Hydrocarbon Agency opened a tender procedure in September 2019 in order to select the most favourable bidder for exploration of the geothermal waters for energy purposes in the exploration area of "Križevci" and in the beginning of 2020 gave concession for exploration of the geothermal waters for energy purposes. Hopefully, this new exploration will provide a more conclusive answer about the water temperature.
There are also several locations (Pavljani, Letičani, Grubišno Polje, Velika Trnovitica and Dežanovac) where numerous geothermal aquifers were detected in BM crystalline/carbonate rocks, with water temperatures in the range of 68 to 170°C (Fig. 1).

The basin fill reservoirs
The BF reservoirs are a consequence of geological changes from the Mid Miocene when rifting occurred in the PBS (ROYDEN, 1988;HORVÁTH 1993;PAVELIĆ, 2001;MALVIĆ & VELIĆ, 2010). A deep basin was formed and offshore pelagic sediments developed in the basin interior, while coarse-grained brackish littoral sediments were deposited on the margins. As a result, a great lake was formed (Lake Pannon) which was gradually infilled with sediments during the Late Miocene and Early Pliocene (ROY-DEN, 1988;HORVÁTH 1993;PAVELIĆ, 2001;SAFTIĆ et al. 2003;MALVIĆ & VELIĆ, 2010;SEBE et al. 2020). These sediments are today's geothermal aquifers, known as the 'Pannonian' sedimentary succession -'Lower Pannonian' (LP) and 'Upper Pannonian' (UP) sediments (RMAN et al., 2016a;RMAN et al., 2016b;ROTÁR-SZALKAI et al., 2017;NÁDOR et al., 2020). The LP sediments consist of coarse material on the basin floor deposited by turbidity currents with subsequent sedimentation of silt and argillaceous marl on the slopes. As infilling of the lake progressed, Figure 3. Schematic cross section A-A' of the Križevci basin (according to Fig. 2). a shelf was formed and UP sediments consisting of fine-grained materials were deposited, overlying the LP sediments. Finally, alluvial plain sedimentation completed the infill of the lake basin (PAVELIĆ, 2001;MALVIĆ & VELIĆ, 2010;MAGYAR et al., 2013) (Figs. 2-5). Unfortunately, this type of aquifer in the study area is not utilized and not explored so much as others. However, in the village of Krečaves there is one old borehole that was drilled in the 1960's during oil and gas exploration that detected thermal water. The favorable geothermal aquifer was reached at depth 638 m bsl. The borehole ended in the base of Tertiary sediments, which consists of Palaeozoic sediments, represented by slate and marble limestone. The measured temperature at the bottom after drilling was 56°C (INA, 1966). Unfortunately, a DST test showed that there is no water in these sediments. In contrast, the BF sequence was productive and a DST test was conducted. The thickness of the geothermal aquifer in the borehole is 150 m. During the test, the outflowing yield was 1.5 L/s, the water temperature was 52°C, porosity ranged from 28 to 34%, and the average horizontal permeability was 380 mD (INA, 1966). Today's outflow water temperature is around 45°C, and is dependent on water outflow duration, which in this case is the time of consumption by local people for showering. Sometimes during the sampling campaign, the water temperature reached 47°C.
The second unused borehole is in the village of Veliko Korenovo near the town of Bjelovar, where the water temperature is 30°C (MALJKOVIĆ & GUDMUNDSSON, 2017).

HYDROCHEMICAL CHARACTERISTICS OF THE GEOTHERMAL WATERS
The mean values of the measured hydrochemical parameters in the thermal waters of monitored springs and boreholes are given in Table 1. The highest EC value is the thermal water from Velika-1, Ciglena (23816 µS/cm), followed by water from borehole Kža-1 (2500 µS/cm). EC values for thermal water from the Daruvar Spa on all three objects are very similar, and Krečaves has slightly higher values than the Spa waters. Thermal water of Velika-1, Ciglena contains a high amount of dissolved CO 2 (30 m 3 / m 3 ) comprising 99.5% of gas, about 59 ppm H 2 S and 24.4 ppm NH 4 + . Also, the water contains oil (INA, 1990). Krečaves contains 2.4 ppm H 2 S and 8.5 ppm NH 4 + . Thermal waters of the Daruvar Spa have much lower concentrations of H 2 S and NH 4 + . The highest concentrations of Na + , Cl -, K + , and SO 4 2are present in thermal water from Ciglena, followed by water from Križevci, Krečaves and Daruvar Spa. Fis observed in all the sampled waters, while Bris only observed in the waters of Ciglena.
Hydrochemical types, as recognized from water chemistry data, are a consequence of the chemistry of the recharging water and water-aquifer matrix interactions (e.g. cation exchange), as well as the groundwater residence time within the aquifer. Two general processes contribute solutes to thermal waters: carbonate dissolution and silicate weathering (DOMENICO & SCHWARTZ, 1990). From Fig. 7 it can be observed that the water chemistry of Ciglena and Daruvar Spa are primarily controlled by the dissolution of dolomite, because the Mg 2+ /Ca 2+ ratios are above 0.33. However, according to the chemical composition of thermal waters from Velika-1, Ciglena, the prevailing cation is Na + , indicating that other geochemical processes are present in the thermal aquifer or the origin of the water is different. On the bivariate mixing diagram of Na + -normalized Ca 2+ versus Na + -normalized Mg 2+ (Fig. 8a), the sample falls between the silicate weathering and carbonate dissolution domains, but slightly closer to the latter, indicating that carbonate dissolution is dominant. The influence of cation exchange was evaluated by an equivalent bivariate plot of corrected bivalent cations versus corrected Na + (Fig. 8b). Concentrations of bivalent cations (Ca 2+ and Mg 2+ ) that may have been involved in exchange reactions were corrected by subtracting equivalent concentrations of the associated anions (HCO 3and  SO 4 2-) that would derive from other processes (e.g. carbonate dissolution or silicate weathering). Similarly, Na + that may derive from the aquifer matrix can be accounted for by assuming that Na + contributions of meteoric origin would be balanced by equivalent concentrations of Cl - (MCLEAN & JANKOWSKI, 2000). For active cation exchange taking place in the aquifer, the slope of this bivariate plot should be -1 (i.e. y = -x) (JANKOWSKI et al., 1998). Also, on the bivariate plot of (Ca 2+ + Mg 2+ ) -(HCO 3 -+ SO 4 2-) against (Na + -Cl -) (Fig. 8b) samples of both locations are on the line y = -x, which indicates active cation exchange processes within the aquifer. These geochemical processes have been detected before in Croatian carbonate thermal aquifers (MARKOVIĆ et al., 2015). However, Fig. 8a shows that the chemical composition of thermal waters from Krečaves and Križevci are controlled by silicate weathering. Although in the Križevci area, the aquifer is a ''silicate'' type of aquifer as in Krečaves, the chemical composition is different because the aquifer matrix contains a high amount of sulfide and chlorite minerals. In the past, mining of lead, zinc and silver ores occurred on Kalnik Mt., which is the recharge area of the thermal aquifer. The Na/Cl ratio of the thermal waters from Velika-1, Ciglena is <1 and for others is >10, which indicates a marine origin for Ciglena and a non-marine origin for the rest of the waters. The Br/Cl ratio of the thermal waters from Velika-1, Ciglena is higher than that for sea water. This water also contains traces of oil.
In all the observed thermal waters, measured stable isotope ratios of δ 2 H and δ 18 O are distributed around the local meteoric water line (LMWL) Bjelovar (HUNJAK et al., 2013) indicating meteoric origin (Fig. 9). The stable isotope ratios δ 2 H and δ 18 O in thermal waters from Velika-1, Ciglena are slightly shifted to the right from the LMWL at Bjelovar because the values are enriched. The enrichment in δ 18 O is a consequence of geothermal exchange between the fluid and host rock at high pressures and temperatures, a phenomenon that has been observed elsewhere in the high temperature deep groundwater systems in the world (CARTWRIGHT, 2002(CARTWRIGHT, , 2010ARNORSSON 2000;RMAN, 2016). Measured ratios in the thermal waters of Daruvar Spa area lie on the LMWL of Bjelovar (MARKOVIĆ & LARVA, 2012), while values from thermal water Krečaves are slightly above the Bjelovar LMWL indicating depletion. Depletion of stable isotopes was characteristic for colder periods in the Pleistocene and such depleted values were observed in the thermal waters of the Pannonian basin in Slovenia and Hungary (SZŐCS et al., 2013;RMAN, 2016;RMAN et al., 2016a). Thermal water from Ciglena does not contain 14 C, it is inactive, indicating that the water is very old. In contrast, thermal waters from Daruvar Spa still contain 14 C and the thermal water age was calculated using 13 C correction which was 15.000 years (LARVA et al., 2008). In addition, BOROVIĆ et al. (2019) in the numerical model of the Daruvar Spa thermal aquifer established that the recharge area is from the east Papuk Mt.

UTILIZATION AND POTENTIAL OF GEOTHERMAL AQUIFERS
Thermal waters in Daruvar town have been utilized since ancient Roman times. In 1772, count Janković started building around springs, hoping that the town would become a place of healing, leisure and a recreation centre which was rapidly realised. Nowadays, thermal waters from Ivanovo vrelo, Antunovo vrelo and borehole D-1 are used only for recreation, balneotherapy and sanitary water in the complex of the Special Hospital for Medical Rehabilitation and Thermal water park Aquae Balissae. Considering the amount of thermal water which can be sustainably used from two boreholes and three major springs, the available thermal power is 4.6 MWt. At the moment, only 13 % of this available thermal power is utilized ( Table 2).
The geothermal power plant "Velika 1" was recently built as a 35-million-euro investment by Turkish MB Holding, starting full-time operation in March 2019 after initial operational tests. The power plant uses old INA wells as producers and injectors, and utilizes the heat of geothermal water to generate electricity by the Organic Rankine Cycle (ORC) process, using isobutane as the working fluid rotating the turbine. Geothermal water is produced by two production wells (VC-1 and VC-1A) at the wellhead temperature of 166°C and wellhead pressure of 30 bar. By transferring its heat to the working fluid, the water is cooled down to 75°C and is then reinjected into the same geothermal reservoir through two remote injection wells (VC-2 and Ptk-1). At maximum gross power of 16.5 MWe, this is currently the largest ORC geothermal power plant in continental Europe. With constant net delivery of 10 MWe into the medium voltage 35 kV grid, it generates approx. 7 GWh of electrical energy per month, which is the equivalent of the average electricity consumption of 29.000 Croatian households. The available thermal power for this location is 64.9 MWt.
In September 2019, the Hydrocarbon Agency opened a tender procedure in order to select the most favourable bidder for exploration of geothermal waters for energy purposes in the exploration area "Križevci". According to EIHP study (MALJKO VIĆ & GUDMUNDSSON, 2017), the available thermal power from Kža-1 was estimated between 0.75 and 1 MWt (Table 2) and water is feasible for the heating of three schools, a gymnasium with swimming pool and greenhouses as part of the agriculture department of the High school and College of agriculture.

CONCLUSION
Two major types of geothermal reservoirs in the study area were determined using geophysical, geological and borehole data: (i) 'basement -BM' reservoirs which are Palaeozoic and Mesozoic fractured/karstified carbonate sediments and/or crystalline/metamorphic rocks, (ii) 'basin fill -BF' reservoirs which are 'Lower Pannonian' and 'Upper Pannonian' sediments composed of coarse and fine-grained sand, sandstones and marls. The BM geothermal aquifers are the most important ones in the study area, with the largest geothermal potential which has resulted in the first geothermal power plant in Croatia to be built there.
Sampled thermal waters from Velika-1, Ciglena belong to the Na-ClHCO 3 hydrochemical type, as a consequence of mixing different waters (oil waters, brines and palaeo-meteoric water) and water-rock interaction. Thermal waters from Daruvar Spa belong to the CaMg-HCO 3 type and represent the typical type of dolomite/limestone thermal aquifers. Sampled waters from Krečaves belong to the Na-HCO 3 type as a consequence of silicate minerals and water-rock interaction. At the Križevci loca-tion, thermal waters belong to the Na-HCO 3 SO 4 Cl type as a consequence of silicate and sulfide mineral weathering. The stable isotope ratios δ 2 H and δ 18 O indicate that monitored thermal waters have a meteoric origin without recent replenishment.
The study area has great geothermal potential for electricity production, but even more for low enthalpy usage in district heating, industry, etc. At locations with future geothermal potential it will be necessary to conduct additional geophysical, geological, hydrogeological and hydrochemical investigations to determine the spatial distribution, thickness, hydrochemical properties and to determine their available thermal power. Calculated total available thermal power for Križevci, Velika-1 and Daruvar spa is 70.5 MWt and only 28 % of this thermal power is used. Agrogeothermal use of lower temperature waters (up to 80° C) in the Bjelovar region also has significant potential because it could lead to the modern development of agriculture and consequently contribute to increasing the standard of living of the local population. It is noticeable that the interest of investors for geothermal projects in Croatia is growing, and hopefully this utilization of geothermal resources -especially with regards to strong and optimistic targets of the New Green Deal -will blossom in the forthcoming years.

ACKNOWLEDGEMENT
The paper is part of the DARLINGe project, co-funded by the European Regional Development Fund (1612249.99 €) and by the Instrument for Pre-Accession Assistance II (534646.6 €) under Grant Agreement No. DTP1-099-3.2. and MAMEBio project, cofunded by the Croatian Science Foundation (HRZZ) PZS-2019-02-7373. In addition, the knowledge gained through the Geo-Twinn project that has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 809943 and contributed to the chemical and isotopic data interpretation. We would like to thank all users for support during the sampling campaigns.  Figure 9. Relationship between isotopic ratios of d 18 O and d 2 H in the sampled thermal waters.