Introduction

Hungary has favorable geothermal conditions related to the formation of the Pannonian Basin. During the early–middle Miocene basin formation, the lithosphere was thinned out due to extension, resulting in an elevated geothermal gradient in the lithosphere (45 °C km−1), and thus high surface heat flow (Royden et al. 1983; Bada et al. 2007; Horváth et al. 2012). Heat flow in the Pannonian Basin varies between 50 and 130 mW m−2, with a mean value of 100 mW m−2 (Dövényi et al. 2002; Horváth et al. 2015), which is significantly higher than the continental average of 65 mW m−2 (Pollack et al. 1993). Since the basement and the basin-fill sequences can also serve as geothermal reservoirs, groundwater with outflow temperature above 30 °C is available in 70% of the country (Mádl-Szőnyi 2006).

The utilization of thermal water has a long tradition in the Pannonian Basin and Hungary (Szanyi et al. 2009; Szanyi and Kovács 2010; Rman et al. 2015; Rotár-Szalkai et al. 2017; Szőcs et al. 2018). Thermal water is widely used for balneological and heating purposes, especially in the agriculture sector, representing the major part of direct use (Nádor et al. 2019). The study area in South Transdanubia (Southwestern Hungary, Fig. 1) is the cradle of thermal water abstraction, since Hungary’s first thermal well was constructed in this region, in Harkány, by Vilmos Zsigmondy in 1866 (Zsigmondy 1873). Currently there are about 50 wells yielding thermal water above 30 °C which is utilized by spas, agriculture and industry (Fig. 1). Thermal water is exploited for balneological purposes without reinjection and only very few heating systems use this technique (e.g., in Bóly in the studied region). This limits the amount of heat and water that can be utilized from a reservoir causing nonsustainable production (Rybach and Mongillo 2006). Unfavourable effects of excessive thermal water production (e.g., significant drop of hydraulic head, temperature, yield, changes in water chemistry) have already been reported in various parts of the Pannonian Basin (Szanyi and Kovács 2010; Tóth et al. 2016; Rotár-Szalkai and Ó Kovács 2016). The sustainable utilization of the thermal waters and the long-term productivity of the reservoirs require the clarification of the path(s) of the groundwater flow as well as the driving forces (Mádl-Szőnyi and Simon 2016). Reinjection and further production possibilities of thermal water can be evaluated and planned only if the regional hydrodynamics are known.

Fig. 1
figure 1

Location of the study area in Southwestern Hungary (Europe)

Full size image

Besides thermal water, considerable hydrocarbon resources are also connected to the sedimentary basin-fill of the Pannonian Basin. In South Transdanubia there are also several hydrocarbon fields connected to the Drava Basin (Kőrössy 1989; Saftić et al. 2003; Fig. 1). Basin-scale groundwater flow systems affect the petroleum systems as a result of geologic agency of groundwater flow (Tóth 1980, 1999). Despite this fact, traditionally stratigraphical, tectonical, lithological traps and their combination are only considered. However, hydraulic traps could also have an important role (Tóth 1988). In the deeper sub-basins of the Pannonian Basin a so-called overpressured flow system can be observed below the gravitational flow systems (Tóth and Almási 2001). The overpressure is the result of vertical compaction and lateral tectonic compression, and usually drives fluids in an upward direction. Where overpressured upward flow meets with gravitational downward flow of recharge areas, regional minimum fluid-potential zones evolve, which can serve as hydraulic traps for hydrocarbons (Tóth 1988). Therefore, regional groundwater flow system analysis could help hydrocarbon exploration as well, as was demonstrated by Verweij and Simmelink (2002) and Verweij et al. (2012), and also successfully applied in the eastern part of the Pannonian Basin (Czauner and Mádl-Szőnyi 2013).

Marginal areas of sedimentary basins and outcropping carbonates are usually characterized by natural thermal water discharge (Erőss et al. 2008; Goldscheider et al. 2010; Erőss et al. 2012; Mádl-Szőnyi and Tóth 2015; Erhardt et al. 2017). These marginal, tectonically controlled discharge areas are favourable sites for the development of hypogene caves, where the morphology clearly reflects the effects of upward flowing waters and the minerals show the effect of thermal waters (Klimchouk 2007, 2012; Erőss et al. 2012; Mádl-Szőnyi et al. 2017). There are also several caves connected to uplifted basement carbonate blocks in the foreland of the Villány Hills or to the margin of the hills (Fig. 1), where the morphology and the minerals decorating the cave walls indicate the effect of thermal waters (Dezső et al. 2004; Rónaki 2000; Takács-Bolner 1985; Takács-Bolner 2003a; Takács-Bolner 2003b; Vigassy et al. 2010). Since groundwater is considered as a geologic agent (Tóth 1999), karst phenomena including caves, springs, and the karst system itself can be interpreted as manifestations of flowing groundwater. Therefore, the understanding of groundwater flow systems is crucial in karst research and in clarifying the origin and formation of the caves.

The thermal water and hydrocarbon resources, as well as karst phenomena including the thermal-water-related caves, were hitherto investigated separately in South Transdanubia. However, all can be considered as groundwater-related phenomena, and they can be evaluated only if their context is understood, i.e. if their common cause is revealed: the pattern of groundwater flow and its thermal and geochemical characteristics. Basin-scale groundwater hydraulics supported by hydrogeochemistry can help to (1) delineate the flow systems and flow regimes (discharge and recharge zones), (2) determine the origin and driving forces of different fluid components, (3) evaluate the role of structural elements in fluid flow, and (4) identify and characterize recent rock-water interactions and processes, such as cave formation. In contrast with the previous traditionally applied aquifer-focused hydrogeological approaches, this basin-scale hydrodynamic approach treats the flow field (including aquifers and aquitards) as a whole, using the concept of regional hydraulic continuity (Tóth 1995). During the interpretation of fluid potential anomalies, the effect of hydrostratigraphy or geology (permeability) can be directly concluded. This approach and the hydraulic analysis of preproduction archival measured data were first applied in the area by the present study; however, there are similar studies in marginal carbonate and sedimentary basin regions of Hungary where this approach helped in the evaluation of resources and phenomena (Czauner and Mádl-Szőnyi 2013; Erhardt et al. 2017; Mádl-Szőnyi et al. 2018; Mádl-Szőnyi et al. 2019).

