Reducing the ambiguity of karst aquifer models by pattern matching of flow and transport on catchment scale

transport in a highly heterogeneous karst aquifer in south western Germany. Therefore, an interface for the simulation of solute transport in one-dimensional pipes was implemented into the software Comsol Multiphysics and coupled to the three-dimensional solute transport interface for continuum domains. For reducing model ambiguity, the simulation was matched for steady-state conditions to the hydraulic head distribution in 10


Introduction
Karst systems play an important role in water supply worldwide (Ford and Williams, 2007).They are characterized as dual-flow systems where flow occurs in the relatively lowly conductive fissured matrix and in highly conductive karst conduits (Reimann et al., 2011).There are a number of process-based modelling approaches available for simulating karst aquifer behaviour.Overviews on the various types of distributed process and lumped-parameter models are provided by several authors (Teutsch and Sauter, 1991;Jeannin and Sauter, 1998;Kovács and Sauter, 2007;Hartmann et al., 2014).In most cases, lumped-parameter models are applied, since they are less demanding on input data (Geyer et al., 2008;Perrin et al., 2008;Hartmann et al., 2013;Schmidt et al., 2013).These models consider neither the actual flow process nor the heterogeneous spatial distribution of aquifer parameters, but are able to simulate the integral aquifer behaviour, e.g.karst spring responses.The spatial distribution of model parameters and state variables, e.g. the hydraulic head distribution, need to be addressed with distributed numerical models should the necessary field data be available (e.g.Oehlmann et al., 2013;Saller et al., 2013).A distributed modelling approach suited for the simulation of strongly heterogeneous and anisotropic aquifers with limited data availability is the hybrid modelling approach.The approach simulates the fast flow component in the highly conductive karst conduit system in discrete one-dimensional elements and couples it to a two-or three-dimensional continuum representing the fissured matrix of the aquifer (Oehlmann et al., 2013).Hybrid models are rarely applied to real karst systems because they have a high demand of input data (Reimann Published by Copernicus Publications on behalf of the European Geosciences Union.et al., 2011).They are, however, regularly applied in longterm karst genetic simulation scenarios (e.g.Clemens et al., 1996;Bauer et al., 2003;Hubinger and Birk, 2011).In these models not only groundwater flow but also solute transport is coupled in the fissured matrix and in the karst conduits.Aside from karst evolution such coupling enables models to simulate tracer or contaminant transport in the karst conduit system (e.g.Birk et al., 2005).In addition to serving for predictive purposes, such models can be used for deriving information about the groundwater catchment itself (Rehrl and Birk, 2010).
A major problem for characterizing the groundwater system with numerical models is generally model ambiguity.The large number of calibration parameters is usually in conflict with a relatively low number of field observations, e.g.different hydraulic parameter fields and process variables may give a similar fit to the observed data but sometimes very different results for prognostic simulations (Li et al., 2009).Especially the geometric and hydraulic properties of the karst conduit system are usually unknown and difficult to characterize with field experiments for a whole spring catchment (Worthington, 2009).With artificial tracer test data the maximum conduit volume can be estimated but an unknown contribution of fissured matrix water prevents further conclusions on conduit geometry (Birk et al., 2005;Geyer et al., 2008).It is well known that the use of several objective functions, i.e. several independent field observations, can significantly reduce the number of plausible parameter combinations (Ophori, 1999).Especially in hydrology (e.g.Khu et al., 2008;Hunter et al., 2005) and also for groundwater systems (e.g.Ophori, 1999;Hu, 2011;Hartmann et al., 2013), this approach has been successfully applied with a wide range of observation types, e.g.groundwater recharge, hydraulic heads, remote sensing and solute transport.Particularly, the simulation of flow and transport is known to reduce model ambiguity and yield information on karst conduit geometry (e.g.Birk et al., 2005;Covington et al., 2012;Luhmann et al., 2012;Hartmann et al., 2013).Usually, automatic calibration schemes performing a multiobjective calibration for several parameters are used for this purpose (Khu et al., 2008).However, for complex modelling studies calculation times might be large due to the high number of model runs needed (Khu et al., 2008) and a precise conceptual model is essential as basis for the automatic calibration (Madsen, 2003).In general, numerical models of karst aquifers are difficult to build because of their highly developed heterogeneity (Rehrl and Birk, 2010).Thus, automatic calibration procedures are better suited for conceptual and lumped-parameter models, where calibration parameters include effective geometric properties and no spatial representation of the hydraulic parameter field and conduit geometry is necessary.Complex distributed numerical approaches generally require longer simulation times due to the necessary spatial resolution.Long simulation times limit the number of model runs that can reasonably be performed and man-ual calibration based on hydrogeological knowledge is necessary (e.g.Saller et al., 2013).Therefore, applied distributed numerical models in karst systems usually focus on a smaller number of objective functions.They generally cannot simulate the hydraulic head distribution in the area, spring discharge and tracer breakthrough curves simultaneously on catchment scale.Some studies combine groundwater flow with particle tracking for tracer directions (e.g.Worthington, 2009;Saller et al., 2013) without simulating tracer transport.On the other hand there are studies simulating breakthrough curves without calibrating for measured hydraulic heads (e.g.Birk et al., 2005).For developing process-based models which can be used as prognostic tools, e.g. for the delineation of protection zones, the simulation should be able to reproduce groundwater flow and transport within a groundwater catchment.Especially in complex hydrogeological systems, this approach would reduce model ambiguity, which is a prerequisite in predicting groundwater resources and pollution risks.
This study shows how the combination of groundwater flow and transport simulation can be used not only to develop a basis for further prognostic simulations in a heterogeneous karst aquifer with a distributed modelling approach on catchment scale, but also to reduce model ambiguity and draw conclusions on the spatially distributed karst network geometries and the actual karst conduit volume.The approach shows the kind and minimum number of field observations needed for this aim.Furthermore, a systematic calibration strategy is presented to reduce the number of necessary model runs and the simulation time compared to standard multi-objective calibrations.For this purpose a hybrid model was built and a pattern matching procedure was applied for a well-studied karst aquifer system in south-western Germany.The model was calibrated for three major observed parameters: the hydraulic head distribution derived from measurements in 20 boreholes, the spring discharge of six springs and the tracer breakthrough curves of two tracer tests.

Modelling approach
The simulation is based on the mathematical flow model discussed in detail by Oehlmann et al. (2013).The authors set up a three-dimensional hybrid model for groundwater flow with the software COMSOL Multiphysics ® .As described by Oehlmann et al. (2013) the simulation was conducted simultaneously in the three-dimensional fissured matrix, in an individual two-dimensional fault zone and in one-dimensional karst conduit elements to account for the heterogeneity of the system.Results showed that the karst conduits widen towards the springs and therefore, a linear relationship between the conduit radius and the conduit length s [L] was established.Values for s start with zero at the point farthest away from the spring and increase towards the respective karst spring.
In agreement with these results and karst genesis simulations by Liedl et al. (2003), the conduit radius is calculated as where r c [L] is the radius of a conduit branch and m and b are the two parameters defining the conduit size.b [L] is the initial radius of the conduit at the point farthest away from the spring and m [-] is the slope with which the conduit radius increases along the length of the conduit s.
In the following the equations used for groundwater flow and transport are described.The subscript "m" denotes the fissured matrix, "f" the fault zone and "c" the conduits hereby allowing a clear distinction between the respective parameters.Parameters without a subscript are the same for all karst features in the model.

