Groundwater flow system understanding of the lukewarm springs in Kistapolca (South Hungary) and its relevance to hypogene cave formation

Springs have an important role both in groundwater flow system understanding and in maintaining groundwater-related ecosystems. The aim of the research of the lukewarm karst springs in Kistapolca is to understand the origin of the elevated temperature (22–24 °C), i.e., whether it is the result of mixing of cold and thermal waters or it represents the terminal area of intermediate flow systems. This question has also an important relevance to the evaluation of recent cave formation in the area. In case of mixing of different waters dissolution by mixing corrosion could take place. This scenario was proposed for the formation of the Beremend Cave, which is situated in the close vicinity of the springs, and which is the only known cave of the area, where the passages are partly filled by water. Therefore, several methods, including geophysical, geochemical and stable isotope measurements, evaluation of continuous-time series of water level change, temperature, and electrical conductivity data and hydraulic evaluation of vertical flow conditions by pressure–elevation profiles were combined in order to acquire information about the hydrogeological environment and flow conditions in the area. The results show that a higher order flow system derived from the Villány Hills feeds the springs in Kistapolca without significant local influence. In addition, the results related to the Beremend Cave indicate that cave formation by mixing corrosion is not an active process today in the area, the cave is rather influenced by local processes.


Introduction
Springs are the surface manifestations of groundwater flow systems (Tóth 1971). Their appearance and physicochemical parameters represent the flow system as a whole and report on the rock-water interactions along the flow path. Nowadays most of the springs lost their natural character, because of anthropogenic influence, i.e., they are captured and intensely used. Therefore, the rather rare natural springs are of great importance to study the flow system of an area (Groves 2007). Moreover, natural springs usually maintain groundwater-dependent ecosystems in forms of creeks, lakes, wetlands, etc., which also emphasizes their value.
Lukewarm springs are characterized by slightly higher temperatures than the average annual air temperature, but considerably lower than the human body temperature (Meinzer 1923;Pentecost et al. 2003). They represent either a discharge of intermediate flow systems or a mixture of waters with different temperatures derived from different order flow systems. If the latter case takes place in a carbonate rock environment, the mixing process can lead to the "mixing corrosion" phenomenon (Bögli 1964(Bögli , 1980 resulting in carbonate dissolution, i.e., cave formation, because the waters present in this process are also characterized by different chemical composition (CO 2 content) beside the different temperature. Mixing corrosion is a common cave forming process in hypogene karst areas, which are typical for regional discharge areas of karstic groundwater flow systems. To these areas different order flow systems convey waters with different temperatures and geochemical composition. One of the well-known examples for such situation is the hypogene karst system in Budapest (Hungary) (Erőss et al. 2012;Leél-Őssy 2017;Mádl-Szőnyi et al. 2017). Moreover, there are many other areas where recent thermal or lukewarm waters and caves are associated (Bakalowicz et al. 1987;Plan et al. 2006Plan et al. , 2012Spötl et al. 2009;Dublyansky et al. 2017). The cave morphologies and the mineral assemblages often reflect clearly the effects of upward flowing waters with elevated temperatures (Palmer 2007;Klimchouk 2012;Leél-Őssy 2017).
The karst springs in the little village of South Hungary, Kistapolca ( Fig. 1), are captured, but the only naturally outflowing springs in the Villány karst area nowadays. Their temperatures are about 22-24 °C, well above the annual mean temperature (11 °C). Previous studies by Rónaki et al. (1967) and Lorberer and Rónaki (1978) suggested that the lukewarm temperature is the result of mixing between thermal and cold karst waters. To investigate this question is of great importance since only 6 km away from Kistapolca can be found one of the most protected hypogene caves of Hungary, the Beremend Cave (Fig. 1), which still have underwater passages. The morphology and the minerals decorating the cave walls indicate the effect of thermal waters (Takács-Bolner 1985, 2003aVigassy et al. 2010). The question can be raised whether the lukewarm spring water temperature is the result of mixing, and whether this mixing together with the resulting cave formation is active even today, in the recent hydrogeological settings. Since understanding of karstification requires clarifying the groundwater flow conditions in an area (Klimchouk 2007(Klimchouk , 2012, these questions can be answered by studying the only natural lukewarm springs and its hydrogeological background in the vicinity, i.e., in Kistapolca. The springs in Kistapolca also maintain a lake, which has been dried up several times in the last few years. Flow system understanding can also help to explore the reason of these dry periods and constrain the mechanism, as groundwater-dependent ecosystems maintained by the different order flow system's discharge react differently to climatic changes or anthropogenic effects (Havril et al. 2018).
The aim of the study is to explore the hydrogeological characteristics of the Kistapolca area in order to (1) understand the origin of the elevated temperature; (2) to analyze its relevance to hypogene cave formation; and (3) to evaluate the vulnerability of the lake in Kistapolca as a groundwaterdependent ecosystem maintained by the springs.

Geological and hydrogeological background
Kistapolca is located in the southern foreground of the Villány Hills, in the south-western part of the Pannonian Basin in Hungary (Fig. 1). The 30 km long Villány Hills is built up by various Mesozoic carbonates (Vadász 1935;Nagy and Nagy 1976;Török 1998). The carbonate rocks continue in the foreland basin basement, where they are covered by young sediments (Fig. 2). The whole carbonate rock suite forms a thick (up to 1700 m) karst reservoir. The younger cover sediments of Miocene, Pliocene, Pleistocene age consist of sand, clayey silt and loess, which form 0-450 m cover above the carbonate aquifer with alternating aquifer/ aquitard characteristics in the surroundings of Kistapolca (Deák et al. 1969). Due to the Late Cretaceous nappe formation, which created the seven structural imbrications of the Mesozoic carbonates (Császár 2002;Petrik 2009;Konrád et al. 2010) (Fig. 2), some basement carbonate outcrops can be also found in the southern foreland of the Villány Hills on or close to the surface. Kistapolca is situated at one of these situations (Fig. 2).
Such marginal areas of outcropping carbonates and sedimentary basins with the same carbonate rock basement are favorable sites for the development of hypogene caves ). In the 10 km surroundings of Kistapolca, there are two hypogene caves, in Nagyharsány and Beremend (Fig. 2). These caves are characterized by Fig. 2 Schematic geological build-up of the Villány Hills area (modified after Dezső et al. 2004;Rakusz and Strauss 1953 phenomena that can be related to the influence of thermal water (Takács-Bolner 1985;Takács-Bolner 2003a, b;Vigassy et al. 2010). Recently, only the Beremend Cave has connection to the karst aquifer, as it has underwater cave passages. Liebe and Lorberer (1981) described the groundwater flow system in the Villány area as gravity and temperature (density)-difference driven flow system, where on the surface of bare carbonates (i.e., in the Villány Hills) infiltrating meteoric waters circulate and discharge in karst springs with different temperatures. However, they called attention to the possibility of other water components because according to their investigations in the Villány Hills infiltrating meteoric waters can sustain only the lukewarm springs' discharge. Therefore, they suggested connection possibilities of the carbonate reservoir with other deep aquifer units. In the case of the lukewarm springs in Kistapolca, Lorberer and Rónaki (1978) proposed that the cold karst waters originate from direct infiltration in the Villány Hills, and the thermal component is derived from the SW. Csicsák and Gondár (2008) suggested the effect of tectonic compression-and compaction-driven fluid flow from the South, from the sedimentary layers of the Drava Basin (part of the Pannonian Basin). Based on trace element measurements (Scheuer 2012), active tectonics induced mantle fluid contribution may also play a role in the geochemical composition of thermal waters.
The Kistapolca area was subject of intense research in the 1950s and 1960s for water supply purposes of the cement factory ( Fig. 1). In the framework of this research, geoelectric measurements were conducted to target the uplifted carbonate basement blocks (Rónaki et al. 1967;Deák et al. 1969). In addition, two springs were captured with 8-12 l/s discharge (Kessler 1959;Rónaki et al. 1967), which fed the lake in Kistapolca. Lorberer and Rónaki (1978) have already mentioned the possible effect of the cement factory due to planned intense water abstraction on the springs' discharge by this time.

Methods
To explore the hydrogeological characteristics of the Kistapolca area different techniques, including geophysical, geochemical, and stable isotope measurements, evaluation of continuous-time series of water level, temperature, and electrical conductivity data and hydraulic evaluation of vertical flow conditions by pressure-elevation profiles were combined.
Since the appearance of the springs in Kistapolca is thought to be connected to an uplifted carbonate basement block, geophysical methods were applied first to investigate this question in detail. Pressure-elevation profiles were constructed based on archival and recently measured well data to investigate the groundwater flow direction, i.e., prove the upward flow conditions in the area, which may explain the existence of the springs. The main ions and stable isotopes report on the geochemistry and rock-water interactions, the continuous-time series of water level, temperature, and electrical conductivity data of the springwater reflect the dynamics of the flow system. Finally, the possibility of that scenario, in which the Villány Hills infiltrated meteoric water can be heated up and discharge in Kistapolca with elevated temperature, was tested by numerical simulation.

Geophysical investigation
The aim of the geophysical measurements was to map the surface of the uplifted carbonate basement block in detail, identify the thickness of the cover sediments and the structural elements. Previous geoelectrical measurements by Rónaki et al. (1967) and Deák et al. (1969) provided only a rough overview of the area. Vertical electrical soundings (VES) is often used to determine the depth and thickness of different geological units or the fracture distribution. In this study, VES using Schlumberger array was carried out at 10 VES stations around the lake in Kistapolca (Fig. 3), with the current electrode spacing of AB/2 = 70 M. From the apparent resistivity values (ρa) of VES, the subsurface resistivity distribution was calculated with the help of two computer programs: PISE4 (Geophysical Department of the Université VI of Paris; Gabalda and Tabbagh 1994) and VES (Department of Geophysics, Eötvös Loránd University, Salát and Drahos 2004). In this way, the layer thickness and resistivity values were obtained. The measurements were supplemented by stratigraphic information of nearby boreholes. Mádl-Szőnyi and Tóth (2015) proved that the gravity-driven regional groundwater flow concept (GDRGF concept) is also applicable in karst regions based on the principle of hydraulic continuity. Therefore, this concept and its methodology were applied here as a working tool. Accordingly, archival and recently measured hydraulic data of the wells were analyzed by one of the basic regional hydraulic methods of hydrogeology, using pressure-elevation [p(z)] profiles (Tóth 2009). The main goal was to determine the regime characteristic of the Kistapolca area and the direction of groundwater flow, i.e., downward, horizontal or upward. The required data (well coordinates, wellhead elevation, depth of the static water level, the elevation of the screened interval) for this hydraulic data analysis were collected from archival well documentations and from recent measurements and summarized in Table 1.

Hydraulic data evaluation: pressure-elevation profiles
During the hydraulic calculations, measured hydraulic head data were converted to pore-pressure values with the help of the following equation, which assumes constant fluid density along the flow path: where h is the measured hydraulic head; z is the elevation of the measuring point with respect to the Baltic Sea level; p is pore pressure; (ρ × g) is the specific weight of the fluid, which is numerically equal to the hydrostatic vertical pressure gradient of that fluid. Pore-pressure data and the hydrostatic vertical pressure gradient [9.81 MPa km −1 ] indicated on the p(z) profiles were determined by applying the water density of 1000 kg m −3 .
The p(z) profiles allow for the examination of the vertical component of fluid flow directions by comparing the vertical pressure gradient to the ideal hydrostatic condition. The vertical pressure gradient (γ) is the change in pore pressure per unit vertical length in a flow domain. Ideally, in the hydrostatic condition, there is no vertical fluid flow because the buoyancy and the force of gravity on water particles are in balance. In this case, the vertical pressure gradient is hydrostatic, and its average value is γ st = 9.81 MPa km −1 , by applying the water density of 1000 kg m −3 . In the hydrodynamic state, a driving force can shift the system off balance, therefore the vertical flow component occurs too, and the vertical pressure gradient changes to dynamic (γ dyn ). In recharge areas γ dyn < γ st , so there is a downward flow direction and γ dyn is called subhydrostatic, while in discharge areas γ dyn > γ st , thus there is an upward flow direction and γ dyn is called super-hydrostatic. In midline areas, the vertical pressure gradient is hydrostatic (γ dyn = γ st ); therefore, the fluid flow has no vertical component; nevertheless, the horizontal component cannot be excluded. Wells interpreted in one p(z) profile had to have approximately the same land surface elevation, because the position of the reference hydrostatic pressure gradient line depends on it.

Hydrogeochemical investigation including stable isotopes
For hydrogeochemical and stable isotope investigations, 17 water samples were collected from dug wells, operating drilled wells, from the mainspring in Kistapolca and from the Beremend Cave at the same time (28-30 August, 2017). Water samples for stable isotope measurements were repeatedly collected from the Beremend Cave (4 samples) and from the mainspring in Kistapolca (2 samples), however, the temporal variation of the measured values were Fig. 3 Detailed map of the study area showing the bounding areas pressure-elevation profiles, the location of the wells involved in pressure-elevation profiles, the location of the hydrogeochemical and isotope sampling sites, the location of vertical electrical sounding (VES) profiles, the Beremend Cave, and the location of A-B profile in Fig. 4 not evaluated during this study, due to the small number of samples. Temperature and electrical conductivity were recorded on the field by YSI ProPlus. Among the main ions Ca 2+ , Mg 2+ , HCO 3 − were analyzed by acid-base titrimetry, Cl − by complexometric titration, Na + and K + were detected by flame photometry and SO 4 2− by photometry at the Department of Geology of Eötvös Loránd University.
The major components of the water samples were investigated and illustrated on Piper diagram. Piper diagrams represent the main anion and cation compositions in the percentage of the total equivalents per liter. The waters are also classified using the hydrochemical facies concept of Back (1966), which is also illustrated on the Piper diagram. Based on the dominant anion and cation, the Ca + Mg/HCO 3 ratio was calculated in order to evaluate their origin.
Stable hydrogen and oxygen isotope measurements of groundwater were determined at the Institute for Geological and Geochemical Research in Budapest using the liquid water isotope analyzer, manufactured by Los Gatos Research Ltd. (LWIA-24d). The instrument uses an off-axis integrated cavity ring-down spectroscopy to measure the absolute abundances of 2 H 1 H 16 O, 1 H 1 H 18 O, and 1 H 1 H 16 O via laser absorption (Lis et al. 2008). Thus, the δ 2 H and δ 18 O values were determined simultaneously and expressed relative to VSMOW standard in ‰. Standardization was conducted using laboratory standard water samples that were calibrated by the international standards (Czuppon et al. 2018). The precisions are better than ± 1.0‰ and ± 0.15‰ for hydrogen and oxygen, respectively.

Continuously recorded data
A Dataqua continuous measuring device recorded the change in water level, temperature, and electrical conductivity values in every 6 hours in the captured mainspring in Kistapolca (accuracy: ± 0.1 °C, EC ref. 20 °C, accuracy ± 1%). The daily precipitation data were obtained from the meteorological station in Nagyharsány (Fig. 1) and were compared with the water level, temperature and electrical conductivity changes in the mainspring.

Numerical simulation
Numerical simulation of groundwater flow and heat transport were performed to establish a hydrogeologic conceptual model depicting and confirming the groundwater flow directions and processes in the vicinity of the lukewarm spring. The 2D numerical model was built in COMSOL Multiphysics 5.1. applying the fluid flow/porous media and subsurface flow/Darcy's law and heat transfer/heat transfer in porous media modules operating with the finite element numerical approach (Zimmerman 2006 The boundary conditions were determined on the basis of measured hydraulic and temperature data of wells and they are shown in Fig. 3. Hydraulic head values were specified along the boundaries, on the upper boundary hydraulic head changes from point to point according to the observed water table. The upper condition for heat transport is the annual mean air temperature (11 °C), while the lower temperature was 22 °C with respect to measured water temperature. Symmetry type boundary conditions were applied on the side borders meaning that there is no heat flux across the boundary.
The model consists of two main hydrostratigraphic units: surface siliciclastic confining layer and thick carbonates (the bottom of this unit is unknown) and five numerical domains. Regarding the hydraulic conductivity, one order of anisotropy was applied (the horizontal one is 10 times higher than the vertical hydraulic conductivity). The hydraulic (porosity and hydraulic conductivity) properties and the thermal conductivity of the units were determined on the basis of analogs and literature references (Mádl-Szőnyi and Tóth 2015; Freeze and Cherry 1979;Brassington 2017;Beardmore and Cull 2001;Clark 1966).

Geophysical measurements
Based on the VES measurements and stratigraphic data of boreholes, the position of the carbonate basement block was refined (Fig. 4), compared to the results of the geophysical survey in the 1960s (Rónaki et al. 1967). The sedimentary cover layers of carbonates are characterized by 7-51 Ωm apparent resistivity values, whereas the carbonate basement's apparent resistivity values vary between 175 and 700 Ωm. Thus, the basement carbonates can be clearly separated from the cover sediments by this method. The depth of the carbonate basement, i.e., the sediment cover thickness, varies between 1 and 130 m in the Kistapolca area. The basement situated in the highest position, i.e., closest to the surface in the close vicinity of the lake. The carbonate rock surface sharply deepens to the east, southeast from the lake (up to 130 m in borehole Nr 9., Fig. 4); therefore, structural elements (faults) can be hypothesized. To the west, the depth of the carbonate basement increases evenly and reaches 40 m at the edge of the village (borehole Nr 32., Fig. 3).

Pressure-elevation profiles
Two p(z) profiles were constructed in order to evaluate the vertical flow directions (Fig. 5). The area of p(z) profile #1 is located in the village of Kistapolca (Fig. 3). The data of dug and drilled wells represent the 53.29-95.63 m asl elevation (z) interval (Table 1). The vertical pressure gradient is greater than the hydrostatic (γ st = 9.81 MPa km −1 ), i.e., super-hydrostatic (γ dyn = 10.45 MPa km −1 ) and refers to the upward vertical flow direction (Fig. 5). The area of p(z) profile #2 is located in SE vicinity of the settlement (Fig. 3). The profile was compiled based on the wells of the cement factory and represent a wider elevation interval, between (− 74.96)-71.55 m asl (Table 1). A slightly smaller than hydrostatic, i.e., subhydrostatic vertical pressure gradient (γ dyn = 9.74 MPa km −1 ) could be inferred, which refers rather to hydrostatic conditions (Fig. 5). That means that in this area horizontal flow is dominant, the water flow has no vertical component.

Hydrogeochemistry
Using the main cations and anions, the geochemical facies (Back 1966) was determined and with the help of the Piper diagram it was illustrated (Fig. 6). The HCO 3 , Cl + SO 4 anion facies with dominating bicarbonate anion is prevalent in the area. Regarding cation facies the waters can be divided into two groups: Ca + Mg cation facies and Ca + Mg, Na + K cation facies occur in the area. In the majority of the samples, calcium ion is dominating, in two cases magnesium dominance is observed (Nr. 1 and Nr. 5 in Table 2.). According to the dominant Ca + Mg cation and HCO 3 anion their ratio was calculated (Table 2), which is close to one, except the two magnesium dominating sampling sites (Nr. 1 and Nr. 5).
The onsite recorded temperature of the water samples varies between 14 and 26 °C. There is a pattern regarding the spatial distribution of temperature data (Fig. 7). In the eastern part of the village, temperatures are higher than 20 °C, while in the western part values around 16 °C occur. The electrical conductivity values, which refer to the dissolved solid content of the waters, vary between 425 and 1361 µS/cm ( Table 2). The lowest value was measured in the Beremend Cave as 425 µS/cm (Nr. 1 in Table 2). The highest values were measured in dug wells in Kistapolca (Nr. 15,27 in Table 2).

Stable isotopes
The δD and δ 18 O values of the studied wells, springwater, and the waters from Beremend Cave show relative large variation, and they are plotted close to the local meteoric  (Table 2), while the most negative stable isotope composition (δD ≈ − 77.9‰ and δ 18 O ≈ − 10.9‰) were found in the Kistapolca springwater (Nr. 14) and some drilled and dug wells (Nr. 7,8,17,21,30,31) in Kistapolca. The other studied wells show intermediate values between these "endmembers". Interestingly, these dug wells (Nr. 21,30,31) in Kistapolca are constructed in the cover sediments, which show similar negative isotopic composition to those wells producing waters from the carbonate aquifer (Nr. 7,8,17) and also similar to the springwater (Nr. 14). Water sample from the Beremend Cave (Nr. 1), however, shows similarities to those wells which are producing waters from the sedimentary cover of the carbonate aquifer (Nr. 3, 5, 6, 15).

Continuously recorded data
The data of continuous measurement device in the captured main spring show that both the temperature (22.4-23.5 °C) and electrical conductivity values (890-911 µS/cm) are rather stable in the observed half-year period. The observed variation of the temperature is within 1 °C and the observed variation of the electrical conductivity values is maximum 20 µS/cm. Both can be considered rather stable. Effects of precipitation periods, however, can be visually recognized in water level increase (Fig. 9).

Numerical modeling
The outcropping carbonates function as recharge areas, where the groundwater is replenished and intense, mostly downwardoriented flow occurred. The magnitude of groundwater flow (color scale) was 1-2 orders less in the confined carbonate and the confining bed itself Fig. 10. The direction of water movement is represented by normalized uniform flow vectors, which showed a convergence zone and a discharge area under the village of Kistapolca in the carbonate dome. The spring received groundwater both from the Northern and Southern parts. The infiltrating cold meteoric water transported by advection penetrated down to the bottom of the model geometry. However, the lukewarm groundwater could rise and reach the surface. Under and within the siliciclastic confining layer, the groundwater could be heated up and 15-20 °C of temperature could be expected and observed. Fig. 6 Piper diagram of the sampled waters (data in Table 2, the location of the sampled wells in Fig. 3)

Table 2
Temperature, electrical conductivity, main anions and cations, and stable isotope data of the sampled waters

Discussion
Based on the VES measurements and stratigraphic data of boreholes, the position of the carbonate basement block was refined in Kistapolca ( Fig. 4) compared to the results of the geophysical survey in the 1960s (Rónaki et al. 1967).
The results show that the basement situated in the highest position, i.e., closest to the surface in the close vicinity of the lake, which explains the existence of the springs. The pressure-elevation profiles showed greater than the hydrostatic (γ dyn = 10.45 MPa km −1 ), i.e., super-hydrostatic vertical pressure gradient in Kistapolca, which refers to the upward vertical flow direction (Fig. 5). This also supports and explains the existence of springs. This method as part of the gravity-driven regional groundwater flow (GDRGF) concept's methodology (Tóth 2009) was previously commonly applied only in sedimentary basins. However, recent investigations (Mádl-Szőnyi and Tóth 2015;Erhardt et al. 2017) highlighted its applicability in karst areas. Pressure-elevation profiles were applied by Erhardt et al. (2017) in a similar, marginal karst area, where the springs and discharge areas were verified by this method. The results in Kistapolca further strengthen the applicability of flow system analysis methodology in karst areas. The VES measurements, the stratigraphic data of boreholes and the spatial distribution of temperature in wells (Fig. 7), suggest the existence of faults within the uplifted basement carbonate block, which can serve as pathways for the upward flow. This upward flow and the discharge effect of karst waters into the cover sediments are also evident from the isotopic composition of dug wells (Nr. 21,23,29,30,31) and springwater (Nr. 14), that are characterized by low hydrogen and oxygen isotope values (δD≈ − 77.9‰ and δ 18 O≈ − 10.9‰). This is also reflected by the bicarbonate and calcium dominance in these wells. The Ca + Mg/HCO 3 ratio is close to 1 also implying carbonate dissolution-related origin of these waters (Zaidi et al. 2015). Higher (1268-1361 µS/cm) electric conductivity values of dug wells (e.g., Nr 15 and 27 in Table 2.) may also indicate anthropogenic influence, which may originate from livestock breeding in the backyard of the houses.
The hydraulic characteristics of karst aquifers are very different from other aquifer types, as they hold a complex conduit network (Goldscheider and Drew 2007). Precipitation events may cause sudden water level increase in the Fig. 7 Temperature distribution in the sampled wells (data in Table 1 and Table 2) Fig. 8 Isotopic composition of the sampled waters (data in Table 2). GMWL and LMWL indicates the global meteoric water line and the local meteoric water line (Czuppon et al. 2017), respectively. As a reference, the annual amount weighted isotopic composition of precipitation is also indicated (Czuppon et al. 2017) 132 Page 12 of 15 system and discharge peek in the springs accompanied by the decrease of electrical conductivity (Király and Müller 1979;Groves 2007). The visually recognizable connection between precipitation events and the water level increase indicates the dynamic response of the thick carbonate reservoir on precipitation in the form of regional water level increase (Fig. 9). The stability of electrical conductivity and temperature over time in the mainspring, however, indicates that local infiltration, local flow system have no direct, immediate effect on the spring's parameters. This is also in accordance with the lowest measured isotopic values of the spring and the wells in Kistapolca, as it is discussed above. This suggests the influence of higher order flow systems rather than the local recharge conditions in the wells and springs of Kistapolca, which is also supported by the numerical simulation. Similarly, Malík et al. (2019) found that the response time of spring discharge to sudden groundwater table rise occurred later at springs with a higher water temperature. However, to draw more detailed interpretation about the hydraulic response of the spring and the thick karst reservoir to precipitation events, monitoring data of several years are needed, which is currently not available, but the monitoring is in progress at the mainspring and in the Beremend Cave. Moreover, it should be noted that some wells show more positive isotopic values (δD: from − 73.8 to − 65.4‰ and δ 18 O: from − 10.3 to − 9.2‰) than the above-mentioned wells and spring, indicating that in some cases, the influence of local, recent precipitation-related water is more pronounced. The isotopic composition of water samples from the Beremend Cave (Nr. 1) shows the most positive values plotting (Fig. 8) closest to those that is observed in the recent precipitation implying lack or very limited effect of higher order flow systems. The measured δ 18 O value of the water in the cave is in accordance with the measurements (δ 18 O: − 8.8 ‰) of Vigassy et al. (2010). The lowest electrical conductivity value (425 µS/cm, Table 2) together with the isotopic composition closest to the precipitation (Fig. 8) in the water of Beremend Cave support the scenario showing the effects of recharge and local flow conditions. The numerical simulation is also in good agreement with these observations. Similar results were reported from the hypogene caves of Black Hills by Bakalowicz et al (1987), showing the influence of the local infiltrating water recently. It can be concluded that the karst water present in the Beremend Cave nowadays does not show the effect of deep flow systems and thermal waters implying that significant mixing of deep and local waters might Fig. 9 Continuously recorded water level change, electric conductivity, and temperature data of the mainspring (lines) and the daily precipitations in Nagyharsány (bars) already be inactive as today. However, further detailed study, in the Beremend Cave and surroundings, is needed.
The small scale, 5-10 cm water level changes recorded at the mainspring in Kistapolca might be the effect of the water abstraction of the nearby cement factory and the waterworks in the village, which have together yearly about 120 -145 − 190 × 10 6 m 3 abstracted water amount in the last 3 years, showing an increasing trend (data provided by the General Directorate of Water). Therefore, the groundwater-dependent ecosystem of the lake fed by the spring can rather suffer from this anthropogenic effect. Climate change, i.e., reduced recharge can further strengthen this effect, as large groundwater abstraction superimposed on regional water level decrease may result in longer and more common dry periods of the lake, as it was theoretically investigated by Havril et al. (2018).

Conclusion
Several methods were applied in the hydrogeological characterization of the surroundings of Kistapolca, where the only naturally outflowing springs of the region can be found. The combination of geophysical, geochemical, hydraulic methods enabled to deliver explanations to the existence of the springs, and multiple pieces of evidence that the spring is fed by a higher order-possibly intermediate-flow system without or limited influence of local flow system. These findings are also supported by numerical simulation. In addition, the vulnerability of the groundwater-dependent lake ecosystem was evaluated as well. The results related to the Beremend Cave indicate that mixing corrosion is not an active process today, the Fig. 10 Numerical simulation results of groundwater flow and heat transport. a Model geometry indicating parameters of siliciclastic confining bed and carbonate: K-hydraulic conductivity, n-porosity, λthermal conductivity. b The intensity of groundwater flow with the normalized flow vectors and boundary conditions. c Temperature distribution and the applied boundary conditions of heat transport