Process-based karst modelling to relate hydrodynamic and hydrochemical characteristics to system properties

More than 30 % of Europe’s land surface is made up of karst exposures. In some countries, water from karst aquifers constitutes almost half of the drinking water supply. Hydrological simulation models can predict the largescale impact of future environmental change on hydrological variables. However, the information needed to obtain model parameters is not available everywhere and regionalisation methods have to be applied. The responsive behaviour of hydrological systems can be quantified by individual metrics, so-called system signatures. This study explores their value for distinguishing the dominant processes and properties of five different karst systems in Europe and the Middle East. By defining ten system signatures derived from hydrodynamic and hydrochemical observations, a processbased karst model is applied to the five karst systems. In a stepwise model evaluation strategy, optimum parameters and their sensitivity are identified using automatic calibration and global variance-based sensitivity analysis. System signatures and sensitive parameters serve as proxies for dominant processes, and optimised parameters are used to determine system properties. By sensitivity analysis, the set of system signatures was able to distinguish the karst systems from one another by providing separate information about dominant soil, epikarst, and fast and slow groundwater flow processes. Comparing sensitive parameters to the system signatures revealed that annual discharge can serve as a proxy for the recharge area, that the slopes of the high flow parts of the flow duration curves correlate with the fast flow storage constant, and that the dampening of the isotopic signal of the rain as well as the medium flow parts of the flow duration curves have a non-linear relation to the distribution of groundwater storage constants that represent the variability of groundwater flow dynamics. Our approach enabled us to identify dominant processes of the different systems and provided directions for future large-scale simulation of karst areas to predict the impact of future change on karst water resources.


Introduction
Almost one third of Europe's land surface is composed of karst exposures (Williams and Ford, 2006).In some countries up to 50 % of drinking water is obtained from karst aquifers (Zwahlen, 2003).Projected trends of increasing temperatures and decreasing precipitation (Christensen et al., 2007) may affect water security in karst water regions (e.g.Butscher and Huggenberger, 2009).Hydrological simulation models are necessary to predict the large-scale impact of future environmental change on hydrological variables Published by Copernicus Publications on behalf of the European Geosciences Union.(Wagener, 2007).The strong subsurface heterogeneity of karstified rocks (Bakalowicz, 2005) means that the hydrological behaviour of karst systems can be very distinct from other hydrological systems (Goldscheider and Drew, 2007).Therefore, hydrological models containing an adequate representation of specific karst hydrological processes have to be applied.
Process-based karst models can be separated into lumped and distributed modelling approaches.Distributed approaches discretise the entire karst system in two-or threedimensional elements and provide spatial information about groundwater levels in each element.Many similar reviews concerning different subtypes and applications can be found in the literature (Ford and Williams, 2007;Goldscheider and Drew, 2007;Kovacs, 2003;Sauter et al., 2006;etc.).Since parameterisation requires spatial information on karst system properties, distributed approaches were either applied at well-explored test sites (e.g.Doummar et al., 2012;Birk et al., 2005;Hill et al., 2010) or for theoretical calculations to understand the general behaviour of karst hydrology (e.g.Reimann et al., 2011;Birk et al., 2006).Lumped approaches do not require spatial information about system properties.They consider physical processes by a set of equations that transfer input to output at the scale of the entire karst system (Hartmann et al., 2012a).In preceding studies, conceptual modelling approaches considered karst processes such as separate conduit and matrix systems (Fleury et al., 2009;Geyer et al., 2008;Maloszewski et al., 2002), storage and recharge concentration in the soil and epikarst (Hartmann et al., 2012c;Tritz et al., 2011), allogenic contribution by sinking streams (Le Moine et al., 2008;Bailly-Comte et al., 2012) or discharge by various springs (Rimmer and Salingar, 2006;Charlier et al., 2012).
Because of their integrating structure, the parameters of lumped process-based approaches describe the representative properties of the system and are therefore difficult to measure.For that reason, they are usually derived by calibrating the model with time series of discharge observations at the karst spring (Moussu et al., 2011).In order to avoid overparameterisation (Perrin et al., 2001;Beven, 2006), most of the lumped modelling studies mentioned above used rather simple model structures and omitted some karst processes deemed not important at their respective study sites.Due to their simplicity, these models are difficult to transfer to other sites.Even though research has recently made much progress in this field (Anderson and Goulden, 2011;Carrillo et al., 2011;Harman and Sivapalan, 2009;Oudin et al., 2010), this is one of the reasons why studies addressing the transfer of karst models to ungauged basins are rarely found.
In this study, a realistic process-based karst model is calibrated using ten hydrodynamic and hydrochemical karst system signatures, i.e. metrics that express a system's response behaviour and storage characteristics (Wagener et al., 2007), to five study sites around Europe and the Middle East.A stepwise model analysis is used to identify optimum parameters and parameter sensitivity.Assuming that the model adequately represents the karst systems, the sensitiveness of parameters that control the different process dynamics in the model can serve as a proxy for dominant natural processes and optimised parameters as approximations of system properties.We use this analysis (1) to explore the information content of the different karst system signatures concerning different karst processes and properties, and (2) to establish relations between parameter values and system signatures.