Groundwater flow
Groundwater flow was simulated for steady-state conditions.This approach seems appropriate since this work focuses on the simulation of tracer transport in the conduit system during tracer tests, which are ideally conducted under quasisteady-state flow conditions.Therefore, the simulations refer to periods with a small change of spring discharge, e.g.base flow recession, and are not designed to predict conditions during intensive recharge/discharge events.The groundwater flow in the three-dimensional fissured matrix was simulated with the continuity equation and the Darcy equation (Eq.2a und b).
where Q m is the mass source term [M L −3 T −1 ], ρ the density of water [M L −3 ] and u m the Darcy velocity [L T −1 ].In Eq. (2b) K m is the hydraulic conductivity of the fissured matrix [L T −1 ] and H m the hydraulic head [L].Two-dimensional fracture flow in the fault zone was simulated with the COMSOL ® fracture flow interface.The interface only allows for the application of the Darcy equation inside of fractures, so laminar flow in the fault zone was assumed.In order to obtain a process-based conceptualization of flow, the hydraulic fault conductivity K f was calculated by the cubic law (Eq.3): For groundwater flow in the karst conduits, the Manning equation was used (Eq. 4).
where u c is the specific discharge in this case equalling the conduit flow velocity [L T −1 ], n the Manning coefficient [T L −1/3 ], r c /2 the hydraulic radius [L] and dH c /dx the hydraulic gradient [-].The Manning coefficient is an empirical value for the roughness of a pipe with no physical nor measurable meaning.The hydraulic radius is calculated by dividing the cross section by the wetted perimeter, which in this case corresponds to the total perimeter of the pipe (Reimann et al., 2011).The whole conduit network was simulated for turbulent flow conditions.Due to the large conduit diameters (0.01-6 m, Sect.5) this assumption is a good enough approximation.Hereby, strong changes in flow velocities due to the change from laminar to turbulent flow can be avoided.At the same time, the model does not require an estimation of the critical Reynolds number, which is difficult to assess accurately.
The three-dimensional flow in the fissured matrix and the one-dimensional conduit flow were coupled through a linear exchange term that was defined according to Barenblatt et al. (1960) as where q ex is the water exchange between conduit and fissured matrix [L 2 T −1 ] per unit conduit length L [L], H m the hydraulic head in the fissured matrix [L], H c the hydraulic head in the conduit [L] and α the leakage coefficient [L 2 T −1 ].The leakage coefficient was defined as where 2 π r c is the conduit perimeter [L].Other possible influences, e.g. the lower hydraulic conductivity at the solidliquid interface of the pipe and the fact that water is not exchanged along the whole perimeter but only through the fissures are not considered.The exact value of these influences is unknown and the exchange parameter mainly controls the reaction of the karst conduits and the fissured matrix to hydraulic impulses.Since the flow simulation is performed for steady-state conditions this simplification is not expected to exhibit significant influence on the flow field.

