Integrating field and numerical modeling methods for applied urban karst hydrogeology

Infrastructures constructed on unstable geologic formations are prone to subsidence. Data have been collected in the context of an upgrading project for a highway located beside a river dam constructed on gypsum-bearing formations. Surface water infiltrates upstream of the dam, circulates through the gravel deposits and into the weathered bedrock around and beneath the dam, and exfiltrates downstream into the river. As a result, an extended weathering zone within the bedrock and preferential flow paths within voids and conduits developed as part of a rapidly evolving karst system. Enhanced karstification in the soluble units of the gypsum-bearing formations resulted in subsidence of the dam and the highway. Since 2006, changes in the groundwater flow regime have been investigated by different methods that allowed the evaluation of the long-term performance of the infrastructures. Geological (outcrops, lithostratigraphic information from boreholes), hydrometrical (extensive groundwater monitoring, dye tracer tests) and hydrogeophysical (Electrical Resistivity Tomography, ERT) data were integrated into high-resolution 3-D hydrogeological and 2-D karst evolution models. The applied methods are validated and the sensitivity of relevant parameters governing the processes determined. It could be demonstrated that the applied methods for karst aquifer characterization complement each other. Short-term impacts and long-term developments on system-dynamics and the flow regime could be evaluated. This includes the description of the transient character of the flow regime during and after episodic flood events (surface-groundwater Correspondence to: J. Epting (jannis.epting@unibas.ch) interaction, conduit and diffuse model outflow) as well as the evaluation of time scales for karst evolution. Results allow the optimization of investigation methods for similar subsidence problems, ranging from general measurements and monitoring technologies to tools with predictive utility.


Introduction
While the characterization and modeling of flow in heterogeneous and fractured media has been investigated intensively, there are no well-developed long-term hydrogeological research sites for gypsum karst.Additionally, measurement systems for monitoring the evolution of karst phenomena are exceptional.This case study documents the integration of different methods in the context of an engineering project for the upgrade of a subsided highway located beside a river dam.
Surface and groundwater monitoring during engineering projects usually is restricted in order to comply with existing laws and regulations governing water quality issues during construction activities.In the current case study, sporadic measurements revealed that subsidence of the highway and the river dam has increased rapidly over the last ten years.At the beginning of the project, the knowledge was limited to purely conceptual models and sparse accurate monitoring data.As regards the local engineering problem, and in order to plan appropriate remedial measures, it was recommended to set up instruments that allow this complex system to be examined under various hydraulic conditions.This included an assessment of the current flow regime as well as the subsidence rate and its development over time.Such approaches require, along with the installation of surface and groundwater monitoring systems, specific field campaigns and modeling techniques to investigate the relevant processes.Furthermore, the developed tools should have not only a monitoring character; they should also enable long-term reliable predictions to be made concerning the future evolution of the system and subsidence.
Infrastructures that are constructed on soluble geologic formations are prone to subsidence (Gutiérrez, 1996;Lamont-Black et al., 2002).When located close to artificial hydraulic structures such as river dams, karstification can ensue due to: (1) the presence of evaporites; (2) elevated hydraulic gradients and (3) the presence of undersaturated water.Such boundary conditions can lead to increased leakage, subsidence and the failure of the hydraulic structures and nearby infrastructures.Especially when found within gypsum-bearing formations, karst features develop much more rapidly than in carbonate formations.Engineering aspects related to gypsum karst and dam construction are discussed by James and Lupton, 1978;James, 1992;Klimchouk and Andrechuk, 1996;Breznik, 1998;Milanović, 2000Milanović, , 2004;;Pearson, 2002;Jarvis, 2003;Johnson, 2003aJohnson, , b, c, 2005;;and Payton and Hansen, 2003.An example is the catastrophic failure of the Quail Creek Dike in Southwest Utah in 1989 due to flow of water through an undetected karstified gypsum unit beneath the earth-fill embankment (Johnson, 2008).A list of leakages at 42 dams worldwide is given by Milanović (2000).
The site-specific aspects of investigation projects and available data sets determine the complexity of modeling strategies used, e.g., issues related to water resource management of karst aquifers can be resolved by black box or global approaches.For process understanding and the integration of more complex groundwater flow in karst systems (drain network and matrix), numerical models that represent double continuum media typical of karst aquifers are to be applied (e.g.Kovacs, 2003).Klimchouk et al. (2000) discussed the evolution of karst aquifers and Quinn et al. (2006) summarized the existing modeling approaches for simulating flow in karst environments.In the appendix these include: (1) models using equivalent porous medium in which flow is governed by Darcy's law (Anderson and Woessner, 1992); (2) models in which the preferred flow paths are simulated with a very high hydraulic conductivity relative to the surrounding matrix material (double porosity); (e.g.Teutsch, 1989;Mace, 1995;Kiraly, 1998;Eisenlohr et al., 1997;Josnin et al., 2000); (3) "black-box" approaches in which functions are developed to reproduce input and output data (recharge and flow at discharge springs; e.g.Dreiss, 1989a, b), as well as "global" approaches which include the hydrological dynamics of the conduit and the diffuse flow system (Butscher and Huggenberger, 2008); (4) fracture network simulations in which individual factures are mapped and then studied (Long et al., 1982;Long and Brillaux, 1987), and (5) open channel equivalents (Thrailkill et al., 1991).
For a more fundamental understanding of rock-groundwater interactions and the evolution of flow within karst aquifers, modeling techniques are based on fundamental and well-established physical and chemical principles.These models allow important processes, ranging from initial small-scale fracture networks to the mature karst, to be studied.Furthermore, model simulations facilitate to assess time scales for karst evolution and to quantify rates at which solution processes operate.In the present case study, the simulation of the evolution of flow within the gypsum karst aquifer was approached by a simplified 2-D finite difference method where groundwater flow is coupled with equations of dissolutional widening.The evolution of two-dimensional karst aquifers has been intensively studied during the past decades (Groves and Howard, 1994;Clemens et al., 1996;Hanna and Rajaram, 1998;Siemers and Dreybrodt, 1998;Dreybrodt and Siemers, 2000;Gabrovšek et al., 2000;Gabrovšek andDreybrodt, 2000, 2001;Kaufmann and Braun, 2000;Bauer et al., 2003Bauer et al., , 2005;;Romanov et al., 2002, Liedl et al., 2003), also close to dam structures (Dreybrodt, 1992(Dreybrodt, , 1996;;Dreybrodt et al., 2001Dreybrodt et al., , 2005;;Romanov et al., 2003Romanov et al., , 2007)).These previous approaches to model karst evolution focused on hypothetical karst catchments with synthetic data.Few studies have attempted to integrate field data.
An approach is presented which merges high resolution 3-D hydrogeologic modeling (3-D HGM) with 2-D karst evolution modeling (2-D KEM; Fig. 1).The different modeling techniques reproduce different aspects of the hydrologic processes and were employed by independent modeling teams.This allowed the estimated parameters to be cross-checked and the results from both approaches to be evaluated continuously and interpreted separately.The 3-D hydrogeological model (3-D HGM) presented in this paper includes a deterministic finite difference approach which takes into account an equivalent porous medium for weathered and non-weathered rock, and a coupling of the system with drains that represent a generalization of the conduit component of model outflow (mixed flow in karst settings; Quinn and Tomasko, 2000;Quinn et al., 2006).The calibrated 3-D HGM provided the geometric and hydraulic boundary conditions (modeling domain, groundwater and head of the river, hydraulic conductivity and conductance of the river bed) for the 2-D KEM.The calibrated 2-D KEM allowed the karst aquifer evolution to be simulated to its current state and the sensitivity of changing natural and anthropogenic boundary conditions to be investigated.Subsequently, resulting aquifer heterogeneities simulated using the 2-D KEM were transferred to the 3-D HGM.Applied modeling strategies complemented each other, results from the 3-D HGM could be used for the validation of the 2-D KEM, and vice versa.
Modeling results are compared, validated and discussed, together with multiple data sources (lithostratigraphic information from boreholes, extensive groundwater monitoring, Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1163/2009/  dye tracer tests, hydrogeophysics) as well as water budgets through defined cross sections over time in the context of the flow regime.The integration of different types and a varying quality of data sets into the different models, presented a particular challenge.