Study sites
We consider five karst systems in Europe and the Middle East (Fig. 1).They are located in different climatic regions and cover a range of scales of approximately 0.1 to 500 km 2 .Table 1 summarises their general characteristics.The systems are drained by one or several karst springs.The Austrian site (Fig. 1c) is dominated by Norian dolomite (Hauptdolomit), partly overlain by plattenkalk and Jurassic/Cretaceous limestone and marls (Kralik et al., 2009;Kralik and Keimel, 2003).The Israeli site (Fig. 1f) consists of two major karst springs that drain a Jurassic limestone aquifer with thicknesses of more than 2000 m.Preceding studies (Hartmann et al., 2013b;Rimmer and Salingar, 2006) showed that, for the purpose of system signature modelling, the groundwater systems of the two springs are not directly connected to each other.For that reason, they are regarded separately in this study (referred to by Israeli site 1 and 2).The Palestinian site (Fig. 1e) is a large karst spring draining an Eocene calcareous rock aquifer in a semi-karstified area (Ghanem, 2005;Hartmann et al., 2012b).The Spanish site (Fig. 1d) consists of a main spring and an overflow spring that drain a karst aquifer composed of Jurassic limestones and dolostones with a variable thickness that can exceed 500 m.The base of the aquifer is constituted by Triassic clays and evaporites (Barberá and Andreo, 2011).The Swiss site (Fig. 1b) is a small karst spring located on a karst plateau of Swiss Tabular Jura.Its aquifer consists mainly of Oxfordian limestone with a thickness of 40-70 m (Butscher and Huggenberger, 2008).

Available data
Table 2 summarises all available data.At the Palestinian site information is limited to monthly discharge measurements.A complete record of discharge and hydrochemical parameters exists for the Spanish study site.For Switzerland the daily discharge record shows a gap of 4 yr.For all sites, except for Palestine, hydrochemical parameters are mostly in a weekly to monthly resolution.At the Israeli sites the time when the δ 18 O samples were taken falls outside the time span of the discharge record.The way these data could still be included in the analysis will be elaborated on in the following subsections.

Karst system signatures
For a complete description of the karst systems' characteristics.we define ten system signatures that describe a wide range of aspects of their combined hydrodynamic and hydrochemical behaviour.Table 3 provides the description of the karst-specific system signatures and equations for their calculation.To consider the hydrodynamics, we separate the flow duration curves of the springs (FDCs, Fig. 2a) into the slopes of high flows (exceedance probabilities 0.0 to 0.1),  span 2002-2005 1989-1999 1989-1999 2007-2010   median flows (0.1 to 0.9) and low flows (0.9 to 1.0).In addition, we consider the autocorrelation of discharge time series (Fig. 2b), which is classically used to determine the memory effect of karst systems (Mangin, 1984).It already proved itself to contribute more data to the calibration of karst models (Moussu et al., 2011).Similar to Laroque et al. (1998), we use cross-correlation to characterise the delayed response of NO 3 compared to discharge (Fig. 2c).A linear regression in the log-log space describes the correlation of SO 4 and discharge (Fig. 2d), whereby the regression slope addresses its dynamics and its offset is related to the SO 4 mass balance.
Information inherent in the δ 18 O signal is expressed by the ratio of its variability in discharge and precipitation (Fig. 2e).
Water balance and inter-annual memory of the systems are considered by the annual discharges (Fig. 2f) and the streamflow elasticity (Sawicz et al., 2011).Since these measures provide time-independent descriptions of the karst systems' characteristics, it is possible to include data that were collected during different time periods than the discharge (such as δ 18 O for the Israeli site).

