Influence of aquifer heterogeneity on karst hydraulics and catchment delineation employing distributive modeling approaches

Introduction Conclusions References


Introduction
Karst aquifers are strongly heterogeneous systems due to a local development of large-scale discontinuities such as conduit systems.This heterogeneity also causes a large anisotropy in the hydraulic parameter field.Conceptually, karst aquifers can be described as dual-flow systems consisting of a fissured matrix with a relatively low hydraulic conductivity and highly conductive karst conduits (Liedl et al., 2003).A characteristic attribute of many karst aquifers is their high discharge focused to large springs.This makes them especially interesting as drinking water resources.However, the delineation of catchment areas of karst springs is still a challenge because of the usually unknown location of large-scale heterogeneities, such as karst conduits, within the aquifer.Common approaches for catchment delineation in porous aquifers like the mapping of geomorphological and topographical features and water balance approaches (Goldscheider and Drew, 2007) are only of limited use in karst systems.Delineating catchment areas from hydraulic head contour lines requires an observation well network, which covers the highly conductive conduit system.On groundwater catchment scale these data are scarce in carbonate areas (Sauter, 1992).Artificial tracer tests provide information about point-to-point connections, but the practical restrictions of tracer investigations prevent using them for completely defining the catchment area.In addition, catchment areas may change under different hydrological conditions further complicating the issue.
Numerical groundwater flow simulations are processbased tools that can be used for combining results from different investigation methods (Geyer et al., 2013) and for augmenting them with physical equations (Birk et al., 2005).There are numerous simulation approaches, which are applicable for karst aquifers.Single continuum models assume the aquifer to be a porous medium that can be divided into representative elementary volumes (REV) (Bachmat and Bear, 1986).The dual flow characteristics of karst aquifers are Published by Copernicus Publications on behalf of the European Geosciences Union.

S. Oehlmann et al.: Influence of aquifer heterogeneity on karst hydraulics and catchment delineation
directly addressed by hybrid or double continuum modeling approaches.Double continuum models simulate groundwater flow in two separate overlapping continua: a matrix continuum and a conduit continuum, linked via a linear exchange term (Teutsch, 1989;Mohrlok and Sauter, 1997).Hybrid models include the spatial distribution of local discrete pipe elements representing the major karst conduits coupled to a matrix continuum which represents the properties of the low permeability fissured matrix blocks (Liedl et al., 2003;Birk et al., 2005).Due to the required detailed information and the relatively high numerical effort, the application of hybrid modeling approaches to real karst systems is rare (Reimann et al., 2011a).The highest accuracy regarding the description of aquifer heterogeneities is achieved by discrete multiple fracture set models which represent the fissured system as well as the conduit system as a set of discrete fissures.Due to the intense investigation effort required for characterizing the discrete pathways they are practically not applicable for catchment studies (Teutsch and Sauter, 1991).Thus, the question which degree of complexity within the numerical model is necessary for achieving the aim of the investigation is of primary importance since more complex models require more specific information about the model area and higher numerical effort.
This work analyzes how distributive numerical models can be used to support the delineation of catchment areas of karst springs.The proposed novel approach is illustrated using a karst area in southwestern Germany.It is based on the evaluation of the influence of different types of aquifer heterogeneity on the karst flow system.More specifically, the interdependencies between hydraulic head distribution, hydraulic parameters and spring discharges are examined.For this purpose, a homogeneous continuum model and hybrid modeling approaches for flow simulation of a large-scale karst system were set up employing the finite element simulation software Comsol Multiphysics ® .These two different modeling approaches were chosen since the geometry of the highly conductive conduits was of special interest in this study because of their potential impact on the delineation of the catchment areas.Simulating the conduit geometry with the single continuum approach would have required intense meshing along the karst conduits needing a very flexible mesh and being numerically highly demanding.Steady state flow equations were implemented for both model types.The three-dimensional geometry of the aquifer system was geologically modeled with the software Geological Objects Computer Aided Design ® (GoCAD ® ) and transferred to the Comsol ® software.