The aim of the study is to (1) perform a regional (i.e., basin-) scale groundwater flow system analysis based on measured archival hydraulic data, (2) derive conclusions related to the thermal water and hydrocarbon resources, and facilitate their further exploration and sustainable use, (3) based on the regional flow system context, determine the hydrodynamic conditions responsible for local-scale groundwater-related phenomena such as geochemical characteristics of groundwater and cave formation, focusing particularly on the southern foreland of the Villány Hills (sfVH), and (4) compile a generalized flow field model for this marginal karst area. Additionally, through comparison with the marginal Buda Thermal Karst area (Hungary) generalized conclusions regarding the connections between marginal karst reservoirs and the Pannonian Basin can be drawn, which can be considered in similar geological settings worldwide.

Study area

The 4571-km2 study area is located in the southwestern part of the Pannonian Basin, Hungary, including the Hungarian part of the Drava Basin, the southern foreland of the Mecsek Mountains, and the Villány Hills (Fig. 1). In the geodetic interpretations, the Hungarian ‘EOV’ coordinate system (Uniform National Projection system) was applied, which uses the horizontal YEOV and vertical XEOV axes to express the easting and northing distances in meters, from the system’s origin. Elevations are measured on the vertical ZEOV axis in meters as well and referenced to the Baltic Sea level.

The basement of the study area is built up by fractured and weathered Palaeozoic low- and medium-grade metamorphics (gneiss, mica schists), siliciclastic terrestrial sediments and granite, which can be found mainly in the western and northern part of the study area (Szederkényi 2001; Haas et al. 2010; Szederkényi et al. 2013; Fig. 2). Some of these formations are subaerially exposed in the Mecsek Mountains. Intensely karstified Mesozoic carbonates mainly form the eastern and southeastern basement part (Török 1998; Császár 2002). These deposits were thrusted upon each other forming the NE-vergent imbricated thrust sheet structure of the Villány Hills (Fig. 2). The carbonate bodies also form uplifted basement blocks in its southern foreland (Császár 2002; Petrik 2009; Konrád et al. 2010). The top of the basement gradually deepens towards the southwest from the area between the Mecsek Mountains and the Villány Hills and reaches an elevation of z = −4,500 m above sea level (asl) at the state border (Haas et al. 2010; Haas and Budai 2014; Fig. 2).

Fig. 2
figure 2

Geological map of the Pre-Neogene basement in the study area (modified after Haas et al. 2010)

Full size image

In most of the study area, the Paleozoic and Mesozoic formations are covered by Cenozoic sediments of the Drava Basin, accumulated mostly during the middle and late Miocene (Saftić et al. 2003). Development of the Drava Basin is related to the formation of the Pannonian Basin, which started in the early Miocene. Alpine-Carpathian-Dinaric orogenesis-related formation of the Pannonian Basin can be divided into syn-rift and post-rift phases. During the early Miocene syn-rift phase, development of the basin was characterized by rift-related extensional processes, due to the formation of a back-arc basin caused by subduction along the Carpathian arc. Since the extensional deformation was not uniform in the basin (neither in style nor in scale), a number of sub-basins were formed (Balázs et al. 2016). These processes led to the generation of the NW–SE trending Drava Basin from the end of the early Miocene. In this tectonic period, pull-apart basins were opened in the region due to the strike-slip faults which resulted in the deepening of the basin (Juhász et al. 2007; Horváth 2007; Nagymarosy and Hámor 2013). During this syn-rift period mainly silt, sand and conglomerate deposited (Saftić et al. 2003). The post-rift phase of the basin started in the late Miocene, and this phase was tectonically relaxed and controlled by thermal cooling (Royden and Horváth 1988). During the beginning of the late Miocene, the Carpathian-Pannonian region was isolated from the seas, thus the brackish Lake Pannon formed (Magyar et al. 1999). The deep Drava Basin was eventually filled up by a river system from NW to SE (Uhrin 2011; Magyar et al. 2013). As the shelf edge approached the basin, the sediment input rate exceeded the subsidence rate of the basin and the sediments gradually filled up the basin. The basin fill is composed of profundal marls, turbiditic sandstones, heterolithic slope deposits and stacked deltaic bodies (Saftić et al. 2003). The thickness of these Pannonian sediments can reach up to 3,000 m at the SW part of the study area; however, in the eastern part they are only several hundred meters thick (Saftić et al. 2003).

In the last phase of basin development, the extension was followed by compression (inversion phase). Uplift and exhumation of the Mecsek Mountains and Villány Hills was driven by late Miocene to post-Miocene basin inversion. The Pliocene was dominated by fluvial sedimentation in the Drava Basin, but terrestrial, palustrine sedimentation also took place.

Methods

Hydrostratigraphy

Hydrostratigraphic classification was defined based on literature data. Since the sedimentary fill of the Drava Basin shows considerable similarities to the Great Hungarian Plain, the hydrostratigraphic units of the Great Hungarian Plain defined by Tóth and Almási (2001) and Czauner and Mádl-Szőnyi (2013) were used. The other lithostratigraphic units in the southern foreland of the Mecsek Mountains and the Villány Hills were determined based on reference values of Freeze and Cherry (1979), Brassington (2017), Erhardt et al. (2017) and Tóth (2018) (Fig. 3).