The karst model
In this study we use the process-based VarKarst model introduced by Hartmann et al. (2013a).It consists of subroutines representing the soil, the epikarst, and the groundwater system.The variability of system properties is expressed by distribution functions that consider the variability of soil and epikarst storage capacities, epikarst hydrodynamics, Table 3. Karst-specific system signatures, their equations and description (cov: covariance, var: variance, std: standard deviation, P (X y ): exceedance probability of variable X at the probability interval y, N: number of time steps, log 10 : decadal logarithm, dX: inter-annual change of annual variable X).
Characterises the memory effect of with r Q,100 = 1 Characterises residence time variability Characterises the dynamics of matrix- recharge separation (diffuse/concentrated) and groundwater hydrodynamics (Fig. 3).Similar to other models that consider variability (Hartmann et al., 2012c;Moore, 2007), the Pareto function is used to attribute the variable system properties to a set of N = 15 model compartments.The inclusion of spatial information makes the VarKarst model a hybrid model between lumped and distributed modelling approaches.With such structure it is much more elaborate than classical lumped models (e.g.Fleury et al., 2007;Geyer et al., 2008), but it still has a relatively low number of parameters.
A list of all model parameters is provided in Table 4, and a full description of the VarKarst model is provided in the Appendix.Hartmann et al. (2013a) showed that compared to a classical reservoir model, the VarKarst model provided superior multi-objective performance when hydrochemical information was considered.It was able to consider a time variant recharge area and gave more stable predictions when a split-sample test (Klemeš, 1986) was performed for validation.Snowmelt routines were set on top of the model for the Austrian and Swiss site.They are based on the snow routine of the HBV model (Lindström et al., 1997).A detailed  description of the routine and the selection of the parameters for the Austrian site are provided in Parajka et al. ( 2007) and Hartmann et al. (2012a).The parameters for the Swiss site were adopted from Schulla (1997), who modelled snowmelt at a nearby site.

Model calibration and sensitivity analysis
The model was applied using a modified version of a model evaluation strategy presented by Hartmann et al. (2013b).Our analysis consists of three stages: (1) evaluation of model performance with respect to system signatures, (2) evaluation of parameter identifiability using sensitivity analysis, and (3) combination of the results of stages ( 1) and ( 2) to establish relations between sensitive calibrated model parameters and system signatures.Stage 1 will show whether the model has enough degrees of freedom to reproduce the karst system characteristics that are expressed by the system signatures.
The model is calibrated on each single signature individually by comparing modelled and observed signatures and using automatic calibration by the Shuffled Complex Evolution Metropolis algorithm (Vrugt et al., 2003).That way, 10 × 6 optimum parameter sets are found for each of the ten signatures and the six study sites.If one of the modelled signatures deviates more than 30 % from the observed signature, this signature will not be used for the further analysis.Stage 2 evaluates the information provided by each signature.Sobol sensitivity analysis (Saltelli et al., 2008) is used to evaluate the sensitivity of the model parameters concerning the different signatures.This allows distinguishing informative sensitive) parameters from non-informative (insensitive) parameters.Sobol sensitivity analysis decomposes the model output variance concerning a certain signature into relative contributions from individual parameters and their interactions (van Werkhoven et al., 2008;Saltelli et al., 2008): where D is the total output variance, D i the model output variance due to parameter i, and D ij the output variance due to interactions of parameter i and parameter j .D ij k and D 12...m describe the interactions of three or more parameters; m is the total number of parameters.Information about parameter sensitivity is provided by the total contribution of a parameter to the model output variance T (also referred to as total sensitivity): where T,i is a sensitivity index that represents the single effect of parameter i plus its interactions with the model output variance and D / ∈i the model output variance produced by all model parameters except for parameter i.The contribution of individual parameters to the model output variance is described by the sensitivity index F (also referred to as first-order sensitivity):  parameters as sensitive if they are equal to or larger than 0.2 and 0.1 for T and F , respectively.In stage 3, the calibrated values of sensitive parameters concerning F are compared with values of the system signatures.Doing so for all five study sites, relations between parameters and signatures can be revealed.The parameter interactions T − F are hereby used as a proxy of the uncertainty of sensitive parameters.If parameter interactions are large, the calibrated value of a parameter may still vary when interacting parameters change, and its relation to the system signatures might be biased.The VarKarst model was already evaluated by multi-variate calibration in Hartmann et al. (2013a).Therefore we do not perform a multi-objective calibration during this analysis.Since we only consider sensitive parameters during stage 2 and 3, an interpretation of their values is possible without the problem of equifinality (Beven, 2006).

Identification of dominant processes and karst system properties
Preceding studies already showed ways to separate different processes, usually along the course of an iterative or stepwise calibration (Fleury et al., 2009;Hogue et al., 2006;Jukic and Denic-Jukic, 2009).In this study we use sensitivity analysis on different signatures to explore separately different processes in the karst systems.Similar to Carrillo et al. (2011) we assume that the model is an acceptable representation of the hydrological system.Thus, dominant processes can be identified by considering T for the different signatures and parameters that control the different process dynamics in the model.That way soil storage behaviour (mean soil storage capacity V mean,S and its distribution a SE ), epikarst storage and dynamics (mean epikarst storage capacity V mean,E , its distribution a SE , and its mean storage constant K mean,E ), recharge dynamics (distribution of diffuse and concentrated recharge a fsep ), fast (conduit storage constant K C and the critical volume to activate overflow springs V crit,OF ) and slow groundwater flow dynamics (distribution of groundwater storage constants a GW ), and water balance (A) can be explored.In addition the dissolution dynamics of SO 4 can be revealed (geogenic contributions Geo SO 4 and their distribution a Geo ).Depending on the hydrodynamic or hydrochemical aspect of the system behaviour they consider, the system signatures and the parameter sensitivity concerning them will reveal different processes for the different sites.
All of them together will provide an overall description of the dominant processes of the karst system with the current data availability.

