A new approach to model the spatial and temporal variability of recharge to karst aquifers

In karst systems, near-surface dissolution of carbonate rock results in a high spatial and temporal variability of groundwater recharge. To adequately represent the dominating recharge processes in hydrological models is still a challenge, especially in data scarce regions. In this study, we developed a recharge model that is based on a conceptual model of the epikarst. It represents epikarst heterogeneity as a set of system property distributions to produce not only a single recharge time series, but a variety of time series representing the spatial recharge variability. We tested the new model with a unique set of spatially distributed flow and tracer observations in a karstic cave at Mt. Carmel, Israel. We transformed the spatial variability into statistical variables and apply an iterative calibration strategy in which more and more data was added to the calibration. Thereby, we could show that the model is only able to produce realistic results when the information about the spatial variability of the observations was included into the model calibration. We could also show that tracer information improves the model performance if data about the spatial variability is not included.


Introduction
For a sustainable groundwater management, detailed quantitative knowledge about groundwater recharge is required (Vries and Simmers, 2002).In karst regions, the epikarst, which develops due to higher dissolution activity of the carbonate rock near the surface (Williams, 1983), controls recharge dynamics.Field research (Aquilina et al., 2006;Williams, 1983Williams, , 2008) ) as well as modeling approaches (Kiraly et al., 1995;Perrin et al., 2003) revealed that the epikarst acts as a temporary storage and distribution system for infiltrating water into karst systems.
Physically based approaches, on the one hand, are described by Hughes et al. (2008), Kiraly et al. (1995), Martínez-Santos and Andreu (2010) and Perrin et al. (2003).They require extensive data collection to characterize system properties (Le Moine et al., 2008), but may provide spatial information about recharge rates (Martínez-Santos and Andreu, 2010).Lumped approaches, on the other hand, were used by Fleury et al. (2007), Geyer et al. (2008), Jukic and Denic-Jukic (2009b) and Tritz et al. (2011), mostly as a subroutine of a model of an entire karst system.These approaches are based on a set of equations transferring input to output, conceptually representing physical processes (Hartmann et al., 2012).Because they are easy to implement, they are widely used in karst modeling (Tritz et al., 2011).
Unfortunately, both approaches have deficiencies.Due to the complexity of hydrogeological characteristics, the parameterization of physically based models is usually not possible (Jukic and Denic-Jukic, 2009a).In contrast, lumped approaches include karst-specific processes with strong simplifications and only provide one single recharge time series for the entire system (Scanlon et al., 2002).
In this study we aim to combine the advantages of distributed and lumped recharge modeling approaches.We hypothesize that recharge spatial and temporal variability is due to the physical variability of the epikarst, which can be represented by distribution functions of system properties.Our model produces not only a single recharge time series but also a variety of time series representing the spatial variability of epikarst recharge.To test the model we used measurements of stalactite drips in a karstic cave at Mt. Carmel, Northern Israel (Arbel et al., 2010;Lange et al., 2010).

General conceptual model of the epikarst
The epikarst develops close to the surface due to higher solution activity of infiltrating water with higher carbon dioxide concentrations.It is regarded as a temporary storage and distribution system for infiltrating water into karst systems (Williams, 1983;Aquilina et al., 2006).Temporarily, perched aquifers can develop (Mangin, 1975) to allow lateral flows to the proximate enlarged fissure or conduit (Williams, 2008).Hence, recharge to the lower karst system is (1) slow and diffuse into the fissured porous matrix and (2) fast and concentrated into the conduits (Fig. 1a).The interplay between concentrated and diffuse recharge depends on the variability of system properties, such as lateral and vertical hydraulic conductivities as well as soil and epikarst thickness.