Solute transport
Transient solute transport was simulated based on the steadystate groundwater flow field.COMSOL Multiphysics ® offers a general transport equation with its solute transport interface.This interface was applied for the three-dimensional fissured matrix.In this work saturated, conservative transport was simulated, with an advection-dispersion equation (Eq.7) where θ m is the matrix porosity [-], c m the solute concentration [M L −3 ], D Dm the mechanical dispersion [L 2 T −1 ] and The solute transport interface cannot be applied to onedimensional elements within a three-dimensional model.COMSOL ® offers a so-called coefficient form edge PDE interface to define one-dimensional mathematical equations.There, a partial differential equation is provided (COM-SOL AB, 2012) which can be adapted as needed and leads to Eq. ( 8) in its application for solute transport in karst conduits: where θ c is the conduit porosity which is set equal to 1, f the source term and u c [L T −1 ] the flow velocity inside the conduits, which corresponds to the advective transport component.Flow divergence cannot be neglected, as is often the case in other studies (e.g.Hauns et al., 2001;Birk et al., 2006;Coronado et al., 2007).Different conduit sizes and inand outflow along the conduits lead to significant velocity divergence in the conduit system.This needs to be considered for mass conservation during the simulation.The mechanical conduit dispersion D Dc was calculated with Eq. ( 9) (Hauns et al., 2001).
where ε is the dispersivity in the karst conduits [L].
The source term f [M T −1 L −1 ] in Eq. ( 8) equals in this case the mass flux of solute per unit length L [L] due to matrix-conduit exchange of solute c ex : The first term of the right-hand side of Eq. ( 10) defines the diffusive exchange due to the concentration difference between conduit and fissured matrix.The second term is a conditional term adding the advective exchange of solute due to water exchange.The concentration of the advective exchange c i is defined as When q ex is negative, the hydraulic head in the fissured matrix is higher than in the conduit (Eq.5) and water with the solute concentration of the fissured matrix c m enters the conduit.When it is positive, water with the solute concentration c c of the conduit leaves the conduit and enters the fissured matrix.Since one-dimensional transport is simulated in a three-dimensional environment, the left-hand side of Eq. ( 8) is multiplied with the conduit cross section π r 2 c [L 2 ].These considerations lead to the following transport equation for the karst conduits: 3 Field site and model design The field site is the Gallusquelle spring area on the Swabian Alb in south-western Germany.The size of the model area is approximately 150 km 2 , including the catchment area of the Gallusquelle spring and surrounding smaller spring catchments (Oehlmann et al., 2013).The Gallusquelle spring is the main point outlet with a long-term average annual discharge of 0.5 m 3 s −1 .The model area is constrained by three rivers and no-flow boundaries derived from tracer test information and the dip of the aquifer base (Oehlmann et al., 2013) (Fig. 1).
The aquifer consists of massive and bedded limestone of the stratigraphic units Kimmeridgian 2 and 3 (ki 2/3) (Golwer, 1978;Gwinner, 1993).The marly limestones of the underlying Kimmeridgian 1 (ki 1) mainly act as an aquitard.In the west of the area where they get close to the surface, they are partly karstified and contribute to the aquifer (Sauter, 1992;Villinger, 1993).The Oxfordian 2 (ox 2) that lies beneath the ki 1 consists of layered limestones.It is more soluble than the ki 1 but only slightly karstified because of the protective effect of the overlying geological units.In the catchment areas of the Fehla-Ursprung and the Balinger springs close to the western border (Fig. 1a) the ox 2 partly contributes to the aquifer.For simplicity, only two vertical layers were differentiated in the model: the aquifer and the underlying aquitard.
The geometry of the conduit system was transferred from the COMSOL® model calibrated for flow by Oehlmann et al. (2013).It is based on the occurrence of dry valleys in the investigation area and artificial tracer test information (Gwinner, 1993).The conduit geometry for the Gallusquelle spring was also employed for distributed flow simulations by Doummar et al. (2012) and Mohrlok and Sauter (1997) (Fig. 1).In this work, all highly conductive connections identified by tracer tests in the field were simulated as discrete one-dimensional karst conduit elements.The only exception is a connection in the west of the area that runs perpendicular to the dominant fault direction and reaches the Fehla-Ursprung spring at the northern boundary (Fig. 1).While the element was regarded as a karst conduit by Oehlmann et  2013) it is more likely that the water crosses the graben structure by a transversal cross-fault (Strayle, 1970).Therefore, the one-dimensional conduit element was replaced by a two-dimensional fault element (Fig. 1b).This leads to a small adjustment in the catchment areas compared to the results of Oehlmann et al. (2013) (Fig. 1a).While the discharge data for the Fehla-Ursprung spring are not as extensive as for the other simulated springs, it is approximated to 0.1 m 3 s −1 , the annual average ranging from 0.068 to 0.135 m 3 s −1 .The fault zone aperture was calibrated accordingly (Sect.5).
Due to a large number of studies conducted in the area during the last decades (e.g.Villinger, 1977;Sauter, 1992;Geyer et al., 2008;Kordilla et al., 2012;Mohrlok, 2014) many data for pattern matching are available even though the karst conduit network itself is not accessible.Since the groundwater flow simulation was performed for steady-state conditions, direct recharge, which is believed to play an important role during event discharge (Geyer et al., 2008), was neglected.It is not expected that recharge dynamics exhibit significant influence on the flow field during recession periods.From Sauter (1992) the long-term average annual recharge, ranges  of hydraulic parameters and the average annual hydraulic head distribution derived from 20 observation wells (Fig. 1a) are available.Villinger (1993) and Sauter (1992) provided data on the geometry of the aquifer base.Available literature values for the model parameters are given in Table 1.
The observed hydraulic gradients in the Gallusquelle area are not uniform along the catchment.Figure 2 shows a Sshaped distribution with distance to the Gallusquelle spring.The gradient at each point of the area depends on the combination of the respective transmissivity and total flow.The amount of water flowing through a cross sectional area increases towards the springs due to flow convergence.In the Gallusquelle area, the transmissivity rises in the vicinity of the springs leading to a low hydraulic gradient.In the central part of the area discharge is relatively high while the transmissivities are lower leading to the observed steepening of the gradient starting in a distance of 4000 to 5000 m from the Gallusquelle spring.Towards the boundary of the catchment area in the west the water divide reduces discharge in the direction of the Gallusquelle spring leading to a smoothing of hydraulic gradients.Geyer et al. (2008) calculated the maximum conduit volume for the Gallusquelle spring V c [L 3 ] with information from the tracer test that will be referred to as tracer test 2 in the following.Since the injection point of the tracer test is close to the catchment boundary, it is assumed that it covers the whole length of the conduit system.The authors calculated the maximum volume at 218 000 m 3 .Their approach assumes the volume of the conduit corresponds to the total volume of water discharged during the time between tracer input and tracer arrival neglecting the contribution of the fissured matrix.
The six springs that were monitored and therefore simulated are shown in Fig. 1.Except for the Balinger spring, their discharges were fitted to long-term average annual discharge data.For the Balinger spring discharge calibration was not possible due to lack of data.It was included as a boundary condition because several tracer tests provided a valuable basis for the conduit structure leading to the spring.
Tracer directions were available for 32 tracer tests conducted at 20 different tracer injection locations (Oehlmann et al., 2013).In all, 16 of the tracer tests were registered at the Gallusquelle spring.For this work two of them were chosen for pattern matching of transport parameters.Both of them were assumed to have a good and direct connection to the conduit network.Tracer test 1 (Geyer et al., 2007) has a tracer injection point at a distance of 3 km to the Gallusquelle spring.Tracer test 2 (MV746 in Merkel, 1991;Reiber et al., 2010)  spring (Fig. 1a).Due to the flow conditions (Fig. 1a) it can be assumed that tracer test 2 covers the total length of the conduit network feeding the Gallusquelle spring.The recovered tracer mass was chosen as input for the tracer test simulation.The basic information about the tracer tests is given in Table 2.
Since the tracer tests were not performed at average flow conditions, the model parameters were calibrated first for the long-term average annual recharge of 1 mm d −1 and the longterm average annual discharge of 0.5 m 3 s −1 .For the transport simulations, the recharge was then adapted to produce the respective discharge observed during the tracer experiment (Table 2).

Parameter analysis
An extensive parameter analysis was performed in order to identify parameters determining the hydraulic parameter field in the model area, as well as their relative contributions to the discharge and conduit flow velocities.The fitting parameters include the parameters controlling the respective transmissivities of the fissured matrix and the karst conduit system, i.e. the geometry and roughness of the conduit system, the hydraulic conductivity of the fissured matrix and the fracture aperture for the Fehla-Ursprung spring.Furthermore, the apparent dispersivities for the two artificial tracer tests were calibrated (Table 1).Since all model runs were performed for steady-state conditions parameters controlling the temporal distribution of recharge were not con- sidered.The parameter analysis was performed with COM-SOL Multiphysics ® parametric sweep tool, which sweeps over a given parameter range.Parameter ranges were chosen according to literature values (Table 1).For the conduit geometry parameters, lowest conduit radius b and slope of radius increase m, no literature values are available.Therefore, the ranges were chosen so that conduit volumes ranged below the maximum volume given by Geyer et al. (2008).
In addition to the variation of the fitting parameters, five basic scenarios were compared.They correspond to different conceptual representations of the area and are summarized in Fig. 3 and Table 3.Three objective functions were employed for pattern matching: spring discharge, hydraulic head distribution and flow velocities of the two tracer tests (Sect.3).The average spring discharge of the Gallusquelle spring was set by the difference between simulated and the measured discharge.A difference of 10 L s −1 was considered as acceptable.Param- eter sets, which could not fulfil this criterion, were not considered for parameter analysis.The other low-discharge and less-investigated springs (Sect.3) were used to inspect the flow field and water balance in the modelling area, i.e. they were only considered after parameter fitting to check the plausibility of the deduced parameter set.
The fit of the tracer tests was determined by comparing the arrival times of the highest peak concentration of the simulation with the measured value (peak offset).Since tracer experiments conducted in karst conduits usually display very narrow breakthrough curves, this procedure appears to be justified.The quality of the fit was judged as satisfactory if the peak offset was lower than either the simulation interval or the measurement interval.
The fit of the hydraulic head distribution was determined by calculating the root mean square error (RMSE) between the simulated and the observed values at the respective locations of the observation wells.Since the fit at local points with a large-scale modelling approach generally shows large uncertainties due to low-scale heterogeneities, an overall fit of < 10 m RMSE was accepted.Furthermore, a qualitative comparison with the hydraulic gradients in the area was performed (e.g.Fig. 2) to ensure that the general characteristics of the area were represented instead of only the statistical value.