Fig. 3
figure 3

Hydrostratigraphic units of the study area. Legend: GP-AF – Great Plain Aquifer; A-AT – Algyő Aquitard; SZ-AF – Szolnok Aquifer; E-AT – Endrőd Aquifer; PrePA-AF – Undifferentiated Pre-Pannonian Siliciclastic Aquifer; Undifferentiated PreNg-AF – Pre-Neogene Aquifer

Full size image

Since the main goal of the basin-scale archival hydraulic data analysis (‘basin hydraulics’) presented hereinafter is to reveal the fluid-potential field with the derived hydraulic gradients and flow directions (Tóth 2009), the hydrostratigraphy (hydraulic conductivity, K, values) was not explicitly involved in this analysis. However, during the interpretation of fluid-potential anomalies, one can directly conclude as to the effect of hydrostratigraphy, for instance elevated hydraulic gradients often refer to low-permeability zones or barrier faults (Tóth 2003; Czauner and Mádl-Szőnyi 2011; Mádl-Szőnyi et al. 2019). Consequently, hydrostratigraphic units were represented on the pressure-elevation [p(z)] profiles and hydraulic cross-sections in order to allow for determining their hydraulic function, while the resolution always depended on data availability. For instance, the Neogene succession could not be divided into hydrostratigraphic units due to the lack of data in some p(z) profiles and cross-sections. These are referred to as the undifferentiated Neogene Aquifer-Aquitard Group (Ng-AF-AT) in the figures. The effectiveness of this hydrostratigraphic simplification, which enhances larger-scale heterogeneities such as the higher regional hydraulic conductivity of the carbonate basement (PreNg-AF) and the generally lower regional hydraulic conductivity of the siliciclastic confining layers (Ng-AF-AT) in the southern foreland of the Villány Hills (sfVH), was also demonstrated in other parts of the Pannonian Basin by e.g., Czauner and Mádl-Szőnyi (2013), Erhardt et al. (2017), Mádl-Szőnyi et al. (2019).

Hydraulic data evaluation

Archival measured hydraulic data of wells were analysed by the basic regional-scale methods of ‘basin hydraulics’ (Tóth 2009). Measured hydraulic data were taken from the original documentation of hydrocarbon and water wells in the archives of government institutions (Mining and Geological Survey of Hungary, General Directorate of Water Management). These hydraulic data were measured during or right after well completion, thus they represent the preproduction hydraulic conditions of the aquifers. Hydraulic data of water wells were registered within the well documentation usually without methodological descriptions and determination of precision. Hydraulic data of hydrocarbon wells are formation pressure values which were measured by several methods. During data culling, the methodology of this study took into consideration the reliability of these measurements and used only the most reliable drill stem test (DST) results. Accuracy of these tests was also not registered in the well documentation. However, errors of measurements in a basin-scale study like the present one can be neglected, particularly compared to the methodological problems which were culled during data preparation.

During the calculations of hydraulic heads and pore pressures, one average groundwater density value had to be determined for the whole regional-scale study area, namely 1,000 kg m–3. This is a good approximation, as the effects of salinity and temperature on groundwater density roughly compensate each other in the study area, while the pressure and dissolved gas content have a less significant effect, and the role of free gas content is filtered out by using hydraulic parameters measured in water (Czauner and Mádl-Szőnyi 2013). The database is included as Table S1 of the electronic supplementary material (ESM).

Pressure-elevation profile

Based on the pressure data, pressure vs. elevation [p(z)] profiles were constructed, which sensu stricto allow for the examination of the vertical component of fluid flow direction by comparing the vertical pressure gradients 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 hydrostatic conditions there is no vertical fluid flow, because the pressure potential and the force of gravity on water particles are in balance. In this case, the vertical pressure gradient is hydrostatic, and its value is γst = 9.81 MPa km–1, by applying the water density of 1,000 kg m–3. In a hydrodynamic state, a driving force can shift the system off balance, therefore a vertical flow component also occurs, 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 superhydrostatic. 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.

Consequently, in the p(z) profiles, beyond the data points, the (reference) hydrostatic vertical pressure gradient line (γst = 9.81 MPa km−1), the dynamic vertical pressure gradients (γdyn) fitted on the data points, and the available hydrostratigraphic information are also presented. The ‘margin of error’ of the hydrostatic vertical pressure gradient was taken as ±0.5 MPa km−1 (i.e., the vertical pressure gradient was interpreted as hydrostatic between the values of 9.81 ± 0.5 MPa km−1) based on the results of Czauner and Mádl-Szőnyi (2013) for the Great Hungarian Plain of similar geological build-up. During the interpretation, this margin proved to be applicable in the present study area as well, thus provided reasonable results, which coincide with the interpretation of tomographic fluid-potential maps and hydraulic cross-sections.

Theoretically, only data derived from one well can be interpreted in a p(z) profile. Since most of the wells in the study area provide only one hydraulic datapoint, application of p(z) profiles had to be extended for areas based on two criteria. Namely, wells interpreted in one p(z) profile had to have approximately the same land surface elevation, because this affects the position of the reference hydrostatic pressure gradient line. On the other hand, hydrostratigraphic build-up regarding the units and their thicknesses also had to be approximately concordant, since otherwise phenomena caused by hydrostratigraphic changes could not be jointly interpreted. Beyond these criteria, the areal and depth distribution of hydraulic data were also taken into account. As a result, 44 p(z) profiles were compiled in the study area (Fig. 4).

Fig. 4
figure 4

a Spatial distribution of wells, areas of p(z) profiles, and trace lines of the hydraulic cross-sections, b Distribution of the number of hydraulic data points in relation to the elevation of their measuring points

Full size image

Tomographic fluid-potential map