Study area
The Orenin Cave (Fig. 1b and c) is a karstic cave, which developed in crystalline limestone located at the western escarpment of Mt.Carmel -a triangular-shaped, anticlinal, uplifted block up to 546 m above sea level.It is located close to the northwestern coast of Israel and composed of upper Cretaceous limestone, dolomites, chalks and marls.The area is intensively fractured and jointed, and characterized by various karstic features (Guttman, 1998;Karczs, 1959).The cave is located 28 m below an almost horizontal surface covered with ∼ 48 % rock outcrops and shallow soil pockets of reddish-brown, silty-clay, stony Terra Rossa soil up to 110 cm deep (Arbel et al., 2008;Wittenberg et al., 2007).The regional water table is located about 120 m below the cave within the Upper Cretaceous Judea Group Aquifer of Mt.Carmel.Aquifer recharge by freshwater is indicated by a seasonal rise in water table with a decline in chloride concentrations (Guttman, 1998) occurring after significant rainstorms during winter.
The climate in Mt.Carmel is typically Mediterranean with cool rainy winters, dry hot summers, and average daily potential evapotranspiration rates of 5-6 mm.The mean annual rainfall above the cave is about 550 mm a −1 .The rainy season lasts from October to April, but most rainfall occurs between November and March.Rainfall intensities exceed 30 mm h −1 during short and localized convective rainstorms in autumn, whereas in winter the storms are frontal events, lasting a few days with lower intensities.Significant winter rainstorms have total rainfall amounts between 40 and 180 mm (Wittenberg et al., 2007).Above the cave, vegetation cover is a typical Mediterranean garrigue, i.e. shrubs (Arbel et al., 2010).
Most recent findings about the spatial and temporal recharge variability at the cave can be found in Lange et al. (2010) and Arbel et al. (2010).The former studied drip and tracer responses of several stalactites following a sprinkling experiment (Fig. 2a and b), while the latter observed dripping rates and tracer concentrations at a seasonal time scale (Fig. 2c).Both highlight the presence of large water storages in the soil of the epikarst, which need to become saturated before the drips activate.Hydrochemical analysis indicated that the drip water was largely composed of preevent water with variable but generally low event water fractions.Even though several tracers were used in both studies, we only consider artificially enriched water with high electric conductivity (EC) since it was the only tracer that was applied uniformly over the whole area.

Transforming spatial variability into statistical variables
Instead of considering individual hydrographs and tracer concentration curves for each single drip, integrated hydrodynamic and hydrochemical responses, Q [l h −1 ] and C[-], of all measured drips are calculated: The coefficients of variation CV C [-] and CV Q [-] specify their spatial variability: where Q i (t) [l h −1 ] is the individual drip rate, and C i (t) [-] is the normalized tracer concentration of the drip i, i = 1...N , at time t.Concentrations are normalized by the tracer input concentration C in [µS cm −1 ]: where C i (t) [µS cm −1 ] is the originally observed tracer concentration at drip i.Note that in Eq. (3) Q is divided by N in order to obtain its mean from the sum (Eq.1).Using Eqs.
(1)-( 5) observations at the individual drips were transformed into three time series of integrated responses Q exp , C exp , and Q seas , and three time series of coefficients of variation, CV Q,exp , CV C,exp , and CV Q,seas , for the experiment hydrodynamic and tracer observations, and the seasonal hydrodynamic observations, respectively (Fig. 2).

Model structure
The developed model structure (Fig. 3) follows the general conceptual epikarst model of Williams (1983)

