Regional groundwater flow and karst evolution–theoretical approach and example from Switzerland

In regional scale aquifers in the Rhine Valley and Tabular Jura east of Basel (Switzerland), the groundwater circulation was investigated using regional-scale geological and hydraulic 3D models. The main aquifers in the area comprise the Quaternary aquifer of unconsolidated gravel deposits along the River Rhine and its tributaries, as well as the regional scale karst aquifer within the Upper Muschelkalk. Land subsidence, a process likely associated with salt solution mining, indicates further subordinate groundwater bearing segments and complex groundwater interactions along fault zones. In the aquifer systems we investigated, regional-scale groundwater circulation was simulated and visualized in relation to the geological settings. Lithostratigraphic units and fault structures were parameterized and analyzed, including the sensitivity of hydraulic properties and boundaries. Scenario calculations were used to investigate the sensitivity that the aquifer systems had to hydraulic parameter changes during Quaternary aggradation and degradation in the main valley. Those calculations were also done for base-level changes in the Rivers Rhine and Birs. For this purpose, this study considered probable historic base-levels before river regulation occurred, and before river dams and power plants were constructed. We also focused on scenarios considering increased groundwater recharge rates, e.g. due to exceptional long-lasting precipitation, or heavy rainfall events in the catchment area. Our results indicate that increased groundwater recharge rates in the catchment areas during such events (or periods) are associated with orders of magnitude increases of regional inflow into the Upper Muschelkalk karst aquifer. Furthermore, the groundwater fluctuations and groundwater saturated regions within the karst aquifer shift to places where high densities of sinkholes are documented. When the surface water base-levels adapt to probable historic levels, it leads to increased hydraulic gradients (i.e. local lowering of the groundwater level by up to 7 m). Those increased gradients are associated with increased groundwater flow within some aquifer regions that are particularly prone to karst development.