By definition, a potential or potentiometric map can be constructed for a horizontal and confined aquifer; however, in areas characterized by low data density and/or with dipping and/or faulted strata, a potential map is prepared for a given depth (elevation) interval. In this way, the effects of heterogeneities in hydraulic conductivity can be interpreted, for instance, along strata boundaries and faults. The so-called tomographic fluid-potential maps are constructed for consecutive elevation intervals, thus also allowing for the interpretation of vertical flow directions by comparing the successive maps. Finally, it is worth emphasizing that in the interpreted potential maps, groundwater flow direction (often indicated by arrows) only applies to the horizontal component of the fluid flow driving force, while the flow intensity also depends on the hydraulic conductivity of the flow medium.

Elevation intervals were chosen based on the depth distribution of the data (Fig. 4), as well as on the initial interpretation of the p(z) profiles. As a result, four tomographic fluid-potential maps were made up roughly for the whole study area and for the following elevation intervals: z1 = 300–0 masl, z2 = 0 to −100 masl, z3 = −100 to −220 masl, z4 = −220 to −1,000 masl. Below this depth, data were available from restricted areas within the Drava Basin resulting in only two deeper maps for the following elevation intervals: z5 = −1,000 to −2,000 masl, z6 = −2,000 to −3,700 masl.

Hydraulic cross-section

Traces of the hydraulic cross-sections were laid out by taking into account the data distribution (Fig. 4), as well as the location of Harkány (at profile 26 in Fig. 4), being the focus of local-scale interest. Subsequently, data points were projected into the sections, always perpendicular to the trace line and from within the distance of 3 km of section A and within 1 km of section B (section A starts at profile 9 and ends at profile 25, and section B start at profile 25 and ends near profile 28, see Fig. 2 and subsequent figures). As well as the hydraulic head values of the wells, the intersections of equipotential lines contoured in the potential maps were also represented in the cross-sections as points referring to the middle of the given elevation interval. These data support the contouring along the sections in the zones that lack data.

Data distribution and quantity did not allow the contouring of h values (masl) along the whole sections, while linear interpolation was not possible since distribution of potentials (i.e., density of equipotentials) depends on several factors and it is usually not uniform. Additionally, contouring (i.e. drawing of equipotentials) was made by hand in order to avoid the errors possibly caused by computational interpolation, as well as to provide the chance for applying personal judgement and local knowledge of the given study area. Thus, the completed interpretation is rather conservative, but also more reliable. In order to emphasize the vertical fluid flow directions, vertical exaggeration was applied in the cross-sections (height:width = 4:1 in section A, 2:1 in section B). Furthermore, it is worth emphasizing that arrows in the interpreted sections only show the direction of fluid-flow driving force, while the flow intensity also depends on the hydraulic conductivity of the flow medium. Finally, the hydrostratigraphy and the major fault zones were also demonstrated in the sections in order to support the interpretation.

Hydrochemical data evaluation

One of the main objectives of the present study was to determine the hydrodynamic conditions that are responsible for local-scale groundwater-related phenomena such as geochemical characteristics of groundwater and cave formation, focusing particularly on the southern foreland of the Villány Hills (sfVH). Accordingly, after the regional (i.e., basin-) scale hydraulic data evaluation, targeted hydrochemical data analyses were also carried out in order to support the evaluation of hydrodynamic conditions in the sfVH where hydraulic data were available only from above z = −200 masl. Namely, the aim was to distinguish fresh karst waters and saline basinal fluids, and consequently to confirm the presumed appearance of deep basinal fluids in the sfVH.

Considering the results of hydraulic data evaluation, archival preproduction and recently measured hydrochemical data were also analysed for the areas of the presented p(z) profiles (Table S2 of the ESM). Archival chemical data of the wells were taken from the original well documentation of the Mining and Geological Survey of Hungary, and from the database of the General Directorate of Water Management. The parameters of the elemental analyses of these documentations were measured by accredited laboratories in accordance with the relevant standard water methods. Besides this, selected wells were also sampled during the present study. Among the main ions (accuracy indicated in brackets), Ca2+ (±5 mg L−1), Mg2+(±2 mg L−1), HCO3− (±12 mg L−1) were analyzed by acid-base titrimetry, Cl− (±2 mg L−1) by complexometric titration, Na+(±2 mg L−1) and K+ (±5 mg L−1) were detected by flame photometry, and SO42− (±5 mg L−1) by photometry at the Department of Geology at Eötvös Loránd University in accordance with the relevant standard water methods. The water sample results were investigated using a Piper diagram and classified by the hydrochemical facies concept of Back (1966).

Results

Pressure-elevation profiles

Within the study area, 44 p(z) profiles were constructed, while in this report 10 representative profiles are displayed from different parts of the area. Data were available from the whole study area down to about z = −500 masl. In this zone, in most cases, the vertical pressure gradient values are close to the hydrostatic value (γ = 9.3–10.3 MPa km−1), indicating lateral flow conditions (Fig. 5). Superhydrostatic pressure gradients (γ > 10.3 MPa km−1) prevail in the southern foreland of the Villány Hills (profiles 24, 26, 28, 29 and 31 in Fig. 5) and near the Danube River (profiles 39 and 40 in Fig. 5), which refer to upward vertical flow direction. Subhydrostatic vertical pressure gradients (γ < 9.3 MPa km−1) and downward vertical flow conditions are typical for areas of elevated topography (profiles 2, 3, 4, 6, 15, 16, 18, 35 and 43 in Fig. 5). Data from deeper levels were available only from the western half of the study area and in p(z) profile 44 (Fig. 5). Within this area, in the elevation interval of about z = −500 to 1,500 masl, horizontal and upward flows become dominant.

Fig. 5
figure 5

Vertical flow conditions in the areas of p(z) profiles from the land surface down to z = −500 masl

Full size image

