Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach

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 discharge with the Richards equations and van Genuchten parameters 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 domain. An extensive sensitivity analysis is performed on single parameters as well as parameter combinations. The transient hydraulic response of the karst spring is strongly controlled 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 conduits are draining the matrix. Sensitivities of parameter combinations 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 degree of ambiguity. The application of van Genuchten parameters for simulation of unsaturated flow in karst systems is critically discussed.

storage and a less permeable one with high storage.Therefore, different individual (rapid, slow) flow components with characteristic temporal distributions are induced.Accordingly, the final spring discharge is then a function of the individual flow contributions 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 simulation of coupled saturated and unsaturated flow is still a challenge 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 systems.Therefore, flow in karst systems is often simulated with lumped parameter modeling approaches, which translate precipitation signals to discharge hydrographs by applying simple transfer functions (Dreiss, 1989).Generally, this type of approach is appropriate for situations in which predicted system states are expected to range between already observed events.The simulation of natural karst systems with distributed parameter models is reported only in a few studies (e.g.Jeannin, 2001;Birk et al., 2005;Hill and Polyak, 2010).However, distributive modeling approaches incorporate flow laws and, therefore, are adequate for the process based simulation of karst hydraulics (e.g.Birk et al., 2006;Reimann et al., 2011b).Teutsch and Sauter (1991) demonstrate in 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 limited information about aquifer geometry and still allows 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 diffuse flow in the fissured matrix and the discrete rapid underground drainage by solution conduits in karst systems.Here, we want to assess the relative importance of individual factors and parameter combinations on the discharge behavior of a karst Figures spring without detailed knowledge about the hydraulic parameter field of an aquifer system.This type of information is of major importance to focus characterization efforts in catchment based karst studies.Furthermore, the importance of infiltration dynamics, i.e. the temporal distribution of the rapid and the slow flow component on the discharge dynamics is to be determined.Therefore, an integrated saturated-unsaturated doublecontinuum approach was selected for the study.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.

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 consecutively at the same node and coupled with an exchange term that defines the hydraulic interface and controls the inter-continuum exchange flow.The governing equation in the applied model HydroGeoSphere (Therrien et al., 2006) is the Richards equation (Richards, 1931), which is slightly modified to account for inter-continuum water exchange: and the total porosity is given as The fluxes q m and q c are obtained from where K m and K c denote hydraulic conductivity, ψ m and ψ c are the pressure heads in each continuum and z is the elevation head.
In the unsaturated zone, the relative permeabilities k rm , k rm and k ri (interface) depend on the water saturation which in turn is related to the pressure head according to van Genuchten (1980): ψ ≥ 0 the residual saturations are S wm = S wc = S wi = 1.The relative permeability is given by: with l p being the pore connectivity parameter (equals 0.5 after Mualem, 1976), S e the effective saturation and ν is defined as:

Conclusions References
Tables Figures

Back Close
Full for β > 1.In the saturated zone the storage terms on the right-hand side of Eqs. ( 1) and ( 2) are replaced by: where S sm and S sc are the specific storage coefficients.Water release by compaction of the porous medium is neglected in the unsaturated zone.The term Γ ex in Eq. ( 1) and Eq. ( 2) describes the volumetric fluid exchange rate per unit volume between primary and secondary continuum and is given as: where K i is the hydraulic conductivity of the interface (e.g.sediments) and k ri the relative interface permeability (Barenblatt et al., 1960).The exchange parameter α ex is determined by calibration and defined as (Gerke and Van Genuchten, 1993): where β is a geometry factor (3 for rectangular matrix blocks, 15 for spheres), a is the distance between the center of a matrix block and the adjacent fracture or conduit and γ w is an empirical coefficient usually set to 0.4.Strictly speaking, the van Genuchten approach, adopted in HydroGeoSphere does not apply to fractured and karstified rock materials.The highly heterogeneous flow field and preferential flow paths associated with such media and the consequently greater size of an REV compared to porous media are rendering the parameter determination by laboratory experiments impractical.Still, the van Genuchten parameters reflect properties of an unsaturated porous material and can be considered as an adequate parameter set to describe transient infiltration processes if they are treated as calibration parameters in order to upscale from the Darcy-scale averaging volume to the field scale.Introduction