Water fluxes
Infiltration into the soil originates from precipitation and surface flow arriving from the neighboring reservoir.Water leaves the soil either as actual evaporation or as percolation to the epikarst when its storage volume S soil,i is exceeded.Similar to many other models (HBV, Bergström, 1995;TOP-MODEL, Beven and Kirby, 1979), actual evaporation E act,i [mm h −1 ] is derived from: where E pot (t) [mm h −1 ] is the potential evaporation, and V soil,i (t) [m 3 ] is the water volume stored in soil compartment i at time step t.Inflow to the epikarst either originates from soil percolation or from the neighboring epikarst compartment as a consequence of water level difference h i (t) [m], which is calculated by where V epi,i (t) [m 3 ] is the water volume stored in the epikarst at compartment i.Lateral flow Q lat,i [l h −1 ] is described by the Dupuit-Forchheimer assumption, which proved to be an adequate representation of lateral flows in perched aquifer delivery (Freeze and Cherry, 1979): where T i (t) is the transmissivity, and W i [m] is the width of flow.Since a decrease of the lateral hydraulic conductivity K lat,i [mm h −1 ] with depth can be expected (Perrin et al., 2003), a decay coefficient a lat [-] is introduced: where K lat,max [mm h −1 ] is the lateral hydraulic conductivity at the surface, and z is the depth below surface.As described in more detail in Wigmosta and Lettenmaier (1999), T i (t) is calculated by: Outflow from the saturated zone either occurs by lateral flow to the next reservoir (Eqs.12 to 16) or in the form of vertical recharge which is represented by a simple linear relationship: where R i [l h −1 ] is the vertical recharge, and K vert,i [mm h −1 ] is the vertical hydraulic conductivity.This simple relationship does not take into account that water movement is also gravity driven, it only considers flow due to water pressure.However, trying different equations describing the vertical percolation, this simple relationship was found to perform best.Following the conceptual model of Williams (1983), the variability of vertical hydraulic conductivity is included by: where K vert,max [mm h −1 ] is assumed to be the hydraulic conductivity at the uppermost part of the epikarst (Fig. 2b), and a vert [-] is a coefficient describing the variability of K vert,i .When inflows simultaneously exceed the S soil,i and S epi,i , surface flow is produced.Based on a doline structure (Fig. 1a), a surface gradient towards the doline centre can be assumed.Hence, surface flow is routed from outer compartments with low soil and epikarst thickness to the inner compartments with higher thickness: where Q out,surf,i [l h −1 ] is the surface flow produced at compartment i, and Q in,surf,i+1 is the surface flow reaching compartment i+1 from the neighboring compartment with lower soil and epikarst thickness.

Solute transport
The solute transport in all reservoirs is based on the assumption of complete mixing.Since preceding studies (Arbel et al., 2010;Lange et al., 2010) found evidence for considerable amounts of old water in the cave drips, a maximum old water volume V old,max [mm] and its solute concentration C old [-] are defined as initial conditions.It is assumed that the old water storage is stored below the epikarst (Fig. 3) and that compartments with a lower soil and epikarst thickness have larger old water storage beneath them than regions with higher soil and epikarst depths.Therefore, the same distribution function as for the soil and epikarst thickness is applied, but this time in the opposite direction: With Eqs. ( 6) to ( 20), the model produces N series of recharge rates and tracer concentrations on which Eqs. ( 1) to ( 5) can be applied.

Calibration approach
Including solute transport, the model consists of 13 parameters (Table 1) that have to be determined by calibration, since no field information is available.To identify parameter values, the Shuffled Complex Evolution Metropolis algorithm SCEM (Vrugt et al., 2003) was chosen, which has proven to reliably find optimal parameter sets including information about their uncertainty (e.g. in Feyen et al., 2007;Schoups et al., 2005;Vrugt, 2004).As a measure of model performance, the Nash Sutcliffe Efficiency NS (Nash and Sutcliffe, 1970) was used.In order to show the influence of different information sources, an iterative calibration strategy was applied (Table 2).Four individual efficiencies were calculated and equally weighted.More and more information was iteratively added to the optimization procedure until all available data was included.

Stability of simulation time series
The number of compartments N was arbitrarily set to 15, which is larger than the number of observed drip locations in the cave (in total 9).To investigate whether this number was large enough to provide numerically stable results, we ran the model using calibration step (4) but with different numbers of compartments (N = 3...50).Calculating NS Q,exp , NS Q,seas , NS C,exp , and NS CV for the different N indicated which minimum number for N was necessary to have stable mean values and coefficients of variation.The results were considered stable as long as their relative deviations remained below 20 %.