Scenario 1 -standard scenario
In scenario 1 all features were implemented as described in Sects. 2 and 3.The parameter analysis shows that for each conduit geometry, defined by their smallest conduit radii b and their slopes of radius increase along the conduit length m (Eq.1), only one value of the Manning coefficient n allows a simulated discharge for the Gallusquelle spring of 0.5 m 3 s −1 .The n value correlates well with that for the total conduit volume due to the fact that the spring discharge is predominantly determined by the transmissivity of the karst conduit system.The transmissivity of the conduit system at each point in space is the product of its hydraulic conduc-tivity, which is proportional to 1/n, and the cross sectional area of the conduit A. Thus, to keep the spring discharge at 0.5 m 3 s −1 a higher conduit volume requires a higher calibrated n value (Eq.4).
With scenario 1 it is possible to achieve a hydraulic head fit resulting in a RMSE of 6 m that can be judged as adequate on catchment scale.Regarding the conduit geometry, a good hydraulic head fit can be achieved with small b values independently of the chosen m value (Fig. 2a).The higher the b value, the higher the m value to reproduce the hydraulic gradients of the area (Fig. 2).This implies that the hydraulic head fit is independent of the conduit volume during steady-state conditions but depends on the b/m ratio.The influence of the b/m ratio on the hydraulic head fit depends on the hydraulic conductivity of the fissured matrix K m .For low K m values of ca. 1 × 10 −6 m s −1 the hydraulic head fit is completely independent of the conduit geometry and the RMSE is very high (Fig. 4a).For high K m values of ca. 5 × 10 −4 m s −2 (Fig. 4a) the dependence is also of minor importance and the RMSE is relatively stable at ca. 11 m.Due to the high hydraulic conductivity of the fissured matrix the hydraulic gradients do not steepen in the vicinity of the spring even for high b/m ratios.For K m values between the above values the RMSE significantly rises for b/m ratios above 1000 m.For the range of acceptable errors, i.e. lower than 10 m, it is apparent in Fig. 4a that the best-fit K m value is approximately 1 × 10 −5 m s −1 independent of the conduit geometry.However, no distinct best-fit conduit geometry can be derived.There are several parameter combinations providing a good fit for the Gallusquelle spring discharge and the hydraulic head distribution.
The goodness of the fit of the simulation of the tracer breakthrough is mainly determined by the conduit geometry.The influence of the hydraulic conductivity of the fissured matrix K m on flow velocities inside the karst conduits is comparatively low and decreases even further in the vicinity of the springs (Fig. 4b and c) leading to minor influences on tracer travel times.Instead, the quality of the fit mainly depends on the conduit volume and accordingly on the Manning coefficient n (Fig. 5).It is possible to simulate only one of the two tracer experiments with this scenario (Fig. 5).Given the broad range of geometries for which an adequate hydraulic head fit can be achieved (Figs. 2 and 4) it is possible to simulate one of the two tracer peak velocities and the hydraulic head distribution with the same set of parameters.While the simulation of the breakthrough of tracer test 1 requires relatively high n values, of ca.2.5 s m −1/3 , that of tracer test 2 can only be calibrated with lower values of ca.1.7 s m −1/3 (cf.Fig. 5a and b).For every parameter set, where the travel time of the simulated tracer test 2 is not too long, that of tracer test 1 is too short.For the simulation of tracer test 2, the velocities at the beginning of the conduits must be relatively high.To avoid the flow velocities from getting too high in downgradient direction, the conduit size would have to increase drastically due to the constant additional influx of water from the fissured matrix.In the given geometric range, the conduit system has a dominant influence on spring discharge.Physically, this situation corresponds to the conduit-influenced flow conditions (Kovács et al., 2005).Thus, conduit transmissivity is a limiting factor for conduit-matrix exchange and a positive feedback mechanism is triggered, if the conduit size is increased.A higher conduit size leads to higher groundwater influx from the fissured matrix and spring discharge is overestimated.Therefore, parameter analysis shows that scenario 1 is too strongly simplified to correctly reproduce the complex nature of the aquifer.

Scenario 2 -conduit roughness coefficient K c
In scenario 2 the Manning coefficient n was changed from constant to laterally variable.In the literature, n is generally kept constant throughout the conduit network (e.g.Jeannin, 2001;Reimann et al., 2011) for lack of information on conduit geometry.However, it is assumed that the Gallusquelle spring is not fed by a single large pipe.Rather there is some evidence in the spring area that a bundle of several smallinterconnected pipes feed the spring.Since the number of individual conduits per bundle is unknown and the regional modelling approach limits the resolution of local details, the small diameter conduits, which the bundle consists of, cannot be simulated individually.Therefore, each single pipe in the model represents a bundle of conduits in the field.
It can be assumed that the increase in conduit cross section is at least partly provided by additional conduits added to the bundle rather than a single individual widening conduit.Therefore, while the cross section of the simulated conduit, i.e. the total effective cross section of the conduit bundle, in-creases towards the springs, it is not specified how much of this increase is due to the individual conduits widening and how much is due to additional conduits, not distinguishable in the simulation.If the simulated effective cross sectional area increase is mainly due to additional conduits being included in the bundle, the surface / volume ratio increases with the cross section, contrary to what would be observed, if a single conduit in the model would represent a single conduit in the field.The variation in surface area / volume ratio implicitly leads to a larger roughness in the simulation, even further enhanced by exchange processes between the individual conduits.This effect again leads to an increase in the Manning coefficient n in the downgradient direction towards the spring for a simulated single conduit.Since the number and size of the individual conduits is unknown, it is impossible to calculate the change of n directly from the geometry.Thus, a simple scenario was assumed where the roughness coefficient K c , which is the reciprocal of n, was linearly and negatively coupled to the rising conduit radius (Eq.13).
where r c [L] is the conduit radius and r c,max [L] the maximum conduit radius simulated for the respective spring, which COMSOL ® calculates from Eq. ( 1).m h [L −2/3 T −1 ] and b h [L 1/3 T −1 ] are calibration parameters determining the slope and the lowest value of the roughness coefficient respectively.
For every conduit geometry several combinations of m h and b h lead to the same spring discharge.However, hydraulic head fit and tracer velocities are different for each m h -b h combination even if spring discharge is the same.With the new parameters a higher variation of velocity profiles is possible.This allows for the calibration of the tracer velocities of both tracer tests.The dependence of tracer test 2 on m h is much higher than that of tracer test 1 since it is injected further upgradient towards the beginning of the conduit (Fig. 6).Therefore, tracer test 2 is influenced more strongly by the higher velocities far away from the spring introduced by high m h values and always shows a significant positive correlation with m h (Fig. 6).
Since the slope of K c is negative with respect to the conduit length, the variable K c leads to a slowing down of water towards the springs.As discussed in detail by Oehlmann et al. (2013) a rise of transmissivity towards the springs is observed in the Gallusquelle area.Therefore, adequate hydraulic head fits can only be obtained, if the decrease of K c towards the spring is not too large and compensates the effect of the increase in conduit transmissivity due to the increasing conduit radius.This effect reduces the number of possible and plausible parameter combinations.From these considerations a best-fit model can be deduced capable of reproducing all objective functions within the given error ranges (Fig. 7a).According to the model simulations, karst groundwater discharge and flow velocities significantly depend on the total conduit volume as is to be expected.It can be deduced from the parameter analysis that the conduit volume can be estimated at ca. 100 000 m 3 for the different parameters to match equally well (Fig. 7a).

