2022 | 75/2 | 223–233 | 9 Figs. | 2 Tabs. | www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society 1. INTRODUCTION The North-western part of Croatia represents the southwestern margin of the Pannonian Basin System (PBS), and it is charac- terized by a high average geothermal gradient (49°C/km) and surface heat flow (76 mW/m2) (BOŠNJAK et al., 1998). Favoura- ble geothermal conditions in this area are a consequence of the very shallow Moho (Mohorovičić seismic discontinuity), at be- tween 25–28 km depth (ALJINOVIĆ & BLAŠKOVIĆ, 1984; ŠUMANOVAC et al., 2009; ŠUMANOVAC et al., 2016). Due to these favourable geological conditions, thermal springs occur in several locations, allowing the people who have inhabited this area throughout history, to use them for bathing, cooking and heating (ŠIMUNIĆ, 2008). Increased utilization of thermal aqui- fers has raised awareness of their protection and sustainability, and feasibility studies appeared in the early 19th century in which protection zones and rules on how to behave with thermal waters were defined (BANCALARI, 1868; RAK, 1876; GRUDELJ, 1890; GORJANOVIĆ-KRAMBERGER, 1896). Today, utilization of thermal aquifers varies from recreation, balneotherapy, space and district heating, greenhouse heating, fish farming, and directly as sanitary water, public water supply and electricity production. Enhanced and sustainable utilization of this resource could significantly contribute to the forecasted geothermal development and achievement of Croatia’s National Renewable Energy Action Plan (NREAP) target. NREAP pre- Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression Tamara Marković1,*, Željka Sladović2, Dragutin Domitrović3, Igor Karlović1 and Ozren Larva1 1 Croatian Geological Survey, Department of Hydrogeology and Engineering Geology, 10 000 Zagreb, Croatia; (*corresponding author: tmarkovic@hgi-cgs.hr) 2 Geoda, 10 000 Zagreb, Croatia 3 Calida Aqua, 10360 Soblinec, Croatia doi: 10.4154/gc.2022.21 Abstract 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 δ2H) and physi- co-chemical (EC, T, pH, DO, Na+, K+, Mg2+, NH4 +, Ca2+, SO4 2-, Cl-, Br-, F-, SiO2 and H2S) 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 δ2H and δ18O indicate that the monitored thermal waters have a meteoric origin, but without recent replenishment. Monitored waters belong to mixed hydrochemical types, from Na-ClHCO3 to Na-HCO3 types in the deep basin thermal waters and a CaMg-HCO3 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. 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 differ- ent 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 GA- ZETTE 153/09, 63/11, 130/11, 56/13, 14/14). Alternatively, if the water is used for heating and electricity production, the Explora- tion 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 north- western 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 sub- depression. 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-depres- sion and to give an overview on potential areas suitable for their direct or indirect use. Article history: Manuscript received January 03, 2022 Revised manuscript accepted May 04, 2022 Available online June 23, 2022 Keywords: geothermal aquifers, basin fill type, basement type, utilization, Bjelovar sub-depres- sion, Croatia G eo lo gi a C ro at ic a Geologia Croatica 75/2224 2. MATERIALS AND METHODS 2.1. 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 pre- cipitation 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čko- križevačka, and 10% of the Zagrebačka County. Moreover, these counties support tourism, agriculture, small and medium enter- prises, and one of the many important factors which ensures fu- ture economic growth is a sustainable source of energy for con- sumption (heat and electricity). 2.2. Data collection Intensive geological data collection of the Bjelovar sub-depres- sion was carried out as part of the Danube Region Leading Geo- thermal Energy project (DARLINGe, 2017–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 Decem- ber 2019 in the scope of the DARLINGe and MAMEBio (A mul- tiphasic approach for deciphering the microbial ecology and bio- technological 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. There- fore, hydrochemical interpretation was based on available chem- ical 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 (O2) were determined using portable WTW elec- trodes. The total alkalinity was measured by titration with 1.6 N H2SO4, using phenolphthalein and bromocresol green-methyl red as indicators. Dissolved cations (Na+, K+, Mg2+, NH4 +, Ca2+) and anions (SO4 2-, Cl-, Br-, F-) were measured using Ion Chromatog- raphy on Dionex ICS-6000 DC, at the Hydrochemical laboratory of the Croatian Geological Survey. Silica concentrations (SiO2) and H2S were measured spectrophotometrically at the same lab- oratory using HACH DR3900. The water quality analyses was checked for a cation–anion balance by measuring the relative de- viation from charge balance (Δmeq = 100 x (Σmeq+ − Σmeq−)/(Σmeq+ + Σmeq−) < ±5%) (MANDEL & SHIFTAN, 1980; DOMENICO & SCHWARTZ, 1990). The δ18Ο and δ2H were determined using Picarro L2130i (Santa Clara, USA) in the Hydrochemical Labo- ratory of the Croatian Geological Survey. The instrument uses Figure 1. The position of the study area with locations of the utilized (Velika-1 and Daruvar spa) and potential geothermal areas (Krečaves, Križevci, Letičani, Pav- lovac, Bačkovica, Velika Trnovitica, Grubišno Polje and Dežanovac). G eologia C roatica Marković et al.: Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression 225 Cavity Ring-Down Spectroscopy (CRDS) technology (http:// www.picarro.com/technology/cavity_ring_down_spectroscopy). All measurements were checked with Picarro’s standards (De- pleted −29.6 ± 0.2 δ18Ο; −235 ± 1.8 δ2H; Mid −20.6 ± 0.2 δ18Ο; −159 ± 1.3 δ2H; Zero 0.3 ± 0.2 δ18Ο; 1.8 ± 0.9 δ2H), 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 δ2H. 3. GEOLOGICAL AND HYDROGEOLOGICAL SETTINGS OF GEOTHERMAL AQUIFERS Two major types of geothermal reservoirs were determined us- ing the methodology developed in the project DARLINGe (RO- TÁR-SZALKAI et al., 2018): (i) ‘basement – BM’ reservoirs, (ii) ‘basin fill – BF’ reservoirs. 3.1. 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/Cre- taceous age) and/or (ii) fractured/fissured crystalline/metamor- phic rocks (gabbro, schist, etc.). Triassic dolomite is the most fa- vourable 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 (tem- peratures > 150°C). Water temperatures range from 40 to 170°C. At the northwest, west and southeast margins of Bjelovar sub- depression, 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 reser- voir 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 devel- opment 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 loca- tions, 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 dur- ing the 1990’s (Ptk-1, VC-2, VC-1, VC-1A) (Fig. 1). All boreholes reached the carbonate rocks (dolomite, dolomite breccia) of Me- sozoic 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 de- termined (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 dolo- mite 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. G eo lo gi a C ro at ic a Geologia Croatica 75/2226 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 geo- thermal water in the underlying aquifer, which is further en- hanced 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. Be- sides 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 dolo- mite breccia) is over 200 m. According to pumping tests and geo- physical 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 per- meability (it is not highly tectonised) is 10%, but in the very per- meable 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 bore- hole 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 tem- perature 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 wa- ters 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 ex- ploration 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 geo- thermal aquifers were detected in BM crystalline/carbonate rocks, with water temperatures in the range of 68 to 170°C (Fig. 1). 3.2. 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 lit- toral 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 sedi- ments 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 deposi- ted 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). G eologia C roatica Marković et al.: Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression 227 a shelf was formed and UP sediments consisting of fine-grained materials were deposited, overlying the LP sediments. Finally, al- luvial 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 per- meability was 380 mD (INA, 1966). Today’s outflow water tem- perature is around 45°C, and is dependent on water outflow dura- tion, which in this case is the time of consumption by local people for showering. Sometimes during the sampling campaign, the wa- ter temperature reached 47°C. The second unused borehole is in the village of Veliko Ko- renovo near the town of Bjelovar, where the water temperature is 30°C (MALJKOVIĆ & GUDMUNDSSON, 2017). Figure 5. Schematic cross section C-C’ of the Ilova basin (according to Fig. 2). Figure 4. Schematic cross section B-B’ of the east part of Bjelovar sub-depression (according to Fig. 2). G eo lo gi a C ro at ic a Geologia Croatica 75/2228 4. 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 Ve- lika-1, Ciglena (23816 µS/cm), followed by water from borehole Kža-1 (2500 µS/cm). EC values for thermal water from the Da- ruvar Spa on all three objects are very similar, and Krečaves has slightly higher values than the Spa waters. Thermal water of Ve- lika-1, Ciglena contains a high amount of dissolved CO2 (30 m3/ m3) comprising 99.5% of gas, about 59 ppm H2S and 24.4 ppm NH4 +. Also, the water contains oil (INA, 1990). Krečaves con- tains 2.4 ppm H2S and 8.5 ppm NH4 +. Thermal waters of the Da- ruvar Spa have much lower concentrations of H2S and NH4 +. The highest concentrations of Na+, Cl-, K+, and SO4 2- are present in thermal water from Ciglena, followed by water from Križevci, Krečaves and Daruvar Spa. F- is observed in all the sampled wa- ters, while Br- is only observed in the waters of Ciglena. According to the chemical composition, thermal waters of Daruvar Spa belong to CaMg-HCO3 type, Velika-1, Ciglena be- longs to Na-ClHCO3 type, Krečaves to Na-HCO3 type, and bore- hole Kža-1 belongs to Na-HCO3SO4Cl type (Fig. 6). 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 dissolu- tion of dolomite, because the Mg2+/Ca2+ ratios are above 0.33. However, according to the chemical composition of thermal wa- ters from Velika-1, Ciglena, the prevailing cation is Na+, indicat- ing that other geochemical processes are present in the thermal aquifer or the origin of the water is different. On the bivariate mix- ing diagram of Na+-normalized Ca2+ versus Na+-normalized Mg2+ (Fig. 8a), the sample falls between the silicate weathering and carbonate dissolution domains, but slightly closer to the lat- ter, indicating that carbonate dissolution is dominant. The influ- ence of cation exchange was evaluated by an equivalent bivariate plot of corrected bivalent cations versus corrected Na+ (Fig. 8b). Concentrations of bivalent cations (Ca2+ and Mg2+) that may have been involved in exchange reactions were corrected by subtract- ing equivalent concentrations of the associated anions (HCO3 - and Table 1. Mean values of the measured hydrochemical parameters. Location No. of samples EC (mS/cm) T (°C) pH HCO3 – (mg/L) Ca2+ (mg/L) Mg2+ (mg/L) Na+ (mg/L) K+ (mg/L) Cl– (mg/L) SO4 2– (mg/L) SiO2 (mg/L) NH4 + (mg/L) H2S (mg/L) F– (mg/L) Br– (mg/L) d18O (o/oo) d2H (o/oo) Krečaves 4 734 45.0 8.1 590 4.2 0.5 198.6 3.3 5.4 35.7 31.2 8.5 2.4 1.0 <0.1 -11.42 -79.39 Ivanovo vrelo 4 616 43.4 7.2 326 73.0 24.0 13.8 3.7 1.5 52.0 39.7 0.2 0.4 0.5 <0.1 -10.20 -71.30 Antunovo vrelo 4 647 47.5 6.9 326 75.1 23.7 13.6 3.7 1.6 52.1 39.8 0.2 0.4 0.5 <0.1 -10.34 -71.20 Borehole D-1 4 653 44.3 6.9 330 76.1 24.7 13.2 3.5 1.7 50.3 39.7 0.3 0.4 0.5 <0.1 -10.25 -70.76 Velika-1, Ciglena 3 23816 150.7 6.8 2450 280.5 65.0 8436.0 485.3 12496.4 117.0 256.5 24.5 59.0 4.8 55.2 5.05 -15.20 Borehole Kža-1 1 2500 70.0 7.2 470 130.0 30.0 740.0 40.0 260.0 355.0 - - - - - - - Figure 6. Piper diagram of mean values of chemical major components in sampled thermal waters. G eologia C roatica Marković et al.: Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression 229 SO4 2-) that would derive from other processes (e.g. carbonate dis- solution 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 equiv- alent 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 (Ca2+ + Mg2+) - (HCO3 - + SO4 2-) against (Na+ - Cl-) (Fig. 8b) samples of both locations are on the line y = -x, which indicates active cation exchange pro- cesses 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 chemi- cal 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 con- tains 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 δ2H and δ18O are distributed around the local meteoric water line (LMWL) Bjelovar (HUNJAK et al., 2013) indicating meteoric origin (Fig. 9). The stable isotope ratios δ2H and δ18O in thermal waters from Velika-1, Ciglena are slightly shifted to the right from the LMWL at Bjelovar because the values are en- riched. The enrichment in δ18O 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, 2010; ARNORSSON 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 iso- topes 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 14C, it is inactive, indicating that the water is very old. In contrast, thermal waters from Daruvar Spa still con- tain 14C and the thermal water age was calculated using 13C cor- rection 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. 5. 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. Now- adays, thermal waters from Ivanovo vrelo, Antunovo vrelo and borehole D-1 are used only for recreation, balneotherapy and san- itary water in the complex of the Special Hospital for Medical Rehabilitation and Thermal water park Aquae Balissae. Conside- ring 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, Figure 7. Molar ratios of magnesium and calcium ions in the thermal waters of Velika-1, Ciglena and Daruvar Spa. G eo lo gi a C ro at ic a Geologia Croatica 75/2230 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 well- head 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 maxi- mum 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 gene- rates 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 loca- tion is 64.9 MWt. In September 2019, the Hydrocarbon Agency opened a ten- der procedure in order to select the most favourable bidder for exploration of geothermal waters for energy purposes in the ex- ploration area “Križevci”. According to EIHP study (MALJKO VIĆ & GUDMUNDSSON, 2017), the available thermal power from Figure 8. a) Bivariate mixing diagram of Na+-normalized Ca2+ versus Na+-normalized Mg2+ in sampled thermal waters; b) Bivariate plot of (Ca2+ + Mg2+) - (HCO3 - + SO4 2-) against (Na+ - Cl-) in sampled thermal waters. a) b) G eologia C roatica Marković et al.: Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression 231 Kža-1 was estimated between 0.75 and 1 MWt (Table 2) and wa- ter is feasible for the heating of three schools, a gymnasium with swimming pool and greenhouses as part of the agriculture de- partment of the High school and College of agriculture. 6. 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/meta- morphic 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 geo- thermal 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-ClHCO3 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-HCO3 type and represent the typical type of dolomite/limestone thermal aquifers. Sampled waters from Krečaves belong to the Na-HCO3 type as a consequence of sili- cate minerals and water-rock interaction. At the Križevci loca- tion, thermal waters belong to the Na-HCO3SO4Cl type as a con- sequence of silicate and sulfide mineral weathering. The stable isotope ratios δ2H and δ18O indicate that monitored thermal wa- ters 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 heat- ing, 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 availa- ble thermal power for Križevci, Velika-1 and Daruvar spa is 70.5 MWt and only 28 % of this thermal power is used. Agrogeother- mal 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 tar- gets 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, co- funded 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. Table 2. The available and utilized thermal power at Križevci, Velika-1 and Daru- var spa locations. Location Available thermal power (MWt) Utilized thermal power (MWt) Križevci 0.75 – 1* – Velika-1 64.9 19 Daruvar Spa 4.6 0.6 *according MALJKOVIĆ & GUDMUNDSSON (2017) Figure 9. Relationship between isotopic ratios of d18O and d2H in the sampled thermal waters. G eo lo gi a C ro at ic a Geologia Croatica 75/2232 REFERENCES ALJINOVIĆ, B. & BLAŠKOVIĆ, I. (1984): Comparison of sediment basis and the Mohorovičić discontinuity in the coastal area of Yugoslavia.− Nafta, Zagreb, 35, 65–71. ARNORSSON, S. (2000): Isotopic and chemical techniques in geothermal exploration, development and use: sampling methods, data handling, interpretation.– Vienna: IAEA, 351 p. BAĆ, J. & HERAK, M. (1962): Prijedlog za određivanje užih i širih zaštitnih zona ter- momineralnih izvora u Hrvatskoj [Proposal for determination of narrower and wider protection zones of thermomineral springs in Croatia – in Croatian].– Un- publ. Technical report, Archive of the Croatian Geological Survey. BANCALARI, D. (1868): Das mineralbad Krapina Toeplitz in Croatien. Krapina spa and its warm springs and baths.– Vienna (in German), 85 p. BOROVIĆ, S., POLA, M., BAČANI, A. & URUMOVIĆ, K. (2019): Constraining the recharge area of a hydrothermal system in fractured carbonates by numerical mod- elling.− Geothermics, 82, 128–149. doi: 10.1016/j.geothermics.2019.05.017 BOŠNJAK, R., ČUBRIĆ. S., GOLUB, M., GRABOVSKI, K., JELIĆ, K., KOLIN, I., KOŠĆAK, S., KULENOVIĆ, I., MIOČEV, D., PRAVICA, Z., RAJKOVIĆ, D., SALOPEK, B., SEČEN, J. & STANČIĆ, L. (1998): GEOEN – Program korištenja geotermalne energije [Programme of geothermal energy utilization – in Croatian].– Energy Institute Hrvoje Požar, Zagreb, 123 p. CARTWRIGHT, I., WEAVER, T., TWEED, S., AHEARNE, D., COOPER, M., CZAP- NIK, K. & TRANTER, J. (2002): Stable isotope geochemistry of cold CO2-bear- ing mineral spring waters, Daylesford, Victoria, Australia: sources of gas and wa- ter and links with waning volcanism.− Chem. Geol. 185, 71–91. doi: 10.1016/ S0009-2541(01)00397-7 CARTWRIGHT, I., WEAVER, T.R., SIMMONS, C.T., FIFIELD, L.K., LAWRENCE, C.R., CHISARI, R. & VARLEY, S. (2010): Physical hydrogeology and environ- mental isotopes to constrain the age, origins and stability of a low-salinity ground- water lens formed by periodic river recharge: Murray Basin, Australia.− J. Hydrol., 380, 203–221. doi: 10.1016/j.jhydrol.2009.11.001 CROSCO (2009): Izvještaj o bušenju bušotine Dar-1 [Drilling report of the borehole Dar-1 – in Croatian].– Unpublished, Fund of professional documentation of Cros- co, Zagreb. CVETKOVIĆ, M., MATOŠ, B., RUKAVINA, D., KOLENKOVIĆ MOČILAC, I., SAFTIĆ, B., BAKETARIĆ, T., BAKETARIĆ, M., VUIĆ, I., STOPAR, A., JARIĆ, A. & PAŠKOV, T. (2019): Geoenergy potential of the Croatian part of Pannonian Basin: insights from the reconstruction of the pre-Neogene basement unconfor- mity.− J. Maps, 15/2, 651–661. doi: 10.1080/17445647.2019.1645052 DARLINGe (2017-2019): Danube Region Leading Geothermal Energy project, http:// www.interreg-danube.eu/approved-projects/darlinge DOMENICO, P.A. & SCHWARTZ, F.W. (1990): Physical and Chemical Hydrogeolo- gy.– John Willey and Sons, New York, 824 p. GORJANOVIĆ-KRAMBERGER, D. (1896): Geology and hydrogeology of the area of Stubičke Toplice. Zagreb. GRUDELJ, A. (1890): Zaštitne zone termalnih izvora i kupelji Vraždinskih Toplica u Varaždinskoj županiji na putu prema Novom Marofu [Protection zones of the bath- ing area of the thermal springs in Varaždinske Toplice in the Varaždin County to- wards Novi Marof – in Croatian]. Zagreb, 25p. HORVÁTH, F. (1993): Towards a mechanical model for the formation of the Pannonian Basin.– Tectonophysics, 225, 333–358. doi: 10.1016/0040-1951(93)90126-5 HUNJAK, T., LUTZ, H.O. & ROLLER-LUTZ, Z. (2013): Stable isotope composition of the meteoric precipitation in Croatia.– Isot. Environ. Health Stud., 49/3, 336–45. doi: 10.1080/10256016.2013.816697 INA (1966): Geološki i DTS izvještaj bušotine Ze-1 [Geological and DTS reports of well Ze-1 – in Croatian].– Unpublished, Fund of professional documentation of INA, Zagreb provided by Ministry of Economy and Energy, 25 p. INA (1990): Geološki i DTS izvještaj bušotina VC-1 i VC-1A [Geological and DTS re- ports of wells VC-1 and VC-1A – in Croatian].– Unpublished, Fund of profession- al documentation of INA, Zagreb provided by Ministry of Economy and Energy 50 p. JANKOWSKI, J., ACWORTH, R.I & SHEKARFOROUSH, S. (1998): Reverse ion ex- change in deeply weathered porphyritic dacite fractured aquifer system, Yass, New South Wales, Australia.– In: AREHART, G.B. & HULSTON, J.R. (eds.): Proceed- ings of the 9th International Symposium on Water-Rock Interaction, Taupo, New Zealand, 30 March - 3 April, 1998. Rotterdam, A.A. Balkema, 243–246. LARVA, O., MRAZ, V. & MARKOVIĆ, T. (2008): Daruvarske toplice - Elaborat utvrđivanja eksploatacijske izdašnosti Ivanovog vrela i objekta Š-3 [Daruvarske Spa - Study for determining the exploitation yields of Ivanovo vrelo and Š-3 – in Croatian].– Unpubl. Technical report, Archive of the Croatian Geological Survey, 10 p. LARVA, O., MRAZ, V. & MARKOVIĆ, T. (2009): Hidrogeološki nadzor tijekom bušenja, ugradnje, čišćenja i pokusnog crpljenja istražno – eksploatacijke bušotine Dar-1 u Daruvaru [Hydrogeological monitoring during drilling, installation, clean- ing and experimental pumping of the exploration-exploitation well Dar-1 in Da- ruvar – in Croatian].– Unpubl. Technical report, Archive of the Croatian Geolog- ical Survey, 12 p. MAGYAR, I., RADIVOJEVIĆ, D., SZTANÓ, O., SYNAK, R., UJSZÁSZI, K. & PÓC- SIK, M. (2013): Progradation of the paleo-Danube shelf margin across the Pan- nonian Basin during the Late Miocene and Early Pliocene.– Glob. Planet. Change, 103, 168–173. Doi: 10.1016/j.gloplacha.2012.06.007 MALJKOVIĆ, D. & GUDMUNDSSON J.R. (eds.) (2017): Geothermal Energy Utilisa- tion Potential in Croatia by SANJA ŽIVKOVIĆ (EIHP), RUDOLF VAN HEM- ERT (EIHP), DRAŽEN TUMARA (EIHP), KARMEN STUPIN (EIHP), J. RÚNAR MAGNÚSSON (EFLA Consulting Engineers, Iceland), HEIMIR HJAR- TARSON (EFLA Consulting Engineers, Iceland), JILL ROBINSON HAIZLIP (Geologica Geothermal Group), MANON M. STÖVER, (Geologica Geothermal Group), EIHP, Zagreb, 57 p. MALVIĆ, T. & VELIĆ, J. (2010): Relations between effective thickness, gas production and porosity in heterogeneous reservoirs: an example from the Molve Field, Cro- atian Pannonian Basin.– Pet. Geosci., 16, 41–51. doi: 10.1144/1354- 079309-722 MANDEL, S. & SHIFTAN, Z.L. (1980): Groundwater Resources: Investigation and Development.– Academic Press, New York, 269 p. MARKOVIĆ, T. & LARVA, O. (2012): Hydrochemistry of thermal water in Daruvar area.– In: STEVANOVIĆ, Z. (ed.): Proceedings of 14th Serbian symposium of Hydrogeology with international participation, University of Belgrade, Mining and Geology Faculty, Belgrade, 187–191. MARKOVIĆ, T., BOROVIĆ, S. & LARVA, O. (2015): Geochemical characteristics of thermal waters of Hrvatsko zagorje.− Geol. Croat., 68/1, 67–77. doi: 10.4154/ GC.2015.05 MCLEAN, W. & JANKOWSKI, J. (2000): Groundwater quality and sustainability in an alluvial aquifer, Australia.– In: SILILO, A. (ed.): XXX IAH congress on ground- water: past achievements and future challenges. A.A. Balkema, Rotterdam. MRAZ, V. (1983): Izvještaj o hidrogeološkim istražnim radovima na području Daruvar- skih toplica II. faza [Report on hydrogeological exploration works in the area of Daruvar Spa II phase – in Croatian].– Unpubl. Technical report, Archive of the Croatian Geological Survey, 13 p. NÁDOR, A., SEBESS-ZILAHI, L., ROTÁR-SZALKAI, Á., GULYÁS, Á. & MARKOVIĆ, T. (2020): New methods of geothermal potential assessment in the Pannonian basin.− Neth. J. Geosci., 98, e10. doi: 10.1017/njg.2019.7 OFFICIAL GAZETTE 153/09, 63/11, 130/11, 56/13, 14/14.– Water Act. OFFICIAL GAZETTE 52/18, 52/19.– Exploration and Exploitation of Hydrocarbons Act. OFFICIAL GAZETTE 56/13, 14/14.– Mining Act. PAVELIĆ, D. (2001): Tectonostratigraphic model for the North Croatian and North Bos- nian sector of the Miocene Pannonian Basin System.− Basin Res., 13, 359–376. doi: 10.1046/j.0950-091x.2001.00155.x PRELOGOVIĆ, E., SAFTIĆ, B., KUK, V., VELIĆ, J., DRAGAŠ, M. & LUČIĆ, D. (1998): Tectonic activity in the Croatian part of the Pannonian basin.− Tectono- physics, 297, 283–293. doi: 10.1016/S0040-1951(98)00173-5 RAK, A. (1876): Krapinske toplice [Krapina spa – in Croatian].– Vienna 45 p. RMAN, N. (2016): Hydrogeochemical and isotopic tracers for identification of seasonal and long-term over-exploitation of the Pleistocene thermal waters.− Environ. Monit. Assess., 188(4), 242. doi: 10.1007/s10661-016-5250-2 RMAN, N., GÁL, N., MARCIN, D., WEILBOLD, J., SCHUBERT, G., LAPANJE, A., RAJVER, D., BENKOVÁ, K. & NÁDOR, A. (2016a): Potentials of transbound- ary thermal water resources in the western part of the Pannonian basin.− Geother- mics, 55, 88–98. doi: 10.1016/j.geothermics.2015.01.013 RMAN, N., LAPANJE, A., PRESTOR, J. & SULLIVAN, M.J.O. (2016b): Mitigating depletion of a porous geothermal aquifer in the Pannonian sedimentary basin.− Environ. Earth Sci., 75, 723. doi: 10.1007/s12665-016-5634-1 ROTÁR-SZALKAI, Á., MAROS, GY., BERECZKI, L., MARKOS, L., BABINSZKI, E., ZILAHI-SEBESS, L., GULYÁS, Á., KUN, É., SZŐCS, T., KERÉKGYÁRTÓ, T., NÁDOR, A., RAJVER, D., LAPANJE, A., ŠRAM, D., MARKOVIĆ, T., VRANJEŠ, A., FARNOAGA, R., SAMARDŽIĆ, N., HRVATOVIĆ, H., SKOPLJAK, F. & JOLOVIĆ, B. (2018): Identification, ranking and characteriza- tion of potential geothermal reservoirs.− Report of the DARLINGe project, 82 p. Available at www.interreg-danube.eu/approved-projects/darlinge/outputs. ROTÁR-SZALKAI, Á., NÁDOR, A., SZŐCS, T., MAROS, G., GOETZL, G. & ZEKIRI, F. (2017): Outline and joint characterization of transboundary geothermal reservoirs at the western part of the Pannonian basin.− Geothermics, 70, 1–16. doi: 10.1016/j.geothermics.2017.05.005 ROYDEN, L. (1988): Late Cenozoic Tectonics of the Pannonian Basin System.– In: ROYDEN, L. H. & HORVÁTH, F. (eds.): The Pannonian Basin (Memoir 45), American Association of Petroleum Geologists, 27–48. SAFTIĆ, B., VELIĆ, J., SZTANO, O., JUHASZ, G. & IVKOVIĆ, Ž. (2003): Tertiary subsurface facies, source rocks and hydrocarbon reservoirs in the SW part of the Pannonian Basin (northern Croatia and south-western Hungary).− Geol. Croat., 56/1, 101–122. doi: 10.4154/232 SEBE, K., KOVAČIĆ, M., MAGYAR, I., KRIZMANIĆ, K., ŠPELIĆ, M., BIGUNAC, D., SÜTŐ-SZENTAI, M., KOVÁCS, Á., SZUROMI-KORECZ, A., BAKRAČ, K., HAJEK-TADESSE, V., TROSKOT-ČORBIĆ, T. & SZTANÓ, O. (2020): Cor- G eologia C roatica Marković et al.: Current utilization and hydrochemical characteristics of geothermal aquifers in the Bjelovar sub-depression 233 relation of upper Miocene-Pliocene Lake Pannon deposits across the Drava Basin, Croatia and Hungary.− Geol. Croat., 73/3, 177–195. doi: 10.4154/gc.2020.12 SZŐCS, T., RMAN, N., SÜVEGES, M., PALCSU, L., TÓTH, G. & LAPANJE, A. (2013): The application of isotope and chemical analyses in managing transbound- ary groundwater resources.− Appl. Geochemistry, 32, 95–107. doi: 10.1016/j.ap- geochem.2012.10.006 ŠARIN, A. (1990): Tehnički izvještaj o pokusnom crpljenju termomineralne bušotine Kža-1 u Križevcima [Technical report on the pumping test at the well Kža-1 in Križevci – in Croatian].– Technical report, Archive of the Croatian Geological Survey. ŠIMUNIĆ, A. (ed) (2008): Geotermalne i mineralne vode RH [Geothermal and mineral waters of the Republic of Croatia – in Croatian].– Croatian Geological Survey, Zagreb, 343 p. ŠUMANOVAC, F., HEGEDŰS, E., OREŠKOVIĆ, J., KOLAR, S., KOVÁCS, A.C., DUDJAK, D. & KOVÁCS, I.J. (2016): Passive seismic experiment and receiver functions analysis to determine crustal structure at the contact of the northern Di- narides and southwestern Pannonian Basin.− Geophys. J. Int., 205, 3, 1420–1436. doi: 10.1093/gji/ggw101 ŠUMANOVAC, F., OREŠKOVIĆ, J., GRAD, M. & ALP 2002 Working Group (2009): Crustal structure at the contact of the Dinarides and Pannonian basin based on 2-D seismic and gravity interpretation of the Alp07 profile in the ALP 2002 experiment.− Geophys. J. Int., 179, 615–633. doi: 10.1111/j.1365-246X.2009. 04288.x ZANINOVIĆ, K., GAJIĆ-ČAPKA, M. & PERČEC TADIĆ, M. (2008): Klimatski atlas Hrvatske. Climate atlas of Croatia: 1961.–1990, 1971–2000.− Državni hidrometeorološki zavod. Zagreb, Croatia, 200 p.