Results
The model was run with a six-hour time step.For each calibration step, SCEM was performed until convergence was reached (see Vrugt et al., 2003).Parameter ranges were set to physically reasonable values according to preceding studies.Table 1 shows the optimized parameters and the resulting NS efficiencies for all the calibration steps.The parameters were normalized by their range and compared (Fig. 4).At calibration step 1 the modeled Q exp and Q seas closely compared with the observations, which was expressed by high NS Q,exp and N S Q,seas values.However, simulated CV and C exp showed a strong bias from the observations and NS C,exp and NS CV had low values.In calibration step 2, tracer information consequentially improved the simulated tracer concentrations, even though a moderate decrease in NS Q,exp was observed.NS CV improved slightly.In calibration step 3, the information about spatial variability strongly improved NS CV , while there was also a small improvement of tracer simulations evident by an increase of NS C,exp .Only when tracer and variability information were added in calibration step 4, an acceptable simulation was reached for Q exp , C exp , Q seas and CV with all NS values exceeding 0.5.For calibration steps 1 and 2, some parameters were similar, e.g.n soil , a vert and a lat , but some parameters had unrealistic and contradicting values when only information about drip rates and tracer concentration was used (e.g., a depth and d epi,max ).When the information about spatial variability was added in calibration step 3, the majority of the parameters changed and plotted almost at the same location in calibration step 4, e.g.d epi,max , a vert and log K vert .In many cases, calibration step 2 parameters plotted between calibration step 1 and calibration step 3 and 4 parameters.Exceptions were the hydrochemical parameters, V old,max and C old , and the parameters describing the lateral flow, a lat and L. The hydrochemical parameters grouped together in calibration steps 2 and 4 when information about the tracer data was considered.The lateral flow parameters did not show any systematic pattern between the calibration steps.
Figure 5 shows the cumulative parameter distributions obtained by SCEM for each calibration step.If a parameter is sensitive, it would significantly differ from a uniform distribution (grey diagonals).Using this criterion, many of the parameters seemed sensitive already at calibration step 1.However, the lateral flow parameters, log K lat,max , a lat and L, the hydrochemical parameters, V old,max and C old , and the porosities, n soil and n epi revealed low sensitivities.V old,max and C old became sensitive in calibration steps 2 and 4 when tracer data was added.log K lat,max , a lat and L, as well as n soil and n epi , remained non-sensitive among all calibration steps, which means that they could adapt to any value without changing the simulation results.The grouping of parameters that occurred at calibration steps 3 and 4 was also reflected by a shift of the cumulative parameter distributions, which changed when the information about spatial variability of drip rates and tracer concentrations was added, e.g. for log K vert,max or a vert .
The ability of the model to produce a realistic spatial variability of recharge rates and tracer concentrations is visualized in Fig. 6.Some model compartments react delayed, others rapidly to rainfall events.On the other hand, the delayed compartments sustain flow long after water input, while the rapidly reacting compartments fall dry shortly after rainfall events.Regarding the tracer concentrations, the delayed compartments almost do not react to the input.More rapidly reacting compartments change their concentration, while the most rapid totally adapted to the input concentration.
The water balances in Fig. 7 show the simulated sum of internal model fluxes for calibration step 4, with separate graphs for the sprinkling experiment (Fig. 7a) and the seasonal time scale (Fig. 7b).Even though all compartments receive the same amount of precipitation, the remaining water balance components of the individual compartments show a high spatial variability.At both time scales, the first five compartments produce surface flow, while noteworthy amounts of actual evaporation are visible from the fourth compartment at the seasonal scale, increasing towards compartment 15.While the sprinkling experiment produces the largest recharge amounts for the last five compartments, major parts of the water are recharged by compartment four to seven at the seasonal time scale.
Increasing the number of model compartments (Fig. 8) indicates that the simulated drip rates were stable at N = 3, both for the experiment and for the seasonal time scale.The CV stabilized for N ≥ 5. Tracer concentrations stabilized for N > 10 and diverged again for N ≥ 25.The drip rates as well as the CVs remained stable even for N = 50.