Conclusions References
Tables Figures

Back Close
Full

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 calibrated model.The RMSE is defined as (Bamberg et al., 2007): where n is the total number of data points, m i denotes the spring discharge derived by the model and f i is the calibrated model value.For sensitivity analyses, documented parameter ranges were employed.However, for some variables, in particular the van Genuchten parameters 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;Roulier et al., 2006).Consequently, the sensitivity of these parameters was determined systematically to evaluate the degree of ambiguity of the model.Due to the complex flow model it is likely that some parameters do not show a linear correlation and sometimes the simulated discharge curve is only influenced by specific parameter combinations.The final analysis of parameter sensitivity on an idealized example is subdivided into: (1) insensitive 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 range of parameter B, i.e. δRMSE(A 1 ) This case is particularly important if a properly calibrated model can be achieved for two different parameter combinations where at least one parameter is a pure calibration parameter, i.e. its range is difficult to estimate by field observations.Therefore, more than 1000 model runs were performed and the influence of parameter combinations on the simulated discharge curve analyzed.Introduction

Conclusions References
Tables Figures

Back Close
Full 3 Case study

Description of the field site
The HydroGeoSphere model is employed to simulate flow in the catchment area of the Gallusquelle spring from February 1988 to January 1990.The Gallusquelle spring is situated in Southwest Germany on the Swabian Alb, a northeast-southwest striking Jurassic carbonate plateau.The catchment area has been studied extensively by several authors including aquifer characterization (Birk et al., 2005;Sauter, 1992), speleology (Abel et al., 2002;Gwinner, 1976), and flow processes (Geyer et al., 2008).
The size of the catchment is about 45 km 2 .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 impermeable by tunneling works in the 1960's.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 consistent in the northwestern parts of the project area and consists of calcareous marls and occasional limestone intercalations.In the southeastern area no borehole information about the lower boundary of the unit are available.The uppermost catchment is made up of a sequence of Kimmeridgian 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 • .

Geometry of the flow model
Based on the geological model, a vertical two-dimensional model domain of the catchment area was set up (Fig. 3).The length of the domain is 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 m a.s.l., which is an average elevation between ca.910 m a.s.l. in the north-western part of the catchment and ca.640 m a.s.l. in the south-eastern catchment 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 (Fig. 3).

Boundary conditions
The lateral sides of the matrix continuum and of the conduit continuum, as well as the top of the conduit continuum are defined as no flow boundaries (see Fig. 4).A constant head boundary is applied to the right side of the conduit continuum at 634 m a.s.l. to represent the spring and allow discharge.A specified flux boundary is set at the top of the matrix continuum to account for diffuse recharge.Daily data of total recharge was estimated by Geyer (2008) for the simulation 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 bottom of the whole conduit continuum to add rapid recharge in the aquifer.The location of the boundary condition considers that the transit time of the rapid recharge component through the unsaturated zone is below 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. estimated by Sauter (1997) from event analysis using oxygen isotopes in precipitation and Gallusquelle spring water to differentiate between different flow components.

Parameterization
For the model calibration, known parameters are only varied within reasonable ranges that agree with actual field observations (Table 1).Unknown model parameters are investigated by an extensive sensitivity analysis.The specific storage coefficients for matrix and conduits are negligible since the aquifer is unconfined; hence, water released due to compaction in the saturated zone is irrelevant.As there are no documented values for the hydraulic properties of the interface available, the van Genuchten parameters α i , β i , S wri and the interface hydraulic conductivity K i were set to values equal to the surrounding fissured 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 contains 21 adjustable parameters for the fissured matrix and the conduit continuum.