Settings
The project area is located in the Lower Birs Valley southeast of Basel, Switzerland.Over the last 30 years subsidence of a man-made river dam and an adjacent highway has been observed (Figs. 2, 3 and 4).The dam in its current dimension was constructed in the 1890's (Golder, 1984).However, documentation of water engineering measures in this region, as the diversion of water for early manufacturing purposes in Basel, goes back as far as the 11th century (Fechter, 1856).The height difference to the base level downstream of the dam is 7.3 m.This hydraulic gradient is used for the generation of hydropower from a small hydro-electric power plant.As there is sufficient water supplied by the Birs River, the height of the impounded water upstream of the dam is practically constant at 266.2 m a.s.l.Surface-groundwater interaction is dominated by the hydraulic head of the river Note, as evidence for the subsidence of the dam structure, that the crest is only overflowed on its left side (a diagram is given in Fig. 2).
and variations in river bed conductance upstream of the dam during flood events.Upstream of the dam, river water infiltrates into the highly permeable fluvial gravels and the weathered bedrock, follows the hydraulic gradient around and beneath the dam and exfiltrates downstream into the river.These processes and the drilling of several boreholes in the 1990's that likely enhanced the connection between aquifers and led to the karstification in the soluble units of the "Gipskeuper" and resulted in an extended weathering zone within the bedrock as well as in the development of preferential flow within voids and conduits.
To prevent further subsidence, construction measures were carried out in two major project phases in 2006 and 2007.The highway was supported by 166 piles and by a sealing pile wall, consisting of approximately 300 piles (Fig. 2), to prevent infiltrating river water from circulating around the dam and beneath the foundation of the highway.Piles extend down to the non-weathered rock at a depth of 20 to 25 m.Caves encountered when the piles were being constructed were filled with a total of 168.2 m 3 of supplementary grout mixture, in order to plug underground water channels and stabilize the ground beneath.In compliance with existing regulations, an observation network was installed in order to monitor surface and groundwater quality in the vicinity of the construction site and regional drinking-water supplies further downstream.Additionally, the observations allowed the identification and evaluation of the relevant processes.Online monitoring allowed the early detection and documentation of changes in hydraulic constraints in the vicinity of the construction site and the dam.To prevent surface and groundwater pollution during the construction measures, interventions could have been initiated in case of breakthrough events and the mobilization of void fillings and grout mixture injections.

Geology and hydrogeology
The stratigraphic column in Fig. 5 includes the Triassic, Jurassic and Quaternary units.Karst development mainly occurs in the gypsum-bearing parts of the Triassic.Quaternary gravels, silty floodplain deposits, as well as artificial fillings beneath the highway overlie the westward-dipping Triassic and Jurassic units on the right side of the river.For better visualization of the karst-relevant units and their spatial distribution the Quaternary sequence has been removed.The tectonic settings (see below) result in a complex spatial distribution of the main geological formations which consist of marls and clays, dolomites and sandstones, marls and, for most of the investigation area, of Gipskeuper.The map also shows the course of the Birs River in the year 1798 compared to the situation in 1983.The river, which was straightened in the 19th century, cuts into the Triassic bedrock, resulting in a narrow couloir.
The Gipskeuper is made up of gypsum near the surface, anhydrite and intercalations of marls.The lithological term "Gipskeuper" as used in this paper generally includes the mineral "gypsum" and also refers to "anhydrite", which, in the deeper subsurface, is the more common anhydrous form of calcium sulfate.
In the non-weathered state anhydrite dominates and the Gipskeuper is characterized as being rather low permeable.Hydraulic conductivities of Gipskeuper, as tested in borehole and modeling studies in northern Switzerland, vary between 1E-14 ms −1 and 1E-07 ms −1 (Nagra, 2002).However, in the weathered state, Gipskeuper shows features of a heterogeneous (karstified) aquifer (Klimchouk and Andrechuk, 1996).
The investigation area is characterized by the Eastern Rhinegraben Master fault accompanied by an intense tectonic segmentation into fault-bounded blocks (Schmassmann, 1972).Tectonically the lithological units dip at an angle of approximately 45 • to the West.Several NNE-SSW normal faults subdivide the units into a series of blocks with variable hydraulic properties.Fault and fracture zones are associated with rock weakness and can locally increase permeability within sequences, resulting in an enhanced groundwater leakage and the development of paths for preferential flow.
In general, reduced water velocities behind reservoir dams enable the sedimentation of fine material, resulting in a clogging layer that delays the infiltration of surface water into groundwater systems.During episodic floods, river bed permeability increases and processes of surfacegroundwater interaction are enhanced.These processes are highly transient and consequently, changing hydraulic river bed conductance plays a key role in understanding surface-groundwater interaction.It is worth mentioning that the investigation period includes the 300-year flood of 9 August 2007 (peak-discharge 370 m 3 s −1 , average discharge 2007 is 19 m 3 s −1 ).
Borehole data show that most voids and solution cavities, which generally contain debris, clay, gravel and calcite fillings, are concentrated at the base of the weathered Gipskeuper (lixiviation front).During episodic flooding, these sediments can be partially flushed and subsequently, more aggressive water can enter the system.