In the following, p(z) profiles are presented as representative examples for the above described situations. On p(z) profile 40, along the Danube, based on the vertical pressure gradient (γ = 11.49 MPa km−1), upward flow conditions could be determined (Fig. 6a). Data points of this profile represent a shallow elevation range (z = 56–73 masl) within the Neogene Aquifer-Aquitard group. On the other hand, p(z) profile 44 displays a hydrostatic vertical pressure gradient (γ = 9.66 MPa km−1), thus lateral flow conditions dominate down to about z = −1,500 masl between the Mecsek Mountains and the Villány Hills (Fig. 6b).

Fig. 6
figure 6

a–j Representative p(z) profiles from the study area

Full size image

The p(z) profiles from the southern foreland of the Villány Hills can be separated into two groups. In the first group the vertical pressure gradient is higher than hydrostatic (profiles 28, 29, 26 in Fig. 6). The vertical pressure gradient of p(z) profile 28 shows the upper limit of hydrostatic conditions (γ = 10.3 MPa km−1) (Fig. 6c), whereas for profile 29 upward flow conditions could be determined (γ = 10.45 MPa km−1) (Fig. 6d). The p(z) profile 26 is located in the town of Harkány. Data points of the profile show a slightly higher pressure gradient (γ = 10.39 MPa km−1) than the hydrostatic one (γ = 9.81 MPa km−1) (Fig. 6e). In the second group of profiles (profiles 27, 30, 32 in Fig. 6), hydrostatic vertical pressure gradients, thus lateral flow conditions could be determined in every case. The p(z) profiles 30, 27 and 32 represent γ = 9.9 MPa km−1, γ = 9.82 MPa km−1, and γ = 9.66 MPa km−1 vertical pressure gradients, respectively (Fig. 6f–h, respectively). However, in p(z) profile 30, with the addition of one deep data point from the Pre-Neogene Aquifer, the gradient becomes higher than hydrostatic (γ = 9.81 MPa km−1) with γ = 10.64 MPa km−1 referring to upward flow conditions (Fig. 6f).

Data were available from deeper than z = −1,500 masl in p(z) profiles 2, 3, 6, 8, 11 and 13 (Fig. 5). Most of these data represent overpressured conditions and upward vertical flow directions. For example, in p(z) profile 2 (Fig. 6i) the two deeper data points (below z = −1,000 masl) show a superhydrostatic vertical pressure gradient (γ = 15.09 MPa km−1), thus refer to upward flow, which results in the dissipation of 31% overpressure from the Pre-Neogene basement to the top of the Szolnok Aquifer. In p(z) profile 6 (Fig. 6j), the deeper data points (below z = −2,000 masl) show a very significant superhydrostatic vertical pressure gradient (γ = 41.58 MPa km−1; 9–26% overpressure), thus upward flow conditions in the Pre-Neogene Aquifer. Additionally, it is worth emphasizing that shallow data of both profiles (and profile 3 as well) show downward flow conditions within the Great Plain Aquifer (profile 2: γ = 8.85 MPa km−1, profile 6: γ = 8.89 MPa km−1).

Tomographic fluid-potential maps

According to the shallowest maps (M1, M2, M3, M4; Fig. 7), which cover roughly the whole study area, groundwater flow directions reflect the topography since fluids flow from areas of higher to lower elevations. As a result, the direction of regional horizontal flow is southward, while in the eastern part of the area it is S–SE.

Fig. 7
figure 7

Tomographic fluid-potential maps M1–M5. Background: topographical map

Full size image

Due to the data distribution, only the Drava Basin part of the study area could be examined on maps M5 and M6. On map M5 (Fig. 7), groundwater flow directions still reflect the topography and the regional fluid flow direction tends southward. The data distribution of fluid-potential map M6 is limited to a small area within the Drava Basin (Fig. 8). Here groundwater flow directions do not reflect the topography. However, considering a Pre-Neogene basement map (Fig. 8), flow directions tend from a deeper sub-basin towards its margins, roughly towards the northeast, east and southeast. Furthermore, a significant increase (300–2,000 m) can be observed in hydraulic head values as compared to the previous map.

Fig. 8
figure 8

Tomographic fluid-potential map M6. Background: Pre-Neogene basement

Full size image

Regarding vertical flow directions, by comparing the consecutive maps, they represent and extend the results of p(z) profiles to the whole study area down to z = −220 masl. Below this elevation, upward flow can also be detected beneath the shallow recharge areas, whereas below z = −2,000 masl upward flow can be observed.

Hydraulic cross-sections

Hydraulic cross-section A is NW–SE oriented and runs across the Drava Basin (Fig. 4). The upper (shallow) part of the section is dominated by hydrostatic conditions, thus lateral flow directions (Fig. 9). In the eastern part of the section, eastward lateral flow can be observed with some vertically upward flow components. However, around the western termination of the section where the Pre-Neogene basement deepens, in the deeper parts of the Great Plain Aquifer, upward flow conditions can be observed as well. The Pre-Neogene aquifer is represented by only one deeper data point along this section (below z = −3,000 masl), which exceeds h = 2,000 masl.

Fig. 9
figure 9

Hydraulic cross-section A

Full size image

Cross-section B is SSW–NNE oriented and located between the Drava Basin and Villány Hills (Fig. 4). Most of the data are above z = –200 masl elevation; therefore, just the shallow part of the section can be interpreted. In this part of the section, nearby hydrostatic conditions prevail with a regional southward lateral flow direction. In addition, a positive fluid-potential anomaly (a so-called potential plume) can be observed around Harkány along fault zones in the uplifted Mesozoic carbonates of the Pre-Neogene Aquifer. The thermal Spa of Harkány produces its thermal water from this zone where a swamp existed before the start of thermal water production (Zsigmondy 1873). As a further result of this positive anomaly, north- and southward flows evolved. The former meets with southward flows and thus composes a local discharge area between Harkány and Máriagyüd. Additionally, the upward flow around Máriagyüd, in the north of this local discharge area, could be related to the local discharge of the waters infiltrating in the Villány Hills. Moreover, there were two cold springs (10 °C) in Máriagyüd (Vadász 1949; Kessler 1959), which could be interpreted as discharge phenomena of local flow systems.