Model calibration
The calibrated model shows a good fit with most of the specific characteristics of the discharge hydrograph during the period between 16 February 1988 and 20 January 1990 (see Fig. 5).Calibrated values for all varied 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.
During the time period investigated, two strong discharge events occurred, caused by major snowmelts.The discharged water volume agrees well with the simulated Figures

Back Close
Full data during the time period of the first peak.During the second discharge event (secondary peak) the simulated peak height is overestimated.It is not possible to change the relative peak height difference between the primary and secondary peak with the available calibration parameters.The recession curve slopes after 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 numerical 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 (see Fig. 6) with a difference of a few centimeters.Above the water table the matrix saturation drops to 0.35 near the surface.Flow paths in the unsaturated matrix continuum and conduit continuum are nearly vertical towards the groundwater table, whereas flow in the saturated zone is laterally oriented towards the outlet, i.e. the karst spring.
The saturation in the conduit continuum is close to zero and has a very sharp transition along the water table.In this model, unsaturated flow in the conduits is also calibrated by the k rminc parameter (minimum relative permeability of the conduits).Without this parameter, the relative permeability of the conduit continuum is a function of the residual saturation, i.e. setting of k rminc simply overrides Eq. ( 12).This is the case for most of the unsaturated conduit continuum, where saturation declines very quickly (below 0.05) above the water table for the given van Genuchten parameters.Therefore, with the applied van Genuchten parameters only, water flow in the unsaturated conduit continuum is extremely small such that exchange from the matrix into the conduit system is nearly completely prevented and a proper model calibration is impossible.However, Tokunaga and Wan (1997) showed that gravity driven film flow processes occur on unsaturated 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 meaning to the k rminc 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 expected to be fully resolved by the model.Attempts to refine the original Introduction

Conclusions References
Tables Figures

Back Close
Full van Genuchten approach and include hydraulic features of fractures into a continuum model have been made for example by Ross and Smettem (1993), Durner (1994) and Brouy ère ( 2006) by constructing a continuous bi-modal retention curve.An important role for the water exchange in the double continuum approach plays the exchange parameter α ex .It determines 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 example, is higher for a dendritic system than for a single conduit 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.

Single variation of hydraulic parameters for saturated flow conditions
Figure 7 shows the computed spring discharge for several model runs with varying hydraulic conductivity K c and porosity θ c in the conduit continuum.These parameters strongly influence the simulated spring discharge.An increased conduit conductivity K c results in sharp discharge peaks, very steep recession curves and low base flow levels because the water is quickly released from the conduit continuum.Decreased conduit conductivity K c favors a slow recession and flattens the discharge curve.Discharge peaks are broadened and base flow is higher.In case the conduit drains the matrix system an increase of K c enhances the exchange process between matrix and conduits by decreasing the hydraulic gradient in the conduit continuum and consequently increasing the hydraulic gradient between matrix and conduits.Similar effects are observed for a variation of the conduit porosity θ c .A contrasting behavior is observed by a variation in matrix porosity θ m .With an increase in the parameter the water transfer between the continua during recharge events is lower because of the lower head Introduction

Conclusions References
Tables Figures

Back Close
Full difference between conduit and matrix and the discharge curve is smoothened accordingly (see Fig. 7).K m displays a low sensitivity within the given parameter range (see Table 2) 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 matrix such that water will not be restrained within the matrix system where 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 magnitudes.A reduction to 0.001 lowers the peak height of both main peaks while decreasing recession curve slopes and slightly increasing base levels.Further reduction to 0.0001 drastically decreases peak heights and increases the base levels.The resulting slopes of recession curves are very gentle and recharge events show no more pronounced peaks.An exchange 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.