Relations between system signatures and system properties
Assuming again that the model is an acceptable representation of the hydrological system, the parameters of the VarKarst model can be regarded as proxies of system properties.All calibrated parameters that have high first-order sensitivity F can be attributed to the system signatures they were derived from.When we compare pairs of parameter values and system signatures for all study sites, relationships can be established.If the correlation is large enough, these relations can be used to obtain model parameters and hence system properties simply by knowing the value of the respective system signature.

Model performance and parameter sensitivity
Table 5 provides the values of all system signatures for the different study sites.The test of performance in evaluation stage 1 showed that the model is able to reproduce almost all of them (Table 6).Some small deviations occurred for R Q,100 at Israel 1 and S SO 4 at Israel 2. For V δ 18 O , the Swiss and the Israeli 2 sites show stronger deviations.While ∼ 20 % of deviation for the Swiss site was regarded as still acceptable, V δ 18 O was discarded for the Israeli 2 site for the following  analysis; 80 % of deviation clearly indicated deficiencies in the performance of the model in simulating the δ 18 O variability, which would strongly bias the proceeding analysis.
In stage 2 of the evaluation, the total sensitivity index T concerning all available and not discarded signatures (Fig. 4) shows that similar patterns of sensitive parameters for all sites could be found among some of the signatures: the high flows S HF and the autocorrelation of discharges R Q,100 show always high sensitivity indices for the parameters representing the fast groundwater dynamics (conduits system K C and overflow spring V crit,OF ).In addition, V crit,OF always has a high sensitivity index for the streamflow elasticity E Q .The distribution coefficients of soil and epikarst storages a SE and the distribution of groundwater storage constants a GW show always high sensitivity indices for low flows S LF , the Q-NO 3 cross-correlation L NO 3 and the δ 18 O variability V δ 18 O .a GW always has a high sensitivity index for the medium flows S MF .For all sites, the recharge area A shows a high sensitivity index for the water balance B Q .The regression offset B SO 4 and slope S SO 4 of the Q-SO 4 relationship produced high sensitivity indices for the geogenic contribution Geo SO 4 and its variability a GEO .
Differences of T among the sites were found for the soil storage capacity (V mean,S ) that is high either for the Q-NO 3 cross-correlation L NO 3 (Swiss and Spanish sites) or the E Q (Austrian and Palestinian sites).In addition to S LF and V δ 18 O , a SE also has a high T for S MF and R Q,100 (Swiss site) and E Q (Austrian, Palestinian and Israeli 2 sites).The epikarst storage capacity V mean,E shows always a high value of the sensitivity index for L NO 3 , but only for the Swiss, Spanish and Israeli sites; the same is true for epikarst constant K mean,E (only sensitive to L NO 3 at Swiss and Israeli 1sites).The distribution of recharge dynamics a fsep has high sensitivity indices either for L NO 3 (Swiss and Israeli sites), the discharge autocorrelations R Q,100 (Spanish site) or E Q (Spanish and Israeli 2 sites).The conduit system storage constants K C show high T for S MF for the Austrian, Swiss and Spanish sites, while V crit,OF has a high sensitivity for R Q,100 (all sites except the Austrian site), for B Q (Spanish, Palestinian and Israeli sites) and L NO 3 (Swiss site).The Spanish site is the only place where a GW has a low T for S SO 4 and R Q,100 .Only at the Israeli sites it has a high value for B SO 4 , and only at the Austrian site, T is also high for S HF .

Relation between system signatures and calibrated parameters
In stage 3 of the evaluation, only parameters with a high first-order sensitivity index F (≥ 0.1) were considered and related to the system signatures they were obtained from (Fig. 5).In order to recognise a relation to their system signatures, only sets with more than three pairs of high F parameters and system signatures were included in the analysis.
From ten relationships, six showed correlation.The conduit storage constants K C are clearly correlated to the high flows S HF , the distribution of groundwater storage constants a GW to the medium flows S MF and the δ 18 O variability V δ 18 O .In addition, geogenic contributions Geo SO 4 were correlated to the offset of the Q-SO 4 relationship B SO 4 , the distribution of geogenic SO 4 contributions a Geo to the slope of the Q-SO 4 relationship S SO 4 , and the recharge area A to the water balance B Q .For S HF , S SO 4 , B SO 4 and B Q , the relations were linear (expressed by linear correlation coefficients r Lin , Fig. 5); for S MF and V δ 18 O they were non-linear (expressed by the Spearman rank coefficient of correlation r SR , Fig. 5).

