Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ doi:10.5194/hess-16-3909-2012 © Author(s) 2012. CC Attribution 3.0 License. Hydrology and Earth System Sciences Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach J. Kordilla1, M. Sauter1, T. Reimann2, and T. Geyer1 1Geoscientific Centre, University of Go¨ttingen, Go¨ttingen, Germany 2Institute for Groundwater Management, TU Dresden, Dresden, Germany Correspondence to: J. Kordilla (jkordil@gwdg.de) Received: 16 January 2012 – Published in Hydrol. Earth Syst. Sci. Discuss.: 1 February 2012 Revised: 25 September 2012 – Accepted: 5 October 2012 – Published: 30 October 2012 Abstract. The objective of this work is the simulation of saturated and unsaturated flow in a karstified aquifer using a double continuum approach. The HydroGeoSphere code (Therrien et al., 2006) is employed to simulate spring dis- charge with the Richards equations and van Genuchten pa- rameters to represent flow in the (1) fractured matrix and (2) conduit continuum coupled by a linear exchange term. Rapid vertical small-scale flow processes in the unsaturated conduit continuum are accounted for by applying recharge boundary conditions at the bottom of the saturated model do- main. An extensive sensitivity analysis is performed on sin- gle parameters as well as parameter combinations. The tran- sient hydraulic response of the karst spring is strongly con- trolled by the matrix porosity as well as the van Genuchten parameters of the unsaturated matrix, which determine the head dependent inter-continuum water transfer when the con- duits are draining the matrix. Sensitivities of parameter com- binations partially reveal a non-linear dependence over the parameter space. This can be observed for parameters not belonging to the same continuum as well as combinations, which involve the exchange parameter, showing that results of the double continuum model may depict a certain de- gree of ambiguity. The application of van Genuchten pa- rameters for simulation of unsaturated flow in karst systems is critically discussed. 1 Introduction Discharge dynamics in karst aquifers are determined by superposition of several effects: (1) water infiltration into soil, (2) water percolation through the unsaturated zone, (3) groundwater flow in highly conductive karst conduits and interaction with (4) groundwater flow in the low-conductive fissured and fractured rock matrix. These different effects, without even having considered the variability of precipi- tation and evapotranspiration, are a result of the particu- lar properties of the individual compartments: soil-epikarstic zone, vadose zone, and phreatic zone. Each of these com- partments is, in turn, characterized by two coupled flow sys- tems: a highly permeable one with low storage and a less per- meable one with high storage. Therefore, different individ- ual (rapid, slow) flow components with characteristic tem- poral distributions are induced. Accordingly, the final spring discharge is then a function of the individual flow contri- butions of each of these compartments (Smart and Hobbs, 1986), which makes the inverse analysis of spring discharge a major challenge, requiring elaborate modeling tools and a large spectrum of data to constrain the model. The simula- tion of coupled saturated and unsaturated flow is still a chal- lenge in hydrogeology in particular in fractured (Therrien and Sudicky, 1996) and karstified systems (Reimann et al., 2011a). This is predominantly a result of the data scarcity respecting the hydraulic parameter field of real karst sys- tems. Therefore, flow in karst systems is often simulated with lumped parameter modeling approaches, which trans- late precipitation signals to discharge hydrographs by ap- plying simple transfer functions (Dreiss, 1989). Generally, this type of approach is appropriate for situations in which predicted system states are expected to range between al- ready observed events. The simulation of natural karst sys- tems with distributed parameter models is reported only in a few studies (e.g. Jeannin, 2001; Hill and Polyak, 2010). However, distributive modeling approaches incorporate flow Published by Copernicus Publications on behalf of the European Geosciences Union. 3910 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems laws and, therefore, are adequate for the process based sim- ulation of karst hydraulics (e.g. Birk et al., 2006; Reimann et al., 2011b) and transport problems (e.g. Dreybrodt et al., 2005; Birk et al., 2005). Teutsch and Sauter (1991) demon- strate how far the different mathematical model approaches are suitable for different types of problems (flow, transport, regional, local). An approach that takes into account the lim- ited information about aquifer geometry and still allows for the simulation of the dynamics of the karst system at an event basis, i.e. considers the dual flow behavior of karst systems is the double continuum approach (e.g. Teutsch and Sauter, 1991; Sauter et al., 2006). The approach was introduced by Barenblatt et al. (1960) and applied for simulation of karst hydraulics on catchment scale by Teutsch (1988) and Sauter (1992). It yields equations for simulation of slow and dif- fuse flow in the fissured matrix and the discrete rapid un- derground drainage by solution conduits in karst systems. Here, we want to assess the relative importance of individ- ual factors and parameter combinations on the discharge be- havior of a karst spring without detailed knowledge about the hydraulic parameter field of an aquifer system. This type of information is of major importance to focus characteriza- tion efforts in catchment based karst studies. Furthermore, the importance of infiltration dynamics, i.e. the temporal dis- tribution of the rapid and the slow flow component on the discharge dynamics is to be determined. We employ the in- tegrated saturated-unsaturated double-continuum approach HydroGeoSphere (Therrien et al., 2006) to simulate recharge and discharge dynamics in a karst aquifer with a thick un- saturated zone. Flow simulations are based on the Richards equation and respective parameters are described via the van Genuchten parametric model. Given that the van Genuchten parameters have been developed for porous media on REV scales, parameters obtained for the large scale conduit system are pure calibration parameters without physical meaning in nature. Furthermore, the Richards equation is not able to ex- press the highly nonlinear flow regime observed in unsatu- rated fractured rocks, for example, rivulet flow (Dragila et al., 2004; Germann et al., 2007). The application and limitations of the approach for flow simulation in karst systems are dis- cussed. A comprehensive parameter study was conducted in order to elucidate sensitive and important model parameters as well as parameter dependencies, and to reduce the model ambiguity to assist in focused karst characterization. 2 Methods 2.1 Modeling approach The application of the double-continuum approach requires two sets of flow equations, one for the matrix (primary) and one for the conduit (secondary) continuum, solved consecu- tively at the same node and coupled with an exchange term that defines the hydraulic interface and controls the inter- continuum exchange flow. The applied HydroGeoSphere model (Therrien et al., 2006) is a non-commercial code avail- able to the interested user under http://hydrogeosphere.org/. The model has been extensively used for various studies in- volving dual porosities such as Mclaren et al. (2000), Rosen- bom et al. (2009) and Schwartz et al. (2010). The govern- ing equation in the applied model is the Richards equation (Richards, 1931), which is slightly modified to account for inter-continuum water exchange: −∇wm(qm)+0ex ±Rm = wm ∂ ∂t (θsmSwm) (1) −∇wc(qc)+0ex ±Rc = wc ∂ ∂t (θscSwc), (2) where wm and wc are the volumetric fractions of each con- tinuum of the total porosity, such that wm = 1.0−wc. Swm and Swc are the water saturations of the respective continuum andRm andRc denote a volumetric fluid flux per unit volume (source/sink term) for each continuum. The saturated water content of the matrix and conduit system are assumed equal to the the effective matrix porosity θsm and conduit porosity θsc and are related to the water content of the matrix θm and of the conduit θc according to θm = Swmθsm (3) θc = Swcθsc. (4) The conduit and total porosity are given as θtotal = θsm(1−wc)+ θscwc = θsm(wm)+ θscwc, (5) i.e. the whole simulation domain consists of nodes with pri- mary porosity θsm with a volumetric fraction of wm = 1.0− wc and secondary porosity θsc with a volumetric fraction of wc. Given that the local conduit porosity is chosen to be 1.0 and that both continua cover the whole domain the overall conduit porosity can simply be evaluated as: θsc ∧= θsc(local)wc. (6) The fluxes qm and qc are obtained from qm =−Kmkrm∇(ψm + z) (7) qc =−Kckrc∇(ψc + z), (8) where Km and Kc denote hydraulic conductivity, ψm and ψc are the pressure heads in each continuum and z is the eleva- tion head. In the unsaturated zone, the relative permeabilities krm, krc and kri (interface) depend on the water saturation which in turn is related to the pressure head according to van Genuchten (1980): Swm = Swrm + (1− Swrm) [ 1+ |αmψm|nm ]−mm (9) Swc = Swrc + (1− Swrc) [ 1+ |αcψc|nc ]−mc (10) Swi = Swri + (1− Swri) [ 1+ |αiψm|ni ]−mi (11) Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3911 for ψ < 0, where Swrm, Swrc and Swri are the residual satu- rations, αm, αc and αi denote the inverse air-entry pressure head, nm, nc and ni are the pore-size distribution indices of each continuum and the interface. Note that the evaluation of the interface relative conductivity is based on the pres- sure head of the matrix. In the saturated zone where ψ ≥ 0, the saturations are Swm = Swc = Swi = 1. The relative per- meability is given by krm(Swm) = S(lp)em [ 1− ( 1− S1/mmem )mm]2 (12) krc(Swc) = S(lp)ec [ 1− ( 1− S1/mcec )mc]2 (13) kri(Swi) = S(lp)ei [ 1− ( 1− S1/miei )mi]2 (14) with lp being the pore connectivity parameter (equals 0.5 af- ter Mualem, 1976), Se the effective saturation Sem = Swm − Swrm1− Swrm (15) Sec = Swc − Swrc1− Swrc (16) Sei = Swi − Swri1− Swri (17) and m is defined as mm = 1− 1 nm (18) mc = 1− 1 nc (19) mi = 1− 1 ni (20) for n > 1. In the saturated zone, the storage terms on the right-hand side of Eqs. (1) and (2) are replaced by SwmSsm ∂ψm ∂t + θsm ∂Swm ∂t (21) SwcSsc ∂ψc ∂t + θsc ∂Swc ∂t , (22) where Ssm and Ssc are the specific storage coefficients. Water release by compaction of the porous medium is neglected in the unsaturated zone. The term 0ex in Eqs. (1) and (2) de- scribes the volumetric fluid exchange rate per unit volume between primary and secondary continuum and is given as 0ex = αexKikri (ψc −ψm) , (23) where Ki is the hydraulic conductivity of the interface (e.g. sediments) and kri the relative interface permeability (Baren- blatt et al., 1960). The exchange parameter αex is determined by calibration and defined as (Gerke and Van Genuchten, 1993): αex = β α2 γw, (24) where β is a geometry factor (3 for rectangular matrix blocks, 15 for spheres), a is the distance between the center of a ma- trix block and the adjacent fracture or conduit and γw is an empirical coefficient usually set to 0.4. As mentioned in the introduction the van Genuchten approach, adopted in Hydro- GeoSphere does not apply to fractured and karstified rock materials. The highly heterogeneous flow field and prefer- ential flow paths associated with such media and the conse- quently greater size of an REV compared to porous media are rendering the parameter determination by laboratory ex- periments impractical. Still, the van Genuchten parameters reflect properties of an unsaturated porous material and can be considered as an adequate parameter set to describe tran- sient infiltration processes if they are treated as calibration parameters in order to upscale from the Darcy-scale averag- ing volume to the field scale. 2.2 Sensitivity analysis An extensive sensitivity analysis is performed to determine the influence of the calibrated parameters on the computed flow. The root-mean-square error (RMSE) is chosen to rate the accuracy of fit and calculate deviations from the cali- brated model. The RMSE is defined as (Bamberg et al., 2007) RMSE= √√√√1 n n∑ i=1 (mi − fi)2, (25) where n is the total number of data points, mi denotes the simulated spring discharge derived by the sensitivity ana- lyzes and fi is the calibrated model value. For sensitivity ana- lyzes the RMSE has been evaluated using daily data and doc- umented parameter ranges were employed if possible. How- ever, for some variables, in particular the van Genuchten pa- rameters for the unsaturated zone of hard rocks, only few data and estimates can be found in the literature, e.g. Weiss (1987), Sauter (1992), Contractor and Jenson (2000) and Roulier et al. (2006). Consequently, the sensitivity of these parameters was determined systematically to evaluate the de- gree of ambiguity of the model. Depending on field obser- vations parameter ranges were varied on linear or log-scale. Parameter spaces were assigned to cover at least the reported values from field experiments. In order to provide a further quantitative comparison recession coefficients are given for one important recharge event during March 1988 (α1) and the following recession until beginning of April 1988 (α2). Due to the complex flow model it is likely that some pa- rameters do not show a linear correlation and sometimes the simulated discharge curve is only influenced by specific pa- rameter combinations. The final analysis of parameter sensi- tivity on an idealized example is subdivided into (1) insen- sitive parameters, (2) one sensitive parameter and (3) both parameters are sensitive (see Fig. 1, from left to right). In the latter case parameter A may be more sensitive for a certain www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3912 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems range of parameter B. Given a constant parameter Bn and Bn+dn where n denotes the parameter value and dn a change in the parameter, the parameter combination is refered to as non-linear if the change in RMSE over the whole range of parameter A is different for Bn and Bn+dn, i.e. if f (A,Bn) 6= f (A,Bn+dn). (26) This case is particularly important if a properly calibrated model can be achieved for two different parameter combi- nations where at least one parameter is a pure calibration parameter, i.e. its range is difficult to estimate by field ob- servations. Therefore, more than 1000 model runs were per- formed and the influence of parameter combinations on the simulated discharge curve analyzed. 3 Case study 3.1 Description of the field site The HydroGeoSphere model is employed to simulate flow in the catchment area of the Gallusquelle spring from Febru- ary 1988 to January 1990. The Gallusquelle spring is situ- ated in Southwest Germany on the Swabian Alb, a northeast- southwest striking Jurassic carbonate plateau. The catchment area has been studied extensively by several authors includ- ing aquifer characterization (Birk et al., 2005; Sauter, 1992), speleology (Abel et al., 2002; Gwinner, 1976), and flow pro- cesses (Geyer et al., 2008). The size of the catchment is about 45 km2. It is delineated by a water divide in the northwest and the River Fehla in the northeast (see Fig. 2). In the south the catchment is bounded by the northeastern fault zone of the Hohenzollerngraben, which was found to be imperme- able by tunneling works in the 1960s. After Sauter (1992) the base of the aquifer is formed by Kimmeridge marls (ki1). With a 70 m to 90 m thickness, this unit is rather consis- tent in the northwestern parts of the project area and con- sists of calcareous marls and occasional limestone interca- lations. In the southeastern area, no borehole information about the lower boundary of the unit are available. The up- permost catchment is made up of a sequence of Kimmerid- gian limestones (ki2-5) from which the majority is developed as algal-sponge facies and, therefore, belongs to the Unterer Massenkalk or Oberer Massenkalk. The largest part of the catchment comprises limestones belonging to the “Unterer Massenkalk” (ki2 and ki3). The whole Jurassic sequence dips southeast at an angle of 1.2 degrees. 3.2 Geometry of the flow model Based on the geological model, a vertical two-dimensional model domain of the catchment area was set up. The length of the domain is 10 km with a vertical thickness of 225 m. It reflects a cross section parallel to the direction of flow in the Gallusquelle catchment (see Fig. 2, lower figure). The model domain is represented by two continua reflecting flow in the low permeability fissured matrix (matrix continuum) and the highly permeable conduits (conduit continuum). The top of the model domain is set to 775 ma.s.l., which is an average elevation between ca. 910 ma.s.l. in the northwestern part of the catchment and ca. 640 ma.s.l. in the southeastern catch- ment and higher than the maximum groundwater head in the catchment. Every continuum is spatially discretized into 50 columns with a length of 200 m and a width of 1 m and 44 layers with a thickness of 5 m. 3.3 Boundary conditions The lateral sides of the matrix continuum and of the con- duit continuum, as well as the top of the conduit continuum are defined as no flow boundaries (see Fig. 3). A constant head boundary is applied to the right side of the conduit con- tinuum at 634 m a.s.l. to represent the spring and allow dis- charge. A specified flux boundary is set at the top of the ma- trix continuum to account for diffuse recharge. Daily data of total recharge was estimated by Geyer (2008) for the sim- ulation period on the Gallusquelle catchment. The applied water balance approach accounts for canopy storage, snow storage and soil-moisture storage before water entering the bedrock. A second specified flux boundary is set on the bot- tom of the whole conduit continuum to add rapid recharge in the aquifer. The location of the boundary condition con- siders that the transit time of the rapid recharge component through the unsaturated zone is less than one day (Geyer et al., 2008) and, therefore, negligible with regard to the daily time steps. The simulation of rapid water percolation from the top of the conduit continuum to the groundwater surface is physically not possible with the van Genuchten approach, because it does not consider gravity driven flow processes like film and droplet flow. The initial head distribution for transient discharge simulations is computed with a steady state simulation. The applied total recharge for the simulation is 1.5 mm d−1, which corresponds to the average recharge across the catchment area during the year 1988. Ten percent of the total recharge is employed as rapid recharge compo- nent at the bottom of the conduit continuum. The amount is in the range of the rapid recharge component estimated by Sauter (1997) from event analysis using oxygen isotopes in precipitation and Gallusquelle spring water to differentiate between different flow components. 3.4 Parameterization For the model calibration, known parameters are only var- ied within reasonable ranges that agree with actual field observations (Table 1). Unknown model parameters are in- vestigated by an extensive sensitivity analysis. The spe- cific storage coefficients for matrix and conduits are negli- gible since the aquifer is unconfined; hence, water released due to compaction in the saturated zone is irrelevant. As Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3913 R M SE A B A B R M SE A B R M SE A1 A2 A1 A2 A1 A2 ∂RMSE(A1) = ∂RMSE(A2) ≠ 0 ∂B ∂B ∂RMSE(A1) = ∂RMSE(A2) = 0 ∂B ∂B ∂RMSE(A1) ≠ ∂RMSE(A2) ≠ 0 ∂B ∂B Fig. 1. Examples for inter-parameter dependencies. From left to right: (1) both parameters insensitive, (2) both parameters sensitive with linear dependency, (3) both parameters sensitive but non-linear dependency. Table 1. Estimated values for the flow model derived from the model calibration and value ranges reported in the literature. Subscript m and c denote the matrix resp. the conduit continuum. Parameter Unit Value Literature Values Reference Kc (m s−1) 10.0 3.0–10.0 Sauter (1992) Km (m s−1) 2.9× 10−6 1.0× 10−6–1.0× 10−4 Sauter (1992) θsc a (–) 0.00023 0.00016–0.00064 Sauter (1992) θsm (–) 0.042 0.007–0.025, >0.0–0.12 Sauter (1992); Weiss (1987) αex (m−2) 1.0 – – Q (%) 90 – – αm (m−1) 0.0365 0.0365, 0.0328–0.623 Roulier et al. (2006); Contractor and Jenson (2000) nm (–) 1.83 1.83, 0.01–3.0 Roulier et al. (2006); Contractor and Jenson (2000) Swrm (–) 0.05 0.01–0.05 Contractor and Jenson (2000) αc (–) 5.1 5.1 Roulier et al. (2006) nc (–) 2.56 2.56 Roulier et al. (2006) Swrc (–) 0.0 – – krminc (–) 0.05 – – a The local conduit continuum porosity is 1.0 i.e. wcθsc(local) = θtotal −wmθsm implicitly gives the total conduit porosity such that θsc ∧= wcθsc(local). The total porosity is θtotal = 0.0422. there are no documented values for the hydraulic proper- ties of the interface available, the van Genuchten param- eters αi , ni , Swri and the interface hydraulic conductiv- ity Ki were set to values equal to the surrounding fis- sured matrix. Accordingly, inter-continuum water exchange is solely controlled by adjusting the exchange parameter αex. Model calibration is accomplished by fitting the observed and simulated discharge curves. Finally, the flow model con- tains 21 adjustable parameters for the fissured matrix and the conduit continuum. 4 Results and discussion 4.1 Model calibration The calibrated model shows a good fit with most of the spe- cific characteristics of the discharge hydrograph during the period between 16 February 1988 and 20 January 1990 (see Fig. 4). Please note that the discharge has been normalized to the catchment area (45 km2). Calibrated values for all var- ied parameters are comparable to values documented in the literature (Table 1). The observed discharge curve shows less sharp peaks and is smoother than the simulated curve. Sauter (1992) did get a comparable fit with a double continuum model for the saturated zone. The author did apply a function www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3914 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems Hohenzollerngraben River Lauchert River Fehla 67 0 66 0 650 64 0 68 06 90 720 73 0 740 70 0 71 0 Gallusspring Germany Faults Catchment Simulation profile 1000 m 400 500 600 700 800 900 0 1 2 3 4 5 6 7 8 9 10 ox1 marls ox2 bedded limestone ki1 marly limestone ki2/3 massive limestone Aquifer Highly transmissive zone m a.s.l distance to spring (km) spring Fig. 2. Gallusquelle catchment and crosssection used for the simu- lation. for the transfer of water from the soil zone to the groundwater surface, which is not necessary for our model. During the time period investigated, two strong discharge events occurred caused by major snowmelts, which are re- ferred to as first and second peak (see Fig. 4). The discharged water volume agrees well with the simulated data during the time period of the first peak. During the second discharge event (second peak) the simulated peak height is overesti- mated. It is not possible to change the relative peak height difference between the first and second peak with the avail- able calibration parameters. The recession curve slopes af- ter discharge events show a good fit, except during low flow conditions between July and October 1989. This behavior could be attributed to the simplified geometry of the numer- ical model, which does not include the documented slightly inclined aquifer base and the geometry of different karstified zones in the karst system. The hydraulic heads in the matrix continuum and the conduit continuum are nearly identical during the simulation period with a difference of a few cen- timeters. Above the water table the matrix saturation drops to 0.35 near the surface (see Fig. 5). Flow paths in the unsatu- rated matrix continuum and conduit continuum are slightly inclined towards the spring, whereas flow in the saturated zone is laterally oriented towards the outlet, i.e. the karst spring. The flow paths of the unsaturated matrix continuum, which would be expected to be vertical for such a large scale porous medium, are caused by the strong influence of the conduit continuum, which imposes a strong hydraulic gra- dient all over the matrix continuum. This behavior cannot be prevented unless the secondary continuum would be re- stricted to cover only the saturated zone. However, this is not an adequate solution considering the transient behavior of the system, i.e. the variation of water levels within both continua. The saturation in the conduit continuum is close to zero and has a very sharp transition along the water ta- ble. In this model, unsaturated flow in the conduits is also calibrated by the krminc parameter (minimum relative per- meability of the conduits). Without this parameter, the rel- ative permeability of the conduit continuum is a function of the residual saturation, i.e. setting of krminc simply overrides Eq. (13). This is the case for most of the unsaturated con- duit continuum, where saturation declines very quickly (be- low 0.05) above the water table for the given van Genuchten parameters. Therefore, with the applied van Genuchten pa- rameters only, water flow in the unsaturated conduit contin- uum is extremely small such that exchange from the matrix into the conduit system is nearly completely prevented and a proper model calibration is impossible due to numerical in- sufficiencies. However, Tokunaga and Wan (1997) showed that gravity driven film flow processes occur on unsatu- rated fracture walls, which contribute to water percolation along surfaces and may act as an interface from the conduit system to the matrix system, thus giving a physical mean- ing to the krminc parameter. As the original van Genuchten model relies on a unimodal distribution of the pore space the hydraulic response of such flow processes cannot be ex- pected to be resolved by the model. Attempts to refine the original van Genuchten approach and include hydraulic fea- tures of fractures into a continuum model have been made for example by Ross and Smettem (1993), Durner (1994) and Brouye`re (2006) by constructing a continuous bi-modal retention curve. An important role for the water exchange in the double continuum approach is the exchange parameter αex. It deter- mines the ability of water to move in and out of the conduit continuum and lumps geometrical and hydraulic properties of the karst matrix system. The surface-volume ratio, for ex- ample, is higher for a dendritic system than for a single con- duit with the same conduit volume. The exchange parameter in the calibrated model is set to a high value such that it does not act as an additional barrier for water transfer between both continua and water transfer is mainly controlled by the hydraulic properties of the two continua. Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3915 550 600 650 700 750 -10000 -8000 -6000 -4000 -2000 0 Y X Z constant head (Dirichlet) 640 m asl (Ψ=0) specified flux (Neumann) top boundary of matrix distance from spring (m) he ig ht (m ) impermeable boundaries Element Size: 1m200 m 5m dual-node approach: 2 nodes at same physical location fractured matrix (primary continuum)conduits (secondary continuum) specified flux (Neumann) bottom boundary of bonduits Fig. 3. Model grid of the two-dimensional model. A constant head boundary is applied at 640 m in the second continuum. Recharge is applied along the top boundary in both continua. Every node in the primary continuum has its counterpart in the secondary continuum at the same physical location. F M A M J J A S O N D J F M A M J J A S O N D J 0 5 10 15 month (Feb./11/1988 − Jan./15/1990) sp rin g di sc ha rg e (m m/ d) 0 10 20 30 40 50 60 70 80 re ch ar ge (m m/ d) RMSE = 0.59 mm/d simulated observed recharge 1. peak 2. peak Fig. 4. The calibrated curve shows an acceptable fit. The secondary peak height is overestimated by the model. The first recharge event (first peak) reaches the spring too early. Recession curve slopes show a good overall correlation. 5 Sensitivity analysis 5.1 Single variation of hydraulic parameters for saturated flow conditions Table 2 gives an overview of the recession coefficients and RMSE values obtained for the sensitive parameters. A pa- rameter has been discarded as insensitive if the maximum RMSE is below 0.05 mm d−1. The recession coefficients have been measured at the first strong recharge event be- ginning in March 1988 (α1) and during the low flow reces- sion beginning in April 1988 (α2). The calibrated values are α1 = 0.23 and α2 = 0.03, which is close to what has been reported by Sauter (1989) for a conduit dominated reces- sion (α = 0.25) for the same recharge event. Figure 6 (up- per two graphs) shows the computed spring discharge for several model runs with varying hydraulic conductivity Kc and porosity θc in the conduit continuum. These parame- ters strongly influence the simulated spring discharge. Fig- ure 7 (upper two graphs) additionally shows the respective recession coefficients. An increased conduit conductivity Kc results in higher α1 recession coefficients and lower base flow levels indicated by the strong decrease of α2. A de- creased conduit conductivity Kc favors a slow recession and decreases α1 to 0.06, which according to Sauter (1989) al- ready indicates a mixed system (fractured matrix + conduits) response. Discharge peaks are broadened and the base flow is higher. In case the conduit drains the matrix system, an in- crease of Kc enhances the exchange process between matrix and conduits by decreasing the hydraulic gradient in the con- duit continuum and consequently increasing the hydraulic www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3916 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems sp rin g 0.0 0.2 0.4 0.6 0.8 1.0 saturation matrix max. water table (Ψ = 0) after 45 days (Mar./26th/1988) min. water table (Ψ = 0) after 656 days (Nov/27th/1988) he ig ht a sl (m ) 0.0 0.2 0.4 0.6 0.8 1.0 saturation conduit 550 600 650 700 750 550 600 650 700 750 he ig ht a sl (m ) distance from spring (km) 10.0 8.0 6.0 4.0 2.0 0.0 conduitflow velocityx-direction 1E-8 1E-7 1E-6 1E-5 1E-4 1E-3 1E-2 m /d distance from spring (km) 10.0 8.0 6.0 4.0 2.0 0.0 matrix 0.1 1 10 100 1000 flow velocity x-direction m /d Fig. 5. Results of the transient flow model (26 March 1988). Shown water tables apply for both continua as the height is nearly equal and differs 1 cm at most. The van Genuchten parameter αm leads to a strong difference between the continua. The capillary rise is more pronounced in the matrix system, where saturation is ca. 0.35 near the surface. Saturation in the conduit system is lower than 0.1 above the water table and below 0.0001 near the surface. Flow velocities (only x-direction vector, note the different scaling for matrix and conduits) are apparently higher in the conduit continuum. 0 5 10 15 100 10 (calibr.) 1 K c (m s−1) 0 5 10 15 0.00029 0.00023 (calibr.) 0.00014 θ sc (−) 0 5 10 15 0.102 0.042 (calibr.) 0.012 θ sm (−) F M A M J J A S O N D J F M A M J J A S O N D J 0 5 10 15 month (Feb./11/1988 − Jan./15/1990) sp rin g di sc ha rg e (m m/ d) 1.0 (calibr.) 0.001 0.0001 α ex (m−2) Fig. 6. Variation of hydraulic parameters (Kc, θsc, θsm) and the exchange parameter (αex). Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3917 Table 2. RMSE values and recession coeffcients for all sensitive pa- rameters. Parameters with a maximum RMSE below 0.05 mm d−1 have been considered as insensitive. Subscript “m” and “c” denote the matrix resp. the conduit continuum. Bold numbers denote the calibrated value. Parameter Value Recess. coeff. (1/d) RMSE (mm d−1) α1 α2 Kc 100 0.26 0.160 1.329 10 0.23 0.030 – 1 0.06 0.015 0.645 θsc 0.00029 0.29 0.033 0.188 0.00023 0.23 0.030 - 0.00014 0.23 0.026 0.108 θsm 0.102 0.33 0.066 0.318 0.042 0.23 0.030 – 0.012 0.19 0.023 0.317 αex 1.0 0.23 0.033 – 0.001 0.06 0.030 0.231 0.0001 0.36 0.002 0.555 αm 0.365 0.28 0.025 0.329 0.0365 0.23 0.030 – 0.00365 0.36 0.122 1.285 nm 2.23 0.31 0.029 0.064 1.83 0.23 0.300 – 1.23 0.26 0.049 0.563 Swrm 0.6 0.29 0.050 0.54 0.3 0.32 0.034 0.142 0.05 0.23 0.030 – krminc 0.8 0.32 0.040 0.334 0.05 0.23 0.030 – 0.01 0.16 0.031 0.068 gradient between matrix and conduits. The conduit poros- ity θsc follows a similar pattern, i.e. an increase will enhance the exchange process, however, the impact on the discharge curves is far less pronounced within the given ranges and both recession coefficients are all in the same order of mag- nitude indicating a conduit dominated recession. A contrast- ing behavior is observed by a variation in matrix porosity θsm. With an increase in the parameter the water transfer be- tween the continua during recharge events is decreased be- cause of the lower head difference between conduit and ma- trix and the discharge curve is smoothened accordingly (see Fig. 6). The recession coefficient α1 and α2 are consequently slightly lower while for a very low matrix porosity of 1.2% it is apparent that recessions coefficients represent a strongly conduit dominated system α1 = 0.33 and α2 = 0.066. Km displays a low sensitivity within the given parameter range which can be attributed to the high exchange parameter of the calibrated model of αex = 1.0. The high value leads to an immediate equalization of heads between conduit and ma- trix such that water will not be restrained within the matrix system when total heads are slightly higher than within the conduit system. The matrix system always depends on the hydraulic state of the conduit continuum, which discharges water rapidly to the spring. However, in a three-dimensional karst system, flow velocities within the matrix will be little influenced by the conduit system with increasing distance to the conduit. The exchange parameter αex is sensitive only for strong reductions on the order of three to four magni- tudes. A reduction to 0.001 lowers the peak height of both main peaks while decreasing recession curve slopes α1 to 0.06 and slightly increasing base levels. Further reduction to 0.0001 drastically decreases peak heights and increases the base levels. The resulting α2 coeffcients are very low (0.002) and recharge events show no more pronounced peaks. An ex- change reduction to 0.1 or 0.01 has no significant influence on the discharge curves, which indicates a sensitive interval between 0.001 and 0.0001. 5.2 Single variation of unsaturated zone parameters Variations of the sensitive van Genuchten parameters for the vadose zone are shown in Fig. 8 and the corresponding re- cession coefficients in Fig. 9. The decrease of αm and nm results in a strong rise of peak heights and increase of reces- sion slopes (α1 = 0.36 and 0.26 respectively). The influence of the van Genuchten αm parameter on the discharge curve is connected to the inter-continuum water exchange pro- cess. Lowering the parameter increases the capillary rise, i.e. the matrix has higher saturation (and relative permeability) above the water table. Consequently the increased permeabil- ity leads to a stronger and earlier exchange of water from the matrix into the conduit continuum, such that recharge events affect spring discharge a lot earlier (pronounced event peaks). The opposite can be observed for a value of 0.365 where the saturation fringe declines very quickly with lower saturations above the water table. This reduces the main exchange inter- face to a smaller area above the water table. Thus during high recharge events, peak heights are reduced since water will remain longer in the matrix continuum and the α2 recession coefficient becomes slightly lower (0.025) reflecting the de- layed discharge via the conduit system. The van Genuchten parameter nm can be considered insensitive compared to αm. The conduit van Genuchten parameters αc and nc are as well insensitive for the shown simulations. In the range of chosen values, the conduits do not produce a strong capillary rise, i.e. the unsaturated zone above the conduit water table al- ways displays a sharp transition from saturated to strongly unsaturated. As mentioned earlier the application of the van Genuchten parameters to a highly conductive and discrete flow system such as a conduit implies a general abstraction of the physically based van Genuchten parameter set in order to create an upscaled continuum system with a characteris- tic infiltration behavior and travel time distribution as well as www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3918 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 2 4 6 8 10 12 14 α1 =0.26 α1 =0.23 α1 =0.06 α2 =0.160α2 =0.030α2 =0.015 100 10 (calibr.) 1 K c (m s−1) 2 4 6 8 10 12 14 α1 =0.29α1 =0.23α1 =0.23 α2 =0.033α2 =0.030α2 =0.026 0.00029 0.00023 (calibr.) 0.00014 θ sc (−) 2 4 6 8 10 12 14 α1 =0.33 α1 =0.23 α1 =0.19 α2 =0.066α2 =0.030α2 =0.023 0.102 0.042 (calibr.) 0.012 θ sm (−) 0 5 10 15 20 25 30 35 40 45 50 2 4 6 8 10 12 14 α1 =0.23 α1 =0.06α1 =0.36 α2 =0.033α2 =0.030α2 =0.002 days (Feb./21/1988 − Apr./11/1988) sp rin g di sc ha rg e (m m/ d) 1.0 (calibr.) 0.001 0.0001 α ex (m−2) Fig. 7. Variation of hydraulic parameters (Kc, θsc, θsm) and the exchange parameter (αex). an exchange interface in the unsaturated zone. In this work the exchange process in the unsaturated zone can be con- trolled by the αc parameter in order to increase the capillary rise in the conduit continuum and enhance inter-continuum water transfer. However, such an approach also introduces a spatial information (i.e. the thickness of the conduit capil- lary fringe in vertical direction) which is not known in real karst systems. As described before the krminc parameter is used instead to maintain a constant water exchange in the unsaturated zone independent of the hydraulic state of the conduit system if saturations are too low. The residual wa- ter saturation of the matrix Swrm and the minimum relative permeability of the conduits krminc both show a similar be- havior regarding parameter variations. Increasing the param- eters yields an enhanced exchange from the matrix to the conduit continuum due to a higher relative conductivity. Con- sequently, recharge events are transmitted faster to the model outlet, i.e. the spring. 5.3 Combined parameter variations The above presented sensitivity analyzes imply only one sin- gle parameter varied at a time. However, a further important observation is that certain parameter combinations may show non-linear behavior with respect to their sensitivity, i.e. the influence of one parameter on the RMSE is not linear over the whole range of a second parameter. Table 3 shows maxi- mum RMSE obtained for each parameter combination and if a non-linear relationship can be observed (bold RMSE val- ues). For example, the simultaneous variation of the matrix van Genuchten parameter αm and the conduit conductivity Kc displays a pronounced sensitivity for low αm values (Fig. 10). While for the calibrated αm value of 0.0365 m−1, the conduit conductivity Kc is almost insensitive in the range of 1–10 m s−1; a lower αm value of 0.00365 m−1 yields a high RMSE of 1.6 mm d−1 already at a Kc value of 10 m s−1, i.e. ∂RMSE(αm = 0.00365)/∂Kc is much higher. A similar be- havior can be shown for a combination of matrix porosity θsm and the conduit conductivity where lower porosities yield higher RMSE values with an increase in conduit conductiv- ity to 100 m s−1, whereas for rather high matrix porosities of 0.102 the increase in RMSE is less pronounced. The con- duit conductivityKc exhibits a higher sensitivity for the cali- brated exchange parameter αex = 1.0 (see Fig. 10) such that a high conductivity value (100 m s−1) results in RMSE values Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3919 0 5 10 15 0.365 0.0365 (calibr.) 0.00365 α m (m−1) 0 5 10 15 2.23 1.83 (calibr.) 1.23 n m (−) 0 5 10 15 0.6 0.3 0.05 (calibr.) S wrm (−) F M A M J J A S O N D J F M A M J J A S O N D J 0 5 10 15 month (Feb./11/1988 − Jan./15/1990) sp rin g di sc ha rg e (m m/ d) 0.8 0.05 (calibr.) 0.01 k rminc (−) Fig. 8. Variation of the van Genuchten parameters (αm, nm, Swrm, krminc). Table 3. Parameter combinations used for the sensitivity analyzes and highest RMSE (mm d−1) observed. Bold RMSE values denote a non-linear inter-parameter dependence. Range Low High Kc Km θsc θsm αex Q αm nm Swrm αc nc Swrc krminc Kc 1.0 100.0 – Km 2.89× 10−7 2.89× 10−4 1.34 – θsc 1.40× 10−4 2.90× 10−4 1.48 0.19 – θsm 0.022 0.102 2.23 0.83 0.89 – αex 1.00× 10−4 1.0 1.33 0.56 0.61 0.77 – Q 30 95 1.33 0.06 0.23 0.77 0.57 – αm 0.365× 10−3 3.65 3.20 · · · · · · · · · 1.29 · · · – nm 1.23 2.23 · · · · · · · · · · · · 0.58 · · · 1.58 - Swrm 0.05 0.6 · · · · · · · · · · · · · · · · · · 1.87 · · · – αc 3.1 8.1 · · · · · · · · · · · · 0.56 · · · · · · · · · · · · – nc 2.16 3.16 · · · · · · · · · · · · 0.56 · · · · · · 0.56 · · · 0.25 – Swrc 0.0 0.4 · · · · · · · · · · · · 0.56 · · · · · · · · · · · · 0.01 · · · – krminc 0.01 0.8 · · · · · · · · · · · · 0.56 · · · · · · · · · · · · 0.33 · · · · · · – www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3920 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 2 4 6 8 10 12 14 α1 =0.36α1 =0.28 α1 =0.23 α2 =0.122α2 =0.030α2 =0.025 0.365 0.0365 (calibr.) 0.00365 α m (m−1) 2 4 6 8 10 12 14 α1 =0.26 α1 =0.23 α1 =0.31 α2 =0.049α2 =0.029α2 =0.030 2.23 1.83 (calibr.) 1.23 n m (−) 2 4 6 8 10 12 14 α1 =0.29 α1 =0.23 α1 =0.32 α2 =0.050α2 =0.034α2 =0.030 0.6 0.3 0.05 (calibr.) S wrm (−) 0 5 10 15 20 25 30 35 40 45 50 2 4 6 8 10 12 14 α1 =0.32 α1 =0.23 α1 =0.16 α2 =0.040α2 =0.030α2 =0.031 days (Feb./21/1988 − Apr./11/1988) sp rin g di sc ha rg e (m m/ d) 0.8 0.05 (calibr.) 0.01 k rminc (−) Fig. 9. Variation of the van Genuchten parameters (αm, nm, Swrm, krminc). of ca. 1.4 mm d−1 while for a lower exchange parameter the same conductivity yields a deviation of only 0.4 mm d−1. The exchange parameter has a higher sensitivity for high conduc- tivity values of the conduit system, while it is nearly insen- sitive for low values (1 m s−1). In sum, the variation of the exchange parameter influences the discharge curve depend- ing on the combination with other parameters. This behavior can also be observed for the combination of matrix poros- ity θsm and exchange parameter αex. Here the exchange pa- rameter has a higher sensitivity for matrix porosities between 0.032 and 0.102, while at the lower limit (0.012–0.022) this sensitivity vanishes. 6 Conclusions The applied vertical two-dimensional double continuum ap- proach is successfully employed to simulate the discharge behavior of the karst system Gallusquelle for a period of two years. Because of the high amount of model parame- ters of the saturated-unsaturated flow model, a comprehen- sive sensitivity analysis was performed. The analysis shows that the simulated discharge curve displays high sensitivity to a variation of a number of model parameters. The sensitiv- ity study demonstrates that the simulation of karst hydraulics requires strong a-priori knowledge about parameter ranges to reduce ambiguity of the model. However, especially for the conduit continuum, only little about physical parameters is known and further research is needed. Physical parameters of the saturated fissured matrix can be obtained from field experiments on the scale of a representative elementary vol- ume. Furthermore, the analysis shows that the sensitivity of a parameter depends to a large degree on the other calibrated model parameters. Therefore, sensitivity analyzes should si- multaneously take into account parameters of both continua in order to detect deviations from a linear behavior if both pa- rameters are sensitive. It also means that conclusions about parameter sensitivity change from model to model and are not simply transferable. The fissured matrix porosity as well as van Genuchten parameters of the matrix continuum are the most important parameters for an appropriate flow simula- tion. The conduit system drains the fissured matrix and can, due to its high hydraulic conductivity, effectively discharge varying quantities of water transferred from the matrix con- tinuum. The van Genuchten parameters of the fissured matrix Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3921 1 10 100 0.00365 0.0365 0.365 3.65 0 0.5 1 1.5 2 2.5 3 3.5 K c (m/s)αm (m −1) R M SE (m m/ d) 1 10 100 0.022 0.032 0.042 0.052 0.062 0.082 0.102 0 0.5 1 1.5 2 2.5 3 3.5 K c (m/s)θsm (−) R M SE (m m/ d) 1 10 100 0.0001 0.001 0.01 0.1 1 0 0.5 1 1.5 2 2.5 3 3.5 K c (m/s)αex (m −2) R M SE (m m/ d) 0.0001 0.001 0.01 0.1 1 0.022 0.032 0.042 0.052 0.062 0.082 0.102 0 0.5 1 1.5 2 2.5 3 3.5 α ex (m−2)θsm(−) R M SE (m m/ d) Fig. 10. Sensitivity matrices showing nonlinear inter-parameter dependencies. are the most crucial property in terms of sensitivity, uncer- tainty and model limitations. The exchange process between matrix and conduit continuum is mainly controlled by dif- ferences in hydraulic properties. The αex parameter was set to a rather high value during the calibration, i.e. exchange is not limited by a too low exchange coefficient. According to Gerke and Van Genuchten (1993), the parameter is de- fined to express the interface connectivity on a rather small scale, e.g. between a porous medium and macropores. On catchment scale it might implicitly correspond to the type of karst system (i.e. karstified layers vs. karst conduits). The choice of a two-dimensional dual-continuum model for sim- ulation of flow in karst systems leads to certain advantages and limitations. Considering the lack of information for most karst spring catchments the dual-continuum approach pro- vides an efficient way to lump its geometrical features, i.e. only limited informations about the conduit geometry and rapid recharge area are necessary. However, lateral flow dy- namics between the localized conduits and the matrix sys- tem are not considered and water table fluctuations in a continuum only represent the volume-averaged dynamics of the whole system. Vertical two-dimensional modeling ap- proaches are therefore more appropriate for the simulation of karstified layers than for the simulation of local karst con- duit systems embedded in a fissured matrix. Furthermore the vertical two-dimensional approach neglects the lateral geom- etry of a spring catchment, e.g. reduction of cross-sectional area of the matrix contiunuum towards a karst spring, which might affect the discharge behavior. The unified simulation of the saturated and unsaturated matrix and conduit continuum using the Richards equation is limited valid regarding the underlying physical process de- scription. The application of the Richards equation for the matrix continuum can be considered as suitable at large scale (e.g. Kaufmann, 2003). But at this scale, van Genuchten pa- rameters cannot be measured and become calibration param- eters. A clear limitation is the loss of the physical meaning of the Richards equation for the conduit continuum. Flow through localized karst conduits might be highly transient, i.e. the application of open channel flow approaches would be appropriate (Reimann et al., 2011a). However, these ap- proaches increase the number of calibration parameters and require additional knowledge about the geometry of a karst conduit system. Additionally, flow in unsaturated karst fea- tures (e.g. vertical shafts) is not controlled by matrix poten- tial and capillary forces but rather flow processes dominated www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 3922 J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems by gravitational forces such as film flow (Tokunaga and Wan, 1997; Tokunaga et al., 2000), turbulent film flow (Ghezze- hei, 2004), droplet flow (Doe, 2001; Dragila and Weisbrod, 2004) and rivulet flow (Su et al., 2001; Dragila et al., 2004; Su et al., 2004). In order to be able to apply the vertical dual continuum approach, boundary conditions were modi- fied and rapid recharge was directly injected at the bottom of the saturated conduit continuum. This abstraction clearly reduces the predictive capability of the applied approach. Acknowledgements. This work was funded by the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) under grants no. GE 2173/2-2 and SA 501/24-1. Careful reviews of the manuscript by the editor and two anonymous reviewers are gratefully acknowledged. We would also like to thank Thomas Graf and Rob McLaren for their advice. Edited by: I. Neuweiler References Abel, T., Sauter, M., and Hinderer, M.: Karst genesis of the Swabian Alb, south Germany, since the Pliocene, Acta Geol. Pol., 52, 43– 54, 2002. Bamberg, G., Baur, F., and Krapp, M.: Statistik, Oldenbourg, Mu¨nchen, 2007. Barenblatt, G. I., Zheltov, I. P., and Kochina, I. N.: Basic Concepts in the Theory of Seepage of Homogenous Liquids in Fissured Rocks Strata, J. Appl. Math. Mech., 24, 1286–1303, 1960. 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. Brouye`re, S.: Modelling the migration of contaminants through variably saturated dual-porosity, dual-permeability chalk, J. Con- tam. Hydrol., 82, 195–219, 2006. Contractor, D. N. and Jenson, J. W.: Simulated effect of vadose in- filtration on water levels in the Northern Guam Lens Aquifer, J. Hydrol., 229, 232–254, 2000. Doe, T.: What do drops do? Surface wetting and network geometry effects on vadose-zone fracture flow, in: Conceptual models of flow and transport in the fractured vadose zone, Chap. 8, 243– 270, National Academy Press, Washington D.C., 2001. Dragila, M. I. and Weisbrod, N.: Flow in Menisci Corners of Capil- lary Rivulets, Vadose Zone J., 3, 1439–1442, 2004. Dragila, M. I., Weisbrod, N., and Council, N. R.: Fluid motion through an unsaturated fracture junction, Water Resour. Res., 40, 1–11, 2004. Dreiss, S. J.: Regional scale transport in a karst aquifer: 1. Compo- nent separation of spring flow hydrographs, Water Resour. Res., 25, 117–125, 1989. Dreybrodt, W., Gabrovsˇek, F., and Romanov, D.: Processes of Speleogenessis: A Modeling Approach, Karst Research Institute at ZRC SAZU, ZRC Publishing, Ljubljana, 2005. Durner, W.: Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res., 30, 211, doi:10.1029/93WR02676, 1994. Gerke, H. and Van Genuchten, M.: Evaluation of a First-Order Water Transfer Term for Variably Saturated Dual-Porosity Flow Models, Water Resour. Res., 29, 1225–1238, 1993. Germann, P. F., Helbling, A., and Vadilonga, T.: Rivulet Approach to Rates of Preferential Infiltration, Vadose Zone J., 6, 207–220, 2007. Geyer, T.: Process-based characterisation of flow and transport in karst aquifers at catchment scale, Ph.D. thesis, University of Go¨ttingen, 2008. 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. Ghezzehei, T. A.: Constraints for flow regimes on smooth fracture surfaces, Water Resour. Res., 40, 1–14, 2004. Gwinner, M. P.: Origin of the Upper Jurassic limestones of the Swabian Alb (Southwest Germany), E. Schweizerbart, 1976. Hill, C. and Polyak, V.: Karst hydrology of Grand Canyon, Arizona, USA, J. Hydrol., 390, 169–181, 2010. Jeannin, P. Y.: Modeling flow in phreatic and epiphreatic karst con- duits, Water Resour. Res., 37, 191–200, 2001. Kaufmann, G.: Modelling unsaturated ow in an evolving karst aquifer, J. Hydrol., 276, 53–70, doi:10.1016/S0022- 1694(03)00037-4, 2003. Mclaren, R. G., Forsyth, P. A., Sudicky, E. A., Vanderkwaak, J. E., Schwartz, F. W., and Kessler, J. H.: Flow and transport in frac- tured tuff at Yucca Mountain: numerical experiments on fast preferential flow mechanisms, J. Contam. Hydrol., 43, 211–238, 2000. Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, 1976. Reimann, T., Geyer, T., Shoemaker, W. B., Liedl, R., and Sauter, M.: Effects of dynamically variable saturation and matrix-conduit coupling of flow in karst aquifers, Water Resour. Res., 47, 1–19, 2011a. 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, 1–15, 2011b. Richards, L. A.: Capillary Conduction of Liquids Through Porous Mediums, Physics, 1, 318–333, 1931. Rosenbom, A., Therrien, R., and Refsgaard, J.: Numerical analy- sis of water and solute transport in variably-saturated fractured clayey till, J. Contam. Hydrol., 104, 137–152, 2009. Ross, P. J. and Smettem, K. R. J.: Describing soil hydraulic prop- erties with sums of simple functions, Soil Sci. Soc. Am. J., 57, 26–26, 1993. Roulier, S., Baran, N., Mouvet, C., Stenemo, F., Morvan, X., Al- brechtsen, H.-J. R., Clausen, L., and Jarvis, N.: Controls on atrazine leaching through a soil-unsaturated fractured limestone sequence at Bre´villes, France, J. Contam. Hydrol., 84, 81–105, 2006. Sauter, M.: Quantification and Forecasting of Regional Groundwa- ter Flow and Transport in a Karst Aquifer (Gallusquelle, Malm, SW. Germany), Tu¨binger Geowissenschaftliche Arbeiten, 1992. Sauter, M.: Differentiation of flow components in a karst aquifer us- ing the δ18O signature, in: Tracer Hydrology, edited by: Kranjc, A., 435–441, Balkema, 1997. Sauter, M., Geyer, T., Kovacs, A., and Teutsch, G.: Modellierung der Hydraulik von Karstgrundwasserleitern – Eine ¨Ubersicht, Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012 www.hydrol-earth-syst-sci.net/16/3909/2012/ J. Kordilla et al.: Double continuum simulation of saturated and unsaturated flow in karst systems 3923 Grundwasser, 3, 143–156, 2006. Schwartz, F., Sudicky, E., and McLaren, R.: Ambiguous hydraulic heads and C14 activities in transient regional flow, Ground Wa- ter, 48, 366–379, 2010. Smart, P. L. and Hobbs, S. L.: Characterisation of carbonate aquifers: A conceptual base, in: Proceedings of the Environmen- tal Problems in Karst Terranes and Their Solutions Conference, Bowling Green, Ky. Published by NWWA, 17–31, 1986. Su, G. W., Geller, J. T., Pruess, K., and Hunt, J. R.: Solute trans- port along preferential flow paths in unsaturated fractures, Water Resour. Res., 37, 2481–2491, 2001. Su, G. W., Geller, J. T., Hunt, J. R., and Pruess, K.: Small-Scale Features of Gravity-Driven Flow in Unsaturated Fractures, Va- dose Zone J., 3, 592–601, 2004. Teutsch, G.: Grundwassermodelle im Karst: Praktische Ansa¨tze am Beispiel zweier Einzugsgebiete im Tiefen und Seichten Malmkarst der Schwa¨bischen Alb, Ph.D. thesis, Universita¨t Tu¨bingen, 1988. Teutsch, G. and Sauter, M.: Groundwater modeling in karst terranes: Scale effects, data acquisition and field validation, in: Proc. Third Conf. Hydrogeology, Ecology, Monitoring, and Management of Ground Water in Karst Terranes, Nashville, TN, 17–35, 1991. Therrien, R. and Sudicky, E. A.: Three-dimensional analysis of variably-saturated flow and solute transport in discretely- fractured porous media, J. Contam. Hydrol., 3542, 1–44, 1996. Therrien, R., McLaren, R., Sudicky, E., and Panday, S.: HydroGeo- Sphere: A three-dimensional numerical model describing fully- integrated subsurface and surface flow and solute transport, Man- ual (Draft), HydroGeoLogic Inc., Herndon, VA, 2006. Tokunaga, T. K. and Wan, J.: Water film flow along frac- ture surfaces of porous rock, Water Resour. Res., 33, 1287, doi:10.1029/97WR00473, 1997. Tokunaga, T. K., Wan, J., and Sutton, S. R.: Transient film flow on rough fracture surfaces, Water Resour. Res., 36, 1737–1746, 2000. van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898, 1980. Weiss, E. G.: Porosita¨ten, Permeabilita¨ten und Verkarstungserschei- nungen im mittleren und oberen Malm der su¨dlichen Frankenalb, Ph.D. thesis, Universita¨t Erlangen, 1987. www.hydrol-earth-syst-sci.net/16/3909/2012/ Hydrol. Earth Syst. Sci., 16, 3909–3923, 2012