Model performance along the iterative calibration
Considering only integrated recharge rates, calibration step 1 provided acceptable predictions of mean recharge rates, but unacceptable results for the spatial variability and hence internal process dynamics (Table 1, Fig. 4).When tracer data was added in calibration step 2, the tracer and also the predictions of spatial recharge variability improved slightly.However, only information about the spatial variability in calibration step 3 adequately reproduced the observed variability.Using both tracer and variability data in calibration step 4, no further significant improvement could be reached.Thus, the information about spatial variability of recharge rates had the strongest impact on the optimized parameters and on the resulting simulations, while tracer data alone had a slight positive impact on the representation of variability and, in addition, gave evidence of proper process representation and parameter choice.Including information about spatial variability resulted in small reductions of NS Q,exp and NS Q,seas in favor of an overall acceptable multi-objective fit and more realistic results.Similar conclusions can also be found in Kuczera and Mroczkowski (1998), Seibert and McDonnell (2002), and Son and Sivapalan (2007).

Water balance
Figure 7 indicates that there is almost no lateral subsurface flow during the sprinkling experiment as well as at the seasonal time scale.This suggests that lateral flow processes are not important at the study site.This is supported by tracers in preceding field research that showed that even though lateral flow processes exist, only minimum lateral flow concentration occurs during infiltration and percolation (Arbel et al., 2010;Lange et al., 2010), majorly due to a previously vertical fissure orientation (Karczs, 1959).There are also no significant amounts of actual evaporation during the sprinkling experiment, because the observation time was too short to evaporate the newly stored sprinkled water.At the seasonal time scale, actual evaporation constitutes a large part of the outflows, which is comparable to other studies (e.g.Andreo et al., 2008).
Due to the variability of soil and epikarst thickness, surface flow is only produced at the first six compartments during the experiment at the seasonal time scale (Fig. 7), which is 40 % of the modeled area.This fraction closely compares with 48 % of rock outcrops (Arbel et al., 2008;Wittenberg et al., 2007) that are found at the surface.During the sprinkling experiment, the largest recharge amounts can be found in the last five model compartments, whereas at the seasonal scale, compartments five to seven produce the largest amounts of recharge.That indicates that major parts of the recharge are rather slow.Only during strong events, as artificially produced by the sprinkling experiment, fast flow paths activate and produce much recharge in a very short time.Using water isotopes and lumped parameters models, Maloszewski et al. (2002) arrived at similar results for a karst system in the Austrian Alps.