Hydrochemical analysis

The geochemical characterization of groundwater in the study area is illustrated with the help of a Piper diagram (Fig. 11). Most of the wells have similar chemical characteristics, with calcium, magnesium and bicarbonate dominance. As compared with the most characteristic facies, the deep wells of p(z) profiles 26 (Harkány) and 44 show Ca-Na, HCO3-Cl-SO4 facies. In the case of these p(z) profile areas, the dominant cations are sodium and calcium, whereas most of the other wells are characterized by only calcium dominance. Regarding anions, although bicarbonate is still dominant, these wells are separate from the others and shifted towards chloride indicating higher chloride content (Table S2 of the ESM). The deep basin profiles 2 and 6 show Na-K, Cl-SO4-HCO3 and Na-K, Cl-SO4 facies with Na + K and Cl dominance, which are characteristically different from the other wells.

Conclusions and discussion

Regional groundwater flow systems and related phenomena

The presented regional-scale hydraulic assessment of the study area applied archival preproduction measured data and the concept of regional hydraulic continuity. Based on the combined results of p(z) profiles, tomographic fluid-potential maps and hydraulic cross-section analyses, the following conclusions could be drawn. In the shallower regions of the study area down to about z = −500 masl, vertical flow directions usually reflect the topography. Namely, downward flow conditions can be observed in the topographically higher areas, and horizontal and upward flow conditions in the topographically lower areas. Consequently, gravitational flow systems could be identified, including recharge (downward flow), midline (horizontal flow) and discharge (upward flow) areas (Fig. 5). As suggested by the available data, the penetration depth of downward flow (in recharge areas) is around z = −400 to −500 masl. Below this elevation, a dominantly lateral flow with some vertically upward flow component could be interpreted beneath the shallow recharge and discharge areas as well around z = −500 to −1,000 masl. These upward flows are still gravitational flows considering the magnitude of h values and the lack of overpressure. On the whole, horizontal flow conditions are dominant within the gravitational flow system of the study area, while the regional flow direction reflects the topography and tends toward S–SE, i.e. from the mountainous and hilly regions towards the regional base levels of the Drava and Danube rivers. However, it is worth mentioning that discharge areas cannot be found along the Drava River, which is the opposite case to the Danube River (Fig. 5). On the other hand, in the deeper regions of the study area, below about z = −1,500 to −2,000 masl, which are represented by data only in the Drava Basin, an overpressured flow system could be identified with significant overpressures (max. 66%) and vertical pressure gradients (max. 41.58 MPa km−1). Its vertical flow component is dominantly upward, whilst fluids are driven laterally from the deeper sub-basins towards their margins, i.e. northeast, east and southeast, which do not reflect the topography. A transition zone of the gravitational and overpressured flow systems could be determined between about z = −1,500 and –2,000 masl where the overpressure seems to be totally dissipated.

In the eastern half of the study area, data were available only from shallow depths (above z = −200 masl) where gravity-driven upward flows could be identified, particularly in the southern foreland of Villány Hills. However, existence of the overpressured flow system is strongly probable below about z = −2,000 masl in the eastern half of the study area as well, particularly in parts of the sedimentary basin, i.e. in the southeastern continuation of the Drava Basin (in the south of Villány Hills) and in the southwestern corner of the Great Hungarian Plain (in the east of Mecsek Mountains and Villány Hills). This theory can be supported by the proven presence of overpressured flow systems within the northwestern part of the Drava Basin (described by the present report) and in the Great Hungarian Plain’s sub-basins (Tóth and Almási 2001; Mádl-Szőnyi and Tóth 2009; Czauner and Mádl-Szőnyi 2013; etc.) based on the geological similarity of these areas.

Based on the regional-scale evaluation of fluid flow systems, considerable conclusions could be drawn regarding the geothermal and hydrocarbon potential of the area. From a hydrodynamic point of view, thermal water or geothermal energy production need less energy investment in the zones of lateral and particularly upward flows. Resources of the overpressured flow system are nonrenewable and can be produced only with reinjection at high pressure. However, production from gravity-driven flow systems can be sustainable, though long-term productivity usually requires reinjection even into this system with less energy input in the zones of lateral and particularly downward flows (Mádl-Szőnyi and Simon 2016). Based on the present study, favourable areas and depth zones can be determined for production and reinjection as well. The p(z) profile 44 (Fig. 6b) demonstrates the pressure conditions of Bóly town where a geothermal cascade system based on a 1,500-m deep production well and an injection well in the same reservoir supplies the district heating system. Pressure data represent through-flow conditions (i.e., midline area) down to z = −1,500 masl, which can be explained by the topographically transient position of the area between the Mecsek Mountains. (as recharge area) and the lowlands of the Danube and Drava rivers (as discharge areas). On the other hand, the lateral flow conditions allow thermal water production with relatively high yield and reinjection with low pressure. Based on the example of Bóly, areas of through-flow conditions are hydraulically favorable for the installation of neighboring production and injection wells, for instance the areas of p(z) profiles 9, 10, 12, and 14 where through-flow conditions are proven down to z = −1,100, (−800), −1,600, and − 1,600 masl, respectively.