Methods and approach
Comsol Multiphysics ® is a software that conducts multiphysical simulations using the finite element method (FEM).The different physical properties and equations are stored in dif-ferent modules, which can be coupled and adapted as required.The interfaces used in this work belong to the Subsurface Flow Module, which provides equations for modeling flow in porous media, and to the basic module.The basic module includes interfaces, where mathematical equations can be defined by the user and employed for any physical application.This concept is described in more detail for scenario 3 (Sect.2.3).All simulations were performed in the stationary mode, thus neglecting storage effects.Simulations were performed three-dimensionally.To examine the effects of different types of heterogeneity several scenarios were set up including more and more characteristic features of karst catchments.Figure 1 schematically shows the simulated scenarios.Catchment areas were derived by importing the simulated water tables from Comsol ® to ArcGIS ® 10.0 and using the default hydrology tools.Generally, those are used for deriving catchment areas from topographic lines.Since the concept of water flowing towards the lower potential is true for groundwater as well as for surface water, they can be likewise used for delineating groundwater catchments from groundwater contour maps.

Scenario 1
Scenario 1 simulates a completely homogenous case.It takes into account the thickness of the aquifer and boundary conditions given by rivers and surface water divides.Recharge and hydraulic conductivity were kept constant throughout the area.For the flow simulation the Darcy's Law Interface of the Subsurface Flow Module was used.It calculates the fluid pressure p (ML −1 T −2 ) within the model domain with the Darcy equation (Eq.1a and b).
In these equations Q m is the mass source term (ML −3 T −1 ), ρ is the density of the fluid (ML −3 ), K m is the hydraulic conductivity of the matrix (LT −1 ) and u the Darcy velocity (LT −1 ).g is the magnitude of gravitational acceleration (LT −2 ) and D is a unit vector in the direction over which the gravity acts.The hydraulic conductivity K m is the only calibration parameter in this scenario.

Scenario 2
Scenario 2 includes a highly conductive fracture simulated as a discrete vertical 2-D element embedded in the threedimensional continuum model.The module requires the definition of the fracture aperture d f (L) and hydraulic conductivity K f (LT −1 ) inside the fracture.Comsol ® assumes that flow processes in the fracture are basically the same as in the surrounding matrix and calculates flow along the fracture with the tangential version of the Darcy equation.The fracture flow module does not allow the application of different flow laws in the two regions.To simulate two-dimensional fracture flow the term for the fracture aperture is multiplied with both sides of Eq. ( 1): with Q f being the mass source term for the fracture (ML −3 T −1 ) and ∇ T the tangential gradient operator.The hydraulic conductivity of the fracture K f is the second calibration parameter beside the matrix conductivity K m (Eq.1b) in scenario 2.

Scenario 3
In scenario 3, highly conductive conduits were included along the positions of dry valleys, which are believed to be former riverbeds that have dried up during karstification.For these, 1-D structures are the most fitting representation.
Since the Subsurface Flow Module does not offer a similar functionality as fracture flow for 1-D elements in 3-D domains, a hybrid model was set up employing Comsol's PDE Interfaces for simulation of one-dimensional pipes.The interface chosen is called Coefficient Form Edge PDE because it allows calculations along the edges (1-D elements) of a 3-D model.The interface offers a partial differential equation (PDE) (Eq. 3) for which coefficients have to be defined.
In Eq. ( 3), c is defined as the diffusion coefficient, γ as the conservative flux source and f as the source term.By default, the source term is dimensionless.Its unit can be defined in the interface and the units of the coefficients are then calculated accordingly.v is the dependent variable in this equation.In the application using Darcy Flow, v corresponds to the pressure p (ML −1 T −2 ).The source term f equals the mass source term Q m of the Darcy equation (Eq.1a).The first of the remaining terms describes the effect of water pressure gradients, the other the effect of gravitation (compare Eq. 1b).In this case the diffusion coefficient c depends on the hydraulic conduit conductivity K c which is normalized for a unit cross-sectional area.Thus, after multiplying with the conduit area π r 2 Eq.(3) translates to Eq. ( 4).The conduit area term replaces the two missing dimensions while performing simulations in 1-D elements in a 3-D domain.
The source term multiplied with the conduit area π r 2 × Q m is equal to the mass exchange of water per unit length between the matrix and the conduit (ML −1 T −1 ).Reimann et al. (2011b) define the exchange term between a karst conduit and the rock matrix as q ex is the exchange flow per unit length (L 2 T −1 ), h ex is the difference between the hydraulic head in the matrix and the hydraulic head in the conduit (L), P ex the exchange perimeter (L) and K /b the leakage coefficient (T −1 ).For this simulation the equation was simplified by assuming the exchange perimeter equal to the pipe perimeter.Assuming there is no barrier between the conduit and the matrix, the leakage coefficient is equal to the hydraulic conductivity of the matrix divided by the theoretical distance b (L) over which the hydraulic head difference is calculated.b is kept at unit length throughout the simulation.The equation given by Reimann et al. (2011b) is multiplied by the density for obtaining the mass exchange term.The resulting exchange equation is defined in Eq. ( 6): with H c being the hydraulic head in the conduit and H m being the hydraulic head in the matrix (L).2π r is the perimeter of the pipe (L).The exchange term is used as mass flux for the matrix and as mass source for the conduits with a changed algebraic sign.Dirichlet conditions were set as boundary conditions at the springs.