Optimized parameters and water balance
The final parameters derived from calibration step 4 can only be interpreted when their sensitivity is considered.The parameter distributions in Fig. 5 show that all parameters were sensitive when all available information was used for the calibration procedure, except for lateral flow parameters, log K lat,max , a lat and L, and the porosities, n soil and n epi .This can be interpreted as an indicator for over-parameterization (Perrin et al., 2001).However, since lateral subsurface flow was found to be not significant at the study site, a nonsensitivity of these parameters is rather due to the nonimportance of the lateral flow processes.The non-sensitivity of the porosity parameters derives from Eqs. ( 8) and ( 9).The same storage volume can be generated by different combinations of area, depth, and porosity.Since the ranges of area and depth are much wider than those of porosities (see discussion below), the porosities appear non-sensitive because their variations can be totally compensated by variations of area and depth.Similar results were obtained with a splitsample test and multiple bootstrapping of subsets of the observation period (see the Supplement).
The sensitive parameters of calibration step 4 suggested that the maximum soil depth d soil,max was quite low (14 cm) compared to measurements in the field (up to 110 cm, Arbel et al., 2008;Wittenberg et al., 2007).The soil porosity n soil was set to 35 % to 45 % and did not consider the stoniness.Stoniness may drastically reduce overall porosity and would therefore result in a larger calibrated d soil,max .However, the correctness of the final soil volumes S soil,i was corroborated by the realistic amounts of actual evaporation.A part of the maximum measured soil thickness might also be attributed to the epikarst.In the model, its thickness d epi,max was close to the rock thickness above the cave (27 m).Due to the depth variability coefficient a depth , very small storages S soil,i and S epi,i can be found at the first compartments, which resulted in noteworthy rates of surface flow (see Sect. 5.2).The epikarst porosity n epi , although non-sensitive, was calibrated to a realistic value (19 %) compared to other studies (e.g.Williams, 2008).The hydrologically small contributing area A of ∼ 9 m corresponded to the findings of the sprinkling experiment (Lange et al., 2010), which showed that only small parts of the sprinkling water finally reached the cave.The maximum vertical hydraulic conductivity K vert,max was calibrated to 55 mm h −1 (1.46 × 10 −5 m s −1 ), which corresponds with other studies (Williams, 1985;Perrin et al., 2003).
The calibrated old water concentration C old was found to remain at the uppermost part of the predefined range (Fig. 4).A widening of the C old calibration ranges resulted in slightly better model results for NS C,exp , but since the original ranges were set according to the old water concentrations determined for all the drips prior to the sprinkling (Lange et al., 2010), physical realism was prefered instead of an optimal fit.The maximum old water volume V old,max of ∼ 1000 mm is in accordance with Lange et al. (2010) and Arbel et al. (2010), who discovered large old water contributions in the rock above the cave.

Representation of spatial variability
Our model reproduces the spatial variability of recharge rates every time step because distribution functions are included.They control the variability of soil and epikarst depth, vertical hydraulic conductivity and old water storage.The assumption that the observed drip rates in the cave show the real spatial variability of recharge rates raises the question whether one drip really represents one individual flow path.Since the calcite precipitation, which formed the stalactite, produces a less permeable layer on the cave ceiling, it is reasonable to expect an accumulation of several flow paths before the water finally reaches the cave.For that reason we used integrated responses and coefficients of variation instead of directly comparing individual model compartments with individual drip observations, and we chose a number of model compartments (15) that is considerably larger than that of observed drips (9).Varying the number of model compartments (Fig. 8) showed that N = 15 model compartments were enough for stable results for both integrated flow rates and their coefficients of variation.Only hydrochemical predictions diverged for compartment numbers > 25, which was most probably due to their dependence on the flow predictions.Small changes in Q exp might result in very large changes of C exp .