Scenario 3 -extent of conduit network
In scenario 3, a laterally further extended conduit system was employed, assuming the same maximum conduit volume as in scenarios 1 and 2 but with different spatial distribution along the different total conduit lengths.The original conduit length for the Gallusquelle spring in scenarios 1 and 2 is 39 410 m, for scenario 3 it is 63 490 m; therefore, the total length was assumed to be larger by ca.50 % (Fig. 8).The geometry of the original network was mainly constructed along dry valleys where point-to-point connections are observed based on qualitative evaluation from artificial tracer tests.Of the dry valleys without tracer tests, only the larger ones were included, where the assumption of a high karstification is backed up by the occurrence of sinkholes (Mohrlok and Sauter, 1997).Therefore, it represents the minimal extent of the conduit network.For scenario 3 the network was extended along all dry valleys within the catchment, where no tracer tests were conducted.
The results of the parameter variations are comparable to those of scenario 2 (cf.Fig. 7a and b).While the hydraulic head contour lines are smoother than for the original conduit For scenario 4 (c) the root mean square error of the hydraulic heads is given for two different conduit geometries in relation to the hydraulic conductivity of the fissured matrix K m .For the version with laterally variable matrix conductivity the axis shows as an example the hydraulic conductivity of the north-western part.The parameters for the two geometries are given in Table 3.
length the general hydraulic head fit is the same (Fig. 7b).It seems possible to obtain a good fit for all model parameters but the scenario is more difficult to handle numerically.Calculation times are up to 10 times larger compared to the other scenarios and goodness of convergence is generally lower.
Since the calibrated parameters are not significantly different from those deduced in scenario 2 it is concluded that the ambiguity introduced by the uncertainty in total conduit length is small if hydraulic conduit parameters and total conduit volumes are the aim of investigation.

Scenario 4 -matrix hydraulic conductivity K m
In scenario 4, the homogeneously chosen hydraulic conductivity of the fissured matrix K m was changed into a laterally variable conductivity based on different types of lithology and the spatial distribution of the groundwater potential.Sauter (1992) found from field measurements that the area can be divided into three parts with different hydraulic conductivities.Oehlmann et al. (2013) discussed that the major influence is the conduit geometry leading to higher hydraulic transmissivities close to the springs in the east of the area.It is also possible that not only the conduit diameters change towards the spring but the hydraulic conductivity of the fissured matrix as well, since the aquifer cuts through three stratigraphic units (Sect.3).These geologic changes are likely to affect the lateral distribution of hydraulic conductivities (Sauter, 1992).Figure 9 shows the division into three different areas.K m values were varied in the range of the values measured by Sauter (1992).
It was expected that a laterally variable K m value has a major influence on the hydraulic head distribution.All variations of scenario 2 that produce good results for both tracer tests and have a high total conduit volume above 100 000 m 3 yield poor results for hydraulic head errors and spatial dis- tributions of the hydraulic heads (Fig. 7a).For scenario 4, two different conduit configurations (geometries) were chosen that achieve good results with respect to conduit flow velocities.Geometry G1 has a conduit volume of 112 000 m 3 .G2 has a higher b value which leads to the maximum conduit volume of ca. 150 000 m 3 .All parameters for the two simulations are given in Table 4.
It was found that while the maximum root mean square error of the hydraulic head fit is similar for both geometries, the minimum RMSE for the hydraulic head is determined by the conduit system.It is not possible to compensate an unsuitable conduit geometry with suitable K m values (Fig. 7c), which assists in the independent conduit network and fissured matrix calibration.This observation increases the confidence in the representation of the conduits and improves the possibility to deduce the conduit geometry from field measurements.For an adequate conduit geometry, laterally variable matrix conductivities do not yield any improvement.The approach introduces additional parameters and uncertainties because the division of the area into three parts is not necessarily obvious without detailed investigation.From the distribution of the exploration and observation wells (Fig. 1a) it is apparent that especially in the south and west the boundaries are not well defined.

Scenario 5 -conduit intersections
In scenario 5, the effect of the conduit diameter change at intersections was investigated.In the first four scenarios the possible increase in cross sectional area at intersecting conduits was neglected.In nature, however, the influx of water from another conduit is likely to influence conduit evolution and therefore its diameter.In general, higher flow rates lead to increased dissolution rates because dissolution products are quickly removed from the reactive interface.If conditions are turbulent the solution is limited by a diffusion dominated layer that gets thinner with increasing flow velocities (Clemens, 1998).Clemens (1998) simulated karst evolution in simple Y-shaped conduit networks and found higher diameters for the downstream conduit even after short simulation times.Preferential conduit widening at intersections could further be enhanced by the process of mixing corrosion (Dreybrodt, 1981).However, Hückinghaus (1998) found during his karst network evolution simulations that the water from other karst conduits has a very high saturation with respect to Ca 2+ compared to water entering the system through direct recharge.Thus, if direct recharge is present, the mixing with nearly saturated water from an intersecting conduit could hamper the preferential evolution of the conduit downstream slowing down the aforementioned processes.In scenario 5 the influence of an increase in diameter at conduit intersections was investigated.Since the amount of preferential widening at intersections is unknown, the cross sections of two intersecting conduits were added and used as starting cross section for the downstream conduit.The new conduit Table 4. Parameters for the two different conduit configurations compared in scenario 4. b is the minimum conduit radius, m the slope of radius increase towards the springs, b h the highest conduit roughness, m h the slope of roughness decrease away from the spring and V the conduit volume.
where r c2 is the conduit radius downstream of the intersection and r c0 and r c1 the conduit radii of the two respective conduits before their intersection.
Results are very similar to those of scenario 2 (cf.Fig. 7a  and d).Both simulations result in nearly the same set of parameters (Table 1).The estimated conduit volume is even a little smaller for scenario 5 since larger cross sections in the last conduit segment near the spring are reached for a lower total conduit volume.The drastic increase of conduit cross sections at the network intersections leads to higher variability in the cross sections along the conduit segments.The differences between the peak offsets of both tracer tests are higher compared to those of scenario 2. While the peak time of tracer test 2 can be calibrated for large conduit volumes, i.e. conduit volumes above 120 000 m 3 (Fig. 7d), the peak time of tracer test 1 is too late for large conduit volumes.This is due to the fact that the injection point for tracer test 1 is much closer to the spring than that for tracer test 2.In scenario 5 the conduit volume is spatially differently distributed from that of scenario 2 for the identical total conduit volume.The drastic increase in conduit diameters downgradient of conduit intersections leads to rather high conduit diameters in the vicinity of the spring.Therefore, while tracer transport in tracer test 2 occurs in relatively small conduits with high flow velocities and larger conduits with lower velocities, the tracer in tracer test 1 is only transported through the larger conduits whose flow velocities are restricted by the spring discharge.In Fig. 7d the parameter values for the best fit would lie well below the lower boundary of the diagram at negative values below −10 h.However, since the fit for conduit volumes around 100 000 m 3 is similar to that of scenario 2, the two scenarios can in this case not be distinguished based on field observations.