Scenario 4
Scenario 4 was based on the same structure of the conduit system as scenario 3 but differed in the assumption for the conduit radius.While for scenario 3 the radius is constant within the entire conduit system, for scenario 4 a change in conduit radius towards the spring was introduced.Liedl et al. (2003) showed with their karst genesis simulations that for a conduit derived from solution processes a change in diameter is likely to occur along its extent.They introduced several simulations with different boundary conditions and derived different types of solutional widening and resulting conduit shapes.
For situations where diffuse recharge prevails, Liedl et al. (2003) showed a nearly linear increase in conduit diameters towards a karst spring.Thus, in scenario 4 a linear widening function was applied to each conduit along its arc length.
At each intersection the radii of both branches were added to account for the larger volume of water flowing there.The largest simulated radius is 4.6 m at the main karst spring.

Field site
Simulations were performed for several karst springs located at the Swabian Alb in southwestern Germany (Fig. 2).The Gallusquelle spring is the largest of the springs located within the investigation area of approximately 150 km 2 (Fig. 3).The size of its catchment area is estimated to be 45 km 2 based on a water balance approach and artificial tracer tests (Sauter, 1992) (Fig. 3).The spring is used for drinking water supply of approximately 40 000 people and has an average annual discharge of 0.5 m 3 s −1 .It is a suitable location for distributive karst modeling due to the extensive studies that have been conducted in the area before (e.g., Sauter, 1992;Geyer et al., 2007;Hillebrand et al., 2012).
Geologically the area consists of Upper Jurassic limestone and marlstone.The main aquifer is composed primarily of massive and layered limestone of the Kimmeridgian 2 and 3 (ki2/3), which are dominated by an algal sponge bioherm facies (Sauter, 1992).Beneath those rocks there are marly limestones and marlstones of the Kimmeridgian 1 (ki1) which mainly act as aquitards due to their lower hydraulic conductivity.The whole sequence dips with approximately 1.2 • to the southeast (Sauter, 1992).
Two major fault zones cross the model area.The Hohenzollerngraben strikes northwest to southeast, the Lauchertgraben crosses the area in the east striking north-south (Fig. 2).While there is no information about the hydraulic conductivity of the Lauchertgraben fault zones, the Hohenzollerngraben was crossed by tunneling work related to the construction of a regional water pipeline (Albstollen, Bodensee-Wasserversorgung).The northern boundary fault was found to be highly conductive from the significant amount of water entering the tunnel while crossing it (Gwinner et al., 1993).A high hydraulic conductivity of this zone can further be assumed from the fact that the Gallusquelle spring lies exactly at the extension of this fault where it meets the river Lauchert (Fig. 2).Parts of the area show intense fracturing.There are two main fracture directions, one with a strike of 0-30 • and one with a strike of 100-140 • parallel to the Hohenzollerngraben (Sauter, 1992).
The average hydraulic heads in the area were derived by Sauter (1992) for the period 1965-1990.The total range of hydraulic head variations during this time differs between 6 m and 20 m depending on the observation well (Sauter, 1992).The monthly rainfall varied from less than 10 mm to more than 180 mm and the annual rainfall from about 600 mm a −1 to 1200 mm a −1 .Even though these variations are high, Villinger (1977) deduced, that the boundaries of the catchment area for the Gallusquelle spring do not change significantly throughout the year.His analysis is based on equipotential maps constructed from hydraulic head measurements for high and low water levels in the area.Furthermore, several artificial tracer tests especially in the west of the area were repeated under different flow conditions and showed little to no alteration in flow directions.

