Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ doi:10.5194/hess-19-893-2015 © Author(s) 2015. CC Attribution 3.0 License. Reducing the ambiguity of karst aquifer models by pattern matching of flow and transport on catchment scale S. Oehlmann1, T. Geyer1,2, T. Licha1, and M. Sauter1 1Geoscience Center, University of Göttingen, Göttingen, Germany 2Landesamt für Geologie, Rohstoffe und Bergbau, Regierungspräsidium Freiburg, Freiburg, Germany Correspondence to: S. Oehlmann (sandra.oehlmann@geo.uni-goettingen.de) Received: 7 July 2014 – Published in Hydrol. Earth Syst. Sci. Discuss.: 4 August 2014 Revised: 17 December 2014 – Accepted: 17 January 2015 – Published: 12 February 2015 Abstract. Assessing the hydraulic parameters of karst aquifers is a challenge due to their high degree of heterogene- ity. The unknown parameter field generally leads to a high ambiguity for flow and transport calibration in numerical models of karst aquifers. In this study, a distributed numerical model was built for the simulation of groundwater flow and solute transport in a highly heterogeneous karst aquifer in south-western Germany. Therefore, an interface for the sim- ulation of solute transport in one-dimensional pipes was im- plemented 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 hy- draulic head distribution in the model area, the spring dis- charge of several springs and the transport velocities of two tracer tests. Furthermore, other measured parameters such as the hydraulic conductivity of the fissured matrix and the maximal karst conduit volume were available for model cal- ibration. Parameter studies were performed for several karst conduit geometries to analyse the influence of the respective geometric and hydraulic parameters and develop a calibra- tion approach in a large-scale heterogeneous karst system. Results show that it is possible not only to derive a con- sistent flow and transport model for a 150 km2 karst area but also to combine the use of groundwater flow and transport parameters thereby greatly reducing model ambiguity. The approach provides basic information about the conduit net- work not accessible for direct geometric measurements. The conduit network volume for the main karst spring in the study area could be narrowed down to approximately 100 000 m3. 1 Introduction Karst systems play an important role in water supply world- wide (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 simulat- ing karst aquifer behaviour. Overviews on the various types of distributed process and lumped-parameter models are pro- vided 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 be- haviour, 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 nu- merical 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 het- erogeneous and anisotropic aquifers with limited data avail- ability is the hybrid modelling approach. The approach sim- ulates 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 repre- senting 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. 894 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching et al., 2011). They are, however, regularly applied in long- term 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 pre- dictive purposes, such models can be used for deriving in- formation about the groundwater catchment itself (Rehrl and Birk, 2010). A major problem for characterizing the groundwater sys- tem 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 vari- ables may give a similar fit to the observed data but some- times very different results for prognostic simulations (Li et al., 2009). Especially the geometric and hydraulic prop- erties 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 observa- tions, can significantly reduce the number of plausible pa- rameter combinations (Ophori, 1999). Especially in hydrol- ogy (e.g. Khu et al., 2008; Hunter et al., 2005) and also for groundwater systems (e.g. Ophori, 1999; Hu, 2011; Hart- mann et al., 2013), this approach has been successfully ap- plied with a wide range of observation types, e.g. ground- water 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 multi- objective 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 num- ber of model runs needed (Khu et al., 2008) and a precise conceptual model is essential as basis for the automatic cal- ibration (Madsen, 2003). In general, numerical models of karst aquifers are difficult to build because of their highly developed heterogeneity (Rehrl and Birk, 2010). Thus, auto- matic calibration procedures are better suited for conceptual and lumped-parameter models, where calibration parameters include effective geometric properties and no spatial repre- sentation of the hydraulic parameter field and conduit geom- etry is necessary. Complex distributed numerical approaches generally require longer simulation times due to the neces- sary spatial resolution. Long simulation times limit the num- ber of model runs that can reasonably be performed and man- ual calibration based on hydrogeological knowledge is nec- essary (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 simu- late the hydraulic head distribution in the area, spring dis- charge and tracer breakthrough curves simultaneously on catchment scale. Some studies combine groundwater flow with particle tracking for tracer directions (e.g. Worthing- ton, 2009; Saller et al., 2013) without simulating tracer trans- port. On the other hand there are studies simulating break- through 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 hydrogeolog- ical 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 de- velop a basis for further prognostic simulations in a het- erogeneous karst aquifer with a distributed modelling ap- proach on catchment scale, but also to reduce model am- biguity and draw conclusions on the spatially distributed karst network geometries and the actual karst conduit vol- ume. The approach shows the kind and minimum number of field observations needed for this aim. Furthermore, a sys- tematic calibration strategy is presented to reduce the num- ber of necessary model runs and the simulation time com- pared to standard multi-objective calibrations. For this pur- pose a hybrid model was built and a pattern matching proce- dure 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 dis- charge of six springs and the tracer breakthrough curves of two tracer tests. 2 Modelling approach The simulation is based on the mathematical flow model dis- cussed 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 simul- taneously in the three-dimensional fissured matrix, in an in- dividual 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. Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 895 In agreement with these results and karst genesis simulations by Liedl et al. (2003), the conduit radius is calculated as rc =ms+ b, (1) where rc [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 parame- ters. Parameters without a subscript are the same for all karst features in the model. 2.1 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 dur- ing tracer tests, which are ideally conducted under quasi- steady-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). Qm =∇ (ρum) , (2a) um =−Km∇Hm, (2b) where Qm is the mass source term [M L−3 T−1], ρ the den- sity of water [M L−3] and um the Darcy velocity [L T−1]. In Eq. (2b) Km is the hydraulic conductivity of the fissured matrix [L T−1] and Hm the hydraulic head [L]. Two-dimensional fracture flow in the fault zone was sim- ulated with the COMSOL® fracture flow interface. The in- terface only allows for the application of the Darcy equation inside of fractures, so laminar flow in the fault zone was as- sumed. In order to obtain a process-based conceptualization of flow, the hydraulic fault conductivityKf was calculated by the cubic law (Eq. 3): Kf = d 2 f ρg 12µ , (3) where df is the fault aperture [L], ρ the density of water [M L−3], g the gravity acceleration [L T−2] and µ the dy- namic viscosity of water [M T−1 L−1]. For groundwater flow in the karst conduits, the Manning equation was used (Eq. 4). uc = 1 n ( rc 2 ) 2 3 √ dHc dx , (4) where uc is the specific discharge in this case equalling the conduit flow velocity [L T−1], n the Manning coefficient [T L−1/3], rc/2 the hydraulic radius [L] and dHc/dx the hy- draulic gradient [–]. The Manning coefficient is an empirical value for the roughness of a pipe with no physical nor mea- surable meaning. The hydraulic radius is calculated by divid- ing 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 approxima- tion. 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 ac- curately. 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 qex = α L (Hc−Hm) , (5) where qex is the water exchange between conduit and fissured matrix [L2 T−1] per unit conduit length L [L], Hm the hy- draulic head in the fissured matrix [L],Hc the hydraulic head in the conduit [L] and α the leakage coefficient [L2 T−1]. The leakage coefficient was defined as α = 2pircKm, (6) where 2pi rc is the conduit perimeter [L]. Other possible in- fluences, e.g. the lower hydraulic conductivity at the solid– liquid interface of the pipe and the fact that water is not ex- changed along the whole perimeter but only through the fis- sures 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 hy- draulic impulses. Since the flow simulation is performed for steady-state conditions this simplification is not expected to exhibit significant influence on the flow field. 2.2 Solute transport Transient solute transport was simulated based on the steady- state groundwater flow field. COMSOL Multiphysics® offers a general transport equation with its solute transport inter- face. This interface was applied for the three-dimensional www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 896 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching fissured matrix. In this work saturated, conservative trans- port was simulated, with an advection–dispersion equation (Eq. 7) ∂ ∂t (θmcm)+∇ (umcm)=∇ [(DDm+De)∇cm]+ Sm, (7) where θm is the matrix porosity [–], cm the solute concentra- tion [M L−3], DDm the mechanical dispersion [L2 T−1] and De the molecular diffusion [L2 T−1]. Sm is the source term [L3 T−1]. The solute transport interface cannot be applied to one- dimensional elements within a three-dimensional model. COMSOL® offers a so-called coefficient form edge PDE in- terface 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 con- duits: θc ∂cc ∂t +∇ (−Dc∇cc+uccc)= f, (8) where θc is the conduit porosity which is set equal to 1, Dc [L2 T−1] the diffusive/dispersive term Dc= (DDc+De), f the source term and uc [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 in- and outflow along the conduits lead to significant velocity di- vergence in the conduit system. This needs to be considered for mass conservation during the simulation. The mechanical conduit dispersion DDc was calculated with Eq. (9) (Hauns et al., 2001). DDc = εuc, (9) 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 cex: f = cex =−De 2pirc L (cm− cc)− qexci . (10) The first term of the right-hand side of Eq. (10) defines the diffusive exchange due to the concentration difference be- tween conduit and fissured matrix. The second term is a con- ditional term adding the advective exchange of solute due to water exchange. The concentration of the advective exchange ci is defined as ci = { cc if qex > 0 cm if qex ≤ 0 . (11) When qex is negative, the hydraulic head in the fissured ma- trix is higher than in the conduit (Eq. 5) and water with the solute concentration of the fissured matrix cm enters the con- duit. When it is positive, water with the solute concentration cc 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 pi r2c [L2]. These considerations lead to the following transport equation for the karst conduits: pir2c ∂cc ∂t +pir2c∇ (−Dc∇cc+uccc) =−De 2pirc L (cm− cc)− qexci . (12) 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 km2, including the catchment area of the Gallusquelle spring and surrounding smaller spring catch- ments (Oehlmann et al., 2013). The Gallusquelle spring is the main point outlet with a long-term average annual dis- charge of 0.5 m3 s−1. The model area is constrained by three rivers and no-flow boundaries derived from tracer test in- formation 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) (Gol- wer, 1978; Gwinner, 1993). The marly limestones of the un- derlying Kimmeridgian 1 (ki 1) mainly act as an aquitard. In the west of the area where they get close to the sur- face, 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 iden- tified 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 perpendic- ular 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 Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 897 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 cali- bration. 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. al. (2013) it is more likely that the water crosses the graben structure by a transversal cross-fault (Strayle, 1970). There- fore, 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 re- sults 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 m3 s−1, the annual average ranging from 0.068 to 0.135 m3 s−1. The fault zone aperture was calibrated accordingly (Sect. 5). Due to a large number of studies conducted in the area dur- ing 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 con- duit 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 www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 898 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching Table 1. Calibrated and simulated parameters for the best-fit simu- lations. Literature values are given if available. TT 1 and TT 2 refer to the two tracer tests. Parameter Simulated Simulated Literature values values values scenario 2 scenario 5 Km (m s−1) 8× 10−6 1.5× 10−5 1× 10−6–2× 10−5 (local scale)e 2× 10−5–1× 10−4 (regional scale)e mh (m−2/3 s−1) 0.3 0.3 – bh (m1/3 s−1) 0.22 0.18 – n (s m−1/3) 1.04–4.55 1.05–5.56 0.03–1.07a b (m) 0.01 0.01 – m (–) 2.04× 10−4 1.42× 10−4 – ε1 (m) for TT 1 7.15 7.5 4.4–6.9f, 10e ε2 (m) for TT 2 30 23 20g Ah (m2) 11.9 13.4 13.9f V (m3) 109 351 89 286 ≤ 200 000b RMSE H (m) 5.61 5.91 – Peak offset TT 1 (h) −0.28c −0.28c – Peak offset TT 2 (h) 2.5d −1.39d – a Jeannin (2001); b Geyer et al. (2008); c measurement interval 1 min, simulation interval 2.7 h; d measurement interval 6 h, simulation interval 2.7 h; e Sauter (1992); f Birk et al. (2005); g Merkel (1991); h average for the interval between tracer test 1 and the spring. 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 S- shaped distribution with distance to the Gallusquelle spring. The gradient at each point of the area depends on the com- bination of the respective transmissivity and total flow. The amount of water flowing through a cross sectional area in- creases 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 trans- missivities 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 di- rection of the Gallusquelle spring leading to a smoothing of hydraulic gradients. Geyer et al. (2008) calculated the maximum conduit vol- ume for the Gallusquelle spring Vc [L3] 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 calcu- lated the maximum volume at 218 000 m3. Their approach assumes the volume of the conduit corresponds to the total volume of water discharged during the time between tracer 2000 4000 6000 8000 10000 Distance to the Gallusquelle spring (m) 640 660 680 700 720 740 H yd ra ul ic h ea d (m ) Legend measured values measured trend m = 1x10–5 m = 1x10–4 Hydraulic head values for the simulation with b = 0.01 2000 4000 6000 8000 10000 Distance to the Gallusquelle spring (m) 640 660 680 700 720 740 H yd ra ul ic h ea d (m ) Hydraulic head values for the simulation with b = 0.05 2000 4000 6000 8000 10000 Distance to the Gallusquelle spring (m) 640 660 680 700 720 740 H yd ra ul ic h ea d (m ) Hydraulic head values for the simulation with b = 0.3 2000 4000 6000 8000 10000 Distance to the Gallusquelle spring (m) 640 660 680 700 720 740 H yd ra ul ic h ea d (m ) Hydraulic head values for the simulation with b = 0.5 a) b) c) d) Figure 2. Hydraulic head distributions for different combinations of geometric conduit parameters for scenario 1. b is the lowest conduit radius andm 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. input and tracer arrival neglecting the contribution of the fis- sured matrix. The six springs that were monitored and therefore simu- lated 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 ba- sis for the conduit structure leading to the spring. Tracer directions were available for 32 tracer tests con- ducted 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) was conducted at 10 km distance to the Gallusquelle Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 899 Figure 3. Conceptual overview of the simulated scenarios. The conduit geometry and the varying parameters are shown. 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 recov- ered tracer mass was chosen as input for the tracer test sim- ulation. 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 long- term average annual discharge of 0.5 m3 s−1. For the trans- port simulations, the recharge was then adapted to produce the respective discharge observed during the tracer experi- ment (Table 2). 4 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 contribu- tions to the discharge and conduit flow velocities. The fit- ting parameters include the parameters controlling the re- spective transmissivities of the fissured matrix and the karst conduit system, i.e. the geometry and roughness of the con- duit system, the hydraulic conductivity of the fissured ma- trix 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 con- trolling the temporal distribution of recharge were not con- Table 2. Field data of the simulated tracer tests. Tracer Tracer test 1 test 2 Input mass (kg) 0.75 10 Recovery (%) 72 50 Distance to spring (km) 3 10 Spring discharge (m3 s−1) 0.375 0.76 Sampling interval 1 min 6 h Peak time (h) 47 79.5 sidered. The parameter analysis was performed with COM- SOL Multiphysics® parametric sweep tool, which sweeps over a given parameter range. Parameter ranges were cho- sen 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. There- fore, 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 ba- sic 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- www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 900 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 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. Parameter Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5 Kc constant linear increase linear increase linear increase linear increase Lateral network minimal minimal extended minimal minimal Km constant constant constant variable constant Intersection radius rc2 rc0 rc0 rc0 rc0 √ r2c0 + r2c1 Main results Hydraulic head fit + + + + + Fit of breakthrough – + + + + Model applicability + o – – o eter sets, which could not fulfil this criterion, were not con- sidered 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 simula- tion with the measured value (peak offset). Since tracer ex- periments 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 lo- cations 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 per- formed (e.g. Fig. 2) to ensure that the general characteristics of the area were represented instead of only the statistical value. 4.1 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 lengthm (Eq. 1), only one value of the Manning coefficient n al- lows a simulated discharge for the Gallusquelle spring of 0.5 m3 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 m3 s−1 a higher conduit volume requires a higher cali- brated 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 ade- quate 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 hy- draulic gradients of the area (Fig. 2). This implies that the hydraulic head fit is independent of the conduit volume dur- ing steady-state conditions but depends on the b/m ratio. The influence of the b/m ratio on the hydraulic head fit de- pends on the hydraulic conductivity of the fissured matrix Km. For low Km 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 Km 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 ma- trix the hydraulic gradients do not steepen in the vicinity of the spring even for high b/m ratios. For Km 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 Km 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 provid- ing 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 geome- try. The influence of the hydraulic conductivity of the fis- sured matrix Km on flow velocities inside the karst conduits is comparatively low and decreases even further in the vicin- ity of the springs (Fig. 4b and c) leading to minor influ- ences on tracer travel times. Instead, the quality of the fit Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 901 100 1000 10000 100000 Conduit geometry parameter b/m (m) 10 100 R oo t m ea n sq ua re e rro r o f t he h yd ra ul ic h ea d di st rib ut io n (m ) Legend Km = 1x10-6 Km = 5x10-6 Km = 1x10-5 Km = 5x10-5 Km = 1x10-4 Km = 5x10-4 Hydraulic head error 0 1000 2000 3000 Distance to the Gallusquelle spring (m) 0.01 0.02 0.03 0.04 0.05 Fl ow v el oc ity (m s -1 ) Legend Km = 1x10–5 Km = 5x10–5 Km = 1x10–4 Conduit flow velocity for tracer test 1 0 2000 4000 6000 8000 10000 Distance to the Gallusquelle spring (m) 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Fl ow v el oc ity (m s -1 ) Legend Km = 1x10–5 Km = 5x10–5 Km = 1x10-4 Conduit flow velocity for tracer test 2 Objective functions in relation to the hydraulic conductivity of the fissured matrix Km b) c)a) 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. mainly depends on the conduit volume and accordingly on the Manning coefficient n (Fig. 5). It is possible to simu- late 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 pa- rameters. 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 simu- lation 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 cor- responds 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 mech- anism is triggered, if the conduit size is increased. A higher conduit size leads to higher groundwater influx from the fis- sured matrix and spring discharge is overestimated. There- fore, parameter analysis shows that scenario 1 is too strongly simplified to correctly reproduce the complex nature of the aquifer. 4.2 Scenario 2 – conduit roughness coefficient Kc In scenario 2 the Manning coefficient n was changed from constant to laterally variable. In the literature, n is generally 0 1 2 3 4 5 Manning coefficient n (s m–1/3) -20 0 20 40 60 Pe ak -o ffs et ti m e (h ) Peak-offset time for TT 1 in relation to the Manning coefficient n 0 1 2 3 4 5 Manning coefficient n (s m–1/3) -80 -40 0 40 80 Pe ak -o ffs et ti m e (h ) Peak-offset time for TT 2 in relation to the Manning coefficient n a) b) 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. kept constant throughout the conduit network (e.g. Jeannin, 2001; Reimann et al., 2011) for lack of information on con- duit 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 small- interconnected 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, can- not 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 sec- tion 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- www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 902 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 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 in- cluded 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 im- plicitly leads to a larger roughness in the simulation, even further enhanced by exchange processes between the indi- vidual 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 impossi- ble to calculate the change of n directly from the geometry. Thus, a simple scenario was assumed where the roughness coefficient Kc, which is the reciprocal of n, was linearly and negatively coupled to the rising conduit radius (Eq. 13). Kc = 1 n =−mhrc+mhrc,max+ bh, (13) where rc [L] is the conduit radius and rc,max [L] the max- imum conduit radius simulated for the respective spring, which COMSOL® calculates from Eq. (1). mh [L−2/3 T−1] and bh [L1/3 T−1] are calibration parameters determining the slope and the lowest value of the roughness coefficient re- spectively. For every conduit geometry several combinations of mh and bh lead to the same spring discharge. However, hydraulic head fit and tracer velocities are different for each mh–bh combination even if spring discharge is the same. With the new parameters a higher variation of velocity profiles is pos- sible. This allows for the calibration of the tracer velocities of both tracer tests. The dependence of tracer test 2 on mh is much higher than that of tracer test 1 since it is injected fur- ther 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 mh values and always shows a significant positive correlation with mh (Fig. 6). Since the slope of Kc is negative with respect to the con- duit length, the variable Kc leads to a slowing down of wa- ter 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 hy- draulic head fits can only be obtained, if the decrease of Kc towards the spring is not too large and compensates the effect of the increase in conduit transmissivity due to the increas- ing conduit radius. This effect reduces the number of possi- ble and plausible parameter combinations. From these con- siderations a best-fit model can be deduced capable of repro- ducing all objective functions within the given error ranges (Fig. 7a). According to the model simulations, karst ground- water discharge and flow velocities significantly depend on 0.1 1 10 Slope of the roughness coefficient of the karst conduit system mh -20 0 20 40 R oo t m ea n sq ua re e rro r o f t he h yd ra ul ic h ea d di st rib ut io n (m ) -20 0 20 40 Pe ak -o ffs et ti m e (h ) Legend RMSE of the hydraulic head distribution peak-offset time of TT 1 peak-offset time of TT 2 Fitting values in relation to mh Figure 6. Hydraulic head errors and differences between peak con- centration 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 ofm= 2× 10−4. Eachmh (m−2/3 s−1) value corresponds to a respective value of the highest conduit roughness bh (m1/3 s−1) and each combination results in the same spring dis- charge. the total conduit volume as is to be expected. It can be de- duced from the parameter analysis that the conduit volume can be estimated at ca. 100 000 m3 for the different parame- ters to match equally well (Fig. 7a). 4.3 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 con- duit 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 ob- served 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 karstifi- cation is backed up by the occurrence of sinkholes (Mohrlok and Sauter, 1997). Therefore, it represents the minimal ex- tent 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 Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 903 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 matrixKm. 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. Cal- culation 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 am- biguity introduced by the uncertainty in total conduit length is small if hydraulic conduit parameters and total conduit vol- umes are the aim of investigation. 4.4 Scenario 4 – matrix hydraulic conductivity Km In scenario 4, the homogeneously chosen hydraulic conduc- tivity of the fissured matrix Km was changed into a later- ally variable conductivity based on different types of lithol- ogy 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 con- ductivities. Oehlmann et al. (2013) discussed that the ma- jor influence is the conduit geometry leading to higher hy- draulic 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 con- ductivities (Sauter, 1992). Figure 9 shows the division into three different areas. Km values were varied in the range of the values measured by Sauter (1992). It was expected that a laterally variable Km value has a major influence on the hydraulic head distribution. All vari- ations of scenario 2 that produce good results for both tracer tests and have a high total conduit volume above 100 000 m3 yield poor results for hydraulic head errors and spatial dis- www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 904 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching Bitz Burladingen Winterlingen Neufra Gammertingen Albstadt/Ebingen Veringenstadt ±0 2 4Kilometers Fehla Lauchert Schmiecha Gallusquelle Legend Gallusquelle spring extent scenarios 1 and 2 scenario 3 dry valley settlement river fault Gallusquelle catchment area Highly conductive conduit network Figure 8. Extended conduit system for scenario 3. The conduit configuration (extent) that is used for the other scenarios is marked in red. Fehla-Ursprung Balinger Quelle B9 B8 B7 B4 B2 B25 B24 B22 B21 B19 B18 B17 B16 B15 B14 B13 B12B11 B10 Abendrain Bitz Burladingen Winterlingen GammertingenNeufra Albstadt/Ebingen Veringenstadt ±0 2 4Kilometers Legend springsdischarge [m3 s–1 ] highly conductive pipe network river fault catchment area Gallusquelle 0 - 0.1 0.1 - 0.2 0.2 - 0.3 0.3 - 0.4 0.4 - 0.5 settlement Fehla Lauchert Schmiecha Gallusquelle Ahlenbergquelle Büttnau-quellen Schlossberg-quelle Bronnen Königsgassenquelle observation well south-eastern part:high conductivity central part:low conductivity north-western part:medium conductivity Figure 9. Model catchment with spatially distributed hydraulic conductivities. The model area is divided into three parts after geologic aspects. For each segment different values of the hydraulic conductivity were examined during parameter analysis in scenario 4. Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 905 tributions of the hydraulic heads (Fig. 7a). For scenario 4, two different conduit configurations (geometries) were cho- sen that achieve good results with respect to conduit flow ve- locities. Geometry G1 has a conduit volume of 112 000 m3. G2 has a higher b value which leads to the maximum con- duit volume of ca. 150 000 m3. All parameters for the two simulations are given in Table 4. It was found that while the maximum root mean square er- ror 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 Km values (Fig. 7c), which assists in the independent conduit network and fissured ma- trix calibration. This observation increases the confidence in the representation of the conduits and improves the possibil- ity 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 ob- vious 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. 4.5 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 con- duits 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 condi- tions are turbulent the solution is limited by a diffusion dom- inated layer that gets thinner with increasing flow velocities (Clemens, 1998). Clemens (1998) simulated karst evolution in simple Y-shaped conduit networks and found higher di- ameters for the downstream conduit even after short simu- lation times. Preferential conduit widening at intersections could further be enhanced by the process of mixing corro- sion (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 re- spect to Ca2+ compared to water entering the system through direct recharge. Thus, if direct recharge is present, the mix- ing with nearly saturated water from an intersecting conduit could hamper the preferential evolution of the conduit down- stream slowing down the aforementioned processes. In sce- nario 5 the influence of an increase in diameter at conduit intersections was investigated. Since the amount of preferen- tial 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, bh the highest con- duit roughness, mh the slope of roughness decrease away from the spring and V the conduit volume. Geometry 1 Geometry 2 b (m) 0.01 0.5 m (–) 2.07× 10−4 1.5× 10−4 bh (m1/3 s−1) 0.17 0.15 mh (m−2/3 s−1) 0.4 0.6 V (m3) 112 564 153 435 radius was then calculated according to Eq. (14) at each in- tersection. rc2 = √ r2c0+ r2c1, (14) where rc2 is the conduit radius downstream of the intersec- tion and rc0 and rc1 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 pa- rameters (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 vari- ability 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 m3 (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 sce- nario 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 trans- port in tracer test 2 occurs in relatively small conduits with high flow velocities and larger conduits with lower veloc- ities, 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 dia- gram at negative values below −10 h. However, since the fit for conduit volumes around 100 000 m3 is similar to that of scenario 2, the two scenarios can in this case not be distin- guished based on field observations. www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 906 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching Figure 10. Comparison of the best-fit simulations with field data for scenarios 2 and 5. (a) Breakthrough curve of tracer test 1, (b) break- through curve of tracer test 2 and (c) spring discharge. 5 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 dis- charges (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 underesti- mates 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 Km values being equal, which basically results in the parameter set of scenario 2. This implies that the ma- jor 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 conduc- tivity 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 pa- rameters 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 dur- ing 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. 6 Discussion 6.1 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 parame- ters and to check the plausibility of the model set-up. Pos- sible 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 ob- served in the field. As seen in Table 1 most of the calibrated parameters range well within values provided in the litera- ture. The calibrated Manning coefficients are relatively high compared to other karst systems. Jeannin (2001) lists effec- tive 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 Jean- nin (2001). This is not surprising since the calibrated n value reflects the total roughness of the conduit bundles and there- fore includes geometric conduit properties in addition to the Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 907 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 de- rived from scenarios 2 and 5 is only 50 % of that estimated with traditional methods (Geyer et al., 2008). Since the con- duit 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 matrix– conduit 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 er- ror incurred by deducing the conduit volume without tak- ing 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 m2. For scenario 2 the simulated average cross sectional area is 11.9 m2 and for scenario 5 13.4 m2, which compares very well with the results of Birk et al. (2005). It was not possible to match the shape of both break- through curves with the same dispersivity. The apparent dis- persion in the tracer test 2 breakthrough is much higher com- pared 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 ma- trix 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 quan- tify 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 in- fluencing the observed retardation and dispersion effects is the input mechanism. The simulation assumes that all in- troduced tracers immediately and completely enter the con- duit 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 disper- sivity 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 sim- ulations. The slight mass difference results from a combi- nation 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 m3 s−1and is not expected to signifi- cantly influence the general flow conditions. It probably re- sults 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 overesti- mates 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 intersec- tions (Fig. 1). Therefore, the conduit volume of the Ahlen- berg and Büttnauquellen springs is 134 568 m3 in scenario 2 and only 75 085 m3 in scenario 5 leading to the different dis- charge values. It is reasonable to assume that a better fit for the spring group can be achieved, if more variations of con- duit intersections are tested. An adequate fit for the Fehla- Ursprung spring of 0.1 m3 s−1 was achieved for both scenar- ios with a fault aperture of 0.005 m. 6.2 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 ap- proach was chosen because tracer tests are supposed to be conducted during quasi-steady-state flow conditions. How- ever, 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 mod- els 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 uni- form 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 dispersiv- ities do not only represent geometrical heterogeneities but also temporal effects as is the case for all standard evalua- tions of dispersion from tracer breakthrough curves. www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 908 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 0 2000 4000 6000 8000 Distance to the Gallusquelle spring (m) 0 0.02 0.04 0.06 0.08 0.1 Fl ow v el oc ity (m s –1 ) Legend no direct recharge 10% direct recharge Scenario 2 0 2000 4000 6000 8000 Distance to the Gallusquelle spring (m) 0 0.02 0.04 0.06 0.08 0.1 Fl ow v el oc ity (m s –1 ) Legend no direct recharge 10% direct recharge Scenario 5 Flow velocities inside the karst conduits with and without a direct recharge component a) b) Figure 11. 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. 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 max- imum percentage of direct recharge of 10 % estimated by Sauter (1992) and Geyer et al. (2008) was used for this sim- ulation. Neither for scenario 2 nor for scenario 5 a gradient reversal between conduit and matrix occurred and the influ- ence 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 sur- face / 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 con- siderably, this fact is believed not to influence the overall re- sults and conclusions of this study, i.e. the relative signifi- cance of the parameters deduced from parameter analysis and the deduced conduit volume, especially since flow is simu- lated for steady-state conditions. For all distributed numerical karst simulations, uncertain- ties regarding the exact positions and interconnectivities of the conduit branches still remain. Due to the extensive inves- tigations 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 in- tegration of any future information that might enhance the numerical model further. 6.3 Calibration strategy For a successful calibration of a distributed groundwater flow and transport model for a karst area on catchment scale cer- tain 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 func- tions 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 mea- surements 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 de- creases 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 ob- served: for a representative calibration, the dye should be in- jected at as large a distance to each other as possible with one of them including the length of the whole conduit sys- tem. 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 Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 909 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 di- vides are majorly determined by the hydraulic conductivity of the fissured matrix Km, 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 Km values. In this study, high Km values above ca. 3× 10−5 m s−1 made the simulation of the spring dis- charge of the Fehla-Ursprung spring (Fig. 1) impossible be- cause 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 catch- ment area for the Gallusquelle spring. There are eight parameters available for model calibration in this study. Two of these parameters define the conduit ge- ometry: b is the lowest conduit radius and m the slope with which the conduit radius increases. One parameter, df, de- fines the aperture of the fault zone. The hydraulic conduc- tivity of the fissured matrix is represented by the parameter Km and the roughness of the conduit system by two parame- ters: bh represents the highest roughness and mh 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 experi- ments (Table 1). For efficiency reasons it is important to know which of these parameters can be calibrated independently. The ap- parent transport dispersivities ε1 and ε2 are pure transport pa- rameters, which influence only the shape of the breakthrough curves and not the flow field. The hydraulic model parame- ters 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 df 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 cali- brated at the beginning of the simulation procedure to repro- duce the catchment and the discharge of the Fehla-Ursprung spring adequately and kept constant throughout all the sim- ulations. In the final calibrated models it was rechecked, but the calibrated value was still acceptable. The hydraulic conductivity of the fissured matrix Km can be calibrated independently in principle as well. The in- fluence on spring discharge is relatively small. The best-fit Km value depends on the conduit parameters, i.e. geome- try 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 be- tween the simulated conduit geometries, the more likely is a slight shift of the best-fit Km value. Therefore, it is ad- visable 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 Km 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 sensi- tive 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 ve- locities, e.g. by use of environmental tracers, would likely be sensitive to the Km values as well. Therefore, the choice of the flow regime and the objective functions determines the strength of the interdependencies between fissured ma- trix and conduit system parameters and therefore whether Km can be calibrated independently. The conduit parameters for geometry and roughness, here four parameters (lowest conduit radius b, slope of radius in- crease m, highest roughness bh and slope of roughness de- crease mh), 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 bh– mh combinations that result in the observed spring discharge. In general, the slowest transport velocities are achieved with a mh value of zero. So, to deduce the range of geometric parameters that reproduce the objective functions, it is ad- visable to check the minimum conduit volume for which the tracer tests are not too fast for a value ofmh equal to zero. For the Gallusquelle area, transmissivities significantly increase towards the springs, which is characteristic for most karst catchments. Therefore, low bh values oppose the general hy- draulic head trend: they increase the conduit roughness at the spring leading to slower flow and higher gradients. The higher the conduit volume, the higher bh is required to repro- duce 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 m3 for the Gal- lusquelle 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- www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 910 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. 7 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 pa- rameters provides new insight into the hydraulics of the Gal- lusquelle conduit system. The model ambiguity was signif- icantly reduced to the point where an estimation of the ac- tual 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 dis- charge, the hydraulic head distribution and two tracer tests. The model allows for the identification of the relevant pa- rameters affecting karst groundwater discharge and transport in karst conduits and the examination of the respective over- all importance in a well-investigated karst groundwater basin for steady-state flow conditions. While a differentiated rep- resentation of the roughness values in the karst conduits is substantial for buffering the lack of knowledge of the ex- act conduit geometry, e.g. local variations in cross section and the number of interacting conduits, variable matrix hy- draulic 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 over- all 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 practi- cally 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 simula- tion is the neglect of flow dynamics, which limits the appli- cability 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 es- sential field conditions and lead to further conclusions re- garding the important karst system parameters, their influ- ences on karst hydraulics and their interdependencies. It can be expected that some parameters, which are of minor im- portance 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. Acknowledgements. The presented study was funded by the German Federal Ministry of Education and Research (promotional reference no. 02WRS1277A, AGRO, Risikomanagement von Spurenstoffen und Krankheitserregern in ländlichen Karsteinzugs- gebieten). This open-access publication is funded by the University of Göttingen. Edited by: M. Giudici References Barenblatt, G. I., Zheltov, I. P., and Kochina, I. N.: Basic concepts in the theory of seepage in fissured rocks (strata), J. Appl. Math. Mech.-USS, 24, 1286–1303, 1960. Bauer, S., Liedl, R., and Sauter, M.: Modeling of karst aquifer gen- esis: Influence of exchange flow, Water Resour. Res., 39, 1285, doi:10.1029/2003WR002218, 2003. Birk, S., Geyer, T., Liedl, R., and Sauter, M.: Process-Based Inter- pretation of Tracer Tests in Carbonate Aquifers, Ground Water, 43, 381–388, 2005. Birk, S., Liedl, R., and Sauter, M.: Karst Spring responses examined by process-based modeling, Ground Water, 44, 832–836, 2006. Clemens, T.: Simulation der Entwicklung von Karstaquiferen, PhD thesis, Eberhard-Karls-Universität zu Tübingen, Tübingen, 1998. Clemens, T., Hückinghaus, D., Sauter, M., Liedl, R., and Teutsch, G.: A combined continuum and discrete network reactive trans- port model for the simulation of karst development, in: Proceed- ings of the ModelCARE 96 Conference, 24–26 September 1996, Golden, Colorado, USA, 309–318, 1996. COMSOL AB: COMSOL Multiphysics® User’s Guide v4.3, 1292 pp., 2012. Coronado, M., Ramírez-Sabag, J., and Valdiviezo-Mijangos, O.: On the boundary conditions in tracer transport models for fractured porous underground formations, Rev. Mex. Fís., 53, 260–269, 2007. Covington, M. D., Luhmann, A. J., Wicks, C. M., and Saar, M. O.: Process length scales and logitudinal damping in karst conduits, J. Geophys. Res., 117, F01025, doi:10.1029/2011JF002212, 2012. Doummar, J., Sauter, M., and Geyer, T.: Simulation of flow pro- cesses in a large scale karst system with an integrated catch- ment model (Mike She) – Identification of relevant param- eters influencing spring discharge, J. Hydrol., 426, 112–123, doi:10.1016/j.jhydrol.2012.01.021, 2012. Dreybrodt, W.: Mixing in CaCO3–CO2–H2O systems and its role in the karstification of limestone areas, Chem. Geol., 32, 221–236, 1981. Ford, D. C. and Williams, P. W.: Karst geomorphology and hydrol- ogy, Wiley, West Sussex, 562 pp., 2007. Geyer, T., Birk, S., Licha, T., Liedl, R., and Sauter, M.: Multi-tracer test approach to characterize reactive transport in karst aquifers, Ground Water, 45, 36–45, 2007. Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/ S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching 911 Geyer, T., Birk, S., Liedl, R., and Sauter, M.: Quantification of tem- poral distribution of recharge in karst systems from spring hy- drographs, J. Hydrol., 348, 452–463, 2008. Golwer, A.: Erläuterungen zu Blatt 7821 Veringenstadt, Geologische Karte 1 : 25 000 von Baden-Württemberg, Geologisches Landesamt Baden-Württemberg, Stuttgart, 151 pp., 1978. Gwinner, M. P.: Erläuterungen zu Blatt 7721 Gammertin- gen, Geologische Karte 1 : 25 000 von Baden-Württemberg, Geologisches Landesamt Baden-Württemberg, Freiburg, Stuttgart, 78 pp., 1993. Hartmann, A., Weiler, M., Wagener, T., Lange, J., Kralik, M., Humer, F., Mizyed, N., Rimmer, A., Barberá, J. A., Andreo, B., Butscher, C., and Huggenberger, P.: Process-based karst mod- elling to relate hydrodynamic and hydrochemical characteristics to system properties, Hydrol. Earth Syst. Sci., 17, 3305–3321, doi:10.5194/hess-17-3305-2013, 2013. Hartmann, A., Goldscheider, N., Wagener, T., Lange, J., and Weiler, M.: Karst water resources in a changing world: Review of hydrological modeling approaches, Rev. Geophys., 52, 1–25, doi:10.1002/2013RG000443, 2014. Hauns, M., Jeannin, P.-Y., and Atteia, O.: Dispersion, retardation and scale effect in tracer breakthrough curves in karst conduits, J. Hydrol., 241, 177–193, 2001. Hu, R.: Hydraulic tomography: A new approach coupling hydraulic travel time, attenuation and steady shape inversions for high- spatial resolution aquifer characterization, PhD thesis, University of Göttingen, Göttingen, 116 pp., 2011. Hubinger, B. and Birk, S.: Influence of initial heterogeneities and recharge limitations on the evolution of aperture distributions in carbonate aquifers, Hydrol. Earth Syst. Sci., 15, 3715–3729, doi:10.5194/hess-15-3715-2011, 2011. Hückinghaus, D.: Simulation der Aquifergenese und des Wärmetransports in Karstaquiferen, C42, Tübinger Ge- owissenschaftliche Arbeiten, Tübingen, 1998. Hunter, N. M., Bates, P. D., Horritt, M. S., De Roo, A. P. J., and Werner, M. G. F.: Utility of different data types for calibrat- ing flood inundation models within a GLUE framework, Hy- drol. Earth Syst. Sci., 9, 412–430, doi:10.5194/hess-9-412-2005, 2005. Jeannin, P.-Y.: Modeling flow in phreatic and epiphreatic karst con- duits in the Hölloch cave (Muotatal, Switzerland), Water Resour. Res., 37, 191–200, 2001. Jeannin, P.-Y. and Sauter, M.: Modelling in karst systems, Bulletin d’Hydrogéologie, 16, Université de Neuchâtel, Neuchâtel, 1998. Khu, S.-T., Madsen, H., and di Pierro, F.: Incorporating multiple observations for distributed hydrologic model calibration: An ap- proach using a multi-objective evolutionary algorithm and clus- tering, Adv. Water Resour., 31, 1387–1398, 2008. Kordilla, J., Sauter, M., Reimann, T., and Geyer, T.: Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach, Hydrol. Earth Syst. Sci., 16, 3909–3923, doi:10.5194/hess-16-3909-2012, 2012. Kovács, A. and Sauter, M.: Modelling karst hydrodynamics, in: Methods in karst hydrogeology, edited by: Goldscheider, N. and Drew, D., Taylor and Fracis, London, 201–222, 2007. Kovács, A., Perrochet, P., Király, L., and Jeannin, P.-Y.: A quanti- tative method for the characterisation of karst aquifers based on spring hydrograph analysis, J. Hydrol., 303, 152–164, 2005. Li, H. T., Brunner, P., Kinzelbach, W., Li, W. P., and Dong, X. G.: Calibration of a groundwater model using pattern infor- mation from remote sensing data, J. Hydrol., 377, 120–130, doi:10.1016/j.jhydrol.2009.08.012, 2009. Liedl, R., Sauter, M., Hückinghaus, D., Clemens, T., and Teutsch, G.: Simulation of the development of karst aquifers using a cou- pled continuum pipe flow model, Water Resour. Res., 39, 1057, doi:10.1029/2001WR001206, 2003. Luhmann, A. J., Covington, M. D., Alexander, S. C., Chai, S. Y., Schwartz, B. F., Groten, J. T., and Alexander, E. C.: Comparing conservative and nonconservative tracers in karst and using them to estimate flow path geometry, J. Hydrol., 448–449, 201–211, doi:10.1016/j.jhydrol.2012.04.044, 2012. Madsen, H.: Parameter estimation in distributed hydrological catch- ment modelling using automatic calibration with multiple objec- tives, Adv. Water Resour., 26, 205–216, 2003. Merkel, P.: Karsthydrologische Untersuchungen im Lauchertgebiet (westl. Schwäbische Alb), Diplom thesis, University of Tübin- gen, Tübingen, 108 pp., 1991. Mohrlok, U.: Numerische Modellierung der Grundwasserströ- mung im Einzugsgebiet der Gallusquelle unter Festlegung eines Drainagesystems, Grundwasser, 19, 73–85, doi:10.1007/s00767- 013-0249-x, 2014. Mohrlok, U. and Sauter, M.: Modelling groundwater flow in a karst terraine using discrete and double-continuum approaches: im- portance of spatial and temporal distribution of recharge, in: Proceedings of the 12th International Congress of Speology, 2/6th Conference on Limestone Hydrology and Fissured Media, 10–17 August 1997, La Chaux-de-Fonds, Switzerland, 167–170, 1997. Oehlmann, S., Geyer, T., Licha, T., and Birk, S.: Influence of aquifer heterogeneity on karst hydraulics and catchment delineation em- ploying distributive modeling approaches, Hydrol. Earth Syst. Sci., 17, 4729–4742, doi:10.5194/hess-17-4729-2013, 2013. Ophori, D. U.: Constraining permeabilities in a large-scale ground- water system through model calibration, J. Hydrol., 224, 1–20, 1999. Perrin, C., Andréassian, V., Serna, C. R., Mathevet, T., and Le Moine, N.: Discrete parameterization of hydrologi- cal models: Evaluating the use of parameter sets libraries over 900 catchments, Water Resour. Res., 44, W08447, doi:10.1029/2007WR006579, 2008. Rehrl, C. and Birk, S.: Hydrogeological characterisation and mod- elling of spring catchments in a changing environment, Aust. J. Earth Sci., 103/2, 106–117, 2010. Reiber, H., Klein, F., Selg, M., and Heidland, S.: Hydrogeologische Erkundung Baden-Württemberg – Mittlere Alb 4 – Markierungsversuche, Abwassereinleitungen, Landesamt für Umwelt, Messungen und Naturschutz Baden-Württemberg, Tübingen, 71 pp., 2010. Reimann, T., Rehrl, C., Shoemaker, W. B., Geyer, T., and Birk, S.: The significance of turbulent flow representation in single-continuum models, Water Resour. Res., 47, W09503, doi:10.1029/2010WR010133, 2011. Saller, S. P., Ronayne, M. J., and Long, A. J.: Comparison of a karst groundwater model with and without discrete conduit flow, Hydrogeol. J., 21, 1555–1566, doi:10.1007/s10040-013-1036-6, 2013. www.hydrol-earth-syst-sci.net/19/893/2015/ Hydrol. Earth Syst. Sci., 19, 893–912, 2015 912 S. Oehlmann et al.: Reducing the ambiguity of karst aquifer models by pattern matching Sauter, M.: Quantification and Forecasting of Regional Groundwa- ter Flow and Transport in a Karst Aquifer (Gallusquelle, Malm, SW Germany), C13, Tübinger Geowissenschaftliche Arbeiten, Tübingen, 1992. Schmidt, S., Geyer, T., Marei, A., Guttman, J., and Sauter, M.: Quantification of long-term wastewater impacts on karst ground- water resources in a semi-arid environment by chloride mass bal- ance methods, J. Hydrol., 502, 177–190, 2013. Strayle, G.: Karsthydrologische Untersuchungen auf der Ebinger Alb (Schwäbischer Jura), in: Jahreshefte des Geologischen Lan- desamtes Baden-Württemberg, Freiburg im Breisgau, 109–206, 1970. Teutsch, G. and Sauter, M.: Groundwater Modeling in karst ter- ranes: scale effects, data aquisition and field validation, in: Pro- ceedings of the 3rd Conference on Hydrogeology, Ecology, Mon- itoring and Management of Ground Water in Karst Terranes, 4– 6 December 1991, Nashville, USA, 17–34, 1991. Villinger, E.: Über Potentialverteilung und Strömungssysteme im Karstwasser der Schwäbischen Alb (Oberer Jura, SW- Deutschland), Geologisches Jahrbuch, C18, Bundesanstalt für Geowissenschaften und Rohstoffe und Geologische Landesämter der Bundesrepublik Deutschland, Hannover, 1977. Villinger, E.: Hydrogeologie, in: Erläuterungen zu Blatt 7721 Gammertingen, Geologische Karte 1 : 25 000 von Baden- Württemberg, edited by: Gwinner, M. P., Geologisches Lan- desamt Baden-Württemberg, Freiburg, Stuttgart, 30–57, 1993. Worthington, S. R. H.: Diagnostic hydrogeologic characteristics of a karst aquifer (Kentucky, USA), Hydrogeol. J., 17, 1665–1678, doi:10.1007/s10040-009-0489-0, 2009. Hydrol. Earth Syst. Sci., 19, 893–912, 2015 www.hydrol-earth-syst-sci.net/19/893/2015/