Single variation of unsaturated zone parameters
Variations of the sensitive van Genuchten parameters for the vadose zone are shown in Fig. 8.The decrease of α m and β m results in a strong rise of peak heights and increase of recession slopes.The residual water saturation of the matrix system (S wrm ) is less sensitive with respect to the modeled hydrograph.The influence of the van Genuchten α m parameter on the discharge curve is connected to the inter-continuum water exchange process.Lowering the parameter increases the capillary rise, i.e. the matrix has higher saturation (and relative permeability) above the water table .Consequently the increased permeability 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 3.65 where the saturation fringe declines very quickly with lower saturations above the water table.This reduces the main exchange interface 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.The van Genuchten parameter β m can be considered insensitive compared to α m .The conduit van Genuchten parameters α c and β c are 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 always 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 characteristic infiltration behavior and travel time distribution as well as an exchange interface in the unsaturated zone.In this work the exchange process in the unsaturated zone can be controlled 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 capillary fringe in vertical direction) which is not known in real karst systems.As described before the k rminc 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 water saturation of the matrix S wrm and the minimum relative permeability of the conduits k rminc both show a similar behavior regarding parameter variations.Increasing the parameters yields an enhanced exchange from the matrix to the conduit continuum due to a higher relative conductivity.Consequently recharge events are transmitted faster to the model outlet, i.e. the spring.

Combined parameter variations
The above presented sensitivity analyses imply only one single parameter varied at a time.However, a further important observation is that certain parameter Figures

Back Close
Full 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.For example, the simultaneous variation of the matrix van Genuchten parameter α m and the conduit conductivity K c displays a pronounced sensitivity for low α m values (Fig. 9).While for the calibrated α m value of 0.0365 m −1 the conduit conductivity K c 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 K c value of 10 m s −1 , i.e. δRMSE(α m = 0.00365)/δK c is much higher.A similar behavior can be shown for a combination of matrix porosity θ m and the conduit conductivity where lower porosities yield higher RMSE values with an increase in conduit conductivity to 100 m s −1 whereas for rather high matrix porosities of 0.102 the increase in RMSE is less pronounced.The conduit conductivity K c exhibits a higher sensitivity for the calibrated exchange parameter α ex = 1.0 (see Fig. 9) such that a high conductivity value (100 m s −1 ) results in RMSE values 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 conductivity values of the conduit system while it is nearly insensitive for low values (1 m s −1 ).In sum, the variation of the exchange parameter influences the discharge curve depending on the combination with other parameters.This behavior can also be observed for the combination of matrix porosity θ m and exchange parameter α ex .Here the exchange parameter 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.

Conclusions
The applied double continuum model lumps the dual flow behavior of a karst aquifer, where conduits represent local features of the aquifer system.The advantage of the approach is that only limited information about the geometry of the aquifer system is Figures

Back Close
Full necessary.Due to their large volume, vertical conduits in a karst unsaturated system would act as flow barriers if simulated by the Richard's equation.However, flow in vertical shafts is not controlled by matrix potential and capillary forces but rather flow processes dominated by gravitational forces such as film flow (Tokunaga and Wan, 1997;Tokunaga et al., 2000), turbulent film flow (Ghezzehei, 2004), droplet flow (Doe, 2001;Dragila and Weisbrod, 2004) and rivulet flow (Su et al., 2001(Su et al., , 2004;;Dragila et al., 2004).In order to be able to use a consistent modeling approach, boundary conditions were modified and conduit recharge was directly injected at the bottom of the saturated conduit system.This procedure allows the simulation of rapid recharge with the given modeling code.The approach is successfully employed to simulate the discharge curve of the karst system Gallusquelle for a period of two years.Because of the high amount of calibration parameters of the saturated-unsaturated flow model, a comprehensive 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 sensitivity study demonstrates that the simulation of karst hydraulics requires a-priori knowledge about parameter ranges of model variables to reduce ambiguity of the model.However, especially for unsaturated zone parameters in double continuum karst systems, only little information about the parameter ranges is documented and further research is needed.Furthermore, the analysis shows that the sensitivity of a parameter depends to a large degree on the other calibrated model parameters.
Therefore, sensitivity analyses should simultaneously take into account parameters of both continua in order to detect deviations from a linear behavior if both parameters 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 simulation.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 continuum.It should be noted that the double-continuum approach assumes Darcian flow for the matrix as well as the conduits.Considering the Introduction