Model design and calibration
The model area is constrained by fixed head boundaries at the rivers Lauchert, Fehla and Schmiecha (Dirichlet boundaries).No flow boundaries are derived from the dip of the aquifer base and artificial tracer test information (Fig. 3).The size of the model area is about 150 km 2 .The assumed catchment area of the Gallusquelle spring lies completely within the model area (Fig. 2).The positions of dry valleys were adapted after Gwinner et al. (1993).Highly conductive pipes connected to the Gallusquelle spring were implemented according to Mohrlok and Sauter (1997) and Doummar et al. (2012).The lateral positions of model boundaries, highly conductive faults and the pipe network along dry valleys were constructed in ArcGIS ® 10.0 and imported to Comsol ® as 2-D dxf-files or interpolation curves.Vertically, the highly conductive conduits were positioned approximately at the elevation of the water table simulated in scenario 1. Lacking other information, it was assumed that the homogeneously simulated water table roughly represents the one existing during the onset of karstification.Therefore, the conduits lie between 710 m and 600 m a.s.l. with a dip towards the springs.The highly conductive 2-D fracture for scenario 2 was positioned along the northern fault of the Hohenzollerngraben.The documented fault was linearly extended to the east to cross the river Lauchert at the position of the Gallusquelle spring (compare Fig. 5a and c).
Vertically the model consists of two layers.The upper one represents the aquifer.In the east it stretches from ground surface to the base of the Kimmeridgian 2 (ki2).The formation is tapering out in the west of the area but reaches a thickness of over 200 m in the east where the Gallusquelle spring is located.In the west the underlying Kimmeridgian 1 (ki1) approaches the surface until it crops out.In that region it shows karstification and thus is part of the aquifer.The depth of the karstification was derived from drilling cores.The unkarstified ki1 acts as aquitard and composes the second vertical layer of the model.It was simulated down to a horizontal depth of 300 m a.s.l.since its lower boundary is not expected to influence the simulation.The ground surface is defined by a digital elevation model (DEM) with a cell size of 40 m.The position of the ki2 base was derived from boreholes and a base map provided in Sauter (1992).Two cross sections were constructed through the model area for illustrating the geology (Fig. 4).Their positions are illustrated in Fig. 2.
Current Comsol ® software has major difficulties interpolating irregular surfaces that cannot be described by  The exact values for all model parameters are provided in Table 1.
The model was calibrated employing Comsol Multiphysics ® Parametric Sweep option, which calculates several model runs considering different parameter combinations.The focus of the calibration lay on the hydraulic head distribution.The measured hydraulic head values are long-term averages derived from twenty exploration or observation wells that were drilled within the model area (Fig. 2).
For the calibration of spring discharges five smaller springs were included in the model besides the Gallusquelle spring.Other springs within the investigation area are either very small or have not been measured on a regular basis for reliably estimating their average annual discharges.The  Gallusquelle spring and three of the other springs considered in the model calibration, the Bronnen spring, the Ahlenbergquelle spring and the Königsgassenquelle spring, are located at the river Lauchert; the Schlossbergquelle spring is situated at the river Fehla; a group of springs called the Büttnauquellen spring is located at a dry valley (Gwinner et al., 1993;Golwer et al., 1978) (Fig. 2).The Büttnauquellen springs and the Ahlenbergquelle spring probably share most of their catchment area and are likely to be fed by the same karst conduit network (Fig. 2).Localized discharge was also simulated into the rivers Fehla and Schmiecha in the west of the area, where several springs exist (Fig. 3).The highly conductive karst conduits used in the simulation connect points in the proximity of the Hohenzollerngraben with the Fehla-Ursprung spring at the Fehla and the Balinger Quelle spring at the Schmiecha.The karst conduits were identified by tracer tests (Fig. 3).However, there is not enough data for the discharges of the Fehla-Ursprung spring and the   Balinger Quelle spring to calibrate the model in this area.
Since the Gallusquelle spring is the most intensively investigated spring in the area and thus not only has the most discharge measurements but the most tracer tests as well, the main weight during calibration was laid on this spring.The simulation had to fit the Gallusquelle spring discharge within a range of 10 L s −1 , if this could be achieved with a reasonable fit for the hydraulic head distribution.
The radii of the highly conductive conduits were calibrated for a conduit volume of 200 000 m 3 for the Gallusquelle spring catchment that was deduced from an artificial tracer test (Geyer et al., 2008).For the other springs in the model area, there was no such information.For scenario 3 a systematic approach for relating the cross-sectional areas of the conduits connected to each spring to the one of the Gallusquelle spring was employed.The conduit area for each spring was defined as the area for the Gallusquelle spring multiplied by the ratio of the spring discharge to the discharge of the Gallusquelle spring.For scenario 4 where a linear relationship between the arc length and the conduit diameter was defined, it was assumed that the shorter conduits of the smaller springs lead to accordingly smaller cross-sectional areas without any further adjustments.At the springs, fixed head boundary conditions were set at the conduits.