Model performance and process sensitivity
The test of performance in stage 1 of our analysis showed whether the VarKarst model is flexible enough to reproduce the observations expressed by the different system signatures at the different sites (Wagener et al., 2001).Except for the δ 18 O variability V δ 18 O , the model performed well (Table 6).
The deficiencies for V δ 18 O occurred at the site with a very strong dampening effect of the isotopic signal of the rain (Israeli 2 site, Table 5).The discrepancy between observed and modelled V δ 18 O may result from differences in the temporal resolution of observation and simulations.While the model provides daily values, the observations for the two sites are in a 2-weekly or even larger resolution (Table 2).Another reason could be the timing of sampling.Due to the coarse sampling resolution, parts of the isotopic variability caused by short rainfall discharge events might have been lost.Hence, errors in the representation of δ 18 O information may be the most probable cause for the model failure.Therefore, instead of discarding the whole model (as in Hartmann et al., 2013b), only V δ 18 O was not considered in the further analysis of the Israeli 2 study site.For the Swiss site, ∼ 22.6 % of deviation was regarded as still acceptable.Similar to van Werkhoven et al. (2009), the results of stage 2 of the evaluation showed that a large number of the system signatures provided information about the same processes for all sites: the fast groundwater dynamics (K C and V crit,OF ) were described by the high flows and the memory effects of the karst systems (S HF and R Q,100 ).Low flows, the interplay of discharge and NO 3 , as well as the dampening of the atmospheric δ 18 O signal (S LF , L NO 3 and V δ 18 O ), provide information about the distributions of epikarst and soil storage capacities (a SE ).The same signatures plus the medium flows (S MF ) describe the distribution of groundwater storage constants (a GW ).The recharge area (A) is described by the water balance (B Q ) and the SO 4 dissolution dynamics (Geo SO 4 and a Geo ) by the interplay of discharge and SO 4 (S SO 4 and B SO 4 ).
The selection of the system signatures was done by combining well-known metrics to describe hydrological (karst) systems and by trial and error.The final large number of sensitive parameters at all sites at stage 2 indicated that the ten signatures that we elected for this study provide in total enough information to describe the different karst systems, even though a wide range of other system signatures would have been available (Yadav et al., 2007).The simultaneous sensitivity of the same parameters to different signatures also indicates that there is an overlapping of information content (i.e.some of the signatures are correlated).In addition, testing the stability of the system signatures by a split-sample test showed that L NO 3 and, for some sites, also the slopes of the flow duration curves varied when only a part of the available data was used for their calculation (see Supplement).Reasons for that are the coverage of extraordinary wet years and the resolution and length of the observation time series.