Hydrochemistry
Chemical analyses of samples taken at specific hydraulic boundary conditions allowed hydrochemical boundary conditions to be defined.SO 4 -Ca content increases in the circulating groundwater within the gypsum formation.Solution rates depend on water chemistry and dynamics of the flow regime.Water entering the system is assumed to be dominated by the chemical composition of infiltrating river water.Calcium and Sulfate concentrations in the river indicate the lowest observed concentration values, ranging from 88.8 to 94.7 mgl −1 and from 7.8 to 15.8 mgl −1 , respectively.Highest concentrations were observed in groundwater samples taken at two different depths from observation well OW4 (Fig. 2), ranging from 256.4 to 277.3 mgl −1 for Calcium and from 118.7 to 116.8 mgl −1 for Sulfate.Water samples taken from the groundwater outlets downstream of the dam show intermediate concentrations, ranging from 135.9 to 204.4 mgl −1 for Calcium and from 53.3 to 92.4 mgl −1 for Sulfate.
The geochemical composition of the Gipskeuper formation is strongly heterogeneous, as indicated by mineralogical analysis.Therefore, the transition zone between weathered and non-weathered Gipskeuper should not be considered as a sharp boundary.
Dissolution kinetics and dissolution rates R of pure gypsum during water-rock interaction in laminar or turbulent flow were determined by Dreybrodt (1987Dreybrodt ( , 1988Dreybrodt ( , 1990) ) and Jeschke et al. (2001), respectively.R can be described by a rate law: where c is the concentration of Calcium in the water and c eq is equilibrium concentration with regard to gypsum.c s is a switch concentration, where the dissolution rates (molcm −2 s −1 ) switch from a linear rate law to a non-linear one with the order n.The values of c s , k n and n are characteristic for the mineral because they are entirely dependent on surface reactions.The characteristics of gypsum dissolution imply that concentrations close to equilibrium concentrations are reached rapidly.However, water infiltrating into the karst system will have a remaining potential for solution over long time periods and after long percolation pathways.

Conceptual approach
The conceptual approach and the applied methodologies are based on a series of questions that generally arise in the context of urban hydrogeology (e.g.Epting et al., 2008): (1) are the existing data sufficient to answer the relevant questions, and which hydrogeological concepts can be set up using the data?(2) Which additional data and experiments could improve predictions and allow hypotheses to be tested?
(3) How could data acquisition and experiment design be optimized?Consequently, the investigation process is an iterative procedure of consecutive data acquisition deriving from specific experiments.Hydraulic conductivities in karst aquifers are extremely heterogeneous, ranging from 10E-08 to 10E-05 ms −1 in the fracture systems, and up to 1 ms −1 in conduits (White, 1988;Ford and Williams, 1989).Structural features such as folds and faults can have significant influence within karstified areas.As a result, modeling groundwater flow in site-specific karst environments is extremely challenging.Modeling results often are highly uncertain because of the lack of site-specific information on heterogeneous subsurface structures and the resulting complexity of flow paths.
A stepwise approach is presented where results from different methods are merged.
In a first step, the investigated area was delineated, comprising an inventory of all relevant boundaries defining the flow regime, including an evaluation and parameterization of the fundamental processes governing the system.A core element of the present approach is the merging of high resolution 3-D HGM with 2-D KEM (Fig. 1).The calibrated 3-D HGM provides the geometric and hydraulic boundary conditions for the 2-D KEM.The calibrated 2-D KEM allowed karst aquifer evolution to its current state to be simulated.Resulting aquifer heterogeneities simulated using the 2-D KEM were transferred to the 3-D HGM.The applied modeling strategies complement each other.The results of the 3-D HGM could be used for the validation of the 2-D KEM and vice versa.
The data required for setting up the 3-D HGM and 2-D KEM with accurate boundary conditions can be of quite different type and varying quality which can be termed "hard" or "soft" data (cf.Regli et al., 2002).The most reliable hard data derive from outcrop and laboratory investigations.The information on several outcrops was incorporated into the geological map (Fig. 5), which formed the basis for delineating the models.Drill core data provide limited information on the spatial distribution of subsurface properties, especially in karst environments where the probability of encountering voids is quite low and relies on a hit-or-miss approach.The quality of individual drill core descriptions varies considerably, depending on the geotechnical approach used, thus permitting limited and speculative conclusions.
The same is true for hydrogeophysical data which allow zones with different properties and behaviors to be delineated over time.Consequently, drill core and hydrogeophysical data are regarded as soft data.The terms of "hard" and "soft" data can also be applied to hydrometric data derived from hydraulic measurements and from tracer tests.Whereas hydraulic measurements from the head of the river can be considered as hard data, data from groundwater observation wells are hard data only if they independently sample one aquifer.In the case where observation wells connect aquifers, data interpretation is not distinct and consequently such measurements should be considered as soft data.Tracer tests can undoubtedly confirm hydraulic links between injection and observation locations; this information corresponds to hard data.The path of preferential flow between these locations, however, is ambiguous and relies on further interpretation, resulting in soft data.Additional data that are documented during construction measures, as for example, lithological information derived during the installation of piles, or information on locations and quantities of supplementary grout mixture injections, generally lack accuracy and should also be considered as soft data.Nevertheless, this kind of information can be indispensable in the process of setting up hydrogeological models.

Data sources
Multiple data sources were available: (1) soft data from lithostratigraphic information from borehole logs and the national geological map (1:25 000; Bitterli-Brunner et al., 1984); (2) soft and hard data from continuous groundwater measurements within the Quaternary and Gipskeuper aquifers; (3) soft data from general geological descriptions of piling works and locations of supplementary grout mixture injection; (4) soft and hard data from dye tracer tests, and (5) soft data from hydrogeophysical investigations.A total of 24 vertical boreholes were drilled in several investigation phases from 1993 to 2007 (Fig. 2).Most boreholes were developed as observation wells for groundwater or subsidence measurements.In total, 12 observation wells were fitted with automatic data loggers for monitoring physical parameters (hydraulic head, temperature and electric conductivity).Additional lithostratigraphic information could be derived from reports made during the construction of the piles.Hydraulic links and flow velocities within the study area were investigated by Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1163/2009/In particular, the interaction of the river with the groundwater system at different flow stages plays a major role in the initiation of movements within the void and conduit system.Surface and underwater hydrogeophysical measurements were carried out under different hydrologic boundary conditions, both before and after the construction activities.This interpretation focuses on the description of distinct geological and hydrogeological features in correlation with river discharge (Fig. 2, Epting et al., 2009).
4 Modeling approach

3-D Geological and Hydrogeological Model (3-D HGM)
The main goals in setting up the high resolution 3-D HGM were: (1) the description of the spatial distribution of the geological formations and the delineation of weathered units; (2) the simulation of the current flow regime and changes to it in the context of the construction measures and future evolution of the system and (3) the provision of model geometries and hydraulic boundary conditions for the 2-D karst evolution modeling.3-D HGM simulations were performed using the Groundwater Modeling System GMS v6.0 (Environmental Modeling Systems Inc., 2006) together with the 3-D finite difference code MODFLOW (Harbaugh et al., 2000).The "solid modeling" approach was employed for constructing the 3-D geological structures (Lemon and Jones, 2003).Solids (volumetric layers representing hydrostratigraphic sequences) were built using the horizons method, based directly on the lithostratigraphic data from 25 boreholes, 166 piles beneath the highway and 273 piles at the pile wall, as well as the introduction of 42 support points.Each lithological formation can be represented by a separate solid.In order to simplify the geological model, the formations of non-weathered Gipskeuper were grouped with the formations of Schilfsandstein and Mergel-Dolomit, resulting in three different materials, including the Quaternary cover, weathered Gipskeuper and non-weathered lithological sequences.Vertically, the model extends from the topographic surface to the non-weathered Gipskeuper (Figs. 5 and 6).
The grid for the 3-D HGM was automatically generated from the solid model geometry (Jones et al., 2002).The horizontal discretisation of the grid is regular (5 by 5 m).To represent the locations of drain features more accurately, a 5-layer approach was chosen.
Hydraulic boundary conditions were defined as follows (Fig. 6): (1) the northern and southern boundaries were defined as specified head, corresponding to available groundwater head measurements of the regional flow regime; (2) the eastern boundary was defined as specified flow from the adjacent catchment; (3) the western boundary was chosen as a no-flow boundary, according to the abundance of comparatively impermeable geological sequences (see above); and (4) the Birs River was simulated as General Head Boundary (GHB), where river infiltration and groundwater exfiltration are calculated in relation to the difference between river level and hydraulic groundwater head, as well as a conductance of the river bed.Modeling the complex flow using a finite-difference approach, with drain networks, representing the conduit component of model outflow, was achieved by using generalized drain features (Quinn et al., 2006).Two drains were introduced in model layer 4 corresponding to information obtained from (a) boreholes, indicating that voids were generally encountered at the bottom of the weathering zone; (b) the 1996 dye tracer test and (c) the location of fracture joints (Figs. 2 and 5).The drain elevation was chosen using values of nearby groundwater head measurements to ensure active flow.
To enhance model certainty, the following procedure was applied: (1) calibration of the 3-D HGM and comparison of observed and calculated heads in numerous groundwater observation wells; (2) inverse modeling, including parameter estimation and sensitivity analysis; (3) groundwater temperature data analysis of riverine observation wells in order to estimate hydraulic conductance values of the GHB upstream of the dam (Appendix A), and (4) scenario development, including the consideration of drains and different extensions of the weathered rock.Hydraulic conductivities of the lithological sequences as well as river bed and drain conductance were determined and their sensitivities evaluated by a combination of manual and automated parameter estimation procedures, which are based on numerical optimization algorithms within the nonlinear regression code PEST (Doherty, 1994).

2-D Karst Evolution Model (2-D KEM)
Distinct time intervals for karstification and the development of the aquifer can be defined for the present case study (Romanov et al., 2009).A first time interval covers the natural karstification process for a time period between several hundred and a few thousand years.The results from this evolution period are used as initial conditions for the second interval, which covers the time period from 1890 to 2007.This period is characterized by anthropogenic alterations to the system, including the construction of the river dam in the 1890's and technical measures to prevent further subsidence of the highway in 2006 and 2007, which considerably changed the evolution of the aquifer.The third time period, covering the evolution of the aquifer after 2007, with regard to the effects of technical measures on boundary conditions, will be the subject of future investigations.This includes the application of the developed modeling tools to forecast aquifer development for the 100 years that follow.
The main goals in setting up the 2-D KEM were: (1) to simulate the spatio-temporal development of karst features within the investigation area over the last 100 years; (2) to determine and evaluate the relationship of the investigated parameters; (3) to determine time scales for future system development, and (4) to provide heterogeneously distributed aquifer properties including distinct high conductive features for the 3-D HGM.
The reason for not modeling karst evolution with existing 3-D codes is the complexity of model geometries and chosen boundary conditions that for the moment cannot be used without unreasonable effort and CPU.However, as 3-D codes are gaining in efficiency, the available data sets will be transferred to 3-D KEM in the near future.
The 2-D KEM includes representative geological horizontal cross sections of the area around the dam structure.With respect to the high information density in the vicinity of the infrastructures, horizontal cross sections were preferred over the generally used vertical model setups.The geometric dimensions as well as the hydraulic and technical boundary conditions (modeling domain, groundwater and head in the river, hydraulic conductivities) of the 2-D KEM correspond to those of the layers 2, 3 and 4 of the 3-D HGM (see above).The modeling domain is 130 m wide (W-E), and 465 m long (S-N).Figure 7 shows the conceptual model setup, the impervious and insoluble dam structure and zones used to assign different hydraulic conductivities and solubilities.Flow within the karst aquifer is driven by infiltrating river water upstream of the dam that is directed towards the river downstream of the dam as a base level.In advanced scenarios, regional groundwater flows from South to North and from the adjacent slope to the East were considered.
Whereas the hydraulic properties for the 3-D HGM are represented with a continuum approach, those of the 2-D KEM are represented by a fracture network.The network of the 2-D KEM is characterized by the average spacing s of the fractures, their aperture widths a 0 and their widths b 0 .This fracture system exhibits a hydraulic conductivity (Lee and Farmer, 1993) of: where ρ is the density of water, η is the viscosity and g is earth's acceleration.Equation (3) was used consecutively for the transfer of hydraulic parameters of the 3-D HGM to the 2-D KEM and vice versa.
The application of gypsum dissolution kinetics (Dreybrodt, 1987(Dreybrodt, , 1988(Dreybrodt, , 1990;;Jeschke et al., 2001) enables the simulation of void and conduit enlargement with time by chemical dissolution and solutional widening within the system.This selective enlargement increases conductivity by several orders of magnitude during the early phase of karstification.Thus flow characteristics change from more uniform to non-uniform flow.An important moment for the stage of karst aquifer evolution is the so-called "breakthrough time", when flow changes from laminar to turbulent and increases by several orders of magnitude in a relatively short interval.(2) subsidence of the highway has been reported since the 1990's (Cantonal Archive Basel, unpublished reports), and (3) previous dam site modeling indicates that under normal conditions, evolution periods of several thousands of years can be shortened to periods of several decades for karst aquifers in the vicinity of hydraulic structures such as dams (Dreybrodt, 1992(Dreybrodt, , 1996)).
A first set of scenarios allowed the determination of the required minimal model resolution and initial hydraulic conductivities for the different zones within the modeling domain.As the hydraulic conductivity depends on the fracture aperture widths and the distances between them, the same hydraulic conductivities are obtained from a network with large but sparsely distributed fractures, as well as for fine fractures with small distances between them.To analyze the sensitivity of this interrelationship adequately, numerous scenarios had to be calculated in order to obtain reasonable values for the fracture aperture widths, while keeping the distance between fractures as large as possible (reduction of CPU).Therefore, calculations have been performed with relatively large time steps of 1 year, and only until turbulent flow after breakthrough occurred in the domain.
As structural initial conditions can influence karst evolution in a significant way (Birk et al., 2005), subsequent scenarios focused on the incorporation of specific subsurface information for the description of aquifer heterogeneity, including: (1) statistical distributions of hydraulic conductivity and solubility zones, and (2) discrete prominent fractures.

Results from borehole data and the dye tracer test
Analysis of drill cores enabled the determination of several permeable zones within the underlying bedrock and the already developed voids.Although the probability of encountering voids is fairly low and relies on a hit-or-miss approach, it was nevertheless possible to detect a total of 7 voids, with diameters ranging from 0.3 to 2.7 m, at depths ranging from 15 to 18 m.This information, considered as soft data, suggests that the vertical extension of the weathered Gipskeuper ranges from 2.2 to 14.8 m.It could be recognized that the drilling of several boreholes in the 1990's resulted in the connection of aquifers.These connections are documented in the drill core records and were confirmed independently by geochemical and hydraulic data (flowmeter measurements).It was possible to derive additional soft data from lithostratigraphic information in the reports made during the installation of the piles, and relatively precise cross sections of the investigation area were constructed based on the lithostratigraphic information from the boreholes (Fig. 5).Hydraulic links within the investigation area were confirmed by the 1996 dye tracer test, representing hard data (Cantonal Archive Basel, unpublished reports).The direction of the hydraulic links corresponds to the main directions of fracture joints within the investigation area (Fig. 5).The dye was injected in OW3 and sampled in OW2, 4 and 5 (Fig. 2).Highest concentrations were measured in OW5, which is nearest to the injection well.A secondary maximum was observed, representing a further preferential pathway.Measurements in the other observation wells resulted in lower concentration values indicating that they are influenced to a certain degree by infiltrating river water.Maximal groundwater flow velocities range from 85 to 111 md −1 , values typical for conduits within well-developed mature karst systems (Birk et al., 2004).Determined dispersivities range between 1 and 3 m.

Model calibration and sensitivity analysis
Table 1 summarizes the initial set of selected parameters, the calibration range and calibrated hydraulic parameters for the extended model with drains.Additionally, normalized composite scaled sensitivities are summarized for the calibrated parameters of all investigated scenarios.The upper value of 1.0E-07 ms −1 resulting from modeling studies in northern Switzerland (Nagra, 2002) was chosen as initial value for the hydraulic conductivity of non-weathered rock.
For the hydraulic conductivity of the weathered Gipskeuper, an initial value of 1.0E-05 ms −1 was chosen, assuming that the hydraulic conductivity of this lithological formation is much higher than the one of the non-weathered rock.The initial hydraulic conductivity value of 5.0E-03 ms −1 for the Quaternary cover represents an empirical value for fluvial deposits.In order to assess initial values for the conductance of the GHB and drains, these parameters were first calibrated manually with the basic model, resulting in 2.0E-06 s −1 and 2.0E-03 ms −1 , respectively.These initial parameters vary by one order of magnitude for the hydraulic conductivity of the Quaternary and the weathered Gipskeuper and by two orders of magnitude for the hydraulic conductivity of the non-weathered rock and the conductance of the GHB and the drains.
Features included in the extended model setups are examples of potentially important processes that can improve the modeling results.Upper boundary values for the hydraulic conductivity within the weathered Gipskeuper and non-weathered rock are attained, indicating that the inclusion of an extended weathered zone beneath the dam and drains is adequate.The reason for the calibrated high values of hydraulic conductance downstream of the dam structure compared to those upstream could be the occurrence of groundwater outlets in the more heterogeneous river bed downstream (high flow velocities and possible turbulent flow in the vicinity of the dam overflow) and the more clogged river bed upstream of the dam structure (area with low flow velocities and sedimentation of fine material).The sensitivity analysis indicates the highest parameter sensitivity for the hydraulic conductivity of the non-weathered and weathered rock and for those of the GHB upstream of the dam.Drain conductance values resulted in rather low sensitivities.However, the arrangement, connectivity and numbers of drains were not covered by the preceding sensitivity analysis.Therefore, scenarios were set up with the 3-D HGM considering (1) only the drain near to the river; (2) only the drain more distant from the river; (3) multiple drains, and (4) the connection of drains.Results show that, with the chosen boundary conditions, total drain-outflow for all simulations is generally identical.Whereas the consideration of multiple drains merely resulted in a further partitioning of outflow among the individual drains, the connection of drains resulted in an altered distribution of outflow within the drains (cf.Fig. 10).The setup of a transport model with the data from the 1996 dye tracer test indicates that the resulting breakthrough with primary and secondary maximum in some of the observation wells can be simulated adequately only with connected drains.As connectivity is determined by the stage of karst evolution, the scenarios with two drains which are unconnected or connected are discussed together with the results from the transient 3-D HGM and 2-D KEM (see below).Results from the temperature data analysis reveal the transient character of river water infiltration and allowed time-dependent conductance values to be provided which could be incorporated into the transient hydrogeologic model (Appendix A).

3-D Hydrogeologic Modeling (3-D HGM)
Figure 8 illustrates the simulated flow regime at average discharge before and after the construction measures.The flow regimes for both model scenarios clearly show the influence of the dam structure.The effect of the drains is striking, as they focus flow paths.The gradient in the non-weathered rock is steeper than in the weathered Gipskeuper and the Quaternary cover.For the simulation after the construction measures, the sealing pile wall was integrated and its hydraulic conductivity calibrated by minimizing the divergence between observed and calculated heads in OW3.As an effect of the sealing pile wall, backwater can be observed towards the river.Water budgets across model boundaries and through defined cross sections are summarized in Table 2.As an effect of the sealing pile wall, flow through Zone 2 (see Fig. 5) after the construction measures is 1/10 of the original flow around the dam.Calibrated residual flow through the sealing pile wall reaches 1.8 ls −1 .
In Fig. 9 observed and calculated groundwater heads of 4 selected observation wells are compared for a 465-day period (23 August 2006 to 30 November 2007).OW12 and 13, both located at the upstream river bank (Fig. 2), illustrate the very good agreement of observed and calculated groundwater heads as a result of the simulated river water infiltration upstream of the dam.This result confirms applying transient conductance values derived from the temperature data analysis (Appendix A).Whereas OW20 is located at the downstream river bank and indicates groundwater exfiltration into the river, OW10 is located centrally within the model domain and is a sign of the more distant flow around the dam.In Fig. 10 water budgets for the GHB up-and downstream of the dam are illustrated together with model outflow through the drains.Model inflow is dominated by the GHB upstream, indicating river water infiltration; model outflow is dominated by the GHB downstream, an indication of groundwater exfiltration (an approximate measure for the diffuse component of flow) and through the drains (approximate measures for the conduit component of model outflow).Model inflow (GHB upstream) generally equals model outflow (GHB downstream and drains).The figure also shows the relative contributions to the flow systems (diffuse and conduit model outflow; Drains 1 and 2) before and after the major flood event.Flow through the drains amounts to approximately 64% of total model outflow before the major flood event, ranging from 12 ls −1 at low to 66 ls −1 at high river discharge.After the major flood event, flow through the drains amounts to approximately 80% of total model outflow, ranging from 30 ls −1 at low to 70 ls −1 at high river discharge.Model outflow by the GHB downstream is comparatively balanced and ranges from 6 ls −1 at low to 12 ls −1 at high river discharge.Whereas the diffuse component amounts to approximately 36% before the major flood event, after the event it only amounts to approximately 20% of total model outflow.The data illustrate that in general during flood events, the relative contributions of the conduit component of model outflow increases.After moderate to medium-scale flood events, the relative contributions of flow components return to ratios observed before the events.However, major episodic flood events, as the 300-year flood of 9 August 2007, can considerably change the relative amounts of flow components also after the event.The results of long-term data monitoring will show if and when the relative amounts of the two model outflow components will return to ratios observed before the major flood.
Figure 10  the scenario with unconnected drains before the major flood event is dominated by Drain 2 amounting to 77% of total outflow through the drains.The dominance of Drain 2 can be explained by its location.It lies within the bedrock troughs formed by the ancient river course resulting in a steep hydraulic gradient from the surrounding aquifer to the drain feature (Fig. 5).Simulated outflow through Drain 1 before the major flood event is comparatively low and amounts to only 23% of total outflow through the drains.However, during flood events outflow through Drain 1 can have the same order of magnitude as outflow through Drain 2. This change in relative contributions of the two drains can be explained by the vicinity of Drain 1 to the infiltrating river.
After the major flood event simulated outflow for both drains rises significantly.The elevated outflow through Drain 1 after the major flood event results in a redistribution of the relative amounts of simulated outflow through the drains (see explanation above).The distribution of model outflow through the drains for the scenarios with connected and unconnected drains shows similar results.However, due to the connection more water is diverted to Drain 2. Before the major flood event outflow through Drains 1 and 2 amounts to 15% and 85%, respectively.After the major flood event drain outflow through Drains 1 and 2 amounts to 27% and 63%, respectively (see explanation above).

2-D Karst Evolution Model (2-D KEM)
The time interval investigated with the 2-D KEM extends from 1890 to 2007.This time period is characterized by anthropogenic alterations to the system, including the construction of the river dam in the 1890's and technical measures to prevent further subsidence on the highway in 2006 and 2007, which considerably changed the evolution of the aquifer.The evolution of the weathered rock zones used as initial conditions for the 2-D KEM can be attributed to the natural karstification beneath the ancient river bed for time periods between several hundreds up to a few thousands of years.
Preliminary scenarios and sensitivity analyses with the 2-D KEM allowed the definition of the sufficient model resolution.These scenarios demonstrated that, for a network having 100 cm spacing, the chosen size of the modeled domain is sufficient and that a time span of several hundreds years is a reasonable time scale for karst evolution within the investigation area and with regard to the hydraulic structure.Whereas calculations with uniform aquifer properties and simple boundary conditions were carried out mainly for the illustration of system behavior and processes, statistical approaches used in more complex scenarios allowed the heterogeneity of the system to be described in more detail.In terms of the heterogeneity of fracture development within the system, statistical distribution of aperture width is more realistic, resulting in a more dendritic and diffuse distribution of fracture development.
In the following, the results of one scenario are discussed, including (1) statistical distributions of hydraulic conductivity and solubility zones, based on information on locations where caves were encountered and supplementary grout mixture was injected, as well as where ERT resulted in low resistivity, and (2) regional hydraulic gradients (Fig. 7).
Figure 11 shows the statistical distributions of hydraulic conductivities for the non-weathered and weathered rock used for the setup.The development of aquifer properties from the non-weathered to the weathered state for the period of natural karstification can be derived from the shift of the statistically distributions for hydraulic conductivities and solubilities.
For the generation of statistical distributed properties of the aquifer (aperture width and solubility), spatial correlation lengths were incorporated which characterize the geometric anisotropy of the sedimentary units.Spatial correlation lengths of 1/10 of the longitudinal (i.e., 46.5 m S-N) and lateral (i.e., 13 m W-E) model domain were chosen.Values for the solubility of the non-weathered and weathered zones were averaged from borehole descriptions, resulting in 15% of soluble fractures for the weathered and 40% of soluble fractures for the non-weathered zone (Fig. 7). Figure 11 also illustrates the development of hydraulic conductivities from the initial state at 0 years (i.e. the 1890's) to 100 years (i.e. the 1990's).Note the logarithmic scale of conductivities on the y-axis.While the gross distribution of conductivities does not change significantly, the resulting curve after 100 years illustrates: (1) a slight shift from lower to higher conductivity values, and (2) the development of a few fractures with very high conductivities.Comparing this to the results from simulated model outflows, it is obvious that the occurrence of a limited number of highly permeable structures and their interconnection can dominate the flow processes.
Figure 11 illustrates the development of leakage around the dam for the 100-year time period from the 1890's to the 1990's.The graph shows a stepwise progression, while modeled outflow increases more rapidly at the beginning of the modeled time period, as remaining soluble zones within the weathered rock are dissolved.Subsequently, time spans between single steps increase.The single steps can be interpreted as local breakthrough events that lead to an abrupt increase of outflow.A steady increase in outflow can be observed between the single steps.As for the 2-D KEM, average boundary conditions are chosen, modeled outflow after 100 years evolution time, resulting in approximately 24 ls −1 , can be compared to the lower values of modeled outflow resulting from the 3-D HGM.Outflow simulated with the transient 3-D HGM resulted in values between 6 and 13 ls −1 for the GHB downstream and between 12 and 77 ls −1 for the Drains (cf.Fig. 10).Hence, modeled outflow with both modeling approaches during low to medium hydrologic boundary conditions are in the same order of magnitude.
The results from the 2-D KEM clearly illustrate that the occurrence of gypsum within the non-weathered and weathered rock determines karstification and the development of connected percolation pathways necessary for breakthrough from infiltration locations to base levels.Moreover, local breakthrough events lead to localized subsidence events as has been observed within the real world.The outflow progression, resulting from the 2-D KEM and heterogeneously distributed solubility, can illustrate patterns of karstification within gypsum rock.

Integration of modeling approaches
An iterative integration and combination of investigative methods together with the analysis of multiple data sets of varying quality considerably improved the description of the gypsum karst system within the investigation area.Associated uncertainties in data interpretation and numeric modeling were approached by (1) classifying data quality (soft and hard data); (2) parameter sensitivity analysis, and (3) scenario development and evaluation.As part of the sensitivity analysis, the most relevant parameters governing the system were determined using the two modeling approaches, which are the hydraulic conductivities of the rock and the conductance of the GHB upstream of the dam.On the basis of calculated scenarios, variable geological, hydrological and geotechnical boundary conditions that influence the flow regime were evaluated.
Boreholes provided soft data and general lithostratigraphic information, including details about the vertical extension of the weathered Gipskeuper and significant permeable zones, as well as already developed voids.Ever since the 1990's, the drilling of several boreholes has left a stratigraphic connection and locally stimulated favored karstification.Initial, coarse cross sections could be developed using the information from the national geological map.Additional soft data from lithostratigraphic information was obtained from the reports made during the installation of the piles, resulting in more precise cross sections of the investigation area.Hydraulic links within the investigation area were confirmed by a dye tracer test with groundwater flow velocities typical for conduit systems.These results indicate that the karst system is already well developed, whereas solution conduits developed along a system of fractures and interconnected joints, suggesting a three-dimensional conduit network.
Results from surface and underwater ERT measurements, taken at different hydrologic and geotechnical boundary conditions, both before and after the construction measures, provided soft data and allowed the description of (1) preferential flow in the shallow subsurface; (2) zones that are related to groundwater flow around the dam, including flow dynamics; (3) zones that are related to groundwater flow beneath the dam; (4) drainage phenomena of karst features such as voids and conduits; (5) the weathering horizon within the Gipskeuper; (6) near-surface faults and fracture zones, (7) buried paleochannels, and (8) weathered zones within the Gipskeuper beneath the river upstream of the dam and below the river sediments.Due to the multiple data sources of varying quality and hydraulic data from high-resolution 3-D HGM, it was possible to partially eliminate ambiguity in data interpretation and to describe the relationship between the different observed features in a spatial context (Epting et al., 2009).
Mass balances allowed the estimation of the amount of gypsum removed from the system over the last 100 years, Figure 12 illustrates how simulated aquifer heterogeneities from the 2-D KEM are integrated into the 3-D HGM.For the model layers 2 to 4 of the 3-D HGM, separate karst evolution models were set up and calculated.In order not to lose connectivity of the simulated developed karst features, the original model resolution of 5 by 5 m had to be refined to 1 by 1 m.Hydraulic boundary conditions of the different layers for both modeling approaches correspond to each other.For model layer 1, a statistical distribution for the non-weathered rock of the GHB downstream and a uniform value for the Quaternary cover was chosen.For model layer 5, a statistical distribution of non-weathered rock was generated.As the modeling domain of the 3-D HGM also comprises weathered zones beneath the river upstream of the dam (see Figs. 6,7,8) that were verified by ERT and that are not covered by the 2-D KEM, data had to be interpolated to these zones.This was achieved by: (1) using the existing evolution patterns generated with the 2-D KEM; (2) delimiting weathered zones beneath the river by the more resistant Schilfsandstein, and (3) assuming that weathering was intensified beneath the abandoned old meandering river bed and within zones of rock weakness (faults, fractures, etc.), which also resulted in the bulge located in the southern modeling domain (see Figs. 7,8).
Figure 13 illustrates the simulated flow regime at average discharge for the model with heterogeneous hydraulic conductivities.As for the model with uniform hydraulic conductivities, the flow regimes clearly show the influence of the dam structure.The gradient in the non-weathered rock is steeper than in the weathered Gipskeuper and the Quaternary cover.Generally the progression of hydraulic heads is comparable to the uniform model (cf.Fig. 8).However, hydraulic heads are more undulating and flow paths are focused to high conductivity zones.Due to the abundance of very low conductivity values and the predetermination of previously calibrated boundary conditions, calculated heads are generally too high compared to the measured ones.However, water budgets through defined cross sections of the 3-D HGM with uniform and heterogeneous hydraulic conductivities show a good agreement (Table 2).While calculated flow budgets through the GHBs and beneath the dam are practically identical, flow budgets around the dam for the heterogeneous model are about half the size compared to those of the uniform model.This discrepancy can be explained again by the abundance of very low conductivity values and the predetermination of previously calibrated boundary conditions within the heterogeneous model.

Conclusions
The applied concepts and methods have significantly contributed to a better understanding of the hydraulics and evolution of the karst system.The approach was illustrated by means of the following procedures: (1) determination of the extent of weathered and non-weathered rock and Comprehensive studies of transient 3-D HGM facilitated the evaluation of the relevant groundwater hydraulics and revealed the dynamic character of the flow regime during low frequency flood events, including river infiltration and diffuse and conduit components of model outflow.The magnitude of the calibrated parameters corresponds to regional hydrogeological investigations and field experiments.This indicates that calculated flow paths and flow budgets through defined zones, and especially their proportions, are plausible.Temperature measurements and the applied heat pulse method (Appendix A) offer an attractive approach for monitoring time-variant infiltration rates through losing stream-reaches, and results can be incorporated in hydrogeological models to describe transient hydraulic conductance.
3-D HGM facilitates current state descriptions of karst systems and provides sufficient information for: (1) estimating the transient composition of water budgets; (2) describing the transient character of the flow regime, and (3) simulating and evaluating short-term impacts on processes such as those which occur during episodic flood events.
However, long-term evolution of karst systems and prognosis could only be achieved by setting up models that account for the change in hydraulic properties with time.The results presented in this paper show how 2-D KEM can be calibrated to describe the current state of karst systems, and can be used for prognosis of system development and subsidence hazard assessment in the near future.2-D KEM results illustrate that the fraction of gypsum within the soluble Gipskeuper can inhibit karstification and the development of connected flow pathways necessary for system breakthrough.Moreover, local breakthrough events can lead to localized subsidence events as can be observed in the field.Time scales for the evolution of the present karst system could be estimated.
The application of different modeling techniques allowed specific aspects of the hydrologic processes to be represented and the required level of model complexity to be defined.Results from the independent approaches enabled the identification of sets of parameters for which the system behavior is satisfactorily described.
Generally, the transferability of the two approaches could be confirmed.The strengths of each individual model could be exploited.
The results of both modeling approaches were evaluated and interpreted continuously.This allowed the identification of sources of associated uncertainties and the determination of relevant parameters governing the processes within the complex geological settings, including varying properties of the formations with regard to hydraulic conductivity and solubility.It is suggested that a significant reduction in the uncertainty of modeling karstic environments can be achieved by an appropriate, complementary combination of modeling approaches viewed as a multi-model ensemble (cf.Makropoulos et al., 2008).
In addition, the development infrastructures (e.g., railways; Gutiérrez, 1996;Guerrero et al., 2008) in areas prone to subsidence requires special sets of rules and regulations (e.g.prohibiting the connection of aquifers) to minimize potential problems from present and future development.The described models require a step-wise approach; can have a predictive character and can be the basis for the development of effective long-term strategies within transient hydrogeological environments.As mentioned above, an important assumption in the heat pulse method is that the rate of thermal flux is controlled by the downward advection of the surface water and is therefore a proxy for the infiltration rate q i .Although temperature patterns strongly suggest advective-dominated heat transfer, the relative importance of conduction and convection in the substrate of the losing reach was assessed based on the calculated Peclet number.Following the approach of Silliman et al. (1995), the dimensionless Peclet number can be determined by: where the thermal diffusivity D is given by: K e is thermal conductivity (1 Jm −1 s −1• K −1 ) and l is the characteristic length and set as 1 m.Determined flow velocities range between 8.7E-06 and 7.2E-05 ms −1 , resulting in Peclet numbers for the losing reach of 1.02 to 8.48.Hence, the Peclet number is greater than 1.0 and advection dominates.Table A1 summarizes input parameters used and the results of the temperature data analysis.q i ranges from 4.7E-07 to 3.9E-06 ms −1 , which is in good agreement with other determined q i of Swiss rivers (Höhn, 2002;Huggenberger et al., 2006).
Since August 2007, generally higher conductivity values have been observed indicating enhanced river water infiltration after the 300-year flood of 9 August 2007.C d obtained from the temperature analysis resulted in 5.8E-07 s −1 before the major flood event, which is in very good agreement with the calibrated C d average value of 2.7E-07 s −1 for the GHB upstream of the dam of the steady state 3-D HGM (situation on 7 February 2006).After the flood event, C d resulted in an average value of 1.1E-06 s −1 .
The primary result of the heat pulse method is the qualitative identification of the gaining and losing reaches of small rivers.Transformations in the hydraulic properties of the streambed, as caused by flooding events, can be captured by applying this method.Furthermore it could be observed that gradients in the maximum temperatures are generally higher than those in the minimum temperatures (cf.Constantz and Thomas, 1996).A correlation between the magnitude of maximum or minimum temperature gradients and low and high flow could not be found.
The ranges in values for hydraulic aquifer properties determined using the heat pulse method is much smaller than the range obtained using Darcy-based methods.

Modeled lithological volumes and drain diameters
Measured values for Calcium and Sulfate at the groundwater outlets represent an average outflow concentration of 256.4 and 148.2 mgl −1 , respectively.Infiltrating river water has Calcium and Sulfate concentrations of 92.2 and 9.6 mgl −1 , respectively.This equals approximately 300 mgl −1 gypsum removed from the system.Modeled outflow at the GHB downstream and in the drains resulted in approximately 20 ls −1 .Consequently, 6 gs −1 gypsum are removed from the system, which means a total of approximately 500 kgd −1 , 200 ta −1 and 20 000 t in 100 years.Assuming a density of the gypsum of 2.3 to 2.4 gcm −3 would result in approximately 8500 m 3 gypsum removed from the system over the last 100 years.This represents a cube with approximately 20 m edge length.However, system outflow increased over the last 100 years, and assuming an outflow of approximately 20 ls −1 over the whole time period is not justified.Incorporating modeled outflow from the 2-D KEM (Fig. 11) would result in approximately 15 000 t or 6400 m 3 gypsum removed from the system over the last 100 years.
The solid representing the weathered gypsum formation of the diffuse flow system within the 3-D HGM has a volume of approximately 200 000 m 3 .If 6400 to 8500 m 3 gypsum was removed from the above volume, the fraction of Calcium-sulfate minerals of the weathered gypsum formation was reduced by 3 to 4% over the last 100 years.The concentration of gypsum removed from the system is assumed to be constant and effects during flood events are not considered.For 5 borehole profiles, the fraction of Calcium-sulfate minerals within the non-weathered and weathered Gipskeuper formation was determined ranging between 30-50% and 5-15%, respectively.Hence, the reduction of Calcium-sulfate minerals lies within the determined values.
Drain diameters can be calculated by applying the transformed Darcy-Weisbach equation (Bobok, 1993): where d is the drain diameter, λ is the friction factor, l is the length of the pipe, h is the head difference within the drain, u=4Q/(π d 2 ) is the average velocity with Q presenting discharge, and g is the earth's gravitational acceleration.
The friction factor depends on the velocity in the drain via the Reynolds number Re=ud/v, with v the kinematic viscosity of water (1E-06 m 2 s −1 ).For Re>2300 turbulent flow is assumed.With regard to the measured flow velocities of around 100 md −1 , turbulent flow would start with drain diameters greater than 1.8 to 2.3 m.Although voids with a maximum height of up to 2.7 m were detected during drilling, it is assumed that there is no connected conduit system exceeding the drain diameters previously cited, that would be necessary for turbulent flow.Additionally, most detected voids were filled.Consequently, laminar flow is assumed and the Hagen-Poiseuille equation can be applied.The friction for laminar flow is calculated as: This expression can be substituted in the equation used for calculating drain diameters, and together with the simulated discharge Q passing through the drain cross section the drain diameter can be calculated as: Equivalent drain diameter along the entire length of the modeled drain was calculated for Drains 1 and 2. For the scenarios with unconnected drains, calculations before the major flood event resulted in diameters for Drains 1 and 2 of 0.11 m and 0.15 m, respectively, and of 0.16 m and 0.18 m, respectively, after the major flood event.This means that the cross section for flow within Drains 1 and 2 increases by 31% and 17%, respectively, after the major flood event.For the scenarios with connected drains, calculations before the major flood event resulted in diameters for Drains 1 and 2 of 0.09 m and 0.16 m, respectively, and of 0.15 m and 0.18 m, respectively, after the major flood event.This means that the cross section for flow within Drains 1 and 2 increases by 40% and 11%, respectively, after the major flood event.
Diameters calculated for both scenarios illustrate the increase in the cross section for flow during flood events, especially for Drain 1.As the determination of h is uncertain (as observed in Oswald and Kinzelbach, 2004;Konz et al., 2009), the influence of lower h was investigated.While it is obvious that lower h result in an increase in calculated drain diameters, relative changes in diameters before and after the major flood event remain unaffected.The increase of the cross section for flow after the major flood event is interpreted as follows: During the flood event the hydraulic head forced water through the cavities and the clay fillings eroded.

Fig. 4 .
Fig. 4. River dam at low river discharge (9 m 3 s −1 ; 16 June 2006).Note, as evidence for the subsidence of the dam structure, that the crest is only overflowed on its left side (a diagram is given in Fig.2).

Fig. 7 .
Fig. 7. Geometry and hydraulic boundary conditions of the horizontal 2-D KEM superimposed on the model geometries of the 3-D HGM.Illustrated also are zones for weathered and non-weathered aquifer properties as well as locations with information from ERT, caves and injections of supplementary grout mixture.

Fig. 8 .
Fig. 8. Visualization of hydraulic heads and particle tracks in model layer 2 (0.1 m resolution).Left: groundwater flow regime at average river discharge before construction measures (7 February 2006).Right: groundwater flow regime at average river discharge after construction measures (15 January 2008).

Fig. 9 .
Fig. 9. Comparison of observed and calculated groundwater heads for 4 selected observation wells (see Fig. 2 for location of observation wells).

Fig. 11 .
Fig. 11.Upper graph: simulation results of karst evolution illustrated by the change of hydraulic conductivity distribution from the initial state to 100 a. Lower graph: simulated model outflow.

Fig. 12 .
Fig. 12. Visualization of heterogeneous hydraulic conductivity distributions transferred from three realizations of the 2-D KEM to 3-D HGM (Layers 2 to 4).Aquifer properties distributions for Layers 1 and 5 are described in the text.

Fig. 13 .
Fig. 13.Visualization of heterogeneous distribution of aquifer parameters derived from the 2-D KEM (left) as well as hydraulic heads (0.1 m resolution) and particle tracks in model layer 2 before the construction measures (right).

Table 1 .
Initial parameters, calibration range and calibrated hydraulic parameters for the extended model with drains as well as normalized composite scaled sensitivities for the calibrated parameters of all investigated scenarios.

Table 2 .
Water budgets (in ls −1 ) across model boundaries and defined zones of the cross section (see Figs.5 and 6for location of boundaries and zones for water budgets).

Table A1 .
Results from one-dimensional heat flux method.