Results and discussion
The four scenarios were evaluated and compared regarding hydraulic head distribution, hydraulic parameters, spring discharges and catchment area delineations.Figure 5 shows the simulated hydraulic head distributions for all scenarios.They are compared to a hydraulic head contour map that Sauter (1992) constructed based on field measurements (Fig. 5a). Figure 6 gives a detailed overview of the measured and simulated hydraulic heads and hydraulic gradients.The calibration parameters can be found in Table 1.Table 2 and Fig. 7 compare the simulated and observed spring discharges.

Hydraulic head distribution
The model can approximate the hydraulic head distribution in all scenarios.However, there is a significant difference of the model fit between scenario 1 with a root mean square error (RMSE) of 15 m and the best fit (scenario 4) with a RMSE of 7.7 m.Scenario 2 and 3 show similar RMSE of about 13 m.The measured hydraulic head values in the observation wells and the difference between measured and simulated head for each scenario are given in Table 3.
The measured hydraulic heads show a lateral change in hydraulic gradients.In accordance with observations in the karst aquifer of Mammoth Cave (Kentucky, USA) reported by Worthington (2009), the Gallusquelle spring catchment shows lower hydraulic gradients in the east towards the spring than in the rest of the area.This is probably caused by the higher hydraulic conductivity due to the higher karstification in the vicinity of the karst spring.After Worthington (2009) this is one of the typical characteristics of karst areas.The observation is also supported by Liedl et al. (2003) who found a widening of karst conduits in spring direction.At the field site, the steepest hydraulic head gradients were observed in the central area.
Scenario 1 cannot reproduce this behavior of the hydraulic gradient (Fig. 5b and Fig. 6a).It shows the opposite of the observed gradient distribution with steeper gradients close to the river Lauchert, where most of the springs are located.This effect usually occurs in homogeneous aquifers with evenly distributed recharge conditions.The highly conductive fracture in scenario 2 crosses the model area completely from west to east.Therefore, it mainly lowers the hydraulic head values in the central and western part, thus opposing the observed gradient distribution.In the west, where the fault starts to drain the area, its very high transmissivity leads to a strong distortion of hydraulic head contour lines (Fig. 5c).
The conduit network in scenario 3 drains the area predominantly in the central part.This results in a much lower hydraulic gradient than actually observed in the field (Fig. 5d and Fig. 6c).This effect is due to the constant and relatively high conduit diameter of 2.56 m for the conduits connected to the Gallusquelle spring.This allows large amounts of water to flow into the conduits in the central part of the catchment.
While the low hydraulic conductivity of the matrix is limiting groundwater flow in this part of the catchment, the ability of the conduits to conduct water becomes limiting close to the Gallusquelle spring and causes water to flow out of the conduits and back into the matrix.According to the classification after Kovács et al. (2005) the flow regime in this part of the model area thus is conduit influenced.Scenario 4 shows a significantly better fit for the hydraulic gradient distribution (Fig. 5e and Fig. 6d).The increase of conduit diameters towards the spring represents the higher degree of karstification and thus higher transmissivity close to the spring.As a consequence, the hydraulic gradient is steeper in the central part of the catchment than close to the spring (Fig. 5e).This corresponds to the matrix-influenced flow regime according to Kovács et al. (2005), where the discharge is controlled by the matrix rather than by the conduits.The effect is not strong enough to completely avoid an overestimation of hydraulic heads in the east and an underestimation in the central part and in the west (Fig. 6d).This leads to the assumption that the change in gradient is not purely derived from the higher karstification but that other, probably geologic factors contribute to the lateral differences in hydraulic conductivity.A more dendritic and farther extended conduit system could also lower the hydraulic head in the east.Due to the gradual widening of the conduits, the troughs in the hydraulic head contour lines are less pronounced in scenario 4 than in scenario 3 and occur further east.