System signatures and dominant processes
Since in stage 2 only the total sensitivity of parameters T was considered, parameter interactions (Saltelli et al., 2008) prohibit a direct quantification of dominant processes by parameter values.However, it allowed determining the critical processes for the different system signatures and how they change among the sites.Stage 1, test of performance, showed that the model is able to reproduce almost all of the observed signatures (Tables 5 and 6).Using the calibrated parameters, it is therefore possible to distinguish responsive from less dynamic systems.Combining this information with the relations between critical processes and system signatures from stage 2, we can identify the dominant processes at the different systems and attribute them to the signatures.
Steep slopes at the high flows S HF of the flow duration curves are found for the Austrian, Swiss, Palestine and Israeli 2 sites (Fig. 2a); at the Spanish and the Israeli 1 sites S HF flows are rather low.The high sensitivity of the K C for all sites (Fig. 4) indicates that the conduit system dynamics are the dominant process controlling the high flow behaviour of all springs.Even though the activation threshold for overflow springs V crit,OF is sensitive for all springs, it does not mean that there are overflow springs at all the systems.Field studies showed that overflow springs can be found at the Spanish site (Barberá and Andreo, 2011) and Austrian site (Kralik et al., 2009).For the Spanish site, the low S HF would indicate the dominance of the overflow spring on the high flow behaviour, but for the Austrian site the overflow behaviour is less pronounced.Indeed preceding studies (Hartmann et al., 2012a) showed that other processes have also a significant control on its discharge behaviour.
For all sites, the distribution of groundwater storage constants a GW is sensitive for medium flows S MF (Fig. 4).That means that, under medium conditions, the hydrodynamic behaviour of all springs is controlled by the variability of aquifer characteristics, such as aquifer geometry, topography and hydraulic properties.For the Austrian, Swiss and Spanish sites, S MF is also sensitive on K C indicating that fast flow processes also contribute to their median flow behaviour.The slope of medium flows S MF is steepest for the Spanish site indicating high hydraulic conductivities.These may be due to the high degree of karstification and the large number of wells in the surroundings that facilitated the drainage of groundwater and resulted in a seemingly more dynamic behaviour of the spring (Barberá and Andreo, 2011).
The slope of the low flows S LF is steepest for Austria and Palestine (Fig. 2a).Again, the distribution of groundwater storage constants a GW controls the flow behaviour.At the Palestine sites an inclination of the bedding plane towards the spring may be the reason for the fast drainage (Ghanem, 1999), while at the Austrian site a dipping of the bedrock stratigraphy towards south-east in the deeper parts of the aquifer results in a preferential flow towards south-east, away from the spring outlet during low flow conditions (Kralik and Keimel, 2003).In terms of autocorrelation of discharges R Q,100 , the Spanish spring shows the lowest memory.Parameter sensitivity indicates that the reason for that is the above-mentioned dominance of fast groundwater flow processes (K C and V crit,OF ).The Israeli site has the springs with the largest memory.For them, the distribution of groundwater storage constants a GW (i.e. also the contribution of slowly reacting parts of the aquifer) is important.Since the numbers we obtained for R Q,100 are also influenced by the climatic variability, they cannot be used directly to understand our karst systems (Jeannin and Sauter, 1998).However, their relation to model processes can be used to infer about the system dynamics that control R Q,100 .
Except for the Austrian site, all systems that were not discarded in stage 1 of the evaluation (test of performance) show a rather strong dampening of the climatic isotope signal V δ 18 O .This is contradictory to the results obtained by the signatures concerning the discharge time series (S HF , S MF , S LF and R Q,100 ).A reason for that may be found in the resolution of the different sources of information, or in the exchange between mobile and stagnant groundwater (e.g.Małoszewski and Zuber, 1985) that is not considered by the model.This may explain why the model failed for the Israeli 2 site at stage 1 of the calibration.NO 3 observations have a higher resolution than δ 18 O.In addition, the model did not show any problems reproducing the interplay of discharge and NO 3 concentration L NO 3 during the test of performance in stage 1.L NO 3 shows rather short lag times between discharge and NO 3 peaks for the Austrian sites, while all the other systems react much more slowly.For all sites, sensitivity indicates that all processes from the surface to the groundwater are relevant for the NO 3 transport through the karst systems.This is no surprise, since NO 3 originates from the surface, either by natural deposition or by anthropogenic origin, and travels through the whole karst system (Perrin et al., 2007).Even though much of the system showed large values for L NO 3 , fast groundwater flow remains a dominant process for its transport, which was also shown by Mahler and Garner (2009).In addition to the groundwater dynamics, L NO 3 clearly shows the importance of soil and epikarst processes (Fig. 4), which was already stated by preceding field studies (e.g.Aquilina et al., 2006;Williams, 1983).However, since the test of stability of the L NO 3 signature showed that it is the most instable signature, these indications have to be considered with care (see Supplement).
Steep slopes of the Q-SO 4 relation S SO 4 are most pronounced at the Israeli 2 site but also abundant at the Spanish site.At both sites this goes along with a higher offset of the Q-SO 4 relation B SO 4 compared to the other sites.Field studies showed that large sources of SO 4 are abundant at the Israeli 2 site (Brielmann, 2008) and the Spanish site (Barberá and Andreo, 2011).Evaporites are a common source of SO 4 in karst systems, and are mostly dissolved from the lower permeability parts of the karst systems (Ford and Williams, 2007).For that reason, in addition to the parameters that control the dissolution of SO 4 in the model (Geo SO 4 and a Geo ), the distribution of groundwater storage constants a GW has also an important impact on the SO 4 dynamics.
The test of performance was also successful for the total water balances B Q , and its values coincide well with annual water balances provided in Table 1.The sensitivity analysis showed that for all sites the most important control on B Q was the recharge area A (Fig. 4), which was already shown in preceding studies (Hartmann et al., 2013a).In addition, the sensitivity analysis indicated that the abundance of overflow springs has an influence on water balance, too (which is plausible since the discharge of overflow springs is not included in the observations).At the Palestinian site, also soil properties (V mean,S and a SE ) have an impact B Q indicating that actual evaporation, which is controlled by the soil depth, plays another important role for water balance.
Finally stream flow elasticity E Q is > 1 for Spain and Israel 2, while it is < 1 for Austria and Israel 1.A streamflow elasticity E Q larger than 1 at the Spanish and Israeli 2 sites indicates a high climate sensitivity (Sankarasubramanian et al., 2001).This agrees with the findings of Hartmann et al. (2013a), who found that at the Spanish site observed annual discharge strongly depends on climatic conditions.At the Israeli sites, several studies (e.g.Hartmann et al., 2012a;Rimmer and Salingar, 2006) showed that the Israeli 2 system was more responsive to climatic variability than the Israeli 1 system.E Q can be regarded as an indicator for the stability of flow due to changes in precipitation (Sawicz et al., 2011).Sensitivity analysis shows that at our karst systems this stability is controlled by the soil and epikarst (V mean,S , V mean,E , a SE or K mean,E ) dynamics or by the recharge dynamics (a fsep ).