Regarding the hydrocarbon potential, where overpressured and superimposed gravitational upward flows meet with gravitational downward flow of recharge areas, regional fluid-potential minimum zones evolve, which can serve as hydraulic traps for hydrocarbons (Tóth 1988). In the study area, a hydraulic entrapment zone was identified in the surroundings of p(z) profiles 2 (Fig. 6i), 3, and 6 (Fig. 6j) around z = −500 masl. In other words, this zone represents the upper boundary of vertical hydrocarbon migration, thus entrapment of hydrocarbons is also more probable in these areas at greater depth. Producing hydrocarbon fields are indeed known from this region.

Local-scale phenomena

Based on the evaluation of regional groundwater flow systems, notable local-scale phenomena could also be explained. Particularly the southern foreland of Villány Hills (sfVH) proved to be worthy for further investigation where hydraulic, karstification, and hydrochemical phenomena could be jointly interpreted. Geologically, the Pre-Neogene Aquifer built-up by carbonates is situated in an uncovered (i.e., unconfined) position within the Villány Hills, and in a shallow covered (i.e., confined) position within the southern foreland area, then gradually deepens towards the Drava Basin (see the geological framework in Fig. 10). Hydraulically, this confined or unconfined situation of carbonates proved to be decisive in the evolution of vertical flow directions. During the regional groundwater flow system analysis, gravity-driven lateral and upward flow conditions were identified in the area (Fig. 5). Lateral flows could be determined where the carbonate basement (PreNg-AF) outcrops or has a very thin sedimentary cover (p(z) profiles 27, 30, 32 and 34 on Fig. 6), whilst upward flows could be identified where the siliciclastic sediments cover the carbonates in greater thickness (p(z) profiles 26, 28, 29 and 31 on Fig. 6). Since the siliciclastic cover usually has lower hydraulic conductivity than that of the basement carbonates (PreNg-AF; Fig. 3), it can be concluded that the thicker cover causes greater vertical pressure gradients thus inducing upward flow within the underlying aquifer.

Fig. 10
figure 10

Hydraulic cross-section B

Full size image