Hydraulic parameters
In heterogeneous aquifers the hydraulic conductivity strongly depends on the scale of investigation of the applied method (Geyer et al., 2013).Sauter (1992) employed several approaches to determine the hydraulic conductivity in the catchment area of the Gallusquelle spring from local to regional scale.Regional methods like the gradient (Darcy) approach or the baseflow recession method average over the whole aquifer system and yielded values between 2 × 10 −5 m s −1 and 2 × 10 −4 m s −1 .Values obtained with local borehole methods such as pumping or slug tests ranged approximately from 1 × 10 −6 m s −1 to 1 × 10 −5 m s −1 .
The simulated K m values for all scenarios are well within the aforementioned ranges.The highest K m value is obtained in scenario 1 with 5.1 × 10 −5 m s −1 .This is due to the fact that K m for the homogeneous case averages the hydraulic conductivities of all structures in the area, since none of the discrete features is considered individually.Therefore, the calibrated K m is within the range given by Sauter (1992) for the regional scale.The highly conductive fracture in scenario 2 allows rapid local flow and therefore lower hydraulic heads can be achieved with a lower value for the matrix conductivity of 3.1 × 10 −5 m s −1 .This trend continues for scenario 3 and 4, where K m drops to 2.3 × 10 −5 m s −1 and 2.6 × 10 −5 m s −1 , respectively.In these scenarios the hydraulic conductivity values approach those obtained by Sauter (1992) with borehole tests, suggesting that most of the highly conductive features in the area are explicitly taken into account.
The fracture conductivity K f is introduced in scenario 2. Despite being in the typical range of literature of 2-10 m s −1 (Sauter, 1992) the obtained value of 2.7 m s −1 probably is too low, because all other karst features, which can drain water from the Gallusquelle spring catchment towards other springs, are neglected.If additional highly conductive features are included, higher fracture conductivities will be necessary to provide the observed average spring discharge of the Gallusquelle spring.This effect is partly responsible for the relatively high conduit conductivity K c of 6.5 m s −1 in scenario 3.Even though the discharge at the Gallusquelle spring is the same as well as the integrated conduit volume, the conduit conductivity of 2 m s −1 obtained for scenario 4 is significantly lower than the value of 6.5 m s −1 obtained for scenario 3.This is because the karst conduit system with constant diameter needs a higher overall transmissivity to transport the same amount of water due to limiting flow capacity of the conduits close to the spring.
The conduit diameter in scenario 3 corresponds to a representative constant diameter for the Gallusquelle spring.Birk et al. (2005) used artificial tracer tests for calculating the representative diameter.The authors calculated a diameter of about 5 m, which is higher than the 2.56 m simulated with scenario 3.This is probably due to the fact that these tracer tests were conducted approximately 3 km northwest of the spring while in the model the conduits extend approximately 10 km to the northwest.Thus, this supports the idea that the diameters of the conduits closer to the spring are higher than those farther away (see Sect. 2.4).

Spring discharge
Scenario 1 fails to simulate the locally increased discharge at the karst springs (Table 2).Since there are no areas of focused flow, there is only diffuse groundwater discharge into the rivers, mainly the Lauchert.In scenario 2 fracture flow along the fault allows the simulation of increased discharge at the Gallusquelle spring (Table 2).The other springs that were not connected to highly conductive elements show no locally increased discharge (Table 2).The slightly raised discharge of the Schlossbergquelle spring compared to scenario 1 results from generally increased water flow into the river Fehla, not from locally raised discharge at the spring location.The local discharges at all springs can only be represented by scenarios 3 and 4. The simulation is satisfactory for both scenarios.The simulated discharge of the scenarios is very similar for the Gallusquelle spring, the Schlossbergquelle spring and the Königsgassenquelle spring (compare Table 2 and Fig. 7).The fit for these springs is good, even though the discharge is slightly overestimated for the Königsgassenquelle spring and underestimated for the Schlossbergquelle spring.Since the Schlossbergquelle spring is the only spring included at the river Fehla and no registration of discharge values of the river itself was conducted, it cannot be distinguished, if the underestimation at the Schlossbergquelle spring is due to an inexact karst conduit network or to an underestimated discharge into the river.For the Bronnen spring, different results can be observed for the two scenarios.While scenario 3 has a very good fit, scenario 4 underestimates the discharge.This suggests that the conduits leading to the spring are assumed too short in the simulation leading to underestimated conduit diameters in scenario 4.
The most pronounced difference between the two simulations occurs at the Büttnauquellen springs and Ahlenbergquelle spring.Both simulations underestimate their discharge with a significantly stronger underestimation in scenario 4 (Fig. 7).This is probably due to the simplified approach of treating them like a single spring and attaching them to the same conduit.While the Ahlenbergquelle spring is perennial, the Büttnauquellen springs are intermittent.This suggests that there are karst conduits in at least two different depths and thus that the representation with a conduit network in a single depth is not adequate.A too short conduit system with too little side branches has a stronger impact on scenario 4 because of the dependence of diameters on the total length and amount of intersections leading to a stronger underestimation of conduit volumes than in scenario 3.