Transferability of the new approach
Even if the introduced model performs well and is physically reasonable under the study site's conditions, the question about the transferability of the approach to other sites has to be addressed.In detail one has to ask whether the model is able to cope with differences in system properties and whether the model can be applied without information about the spatial variability of recharge dynamics.We theoretically consider two different caves (1) the Soreq cave, Israel, and (2) the Vers-chez-le-Brandt cave, Switzerland, to discuss this question.
The Soreq cave is located in a semi-arid climate 10-50 m below the surface at a sloping terrain.Using drip rate observations at several observation points within the cave and environmental tracers, it was found that, depending on the degree of fractures, distance to the surface, and rainfall duration and intensity, fast flow can occur (Even et al., 1986;Ayalon , 1998).However, most recharge water travels several decades in small fissures before it reaches the cave (Kaufman et al., 2003).Our model already proved that it is able to cope with fast as well as with slow flow paths.However, due to the sloping terrain, lateral flow processes that have shown to be not as important at our study site may be of higher significance at the Soreq cave.On such terrain, surface runoff should also be accounted for, as it has been shown on a slope nearby (Lange et al., 2003).
The Vers-chez-le-Brandt cave is located below relatively flat pasture land in a humid environment, 30 m below the surface with 1-2 m deep soil on top.Drip rate observations at one observation point in the cave and artificial and environmental tracers showed that its typical dynamics during a rain event are characterized by five phases (Pronk et al., 2009): saturation of the soil, initiation of pressure pulse by perched water, arrival of first event water at the cave, increased amounts of event water at the cave, and recession period.Since our model includes a soil layer and allows for perched water tables, we are confident that our approach will most probably be able to reproduce these five characteristic phases.Compared to our site, the soils are much thicker and distributed more evenly.Hence, allowing a deeper soil in the model setup will address this difference.The lack of information about the spatial variability of recharge could be reduced by directly measuring soil depth distributions (as for example in Tromp-van Meerveld and McDonnell, 2006).However, information about the distribution of vertical conductivity, water storage, and decrease of lateral conductivity with depth is difficult to obtain.To compensate for this, we see a high value in karst evolution studies and models, which provide aperture width distributions depending on climate and degree of fractures before karstification began (e.g.Bloomfield et al., 2005;Hubinger and Birk, 2011).In this study we hypothesize that the spatial variability of epikarst recharge is due to the variability of system properties, which can be represented by simple distribution functions.Based on our conceptual understanding of the epikarst, we developed a model structure that includes the variability of selected system properties.Our results indicate that our model is able to realistically produce the spatial variability of recharge and tracer fluxes.The proposed iterative calibration strategy reveals that these results are realistic only when information about spatial variability was included into model calibration.
The use of simple parametric distribution functions of model properties and parameters in an otherwise lumped model may allow reproduction of the temporal and spatial variability of karstic recharge if variability is considered in the model structure and calibration.However, before more general statements on the performance of the newly developed model for practical applications can be made, it has to be compared with other already existing approaches.With the information provided by this study, we can only confirm that the model works well under the observed conditions and data.The model should be applied to other sites with different controls of climate, processes and epikarst properties to evaluate its transferability before it can finally be applied for practical purposes, e.g.water resources management.This will be the scope of further research.

Fig. 2 .
Fig. 2. Individual observations, integrated response and coefficient of variation of (a) drip rates during the sprinkling experiment (Lange et al., 2010), (b) drips' normalized EC concentrations during the sprinkling event(Lange et al., 2010), and (c) drip rates at the seasonal time scale(Arbel et al., 2010).

Fig. 3 .
Fig. 3. Model structure and parameters; parameter a lat controls the decay of K lat,max with depth z, and parameters a vert and a depth control the variability of K vert,max , V old, max , d soil,max and d epi,max among the N compartments. h

Fig. 4 .
Fig. 4. (a) Normalized optimized parameter sets for the four calibration steps and (b) the resulting NS efficiencies.

Fig. 6 .
Fig. 6.Simulation results of all individual model compartments for calibration step (4) and the integrated responses and CVs of all calibration steps for (a) drip rates during the sprinkling experiment, (b) drips' concentrations during the sprinkling event, and (c) drip rates at the seasonal time scale.

Fig. 7 .
Fig. 7. Water balance for each model compartment: (a) sprinkling experiment and (b) seasonal observations for calibration step 4 (the change of storage within the compartments is attributed to the output).

Fig. 8 .
Fig. 8. Maximum relative deviation of NS efficiencies with varying number of compartments from the NS efficiencies using N = 15 compartments.
soil,i [m] and d epi,i [m] are the soil and epikarst thicknesses of reservoir i, d max,soil [m] and d max,epi [m] are the maximum soil and epikarst thicknesses, and a depth [-] is the depth variability coefficient.Defining soil and epikarst porosities n soil[-]and n epi[-], respectively, the storage volumes S soil,i [m 3 ] and S epi,i [m 3 ] are defined by:

Table 1 .
Parameter table, calibration ranges, optimized parameter sets and their efficiencies for all four calibration steps.

Table 2 .
Description of calibration steps and weights associated with the Nash Sutcliffe efficiencies concerning the drip rates observed during the experiment, NS Q,exp , the drip rates during the seasonal times scale, NS Q,seas , the drip water concentration during the experiment, NS C,exp , and their coefficients of variation, NS CV .
*In this calibration step only the variability of CV Q,exp and CV Q,seas are considered.