Conclusions of the parameter analysis
Table 3 provides a comparison, i.e. the characteristics for all scenarios.The parameter analysis shows that there is only a limited choice of parameters with which the spring discharges (water balance), the hydraulic head distribution and the tracer velocities can be simulated.Scenario 1 is the only scenario that cannot reproduce the peak travel times observed in both tracer tests simultaneously (Sect.4.1).It underestimates the complexity of the geometry and internal surface characteristics (e.g.roughness) of the conduit system.
Scenario 4 introduces two additional model parameters.The best fit for this scenario is, however, still achieved with all three K m values being equal, which basically results in the parameter set of scenario 2. This implies that the major influence leading to the differences in hydraulic gradients observed throughout the area is the conduit system and not the variability of the fissured matrix hydraulic conductivity.It was also shown that for the Madison aquifer (USA), by Saller et al. (2013), a better representation of the hydraulic head distribution can be achieved by including a discrete conduit system even for reduced variability in the hydraulic conductivity of the fissured matrix.Their conclusion complies very well with the findings for scenario 4.
Scenario 3 simulates the presence of a couple of additional smaller dendritic branches.The deduced parameter values and the fit of the objective functions are similar to those of scenarios 2 and 5.Because of long calculation times without additional advantage for the presented study, scenario 3 is not considered for further analysis.
Scenarios 2 and 5 are both judged as suitable.Their parameters and the quality of the fit are similar.Therefore, it is not possible to decide which one is the better representation of reality.Regarding the different processes interacting during karst evolution (Sect.4.5) it is most likely that the actual geometry ranges somewhat in between these two scenarios.Table 1 summarizes all parameters of both simulations and Fig. 10 shows the simulated tracer breakthrough curves and spring discharges.

Plausibility of the best-fit simulations
The main objective of the model simulation is not only to reproduce the target values but also to provide insight into dominating flow and transport processes, sensitive parameters and to check the plausibility of the model set-up.Possible ambiguities in parameterizations can also be checked, i.e. different combinations of parameters producing identical model output.
For these aims model parameters and aquifer properties simulated with scenarios 2 and 5 are compared to those observed in the field.As seen in Table 1 most of the calibrated parameters range well within values provided in the literature.The calibrated Manning coefficients are relatively high compared to other karst systems.Jeannin (2001) lists effective conductivities for several different karst networks that translate into n values of between 0.03 and 1.07 s m −1/3 , showing that the natural range of n values easily extends across 2 orders of magnitude and the minimum n values of the simulation lie within the natural range.The maximum n values are significantly higher than those given by Jeannin (2001).This is not surprising since the calibrated n value reflects the total roughness of the conduit bundles and therefore includes geometric conduit properties in addition to the wall roughness that it was originally defined for.This effect is specific for the Gallusquelle area but it might be important to consider for other moderately karstified areas as well where identification of conduit geometries is especially difficult.
The total conduit volume of the Gallusquelle spring derived from scenarios 2 and 5 is only 50 % of that estimated with traditional methods (Geyer et al., 2008).Since the conduit transmissivity increases towards the spring water enters the conduits preferably in the vicinity of the spring in the Gallusquelle area.Therefore, the matrix contribution is high.In addition, the travel time at peak concentration of tracer test 2, which was used for the volume estimation by Geyer et al. (2008), is longer than 3 days, during which time matrixconduit water exchange can readily take place.Based on the results of a tracer test conducted in a distance of 3 km to the Gallusquelle spring Birk et al. (2005) estimated the error incurred by deducing the conduit volume without taking conduit-matrix exchange fluxes into account with a very simple numerical model.The authors found a difference in conduit volumes of approximately 50 %.This fits well with the results of the present simulation.Birk et al. (2005) also the simulated equivalent conduit cross sectional area between their tracer injection point and the spring to be 13.9 m 2 .For scenario 2 the simulated average cross sectional area is 11.9 m 2 and for scenario 5 13.4 m 2 , which compares very well with the results of Birk et al. (2005).
It was not possible to match the shape of both breakthrough curves with the same dispersivity.The apparent dispersion in the tracer test 2 breakthrough is much higher compared to that of tracer test 1, while the breakthrough of tracer test 1 shows a more expressed tailing (Fig. 10a and b).This corresponds to the effect observed by Hauns et al. (2001).The authors found scaling effects in karst conduits: the larger the distance between input and observation point, the more mixing occurred.The tailing is generally induced by matrix diffusion or discrete geometric changes such as pools, where the tracer can be held back and released more slowly.Theoretically, every water drop employs medium and slow flow paths if the distance is large enough, leading to a more or less symmetrical, but broader, distribution and therefore a higher apparent dispersion (Hauns et al., 2001).To quantify this effect, exact knowledge of the geometric conduit shape such as the positions and shapes of pools would be necessary.Furthermore, an additional unknown possibly influencing the observed retardation and dispersion effects is the input mechanism.The simulation assumes that all introduced tracers immediately and completely enter the conduit system, which neglects effects of the unsaturated zone on tracer breakthrough curves.In addition, the shape of the breakthrough curve of tracer test 2 is difficult to deduce since the 6 h sampling interval can be considered as rather low leading to a breakthrough peak which is described by only seven measurement points.Therefore, the apparent dispersivity was calibrated for both breakthrough curves separately.Calibrated dispersivity ranges well within those quoted in lit-erature (Table 1).The mass recovery during the simulation was determined to range between 98.4 and 99.9 % in all simulations.The slight mass difference results from a combination of diffusion of the tracer into the fissured matrix and numerical inaccuracies.
The spring discharge of the minor springs in the area (Sect.3) was slightly underestimated in most cases (Fig. 10c).For most springs the models of scenarios 2 and 5 provide similar results.The underestimation of discharge is in the order of < 0.05 m 3 s −1 and is not expected to significantly influence the general flow conditions.It probably results from the unknown conduit geometry in the catchments of the different minor springs.The only case in which the two scenarios give significantly different results is the spring discharge of the spring group consisting of the Ahlenberg and Büttnauquellen springs (Fig. 10c).Scenario 2 overestimates and scenario 5 underestimates the discharge.This is due to the fact that the longest conduit of the Ahlenberg and Büttnauquellen springs is longer than the longest one of the Gallusquelle spring but the conduit network has less intersections (Fig. 1).Therefore, the conduit volume of the Ahlenberg and Büttnauquellen springs is 134 568 m 3 in scenario 2 and only 75 085 m 3 in scenario 5 leading to the different discharge values.It is reasonable to assume that a better fit for the spring group can be achieved, if more variations of conduit intersections are tested.An adequate fit for the Fehla-Ursprung spring of 0.1 m 3 s −1 was achieved for both scenarios with a fault aperture of 0.005 m.