Catchment area delineation
The spring catchment areas were delineated according to the hydraulic heads within the matrix.For the delineation a bending of contour lines towards the springs is required, meaning they can only be generated with localized discharge at the spring positions.Therefore no catchment areas can be delineated in scenario 1.In scenario 2 a catchment area for the Gallusquelle spring can be delineated.It has approximately the size that can be expected from water balance calculations, but does not include all injection locations of tracer tests with recovery at the Gallusquelle spring.Since the hydraulic conductivity of the fault is assumed to be constant, it receives most of the inflow in the west and cannot receive more water close to the spring.Thus, the catchment area mainly includes the western part of the model area (Fig. 5c).
In scenario 3 catchment areas can be simulated for the Gallusquelle spring and for the Büttnauquellen springs and Ahlenbergquelle spring (Fig. 5d).The unusual looking shape of the areas is caused by the filling of the conduits with water in the west of the model domain which prevents drainage of the fissured matrix by the conduit system in the east of the area.Therefore the Gallusquelle spring mainly receives water from the western part of the area, where its conduits drain enormous water volumes due to their relatively large diameter.Due to outflow of water into the matrix in the east, only part of the water from the shown catchment area is transported to the springs.In the west it can be observed that the catchment areas of the Gallusquelle spring and the Büttnauquellen springs and Ahlenbergquelle spring reach across karst conduits leading to other springs (Fig. 5d).In this case the catchment areas of the springs overlap.The catchment areas were constructed in 2-D according to surface values, so that they envision the flow above the smaller conduits in the west.In the east it can be observed that the catchment areas do not include all parts of the respective karst conduit network.In these areas the conduits cannot accommodate more water and outflow occurs.The catchment area for the Gallusquelle spring that was delineated in scenario 3 includes all but one tracer test conducted.The Gallusquelle spring drains nearly all water from the springs at the river Fehla.The hydraulic heads in the west are lowered leading to influent flow conditions along parts of the western Fehla.This contradicts the development of several springs in this area and makes this scenario highly unlikely (compare Fig. 3).
Scenario 4 is the only simulation leading to reasonable results regarding the catchment areas (Fig. 5e).The size of the Gallusquelle spring catchment area is in accordance with water balance calculations and includes all tracer tests conducted in the catchment of the Gallusquelle spring.The size of the catchment area for the Büttnauquellen springs and Ahlenbergquelle spring is probably underestimated due to the underestimation of spring discharge (Table 2).Since the underestimation is more pronounced for scenario 4 than for scenario 3, the catchment area is significantly smaller (compare Fig. 5d and Fig. 5e).A small overlap of catchment areas can still be observed in the west but in scenario 4 the Gallusquelle spring only drains small amounts of water from the western part, so that the western Fehla is completely efflu-ent.Since the simulation was performed stationary, the delineated catchment areas are only valid for the average hydraulic head distribution.As known from literature (Sect.3) they should be representative for the usually observed variations in the Gallusquelle spring area.For reliably simulating possible shifts in the catchment areas during extreme flow conditions, more detailed information on recession behavior of the aquifer and lateral and temporal recharge distribution should be included.This is beyond the scope of this paper.
For the smaller springs, no catchment areas could be generated in either of the scenarios.They produce a very small ratio of the total discharge of the model area (< 5 %) and the resolution of the simulation was not fine enough to reliably draw their catchment boundaries.