Regional groundwater flow system analysis also helps to better understand the cave forming processes. Those caves in the Villány Hills that show the effect of upwelling thermal waters (i.e., hypogene caves) can be found in uncovered carbonates at the margin of the Villány Hills (Nagyharsány, Máriagyüd in Fig. 1) and in uplifted basement blocks (Siklós, Beremend in Fig. 1). According to the hydraulic evaluation, their formation cannot be related to the present-day flow systems since lateral flow conditions dominate within these areas—profiles 27 (Siklós), 30 (Beremend), and 32 (Nagyharsány) in Fig. 6—moreover they are mostly situated above the water table. Formation of the caves is rather connected to the inversion of the Pannonian Basin where the compressional stress field resulted in the uplift of the Mecsek Mountains, Villány Hills and some basement blocks, which were still covered with sediments. As the recharge through the uncovered hilly ranges and carbonates initiated the gravity-driven flow systems, these waters could meet with overpressure-driven upwelling basinal thermal waters at the marginal areas of the Villány Hills and uplifted basement blocks. Fluid mixing in a carbonate rock environment can lead to the “mixing corrosion” phenomenon (Bögli 1964, 1980) resulting in carbonate dissolution, i.e. cave formation, because the waters present in this process are also characterized by different chemical composition (CO2 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. Different ordered flow systems convey waters with different temperatures and geochemical composition to these areas. One of the well-known examples for such a situation is the hypogene karst system in Budapest (Buda Thermal Karst, Hungary) (Erőss et al. 2012, Leél-Őssy 2017; Mádl-Szőnyi et al. 2017). As direct recharge through the uncovered hilly ranges became more and more pronounced, the fraction of basinal fluids decreased, the caves were flooded by “fresh” waters, and the thermal-water-related hypogene cave formation period stopped, also because some of the caves were disconnected from the water table due to ongoing uplift. This scenario is supported by recent isotope measurements from the cave-filling waters of the Beremend Cave, which showed similar values to the present-day precipitation (Erőss et al. 2020).

Furthermore, within the sfVH, there is a special place around Harkány where effects of the presumably existing overpressured flow system, i.e. groundwater flows derived from the Drava Basin, could be identified. The geochemical uniqueness of the thermal waters of Harkány was explained by Lorberer and Rónaki (1978) and Liebe and Lorberer (1981) with fluid contribution from the north-northwest. Later, the conceptual model of Csicsák et al. (2008) already indicated an interplay between the Drava Basin and the carbonate suite of Villány, whilst as indicators of basinal fluids δ18O and δD values, 14C ages, and organic compounds (hydrocarbon traces) were designated in the thermal waters of Harkány. The present study proved both presumed flow directions based on the analyses of measured data. Namely, in p(z) profile 26 (Fig. 6e) an upward flow direction can be observed, whilst on hydraulic cross-section B (Fig. 10) a supposedly fault-related positive potential anomaly (i.e., elevated h values) can be observed. This flow pattern could be determined within the gravitational flow system since data were available only above z = −200 masl; however, support and contribution of the overpressured system can be presumed at least along the fault zone, which approaches the land surface within the elevated carbonate basement (PreNg-AF) block. Further evidence could be that in the sfVH, the Harkány area is located the nearest to the Drava Basin’s margin which forms a bay around Harkány (Fig. 2, z = −1,000 masl basement elevation contour line around p(z) profile 26). The effect of the basin can also be seen in the temperatures of the thermal waters in Harkány (62 °C, see Table S1 of the ESM), which is the highest groundwater temperature within the sfVH (maximum temperature value of the other wells in p(z) profiles 27–30 and 32 is 33 °C (Table S1 of the ESM). In addition, based on the hydrochemical analysis, groundwater of Harkány of Ca-Na, HCO3-Cl-SO4 facies, and their higher chloride and sodium content, represent a transition between the karst waters of Ca-Mg, HCO3 facies and the deep basinal fluids of Na-K, Cl-SO4-HCO3 facies (Fig. 11). In other words, it is a mixed water with a deep basin component as well. Regarding the N–NW origin of groundwater, it can be supported by the flow directions on fluid-potential maps M1–M4 (Fig. 7) down to about z = −1,000 masl. Regarding the karst water component, Villány Hills represent the primary recharge area, but regional-scale flows derived from the Mecsek Mountains could also reach the Harkány area based on the flow directions on fluid-potential maps M1– M4 (Fig. 7). Finally, as a result of the complex analysis of the Harkány area based on measured data, a generalized flow-field model could be established (Fig. 12) based on hydraulic cross-section B (Fig. 10) and related phenomena.

Fig. 11
figure 11

Piper diagram of water samples from the areas of p(z) profiles presented in Fig. 6 (2, 6, 26–30, 32, 40, 44). Full circles: water samples from the Pre-Neogene Aquifer, empty circles: water samples from the Neogene Aquifer-Aquitard group

Full size image
Fig. 12
figure 12

Generalized flow field model of the Harkány area based on the hydraulic analyses of measured data

Full size image

Comparison with the Buda thermal karst and Hungarian Paleogene Basin (Hungary)

The boundary zone of the dominantly unconfined Pre-Neogene carbonates of the Buda Thermal Karst (BTK) (central Hungary) and confined deep carbonates within the basement of the Hungarian Paleogene Basin (HPB) could serve as an analogy for the southern foreland of Villány Hills (sfVH). The BTK area is the focus of research interest because of its thermal water resources and the on-going hypogene karstification processes. The first complex hydraulic evaluation of the area based on measured data was carried out by Erhardt et al. (2017). As a result, gravitational flow systems, hydraulic continuity, and the modifying effects of aquitard units and faults were identified in the karst area similarly to the sfVH. The significant conclusion, according to which the thicker siliciclastic cover causes greater vertical pressure gradients, thus upward flow within the underlying carbonate aquifer, also applies in the BTK and sfVH. The evolution history of the karst area and the caves (uplift-initiated evolution of the topography-driven groundwater flow related to the inversion of the Pannonian Basin) and the demonstrability of the basinal fluid component with hydrocarbon traces, also show similarities in the two areas (Poros et al. 2012; Mádl-Szőnyi et al. 2018). However, a significant difference can also be determined between the confined deep carbonates of the sfVH-Drava Basin system and the BTK-HPB system. Namely, in the former case, normal and overpressured conditions were detected with lateral and upward flows, while in the latter case, normal and underpressured conditions were identified by Mádl-Szőnyi et al. (2019). They explained the underpressure by erosional disequilibrium reflecting a geologically transient flow field, and by the flow-impeding effect of the aquitard units within the siliciclastic confining layers. On the other hand, in the Drava Basin, overpressure could be generated by the interaction of (1) tectonic compression and vertical compaction as fluid-flow driving forces, and (2) the flow-impeding effect of the aquitard units within the siliciclastic confining layers, similarly to the S–SE sub-basins of the Great Hungarian Plain (Tóth and Almási 2001; Czauner and Mádl-Szőnyi 2013). In fact, the effect of tectonic compression on the flow pattern could be the most significant in the Drava Basin compared to the other sub-basins of the Pannonian Basin, since inversion of the basin first started in its southwestern part at the end of the Miocene (Bada et al. 2007). Afterwards, the onset of inversion gradually migrated towards more internal parts of the system resulting in transpression (local shortening) in the western parts of the Great Hungarian Plain, while in the eastern parts transtension (local extension) still dominates. As a consequence of the difference between the pressure conditions, significant distinction can be made regarding the origin of the basinal fluid components of the two areas as well. In case of the BTK-HPB system, underpressure of the confined carbonates drives the dominantly basinal fluids of the confining strata into the underlying carbonate aquifers via vertical leakage. On the other hand, in the sfVH-Drava Basin system, overpressure of the confined carbonates drives the basinal fluids upward and towards the marginal areas of the basin.

As a result, comparison of the study area with the Buda Thermal Karst allows for the definition of the key factors that determine the connection of the marginal carbonate reservoirs with the Pannonian Basin, i.e. (1) tectonic regime and related geographical evolution, and (2) flow-impeding effect of the aquitard units within the siliciclastic confining layers. It can be also concluded that the basin evolution, particularly the inversion, had a crucial role in the initiation of the flow system of these marginal karst reservoirs and their karstification (cave forming) processes.

Summary

The study deduced the regional flow field of South Transdanubia (Hungary) based on archival measured hydraulic data. During the research, 44 p(z) profiles, 5 tomographic fluid-potential maps and 2 hydraulic cross-sections were constructed. Based on these hydraulic interpretations in the shallower regions (down to about z = −500 masl) of the research area gravitational flow systems could be identified with recharge, midline and discharge areas. At the level z = −500 to −1,000 masl, dominantly lateral flow with some vertically upward flow component could be determined. In the deeper regions of the study area (below about z = −1,500 to −2,000 masl), an overpressured flow system could be identified with significant overpressures and vertical pressure gradients. Its vertical flow component is dominantly upward, whilst fluids are driven laterally from the deeper sub-basins towards their margins. These results can be directly applied in the exploration and sustainable utilization of groundwater-related resources, such as thermal waters and hydrocarbons. The presented study emphasizes that knowledge on regional groundwater flow systems is essential to understand local-scale groundwater-related phenomena such as recent cave formation in an area. The caves of the Villány Hills show the effects of the upwelling thermal waters in the uncovered carbonates of the Villány Hills and in covered but uplifted basement blocks. Based on the results of the hydraulic evaluation, their formation cannot be related to the present-day flow systems, since lateral flow conditions dominate within these areas and the caves are mostly situated above the recent water table. Moreover, the study confirmed that in marginal karst areas the evolution of the adjacent sedimentary basin has a crucial role in the evolution of the karst system itself.