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/