Conclusions
The results show that distributive numerical simulation is a useful tool for approaching the complex subject of subsurface catchment delineation in karst aquifers as long as effects of karstification are sufficiently taken into account.Even though the Gallusquelle spring area is significantly less karstified than for example the Mammoth Cave (Kentucky, USA) (Worthington, 2009) and does not show significant troughs in the hydraulic head contour lines, it cannot be simulated with a homogeneous hydraulic parameter field.The geometry of the conduits is of major importance for the simulation.Although the Gallusquelle spring is positioned on the linear extension of the northern fault of the Hohenzollerngraben the hydraulic conditions cannot correctly be simulated without consideration of dry valleys.For catchment delineation, the approach of using conduits with constant geometric parameters is not satisfactory, either.While it is possible to fit spring discharges with a double continuum model (e.g., Kordilla et al., 2012) or a single continuum model with a highly conductive zone with constant hydraulic properties (e.g., Doummar et al., 2012) the hydraulic head distribution and hydraulic conductivities cannot be correctly approximated with these approaches.
Using numerical models for catchment delineation allows for the combination of several methods and observations under consideration of the geological and hydrogeological properties of the area.The model can be used for advanced simulations of transient groundwater flow and transport and can also account for heterogeneous distributions of recharge or aquifer properties.It therefore represents a flexible tool for risk assessment and prediction in heterogeneous flow systems.
The uncertainty of the results depends mainly on the available input data.The modeling approach allows an integrated analysis of data from different sources.Theoretically, the method requires average annual spring discharge and hydraulic head measurements in the catchment.Nonetheless, the measurement of the discharge of several springs in the proximity of the investigated spring catchment is advisable for the simulation of catchment boundaries.In addition, deriving some knowledge about the location and properties of the karst conduit network from natural or artificial tracers, groundwater contour lines, direct investigations or the morphology of the land surface is highly recommended.
To improve simulation results, future work includes the implementation and simulation of solute transport, e.g., simulation of artificial tracer tests.Since the hydraulic head distribution and the spring discharges were found to be strongly dependent on the selected geometry of the highly conductive elements it seems unavoidable to better constrain their positions and sizes in the area.In case of the Gallusquelle spring area the smooth hydraulic gradients do not allow the localization of conduits by troughs in the hydraulic head contour lines like in some other karst areas (e.g., Joodi et al., 2010).Karst genesis simulation would provide process-based information about conduit widening towards a karst spring.Such simulations were employed for instance by Kaufmann and Braun (1999), Liedl et al. (2003), Bauer et al. (2003), andHubinger andBirk (2011).They simulate the temporal evolution of a small fracture or fracture network due to solution with coupled transport and hydraulic models.Under the constraints of recharge conditions and initial geometries they derive the conduit size distribution.A detailed overview of the basic techniques and processes is given by Dreybrodt et al. (2005).The implementation of a karst genesis module would be possible with Comsol Multiphysics ® , given sufficient input data.

Fig. 1 .
Fig. 1.Conceptual geometry of the simulated scenarios.For explanation of the flow equations see scenario description in Sect. 2.

Fig. 2 .
Fig.2.Model area, including the catchment of the Gallusquelle spring and positions of all simulated springs.The highly conductive elements feeding the Gallusquelle spring were modeled afterDoummar et al. (2012) and the ones along the dry valleys afterGwinner et al. (1993).

Fig. 3 .
Fig. 3. Top view of the model area.Tracer tests within the area are illustrated with their major and minor registration points (excluded: uncertain registrations and registration points in rivers) after information from the Landesamt für Geologie, Rohstoffe und Bergbau (LGRB).Dry valleys were simulated with ArcGIS ® 10.0 and counterchecked with field observations of Gwinner et al. (1993).

Fig. 4 .Fig. 5 .
Fig. 4. Cross sections of the study area as constructed in GoCAD ® from northwest to southeast with a vertical exaggeration of 10 : 1.(a) cross section 1 through the Lauchertgraben and the Gallusquelle spring.(b) cross section 2 through the Hohenzollerngraben, the Lauchertgraben and the Königsgassenquelle spring.

R
= groundwater recharge by precipitation, K m = hydraulic conductivity of matrix, K l = hydraulic conductivity of lowly conductive ki1, K f = hydraulic conductivity of fracture, K c = hydraulic conductivity of conduits, dz = fracture depth, dy = fracture aperture, RMSE = root mean square error for the hydraulic head distribution.

Fig. 6 .
Fig. 6.Comparison of the hydraulic head values measured in the observation wells and those simulated at the well positions.(a) after the homogeneous simulation.(b) after the simulation with fracture flow along the northern fault of the Hohenzollerngraben.(c) after the simulation with a 1-D conduit network with constant radius.(d) after the simulation with a 1-D conduit network with increasing radius.

Fig. 7 .
Fig. 7. Spring discharge: measured and simulated values using a conduit network with constant radius (scenario 3) and with linearly increasing radius (scenario 4).

the Fracture Hydrol. Earth Syst. Sci., 17, 4729-4742, 2013 www.hydrol-earth-syst-sci.net/17/4729/2013/
The 2-D element, in this case, represents a large-scale fault zone observed from geological mapping within the area of investigation.The continuum represents the fissured matrix of the karst aquifer.Groundwater flow in the fracture was simulated with

Table 1 .
Input and calibration values of the different scenarios.The root mean square error of the hydraulic head distribution is given as an index for the quality of the model fit.

Table 3 .
Measured hydraulic head values that were used for calibration.For each scenario the difference of the simulated to the measured hydraulic heads is given in meters.The positions of the observation wells are given in Fig.5a.