Introduction
Karst aquifers are characterized by caves, springs, and conduits that form underground drainage systems, which were shaped by dissolution (subrosion), whereas subsidence processes can lead to sinkholes (Ghasemizadeh et al. 2012). Since sinkholes develop from subsidence, they have been recognized as a major geohazard, resulting in huge financial losses to society, for they pose great threats to infrastructure, such as buildings, roads, bridges, and pipelines. The major factor that triggers karst evolution is limestone dissolution associated with heavy rainfall, high rates of downward groundwater leakage, and rapid declines of potentiometric levels due to groundwater over-exploitation (Ford and Williams 2007;Gutierrez et al. 2014). Meanwhile, land subsidence can also be attributed to salt solution mining, because of subsurface dissolution (subrosion) of evaporites such as halite and gypsum.
Many areas in Central Europe are characterized by complex landforms and geological settings that produce multiscale groundwater flow systems (Tóth 1963b). Among the basic requirements for rock dissolution, one of the most important ones is the energy provided from hydrostatic head differences, which causes groundwater flow through the system (Anderson and Kirkland 1980;Cooper 2002;Johnson 1981Johnson , 2005. Gravity drives groundwater flow in a sedimentary basin, and it is thought to be a superposition of several flow systems. Additionally, topography and anisotropy determine the geometry and the penetration of different flow systems (Tóth 1962(Tóth , 1963a(Tóth , 1995. That concept was later extended by Zijl (1999) as well as by Butscher and Huggenberger (2007), Marklund andWorman (2011), Tóth (2009) and Worman et al. (2007).
Surface recharge and precipitation in the entire catchment area drive flow within karst aquifers, and the flow is often directed toward an entrenched river as a base-level (e.g. Bauer et al. (2005); Kaufmann and Braun (2000)). Groundwater recharge and discharge drives the flow through the aquifer, whereas the height of the base-level establishes the pressure difference for flow within the aquifer. Likewise, changes during extensive groundwater recharge periods are the main drivers controlling hydraulic changes of deeper regional groundwater flow systems. The base-level is defined as the lowest level to which running water or groundwater can flow and erode.
Such "catchment-controlled" recharge and discharge boundary conditions imply that the water collected in the catchment is drained by the karst aquifer and discharged through the base-level (Palmer 1991). A karstic base-level in continental systems is usually considered to be very closely related to the floodplain of the main fluvial course to which the rock massif drains (Ford 1998). A fluvial base-level drop would lead to fluvial vertical incision and headward erosion in valleys and channels, which at the same time would produce a drop in the karstic base-level of the karst aquifers that are connected to them (Columbu et al. 2015;Jennings 1985). Therefore, fluvial and karstic base-levels would co-evolve, although the latter base-level would always be dependent on the evolution of the former one (Ahnert and Williams 1997).
Accordingly, the development of topographic levels of receiving waters, i.e. river bed topography, is an important factor in karst evolution as it defines the base-level of the system. Numerical simulation has shown that conduits close to rivers develop perpendicular to the base-level, following the steepest gradient in hydraulic heads (Kaufmann 2002). The hydraulic gradient is strongly influenced by the aquifer base gradient, and the direction of groundwater flow is controlled by the topography and structure of the aquifer base.
Regional hydrological investigations often stress the influence of the actual aquifer structure and include the spatial orientation and morphology of the aquifer base and its displacement at faults (e.g. Butscher and Huggenberger (2007); Herold et al. (2000)). Likewise, karst evolution often relates to a fault and fold structures, as well as fracture systems. However, detailed studies report many of these structures to be either independent of or even obstacles to karstification (Kovacs and Jeannin 2003;Luetscher and Perrin 2005). The hydraulic connection between different aquifer units depends on their shape, arrangement, and internal structure . Thus, regional flow systems can be segmented by faults, which may either act as barriers to groundwater flow or behave as preferential flow paths (Caine et al. 1996;Moya et al. 2014Moya et al. , 2015Moya et al. , 2016Raiber et al. 2015).
Building dams, dikes, or tunnels may change groundwater flow through the system and increase hydraulic gradients (e.g. Epting et al. (2009);Gutiérrez et al. (2003); James and Kirkpatrick (1980)). Similarly, withdrawing large amounts of groundwater increases hydraulic gradients and can increase karst evolution, especially in connection with fracture-networks (Hauber 1971). In settings with complex geology, hydrogeology, and hydrogeological boundary conditions, including the interaction of surface water with groundwater systems as well as numerous anthropogenic impacts, monitoring and modelling tools must be established which allow researchers to study the relevant processes at the appropriate scales.
Geological 3D models contain derivations of a realistic spatial representation of geological structures, often from sparse data and interpretations (Calcagno et al. 2008;Hassen et al. 2016;Pakyuz-Charrier et al. 2018). Principally, the surfaces that 3D geological models represent must be chosen in accordance with the study's main objectives and the given scale of the study area . The modeling objectives can vary from spatially depicting geological features to generating input data for physical process simulations (Calcagno et al. 2008). Several sources of data are usually employed to establish 3D geological models, including borehole logs, geological maps, and geophysical data (Wellmann et al. 2010). In short, 3D geological models are composed of various surfaces (e.g., strata contacts and faults) defining discrete volumes of lithological units for which lithological and hydraulic properties (among others) can be defined .
In this study, we follow the theory of regional-scale groundwater circulation for an unconsolidated gravel and complex karst aquifer system. This approach is based on geological and steady-state hydraulic 3D modeling at the catchment-scale. The main research questions we discuss include sensitivity analysis, which studies how sensitive the groundwater systems are to changes in groundwater recharge rates, historic base-level changes of surface waters, and hydraulic properties of the subsurface. In addition, we investigate how fluctuations of potentiometric heads in regions that are underlain by carbonate aquifers can at least partially explain sites of increased karst evolution, which can be documented by the location and density of sinkholes. Likewise, we study the impact that adjusting surface water baselevels to probable historic levels has on river regulation as well as aggradation and degradation of the river valley. We hypothesize that such base-level changes are associated with changes in hydraulic gradients that increase karst evolution.

Study area
The study area is located east of the city of Basel, Switzerland, in the Rhine Valley within the Tabular Jura. The Rhine Valley is densely populated, and the main land uses are residential, industrial, agricultural and infrastructure. The shallow groundwater resources of the area are mainly used by industry as well as for drinking water production (drinking water supply in the Hardwald; industrial use in Schweizerhalle).
Previous groundwater modeling results have demonstrated that large-scale industrial pumping affects the groundwater flow field in the Upper Muschelkalk aquifer at distances of up to 2 km to the south (Huggenberger et al., 2010;Spottke et al., 2005). The term "Upper Muschelkalk" (UMK) aquifer is commonly used and was retained here. In 1958, increasing water demand-led officials to install an artificial groundwater recharge system that was also designed to induce a hydraulic gradient toward areas of potential risk (e.g. urban areas, landfills, industrial areas, highways, railways and a harbor for the chemical industry). The groundwater flow regime is strongly influenced by this artificial groundwater recharge scheme as well as by abstraction, together with highly transient river-groundwater interaction processes along the Rivers Rhine and Birs (Fig. 1). However, there are only a few indications of the influence that the deeper (> 100 m bgs; Fig. 1) and regional scale groundwater flow systems have been influenced by these processes (Popp et al. 2019).
Within the last 20 years, there has been increased pressure on these systems. This is due to intense groundwater abstraction for water supply and for the construction of subsurface infrastructure as well as salt dissolution mining in this densely populated vulnerable industrial and residential area. Therefore, monitoring programs and modeling instruments have been established to understand the complexity of the interacting processes that influence groundwater flow in the region at various scales (Huggenberger and Epting 2011).

Geological setting and aquifer units
The study area is characterized by a pronounced relief and by a series of aquifer units, which are typical for many complex groundwater systems in front of mountain chains such as the Alpine foreland and the Jura Mountains of Central Europe. The two most important regional aquifers include: (1) the Quaternary aquifer of the Rhine Valley, which consists, of course, highly permeable unconsolidated gravel deposits; and (2) the karst aquifer within the UMK limestone and dolomite sections (Fig. 1). Morphologically, the carbonate units on the edges of the Rhine Valley within the Tabular Jura represent the recharge area of the basin. Although many aspects of the system dynamics at different scales are not yet well understood, several projects, which have focused on groundwater management and quality issues within the karst and shallow alluvial aquifer systems, offer some insight into the complexity of the multi-scale flow systems (Huggenberger et al. 2010;Moeck et al. 2018;Spottke et al. 2005;Zechner et al. 2011).
Across the western part of the study area (in a NNE direction) is the eastern main border fault of the Upper Rhine Graben, the so-called flexure zone (Fig. 1). With a vertical displacement of about 1400 m, it separates the Cenozoic sediments of the graben from the Triassic and Jurassic strata, which cover the Paleozoic basement at the rift shoulder. Although the main part of the study area is positioned on this eastern rift shoulder, it was still influenced by multistage tectonic evolution mainly during the time period ranging from the Eocene to the early Oligocene (Laubscher 2004(Laubscher , 1982Schumacher 2002;Spottke et al. 2005).
The Cenozoic fill of the Upper Rhine Graben in the western part of the study area is characterized by having a low permeability, which is related to the regional scale we analyzed. Because of very steep and vertically arranged units in the Middle to Lower Dogger-group and in the Keuper-group, the flexure zone is considered to have low permeability. The series of NNE-SSW striking horst and graben structures, which are typical for the Tabular Jura, in the central and eastern part of the study area form a complex aquifer and aquitard geometry, which is overlain by unconsolidated Quaternary sediments. In general, the lithostratigraphic layers ( Fig. 1) of the Tabular Jura dip slightly to the southeast and are characterized by varying hydraulic conductivities (Table 1). This part of the Rhine Valley is bounded to the north and south by hills which are part of the Tabular Jura and which develop a shallow basin in between (Fig. 3). The lithostratigraphic sequence of the area extends from the Paleozoic crystalline basement with a hiatus between the Late Jurassic and the Eocene to the Cenozoic periods ( Fig. 1). Due to their high porosity and permeability, the Quaternary cover of fluvial gravels represents a productive shallow aquifer within the basin.
The Paleozoic basement consists of crystalline basement covered by sediments of Carboniferous to Permian ages, and of the early Triassic Buntsandstein-group. Together, they are considered to be a basal aquifer with hydraulic properties that are dominated by the fracture network. The succession continues with the Lower Muschelkalk and the evaporites of the Middle Muschelkalk, which are considered to have a low permeability throughout northern and northwestern Switzerland (Gürler et al. 1987;Nagra 1988Nagra , 2002. They constitute a lower aquitard, which serves as a confining unit for the overlying aquifer. The highly porous and fractured dolomites, from the late Middle Muschelkalk-group and the UMK-group with limestones and dolomites, constitute a regionally important aquifer (Gürler et al. 1987;Nagra 1988Nagra , 2002Schmassmann 1990). The overlying sediments of the Keuper-group, and the Early to Middle Jurassic sediments, consist of marls and clays and are defined as the upper aquitard. They separate the UMK aquifer from the important regional Middle to Late Jurassic aquifer, which is mainly represented by the Hauptrogenstein in the southern part of the investigation area. The thickness of the UMK aquifer ranges from 50 to 80 m, and the average hydraulic conductivity is 1.3 × 10 4 m s −1 (Saladin (2004), compilation of borehole tests). The sediment permeabilities in the aquitards are below 1.0 × 10 7 m s −1 (Nagra 2002). There are only a few indications of the influence that the deeper (> 100 m bgs; Fig. 1) regional-scale groundwater flow systems have on the upper aquifers (Popp et al. 2019).
The shallow Quaternary aquifer is characterized by typical saturated thicknesses of the River Rhine gravel deposits, ranging between 10 and 20 m, and which locally can reach more than 30 m in thickness. Hydraulic conductivity measurements that have been determined from 56 pumping tests in the study area reported an average value of 3.1 × 10 3 m s −1 for these sediments (Saladin 2004).
Intensively fractured zones related to brittle fault systems that cause increased permeability may favor vertical exchange of groundwater across aquitards (Gürler et al. 1987). At the regional scale (several km), the connectivity of aquifers and aquitards is primarily influenced by the structural setting of horst and graben structures (Schmassmann 1990). Likewise, the horst and graben structures are important, not only because of the fractures and their conductivity but also because of the displacement of aquifer units which connects or separates them. The general dip of the Mesozoic strata to the SE in the study area results in direct contact between the Quaternary and the UMK aquifer, mainly in the western central part of the study area (Fig. 1).
Hydraulic exchange processes between the Rivers Rhine and Birs, and the groundwater, are proportional to the difference between river level and hydraulic head in the connected aquifer. The present river level and riverbed along the River Rhine, however, does not vary significantly upstream of the ''Birsfelden'' river dam and power plant (Fig. 2). Consequently, exchange fluxes are governed more by transient hydraulic groundwater heads.
Large-scale groundwater pumping for industrial use results in increased infiltration in the sections of the River Rhine, located in the western central part of the study area, and in places where the river is impounded by dams (Fig. 2). For all other sections along the River Rhine, the exfiltration of groundwater into the river, and especially where artificial recharge along the Hardwald area (Fig. 2) in the northwestern part of the model takes place, is the dominating exchange processes.
The history of aggradation and degradation in the River Rhine valley have considerably changed the altitude of the riverbed, and consequently, the groundwater flow regime from the Quaternary until today. The most recent influences on potentiometric heads in this valley are the construction of river dams for electricity production as well as the development of artificial groundwater recharge and extraction schemes. With the regional 3D hydraulic model, these changes in groundwater circulation could be investigated in relation to karst evolution following the aquifer base gradient approach for regional-scale karst and fractured flow systems (Butscher and Huggenberger 2007).

Karst evolution and impact of subsurface salt dissolution mining
There are complex relationships between the development of underground cavities and groundwater circulation within the various geological formations. Rock solution is formed where weakly mineralized groundwater can dissolve salt, gypsum, or carbonates. Rock salt, gypsum, or carbonates have very different solubilities, so the cavity formation can  In several locations within the study area, Zechner et al. (2011) investigated the observed land subsidence which was due to subsurface dissolution (subrosion) of evaporites, such as halite and gypsum, which could likely be attributed to salt solution mining. Depending on the geological settings, such as the depth of the karst, and the mechanical properties of overlying sediments, subrosion may cause widespread land subsidence. Salt solution mining started around 1850, but the resulting subsidence was only monitored after 1900. The highest rates of up to 100 mm a −1 , which have been observed during the last 10 years, have affected railway tracks for high-speed trains that are part of the main north-south travel connection through Switzerland. Fig. 2 Overview of the regional 3D hydraulic model with the integrated different geological units (Fig. 1) for the parameterization of groundwater recharge as well as surface topography (derived from a 90 m SRTM DEM). Upper left and below: Longitudinal profiles along the Rivers Rhine (R -R') and Birs (B -B'), illustrating the current river head, the base of the unconsolidated rock which coincides with the "maximum" base-level (scenario A), as well as an "intermediate" base-level along the surface waters without river dams (scenario B). (Abbreviations see Table 1 and Fig. 1)

3D geological model
The 3D structure of the aquifer and aquitards in the region were modeled by using GoCAD © (Geological Objects Computer-Aided Design). GoCAD © is an integrated 3D geologic object modeling and visualization software application. The software is based on a Discrete Smooth Interpolation (DSI) algorithm (Mallet 1992(Mallet , 2002. The 3D geological model was developed using existing geological data, including lithology and stratigraphy, structural information such as folds and faults derived from geological profiles, as well as by incorporating data from digital elevation models (DEM) and from the mapped geology of the surface (GeORG-Projektteam 2013).
Approximately 2740 borehole data of differing depths were available and used from GeoData (Kirchhofer 2006). The model includes 32 faults and 7 modeled lithostratigraphic horizons to derive structural maps of the main aquifer-aquitard boundaries for groundwater modeling (Fig. 1). To define the site-specific hydrostratigraphy, relevant hydrogeological units were identified and, if appropriate, units with similar hydraulic properties were merged.

3D hydraulic model
The concept of the steady-state regional-scale numerical 3D hydraulic model is based on the approach of Tóth (2009), which assumes the topography acts as the hydraulic "engine" for the underlying regional groundwater circulation. Following this theoretical approach, local heterogeneities in hydraulic parameters, such as hydraulic conductivities, are overlooked, resulting in a description of a generalized groundwater flow field over a large area that combines the sequence of aquifer units and fracture systems.
Simulation of the regional scale groundwater flow requires an appropriate transfer of the 3D horizons and faults from the 3D geological model. To characterize the main regional aquifers, geological units were combined for the 3D hydraulic model. For the UMK aquifer and other lithological units a single-continuum approach was used in the present study. Such an approach is considered to be appropriate for flow simulation in regional-scale carbonate aquifers with comparable small solution cavities, and when there is an absence of well-developed and regionally-extensive conduit networks (Ghasemizadeh et al. 2012;Kalhor et al. 2019;Kovács and Sauter 2007).
The geometries of the relevant hydraulic units created by the 3D geological model were imported as surfaces into the hydraulic model. In contrast to the commonly layer-based integration of geological units, the respective volumes could be directly considered. For this purpose, the individual areas were exported from the 3D geological model from GoCAD © as ASCII point files and imported into COMSOL © .
The regional-scale 3D hydraulic model that was developed covers an area of approximately 71 km 2 . In total, it consists of 11 hydrogeological units, including the Rivers Rhine and Birs (Fig. 3). The volume of the hydraulic model is 6.8 × 10 3 m 3 , and the resulting mesh has over 7.7 × 10 6 elements (1.4 × 10 6 mesh vertices). The distribution of the hydraulic aquifer properties is based on calibration results of local hydraulic models (Huggenberger et al. 2010;Spottke et al. 2005), as well as on literature references to corresponding lithologies (Nagra 2002). Furthermore, representative ranges of hydraulic conductivity values for the individual geological units were evaluated by performing a sensitivity analysis. The hydraulic conductivities finally selected for the model parameterization are listed in Table 1. To consider the complex structural items that interlock between the flexure zone and the Tabular Jurassic units in the hydraulic model, simplified fractures were inserted. They are considered to have an aperture of 10 cm and hydraulic conductivities of 1.0 × 10 -4 m s −1 . The anisotropy factor of the hydraulic conductivities (k fz /k fx ) is 0.1.
Regarding the hydraulic constraints of the aquifer unit boundaries, the main challenge was to parameterize the model based on limited data. Groundwater recharge by percolating precipitation (seepage water that reaches the groundwater surface) was subdivided for the different geological units (Fig. 2). The groundwater recharge was integrated as a Neumann flux boundary condition. For the unconsolidated gravel aquifer, groundwater-saturated and -unsaturated areas were distinguished, based on average groundwater levels derived from the INTERREG II project (Wagner et al. 2001). Here, groundwater-saturated areas were considered as a Dirichlet fixed-head boundary condition. The Rivers Rhine and Birs were considered as a Cauchy boundary condition; depending on the river and groundwater heads, the exchange of the surface waters with the groundwater was calculated via a defined conductivity value for the riverbed (i.e., a conductance value). For the basic model, the river head corresponds to an average water level and was derived from a DEM with a 2 m grid resolution.

Sensitivity analysis of hydraulic properties of the aquifer units and the river bed
To see if the aquifer properties of individual geological and /or aquifer units, fractures, or altered hydraulic boundary conditions influence the regional groundwater circulation, the aquifer properties were evaluated by using the sensitivity 201 Page 8 of 18 of the hydraulic pressure (residuals) in the UMK aquifer of selected piezometers in the study area. Since the reaction to different hydraulic parameterization in the UMK aquifer is of importance, the selection of observation wells (OW) was limited to this aquifer ( Fig. 2 and Table 2).

Variations of groundwater recharge in the catchment area
For the basic model, groundwater recharge to the UMK aquifer and for the groundwater unsaturated areas of the unconsolidated gravel deposits was represented by a flux of 3.2 × 10 -9 m s −1 . This corresponds to a recharge rate of 0.3 mm d −1 only and is representative of drought periods with little recharge. For extensive groundwater recharge conditions in the catchment area (i.e., a period of several days of heavy precipitation or long-lasting rainfall events), groundwater recharge was simulated by considering a flux of 5.4 × 10 -8 m s −1 , which corresponds to a recharge rate of 140 mm mo −1 . Those are values which are in accordance with groundwater recharge simulations in the region for comparatively "wet" months (Huggenberger and Epting 2011). Groundwater recharge for the remaining, less permeable areas was further reduced by a factor of 2. These are the order of magnitude estimates which are broadly in agreement with values obtained from regional soil water balance models (Gudera and Morhard 2015).

Reduction of surface water base-levels
The River Rhine in this area is largely a receiving watercourse and, therefore, it is the base-level for the regional groundwater circulation. Two scenarios were evaluated, one where the bottom surface of the unconsolidated gravels was considered as fixed potential (scenario A), and one "intermediate" scenario (B) for which an even gradient between the Fig. 3 Overview to the aquifer -aquitard geometry of the regional 3D hydraulic model with respect to the combined geological units (Abbreviations see Table 1 and Fig. 1) eastern inflow and western outflow of the River Rhine, and the southern inflow and northern outflow of the River Birs, respectively, were assumed to be the base-level (Fig. 2). Scenario A can be associated with the "maximum" baselevels, while for the "intermediate" scenario (B), the effect of the river dams is removed. For the River Birs (Fig. 2), this was not possible for the entire length as the line constructed partly intersects with the base of the unconsolidated gravels. Here, it must be noted that in earlier times, the River Birs changed its course within the valley frequently, and the river could flow around this bedrock mound.

Sensitivity analysis
We evaluated representative ranges of hydraulic conductivity values for the individual geological and /or aquifer units and found they only had a minor effect on pressure conditions. Figure 4 shows the results for different hydraulic conductivities of the Quaternary and UMK aquifer as well as faults, and results for the conductivity of the river bed (conductance) of the Rivers Rhine and Birs as well as groundwater recharge.
The sensitivity analysis for the hydraulic conductivity in the Quaternary aquifer was conducted for values ranging from 2.5 × 10 -4 to 8.0 × 10 -3 m s −1 . A reduction of the hydraulic conductivity leads to lower groundwater heads in the UMK aquifer for most OWs, which can be explained by the reduced exchange between the unconsolidated gravel and the UMK aquifer. Exceptions to this are the German OW9 and OW10, where no reaction was observed. The strongest reactions (variations between the scenarios) were observed for OW1 and OW5 (0.28 m), OW2 (0.31 m), and OW9 (0.47 m; Table 2; Fig. 4).
The sensitivity analysis for the hydraulic conductivity in the UMK aquifer was conducted for values ranging from 6.3 × 10 -5 to 4.0 × 10 -3 m s −1 . As expected, the change in piezometry due to reduced hydraulic conductivities in the UMK aquifer shows exactly the opposite reaction compared to the Quaternary aquifer. The lower the hydraulic conductivities are, the higher the groundwater head is. As in the Quaternary aquifer, the strongest reactions (variations between the scenarios) can be observed for OW2 (0.28 m), OW4 (0.41 m), OW1 (0.55 m) and OW9 (0.76 m; Table 2; Fig. 4).
The conductance determines the exchange fluxes between the surface waters and groundwater in relation to river and groundwater heads. Likewise, higher conductance leads to a higher flux in both directions, which again depends on the hydraulic gradient between the surface water and groundwater. Since most stretches the River Rhine receive groundwater discharge, higher conductance values lead to increased fluxes of groundwater into the River Rhine (exfiltrating conditions), thus reducing groundwater heads. However, the sensitivity for this parameter, which was examined for conductance values ranging from 2.5 × 10 -5 to 8.0 × 10 -4 s −1 , can be considered very low, whereas only local reactions (variations between the scenarios) can be observed (e.g., for OW7 (0 m) and OW9 (0.15 m); Table 2; Fig. 4). A sensitivity analysis for the hydraulic conductivity values in the faults (flexure zones within the Tabular Jura) was conducted for values ranging from 1.0E-5 to 5.0E-2 m s −1 . High reactions (variations between the scenarios) were observed in OW1 (0.85 m), but only locally. However, at other locations either no reactions (e.g. OW2, OW10, OW11), or low reactions were observed (e.g. OW4 (0.13 m); Table 2; Fig. 4).
Groundwater recharge is considered as a flux in the regional 3D hydraulic model. The topographically relevant geological units are considered in different ways (Fig. 2).  Table 2), note the different axis-resolution for the plot "groundwater recharge". The box plots show (see also imbedded sketch lower right) the median (horizontal line within the gray shaded box), the quartiles Q1 and Q3 (shaded box), the upper and lower whiskers (horizontal bars outside of the box) as well as extreme outliers beyond the whiskers. Whiskers mark those values which are minimum and maximum unless these values exceed 1.5 times the inter-quartile range (distance between Q1 and Q3) A very important result is that the range of groundwater recharge rates we calculated for the sensitivity analysis shows that groundwater recharge has a very strong influence on piezometry (variations between the scenarios). Results range from no reactions (e.g. OW2, OW10 and OW11), to intermediate reactions (e.g. OW4(1.90 m); OW8(1.06 m)), to strong reactions (e.g. OW1(7.31 m); OW9(3.49 m); Table 2; Fig. 4).
In summary, the sensitivity analysis of different hydraulic conductivities of the Quaternary and UMK aquifers, as well as for faults and the conductance of the riverbeds, showed that they had no significant influence on the UMK aquifer's piezometry. In contrast, groundwater recharge is a dominant factor, as variations in the magnitude of this parameter can locally increase groundwater levels by up to 7 m. Figure 5 shows results of the calculated groundwater flow regime in the UMK aquifer. Near the River Rhine, exfiltration of groundwater into the surface water is dominant, and high flow velocities can be observed. One exception can be observed locally in vicinity of the Augst river dam and power plant in the east. Here, surface water infiltrates into the aquifer up-gradient of the dam, and groundwater exfiltrates into the surface water down-gradient of the dam. In the southern area, infiltrating conditions can be observed where the River Birs is partially in direct contact with the UMK aquifer. In contrast, very low flow velocities can be observed in the trench structures with increased hydraulic potential; here, the superimposed units with lower hydraulic conductivities (Keuper-group) numerically result in pressure increases. The transform faults can represent vertical hydraulic connections with locally distinct influences (Fig. 4). This is especially true for the faults crossing into the flexure zone in the southwest, where the UMK aquifer also connects to the River Birs (see also reaction of OW1).

Scenarios
The groundwater flow lines in Fig. 6 show that the changes in the groundwater flow field in the UMK aquifer result from increased groundwater recharge in the catchment area and lowered base-levels (scenarios A & B). In the northern area ( Fig. 6 left), almost no changes in the groundwater flow field can be observed for the different groundwater recharge rates we simulated. In contrast, the increased inflow from the southern catchment area (Fig. 6 right) changes the groundwater flow paths significantly compared to the basic model, especially in the western part of the study area. Here, the groundwater flow lines now terminate in the Quaternary aquifer, instead of flowing around the Wartenberg Graben structure (Fig. 3) and ending in the River Rhine as before. If, in addition to high groundwater recharge, the base-levels of the Rivers Rhine and Birs are lowered, then the groundwater flow lines in the northern area reach the River Rhine, which previously terminated in the Quaternary aquifer (Fig. 6  left). Also, in the southern area, a part of the groundwater flow lines reach as far as the River Birs, while the major part terminates in the River Rhine (Fig. 6 right). In general, from that result it can be assumed that solution processes are intensified during extreme events and in areas where the rate of groundwater throughflow is greater than elsewhere in the flow field.
In the following, the distribution of documented sinkholes is discussed (for Swiss territory on the southern side of the River Rhine only). Their distribution is discussed in connection with the fluctuation range of the potentiometric surface in the UMK aquifer. For this purpose, the surfaces of the hydraulic potential within the UMK aquifer for low and high groundwater recharge rates were exported and cut with the base of the unconsolidated gravel surface (base of the Quaternary aquifer). This results in a groundwater outcrop line, which shifts to the south for high groundwater recharge rates (Fig. 7). The outcrop lines indicate the transition between vadose and phreatic conditions, as well as changes in confined conditions for different hydraulic conditions.
Obviously, the resulting outcrop line for high groundwater recharge rates coincides with the distribution of documented sinkholes. However, a difference can clearly be observed between the area west and east of the Wartenberg Graben (Fig. 3). While in the western area the outcrop lines for the different recharge scenarios are almost identical, in the eastern area the outcrop line shifts by up to 900 m towards the south. A visualization of the difference between the two hydraulic potentials for medium and high recharge rates shows small differences for the western area, and up to 12 m differences in the eastern area. This can be explained by units with lower hydraulic conductivities (Keuper-group) superimposed on the UMK aquifer for the eastern area, resulting in confined groundwater conditions. Figure 8 illustrates the results of the base-level changes of the Rivers Rhine and Birs. Compared to the basic model, the scenarios of "maximum" (scenario A) and "intermediate" (scenario B) base-levels change the groundwater flow velocities, and extensively lower the groundwater level in the UMK aquifer by up to 22.5 and 7 m, respectively (Fig. 7).

Discussion
Regional-scale geological and hydraulic 3D modeling made it possible to consider hydrogeological changes in relation to karst evolution. Our study has indicated that geological discontinuities, such as stratigraphic horizons and faults, play a key role in controlling groundwater flow in the region. Commonly used graphical user interfaces for groundwater simulation offer only limited possibilities for mapping geological settings with more complex tectonic structures. This is because of four reasons: (1) geological and geophysical data from different sources cannot be adequately imported and compared; (2) 3D structures, such as faults and faulted stratigraphic horizons, are very difficult to model accurately; (3) the consistency of a geological model of faulted horizons cannot be thoroughly reviewed; and (4) automated discretization of hydrogeological properties in such a setting is often not manageable by the graphical user interface (Spottke et al. 2005). An integration of the geological structures of the Tabular Jura and fault structures required the implementation of a layer-independent software (Comsol © ). In this way, changes in hydraulic boundary conditions can be evaluated. In this study, this was found to be particularly the case for evaluating the influence of changes in precipitating on groundwater recharge and on base-level changes of the Rivers Rhine and Birs, and for evaluating the sensitivity of hydraulic aquifer properties of individual geological units.
The results of a sensitivity analysis confirmed that the regional scale 3D hydraulic model could represent changes in regional groundwater circulation very well, even without a detailed knowledge of hydraulic parameters. Furthermore, the sensitivity analysis that was undertaken for the regionalscale model showed that regional groundwater circulation is determined primarily by the geometry and sequence of aquifer units and fracture systems. As the research was primarily focused on assessing how sinkholes develop in the area under investigation, the parameter sensitivity of the deeper basal aquifers, as well as those of the hydraulic characteristics of the flexure zone, were not further investigated in this study.
The modeling indicated that the groundwater outcrop line shifted to the south in the area east of the Wartenberg Graben during periods of exceptionally high groundwater recharge. This suggests that dissolution and karstification processes take place over long time-periods and in recurrent extreme events, which is documented by sinkhole development.
However, the reason why there are no sinkholes documented to the east of the Adlergraben (Fig. 3) should be investigated, even though the geological and hydrogeological settings are similar. It is possible that sinkholes have not yet been detected in this area or have not been the focus of previous research. Likewise, additional fracture structures Fig. 8 Influence of base-level changes of the Rivers Rhine and Birs on groundwater circulation within the UMK aquifer. Basic model (I) and scenarios "maximum" (II) as well as "intermediate" (III) base-levels (black arrows correspond to Darcy flow velocity). Also shown are maps of hydraulic head differences for the scenarios in relation to the basic model which have not yet been considered to have a hydraulic influence could be included in future scenario considerations. Due to the hydraulic connection of the UMK with the Quaternary aquifer to the west of the Wartenberg Graben, the position of the groundwater outcrop line within the UMK aquifer remains stable, and the Wartenberg Graben itself acts as a hydraulic barrier.
With the Triassic layers dipping towards the southeast, the hydraulic pressure within the Keuper-group sequences changes for different hydraulic and confined conditions. Likewise, the outcrop lines (Fig. 7) that were derived indicate the transition between vadose (i.e. unsaturated) and phreatic conditions within the Keuper-group sequence. Karstification under confined settings generates different morphologies than under unconfined settings (e.g., Klimchouk et al. (2006)). Moreover, 2-D or 3-D densely packed maze conduits (conduit loops along same flow direction) are likely under confined conditions, in contrast to conduit-like morphologies, which tend to be broadly spaced and dendritic due to highly competing developments under unconfined conditions. This particularly applies to local effects, which can hardly be resolved in the regional scale model we developed. In contrast, essentially regional effects of discrete hydraulically effective features, which act as an internal drainage system that controls the overall water budget, can be identified and described better in a quantitative way.
Likewise, the results of geochemical investigations of different groundwaters can only be understood by a combination of geochemical data analysis and regional scale hydraulic modeling. Such a combined approach would allow researchers to differentiate the influence of different groundwater components. Due to the impounding of the Rivers Rhine and Birs, the level of the receiving waters has risen by several meters, compared with the situation before the artificial dams were built. The scenarios we simulated illustrate the significantly lower base-level at which karst evolution could have taken place in times before aggradation of the river valleys occurred, and much later before dam construction. It is worth noting that the Lower Terrace (Quaternary) was deposited during two periods (30-15 ka and 13-11 ka; Kock et al. (2009)).
At the same time, groundwater flow directions and velocities for individual areas that are prone to karstification have changed considerably. Likewise, during extreme groundwater recharge events, water volumes and flow rates for selected regions can be increased by up to a factor of 60, which is accompanied by increased dissolution processes.

Conclusions
The 3D modeling approach presented in this study is based on modeling regional-scale groundwater flow regimes. For such regimes, the topography and the hydraulic properties of the geological units, the sequence of aquifer units, and river base-levels are all dominant factors for controlling the groundwater flow regime. This approach makes it possible, even with limited hydrogeological data, to theoretically test different hypotheses and to visualize regional groundwater flow regimes more accurately. Using this approach, we could show that groundwater circulation can change considerably over longer periods, especially during extreme precipitation events with increased groundwater recharge, as well as under changing boundary conditions such as base-level changes.
Our geological and hydraulic 3D models provide an appropriate, yet flexible, foundation for decision-making for local and regional groundwater managers. Likewise, these tools are an excellent platform for evaluating the impact of future infrastructure developments and can be used for any 3D data queries regarding different aspects of an area's geology. They can also be used to qualitatively assess the vulnerability of local drinking water supplies, and to optimize suitable locations for installing groundwater and subsidence monitoring systems. The results of the research presented here form a base for further model development, which will help researchers focus on the causes of the subsurface salt dissolution, and ultimately provide predictions on land subsidence risks in the area.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.