Conclusions References
Tables Figures

Back Close
Full high flow velocities in the conduits it is apparent that strictly Darcian flow will underestimate the heads (no energy loss due to friction) and consequently the exchange from matrix to conduits when the conduit continuum is draining the matrix system.More realistic results may be obtained by evaluating these influences for example by applying turbulent flow in the conduit continuum (Shoemaker, 2008;Reimann et al., 2011b).The van Genuchten parameters of the matrix are the most crucial property in terms of sensitivity, uncertainty and model limitations.The exchange process between matrix and conduit continuum is mainly controlled by differences 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 defined 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 conduit system (i.e.dendritic vs. large single conduits).If the parameter is used to calibrate the model by limiting water transfer between continua, attention should be paid to the non-linear behavior of certain parameter combinations and their resulting sensitivities.The application of van Genuchten parameters to fractured aquifer systems treats them as upscaled calibration parameters.Local scale flow processes, e.g.film and droplet flow along fracture surfaces, are not physically represented.Additionally, the dual-continuum approach lumps the geometrical features of the conduit system and the fissured matrix blocks, respectively, in the saturated and unsaturated zone.Introduction

Conclusions References
Tables Figures

Conclusions References
Tables Figures

Back Close
Full 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 precipitation and evapotranspiration, are a result of the particular properties of the individual compartments: soilepikarstic zone, vadose zone, and phreatic zone.Each of these compartments is, in turn, characterized by two coupled flow systems: a highly permeable one with low Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2) where w m and w c are the volumetric fractions of each continuum of the total porosity, such that w m = 1.0 − w c .S wm and S wc are the water saturations of the respective continuum and R m and R c denote a volumetric fluid flux per unit volume (source/sink Discussion Paper | Discussion Paper | Discussion Paper | term) for each continuum.The effective matrix porosity θ sm and conduit porosity θ sc are related to the water content of the matrix θ m and of the conduit θ c according to θ m = S wm θ sm (3) 10) for ψ < 0, where S wrm , S wrc and S wri are the residual saturations, α m , α c and α i denote the inverse air-entry pressure head, β m , β c and β i 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 pressure head of the matrix.In the saturated zone where Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 10 km with a vertical thickness of 225 m.It reflects a cross section parallel to the direction of flow in the Gallusquelle catchment.The model domain is represented by two continua reflecting flow Discussion Paper | Discussion Paper | Discussion Paper | 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 component at the bottom of the conduit continuum.The amount is in the range of the rapid recharge component Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Su, G. W., Geller, J. T., Hunt, J. R., and Pruess, K.: Small-scale features of gravity-driven flow in unsaturated fractures, Vadose Zone J., 3, 592-601, 2004.1531 Teutsch, G.: Grundwassermodelle im Karst: Praktische Ans ätze am Beispiel zweier Einzugsgebiete im Tiefen und Seichten Malmkarst der Schw äbischen Alb, Ph.D. thesis, Universit ät T übingen, 1988.1517 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.1517 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.1517 Therrien, R., McLaren, R., Sudicky, E., and Panday, S.: HydroGeoSphere: a three-dimensional numerical model describing fully-integrated subsurface and surface flow and solute transport, Manual (Draft), HydroGeoLogic Inc., Herndon, VA, 2006.1516, 1518 Tokunaga, T. K. and Wan, J.: Water film flow along fracture surfaces of porous rock, Water Resour.Res., 33, 1287, doi:10.1029/97WR00473,1997.1526, 1531 Tokunaga, T. K., Wan, J., and Sutton, S. R.: Transient film flow on rough fracture surfaces, Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

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.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 2 .
Parameter combinations used for the sensitivity analyses and highest RMSE (mm d −1 ) observed.Bold RMSE values denote a non-linear inter-parameter dependence.