Uncertainties and limitations
The most important uncertainties regarding the reliability of the simulation include the assumptions that were made prior to modelling.First, flow dynamics were neglected.This approach was chosen because tracer tests are supposed to be conducted during quasi-steady-state flow conditions.However, this is only the ideal case.During both tracer tests spring discharge declined slightly.The influence of transient flow on transport velocities inside the conduits was estimated by a very simple transient flow simulation for the best-fit models in which recharge and storage coefficients were calibrated to reproduce the observed decline in spring discharges.The transient flow only slightly affected peak velocities but lead to a larger spreading of the breakthrough curves and therefore lower calibrated dispersion coefficients.This effect occurred because the decline in flow velocities is not completely uniform inside the conduits and depending on where the tracer is at which time it experiences different flow velocities in the different parts of the conduits, which leads to a broader distribution at the spring.The same breakthrough curves can be simulated under steady-state flow conditions with slightly higher dispersivity coefficients.So, the calibrated dispersivities do not only represent geometrical heterogeneities but also temporal effects as is the case for all standard evaluations of dispersion from tracer breakthrough curves.The influence of rapid recharge is not to considered in the simulation of baseflow conditions.However, there might be an influence on flow velocities during the actual recharge events, i.e. if rapid recharge is intensive and strong enough to lead to a reversal of the flow gradients between conduit and fissured matrix.Therefore, an alternative simulation was performed for tracer test 2, which was conducted during high flow conditions (Table 2) after a recharge event.The maximum percentage of direct recharge of 10 % estimated by Sauter (1992) and Geyer et al. (2008) was used for this simulation.Neither for scenario 2 nor for scenario 5 a gradient reversal between conduit and matrix occurred and the influence on flow velocities was negligible (Fig. 11).
Furthermore, flow in all karst conduits was simulated for turbulent conditions.Turbulent conditions can be generally assumed in karst conduits (Reimann et al., 2011) and also apply to all calibrated model conduit cross sections.Since the conduit cross section presents the total cross section of the conduit bundle, the cross sections of the individual tubes are uncertain, though.The high n values suggest that the surface / volume ratio is relatively high, which implies that the individual conduit cross sections are rather small.Therefore, laminar flow in some conduits is likely.While laminar flow conditions in the conduits influence hydraulic gradients considerably, this fact is believed not to influence the overall results and conclusions of this study, i.e. the relative significance of the parameters deduced from parameter analysis and the deduced conduit volume, especially since flow is simulated for steady-state conditions.
For all distributed numerical karst simulations, uncertainties regarding the exact positions and interconnectivities of the conduit branches still remain.Due to the extensive investigations already performed in previous work (Sect.3) these uncertainties are reduced in the Gallusquelle area and the above scenarios include the most probable ones.However, the flexibility of the modelling approach allows for the integration of any future information that might enhance the numerical model further.

Calibration strategy
For a successful calibration of a distributed groundwater flow and transport model for a karst area on catchment scale certain constraints have to be set a priori.The geometry of the model area, i.e. locations/types of boundary conditions and aquifer base, fixed during calibration, has to be known with sufficient certainty.Furthermore, the objective functions for calibration have to be defined, i.e. the hydraulic response of the system and transport velocities.In a karst groundwater model, these consist of measurable variables, i.e. spring discharges, hydraulic heads in the fissured matrix and two tracer breakthrough curves.The hydraulic head measurements should be distributed across the entire catchment and preferably close to the conduit system, should geometric conduit parameters be calibrated for as well.It is expected that the influence of the conduits on the hydraulic head decreases and the influence of matrix hydraulic conductivities increases with distance to the conduit system.In the design of the tracer experiment, the following criteria should be observed: for a representative calibration, the dye should be injected at as large a distance to each other as possible with one of them including the length of the whole conduit system.Each tracer test gives integrated information about its complete flow path.If the injection points lie close together, no information about the development of conduit geometries from water divide to spring can be obtained.Further, the dye should be injected as directly as possible into the conduit system, e.g.via a flushed sinkhole, to obtain information on the conduit flow regime and to minimize matrix interference.To ease interpretation a constant spring discharge during the tests is desirable.
In this study, the flow field was simulated not only for the catchment area of the Gallusquelle spring, but also for a larger area including the catchment areas of several smaller springs (Fig. 1).This is in general not essential for deducing conduit volumes and setting up a flow and transport model.Simulating several catchments, however, helps to increase the reliability of the simulation.The positions of water divides are majorly determined by the hydraulic conductivity of the fissured matrix K m , so that the simulated catchment areas of the different springs can be used to estimate how realistic the simulated flow field is and decrease the range of likely K m values.In this study, high K m values above ca.3 × 10 −5 m s −1 made the simulation of the spring discharge of the Fehla-Ursprung spring (Fig. 1) impossible because the water divide in the west could not be simulated and most of the water in the area discharged to the east towards the river Lauchert resulting in a very narrow and long catchment area for the Gallusquelle spring.
There are eight parameters available for model calibration in this study.Two of these parameters define the conduit geometry: b is the lowest conduit radius and m the slope with which the conduit radius increases.One parameter, d f , defines the aperture of the fault zone.The hydraulic conductivity of the fissured matrix is represented by the parameter K m and the roughness of the conduit system by two parameters: b h represents the highest roughness and m h the slope of roughness decrease in upgradient direction from the spring.The last two parameters ε 1 and ε 2 are the respective conduit dispersivities obtained from the two artificial tracer experiments (Table 1).
For efficiency reasons it is important to know which of these parameters can be calibrated independently.The apparent transport dispersivities ε 1 and ε 2 are pure transport parameters, which influence only the shape of the breakthrough curves and not the flow field.The hydraulic model parameters influence the shape of the tracer breakthrough curves as well.Therefore, dispersivities ε 1 and ε 2 should be calibrated separately after calibrating the hydraulic model parameters.
Only for hydraulically dominant fault zones knowledge of the fault zone aperture d f is required.For the model area this parameter was required for one fault zone lying in the west of the area feeding the Fehla-Ursprung spring (Fig. 1).Since the Fehla-Ursprung spring has its own catchment area the fault zone has only minor influence on the flow regime in the Gallusquelle catchment.Its hydraulic parameters were calibrated at the beginning of the simulation procedure to reproduce the catchment and the discharge of the Fehla-Ursprung spring adequately and kept constant throughout all the simulations.In the final calibrated models it was rechecked, but the calibrated value was still acceptable.
The hydraulic conductivity of the fissured matrix K m can be calibrated independently in principle as well.The influence on spring discharge is relatively small.The best-fit K m value depends on the conduit parameters, i.e. geometry and roughness, since the hydraulic conductivities of the conduit system and of the fissured matrix define the total transmissivity of the catchment area together.Nonetheless, the best-fit value lies in the same range for different conduit geometries (Figs.4a and 7c).The greater the difference between the simulated conduit geometries, the more likely is a slight shift of the best-fit K m value.Therefore, it is advisable to calibrate it anew for significant model changes, e.g.different scenarios, but to keep it constant during the rest of the calibrations.For the best-fit configuration, potentially used as a prognostic tool, the K m value needs to be checked and adapted if necessary.This observation is, however, only valid for steady-state flow conditions.The dynamics of the hydraulic head and spring discharge might be highly sensitive to the matrix hydraulic conductivity, the conduit-matrix exchange coefficient and the lateral conduit extent.This work focuses on the conduits as highly conductive pathways for e.g.contaminant transport, but the calibration of matrix velocities, e.g. by use of environmental tracers, would likely be sensitive to the K m values as well.Therefore, the choice of the flow regime and the objective functions determines the strength of the interdependencies between fissured matrix and conduit system parameters and therefore whether K m can be calibrated independently.
The conduit parameters for geometry and roughness, here four parameters (lowest conduit radius b, slope of radius increase m, highest roughness b h and slope of roughness decrease m h ), have to be varied simultaneously.All of them have a major influence on spring discharge and cannot be varied separately without introducing discharge errors.For each conduit geometry, there are a number of possible b hm h combinations that result in the observed spring discharge.In general, the slowest transport velocities are achieved with a m h value of zero.So, to deduce the range of geometric parameters that reproduce the objective functions, it is advisable to check the minimum conduit volume for which the tracer tests are not too fast for a value of m h equal to zero.For the Gallusquelle area, transmissivities significantly increase towards the springs, which is characteristic for most karst catchments.Therefore, low b h values oppose the general hydraulic head trend: they increase the conduit roughness at the spring leading to slower flow and higher gradients.The higher the conduit volume, the higher b h is required to reproduce the observed transport velocities.Therefore, the best-fit model likely has the smallest conduit volume for which both tracer tests can be reproduced.In Fig. 7 this condition can be seen to clearly range in the order of 100 000 m 3 for the Gallusquelle area.While the four conduit parameters allow for a good model fit, they are pure calibration parameters.They show that the karst conduit system has a high complexity, which cannot be neglected for distributed velocity and hy-S.Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching draulic head representation.A systematic simulation of the heterogeneities, e.g. with a karst genesis approach, would be a process-based improvement to the current method and give more physical meaning to the parameters.