Calibrated parameters versus system properties and model realism
Sensitive parameters that individually contribute to the model output variance could be identified using the firstorder sensitivity index of the model parameters F .Assuming that the model structure represents the real system and parameter interactions T − F are small, these parameters can serve as proxies of system properties (Carrillo et al., 2011).The parameters identified that way were overflow spring threshold V crit,OF , the conduit storage constant K C , the distribution of groundwater storage constants a GW , the geogenic contribution Geo SO 4 , their variability a Geo and the recharge area A (Fig. 5).The most obvious among them are the relationships that concern the water and solute balances: B Q and A (r Lin = 0.99 in a log-log scale), the slope of the Q-SO 4 relationship S SO 4 and a Geo (r Lin = 0.996), and the offset of the Q-SO 4 relationship B SO 4 and Geo SO 4 (r Lin = 0.97).B Q and A indicate that disregarding effects of evaporation, the recharge area of all considered systems can be derived directly from their mean annual discharge.For the SO 4 balance the correlation indicates that for all considered systems, SO 4 mass balance is not dependent on atmospheric input of SO 4 .Thus, B SO 4 and S SO 4 give an estimate of if and how water gets in contact with evaporites in the system (see Ford and Williams, 2007, dissolution of gypsum and anhydrite).Relations between high flows S HF and K C (r Lin = 0.9), medium flows S MF and a GW (r SR = 1.0), and the δ 18 O variability V δ 18 O and a GW (r SR = 0.8) describe the slow and fast discharge dynamics.Recession analysis is often used to derive the parameters of the slow groundwater system of hydrological models (e.g.Fleury et al., 2007).Our results indicate that the slopes of the flow duration curve during high flows may be used in a same way for the peak flows.Large values of a GW result in very slow groundwater flow dynamics (see Appendix).Accordingly, the established relations show that high values of a GW go along with flat slopes of the flow duration curves for medium flows and dampened isotopic signals.Kovacs et al. (2005) showed that storage constants such as K C can be related to hydraulic properties of a karst system.With our new findings, not only mean hydraulic conductivities, but also their distribution a GW may be approximated when S HF and S MF or V δ 18 O are known.
Having indications that correlations between system signatures and system properties exist, the question arises whether it is possible to transfer system signatures to ungauged karst systems by climatic and topographic information.When we compared the values of the ten system signatures with the climatic and topography descriptors presented in Table 1, we could not find any significant correlations (see Table 7 and Supplement).Hence, the available information was not complete enough to allow a regionalisation of system signatures and, therefore, a model application in ungauged karst basins.However, high uncertainty goes also along with direct regionalisation of model parameters (Wagener and Wheater, 2006) and the measurement of karst system properties in the field (Goldscheider and Drew, 2007).The relations between model parameters and system signatures found in this study encourage following the idea of the approach further to regionalise system signatures instead of model parameters.That way ungauged karst systems can be described without the uncertainties that go along in transferring the parameters of a chosen model.System signatures are easily available, since they often include commonly available data like flow chart characteristics or regionalised flood or low flow indices (Zhang et al., 2008).There are also more possibilities to define new system signatures: for instance, Long and Mahler (2013) suggest using metrics describing the shape of impulse-response functions to characterise and distinguish karst systems.Yadav et al. (2007) propose more than 20 descriptors derived from topography, climate observations and landscape properties which can be used to regionalise system signatures.Unfortunately, due to the usually unknown size and location of the subsurface catchment of karst systems, most of these metrics could not be used in this study.But especially for karst systems, information about general geological properties and degree of karstification may also be quantified, using for instance descriptors of initial porosity, fractures and age of the karst system that can be derived from modelling studies (e.g.Bloomfield et al., 2005;Hubinger and Birk, 2011) or age dating of stalactites (e.g.Vaks et al., 2003;White, 2007).
A main assumption of this approach was an adequate system representation by the model.For this assumption to be correct, certain flexibility in the model is necessary given that the considered karst systems vary in scales, climates, surface and subsurface properties (Table 1).Hartmann et al. (2013a) showed that the VarKarst model includes such flexibility enabling it to consider different aspects of the karst systems' behaviour.Unlike Carrillo et al. (2011) we use automatic calibration and sensitivity analysis on each of the ten signatures and use only parameters with a high sensitivity for interpretation.That way, parameter identification is more objective (Hartmann et al., 2012a).In addition, by a large number of hydrochemical signatures we included more information to improve system understanding (Bishop et al., 2004;Weiler and McDonnell, 2005).