Conclusions
The study presents a large-scale catchment-based distributed hybrid karst groundwater flow model capable of simulating groundwater flow and solute transport.For flow recession conditions this model can be used as a predictive tool for the Gallusquelle area with relative confidence.The approach of simultaneous pattern matching of flow and transport parameters provides new insight into the hydraulics of the Gallusquelle conduit system.The model ambiguity was significantly reduced to the point where an estimation of the actual karst conduit volume for the Gallusquelle spring could be made.This would not have been possible simulating only one or two of the three objective functions, i.e. the spring discharge, the hydraulic head distribution and two tracer tests.
The model allows for the identification of the relevant parameters affecting karst groundwater discharge and transport in karst conduits and the examination of the respective overall importance in a well-investigated karst groundwater basin for steady-state flow conditions.While a differentiated representation of the roughness values in the karst conduits is substantial for buffering the lack of knowledge of the exact conduit geometry, e.g.local variations in cross section and the number of interacting conduits, variable matrix hydraulic conductivities cannot improve the simulation.It was shown that the effect of the unknown exact lateral extent of the conduit system and the change in conduit cross section at conduit intersections is of minor importance for the overall karst groundwater discharge.This is important since these parameters are usually unknown and difficult to measure in the field.
For calibration purposes, this study demonstrates that for a steady-state flow field and the observed objective functions the hydraulic conductivities of the fissured matrix can practically be calibrated independently of the conduit parameters.Furthermore, a strategy for the simultaneous calibration of conduit volumes and conduit roughness in a complex karst catchment was developed.
As discussed in Sect. 5 the major limitation of the simulation is the neglect of flow dynamics, which limits the applicability to certain flow conditions.Therefore, transient flow simulation is the focus of on-going work.This will enhance the applicability of the model as a prognostic tool to all essential field conditions and lead to further conclusions regarding the important karst system parameters, their influences on karst hydraulics and their interdependencies.It can be expected that some parameters, which are of minor importance in a steady-state flow field, e.g. the lateral conduit extent and the percentage of recharge entering the conduits directly, will exhibit significant influence for transient flow conditions.
where d f is the fault aperture [L], ρ the density of water [M L −3 ], g the gravity acceleration [L T −2 ] and µ the dynamic viscosity of water [M T −1 L −1 ].

Figure 1 .
Figure 1.(a) Plan view of the model area.Settlements, fault zones and rivers in the area are plotted, as well as the 20 observation wells used for hydraulic head calibration, the six springs used for spring discharge calibration and the two tracer tests employed for flow velocity calibration.Catchment areas for the Gallusquelle spring and the Ahlenberg and Büttnauquellen springs were simulated according to Oehlmann et al. (2013).(b) Three-dimensional view of the model.The upper boundary is hidden to allow a view of the karst conduit system and the aquifer base.The abbreviation BC stands for boundary condition.At the hidden upper boundary, a constant recharge Neumann BC is applied.

Figure 2 .
Figure 2. Hydraulic head distributions for different combinations of geometric conduit parameters for scenario 1. b is the lowest conduit radius and m the radius increase along the conduit.For comparison, a trend line is fitted to the measured hydraulic head values showing the distribution of hydraulic gradients from the Gallusquelle spring to the western border of its catchment area.

Figure 3 .
Figure 3. Conceptual overview of the simulated scenarios.The conduit geometry and the varying parameters are shown.

Figure 4 .
Figure 4. Influence of the hydraulic conductivity of the fissured matrix on the objective functions.(a) Influence on the root mean square error of the hydraulic head distribution in relation to the conduit geometry.The conduit geometry is represented by the parameter b/m (Eq.1), which is the ratio of the smallest radius to the slope of radius increase along the conduit length.(b) Influence on the conduit flow velocity for tracer test 1.(c) Influence on the conduit flow velocity for tracer test 2.

Figure 5 .
Figure 5. Difference between peak concentration times vs. the Manning n value for scenario 1. High n values correspond to high conduit volumes and high cross sectional areas at the spring (a) for tracer test 1 and (b) for tracer test 2.

Figure 6 .
Figure6.Hydraulic head errors and differences between peak concentration times for both tracer tests for scenario 1.The example is shown for a conduit geometry with a starting value b = 0.01 m and a radius increase of m = 2 × 10 −4 .Each m h (m −2/3 s −1 ) value corresponds to a respective value of the highest conduit roughness b h (m 1/3 s −1 ) and each combination results in the same spring discharge.

Figure 7 .
Figure 7. Calibrated values for the simulated scenarios.For scenarios 2, 3 and 5 (a, b, d) hydraulic head fit and the peak-offset times of both tracer tests (referred to as TT 1 and TT 2) are shown in relation to conduit volume.The thick grey bar marks the target value of zero.For scenario 4 (c) the root mean square error of the hydraulic heads is given for two different conduit geometries in relation to the hydraulic conductivity of the fissured matrix K m .For the version with laterally variable matrix conductivity the axis shows as an example the hydraulic conductivity of the north-western part.The parameters for the two geometries are given in Table3.

Figure 8 .Figure 9 .
Figure 8. Extended conduit system for scenario 3. The conduit configuration (extent) that is used for the other scenarios is marked in red.

Figure 10 .
Figure 10.Comparison of the best-fit simulations with field data for scenarios 2 and 5. (a) Breakthrough curve of tracer test 1, (b) breakthrough curve of tracer test 2 and (c) spring discharge.

Figure 11 .
Figure11.Flow velocities inside the main conduit branch of the Gallusquelle spring during the simulation of tracer test 2. The best-fit simulations for scenarios 2 and 5 are compared to simulations where a direct recharge of 10 % is introduced.

Table 1 .
Calibrated and simulated parameters for the best-fit simulations.Literature values are given if available.TT 1 and TT 2 refer to the two tracer tests.

Table 2 .
Field data of the simulated tracer tests.

Table 3 .
Specifics of the different scenarios.The bold writing indicates the parameter that is analysed in the respective scenario.The results are indicated by comparative markers."+" means good, "o" means average and "-" means bad compared to the other scenarios.Details to the scenarios and results evaluation can be found in Sect. 4.