Conclusions
The main scope of this study was to identify differences between dominant processes and system properties for five karst systems of varying size and in different climatic regions in Europe and the Middle East.Using a set of ten hydrodynamic and hydrochemical system signatures in a processbased karst model, their importance for relating them with karst system properties was explored.During a stepwise analysis the models were calibrated and the parameter sensitivity concerning the signatures was investigated.It was possible to show that sensitivity analysis can be used to identify and distinguish processes for different karst systems.Moreover, relations were found between signatures concerning water and solute balances, and model parameters that express the recharge area and geogenic contributions of hydrochemical compounds (Table 7).It was possible to relate hydrodynamic and hydrochemical karst system signatures to model parameters that represent different properties of the karst systems.The inclusion of hydrochemical information was crucial during all stages of the analysis.While hydrodynamic signatures majorly provided information about the groundwater dynamics, NO 3 described the behaviour of soil and epikarst processes and SO 4 further contributed to the characterisation of the groundwater flow dynamics.Similarly, hydrochemical information contributed to establish relations between climatic and topographic descriptors and system signatures (Table 7).
The stepwise analysis with a process-based karst model including automatic calibration and Sobol sensitivity analysis offered new directions in comparing the process dynamics and properties of karst systems.It allowed (1) investigating the information content of the different hydrodynamic and hydrochemical karst system signatures, (2) explaining the dominant processes that are responsible for the different system signatures at the different karst systems, and (3) establishing relations between system signatures and system properties.Even though the number of these relations was still too small to facilitate a regionalisation of system signatures and model parameters, this study encourages investing more time in the exploration of alternative ways to define system signatures and descriptors of the karst systems.Hereby, descriptors of general geological properties and degree of karstification (e.g.Bloomfield et al., 2005;Hubinger and Birk, 2011)

Fig. 1 .
Fig. 1.(a) Location and (b)-(f) maps of the study sites in Europe and the Middle East.

Fig. 2 .
Fig. 2. Different hydrological and hydrochemical aspects of the system behaviour of the study sites: (a) flow duration curves, (b) autocorrelation of discharge, (c) cross-correlation of discharge and NO 3 , (d) correlation of discharge and SO 4 , (e) ratios of variability (standard deviations) in observed δ 18 O in discharge (Q) and precipitation (P ), and (f ) annual amounts of discharge.

Fig. 3 .
Fig. 3. Sketch of model structure adopted from Hartmann et al. (2013a) illustrating the relevant processes and connections (light blue indicates the unsaturated part of the groundwater aquifer).
where F,I represents the contribution of an individual parameter i to the model output variance D. Due to their definitions in Eqs.(2) and (3), T and F range from 0 to 1.The difference between T and F represents the parameter interactions.Similar toHartmann et al. (2013b), we consider Hydrol.Earth Syst.Sci., 17, 3305-3321, 2013 www.hydrol-earth-syst-sci.net/17/3305/2013/

Fig. 4 .
Fig. 4.T of model parameters concerning the system signatures for all study sites (all parameters with T < 0.2 are considered not sensitive and have been removed).

Fig. 5 .
Fig. 5. Relationships between calibrated parameters with F ≥ 0.1 and system signatures; dot sizes indicate parameter interactions T − F : the smaller the dot, the larger the interactions, the higher the uncertainty of the parameter location (r Lin and r SR are linear correlation coefficient and the Spearman rank correlation coefficient, respectively; p values are only calculated for the linear correlations).

Table 1 .
General characteristics of the study sites.
* Combined discharge of both springs.

Table 2 .
Available data for the analysis.

Table 4 .
Parameters of the VarKarst model, their descriptions, units and calibration ranges for the different study sites.

Table 5 .
Observed karst system signatures obtained by the signature equations provided in Table3.

Table 6 .
Agreement of modelled with observed signatures (Table5) in model evaluation stage 1.

Hartmann et al.: Process-based karst modelling to relate hydrodynamic and hydrochemical characteristics 3317Table 7 .
provide a very promising direction.Summary of correlations between system signatures and model parameters, and between system signatures and climatic and topographic descriptors.