Articles | Volume 11, issue 12
Development and technical paper
06 Dec 2018
Development and technical paper |  | 06 Dec 2018

V2Karst V1.1: a parsimonious large-scale integrated vegetation–recharge model to simulate the impact of climate and land cover change in karst regions

Fanny Sarrazin, Andreas Hartmann, Francesca Pianosi, Rafael Rosolem, and Thorsten Wagener

Karst aquifers are an important source of drinking water in many regions of the world. Karst areas are highly permeable and produce large amounts of groundwater recharge, while surface runoff is often negligible. As a result, recharge in these systems may have a different sensitivity to climate and land cover changes than in other less permeable systems. However, little is known about the combined impact of climate and land cover changes in karst areas at large scales. In particular, the representation of land cover, and its controls on evapotranspiration, has been very limited in previous karst hydrological models. In this study, we address this gap (1) by introducing the first large-scale hydrological model including an explicit representation of both karst and land cover properties, and (2) by providing an in-depth analysis of the model's recharge production behaviour. To achieve these aims, we replace the empirical approach to evapotranspiration estimation of a previous large-scale karst recharge model (VarKarst) with an explicit, mechanistic and parsimonious approach in the new model (V2Karst V1.1). We demonstrate the plausibility of V2Karst simulations at four carbonate rock FLUXNET sites by assessing the model's ability to reproduce observed evapotranspiration and soil moisture patterns and by showing that the controlling modelled processes are in line with expectations. Additional virtual experiments with synthetic input data systematically explore the sensitivities of recharge to precipitation characteristics (overall amount and temporal distribution) and land cover properties. This approach confirms that these sensitivities agree with expectations and provides first insights into the potential impacts of future change. V2Karst is the first model that enables the study of the joint impacts of large-scale land cover and climate changes on groundwater recharge in karst regions.

1 Introduction

Carbonate rocks, from which karst systems typically develop, are estimated to cover 10 %–15 % of the world (Ford and Williams, 2007). More importantly, karst aquifers are a considerable source of drinking water for almost a quarter of the world's population (Ford and Williams, 2007) and play a critical role in sustaining food production because most karst areas contain some form of agricultural activity (Coxon, 2011). In Europe, carbonate rock areas cover 14 %–29 % of the land area, and some European countries such as Austria and Slovenia derive up to 50 % of their total water supply from karst aquifers (Chen et al., 2017; COST, 1995).

Karst systems are characterised by a high spatial variability of bedrock and soil permeability due to the presence of preferential flow pathways (Hartmann et al., 2014). The soluble carbonate bedrock is structured by large dissolution fissures or conduits (Williams, 1983, 2008) and the typically clayey soil often contains cracks (Blume et al., 2010; Lu et al., 2016) where infiltrating water concentrates. Therefore, a large part of the groundwater recharge occurs as concentrated and fast flow in large apertures while the other part moves as diffuse and slow flow through the matrix (Hartmann and Baker, 2017). Preferential flow pathways are particularly developed in karst, but they are also widely found in many other systems, due to root and organism activities, discontinuous subsurface layers, surface depressions, soil desiccation, tectonic processes and physical and chemical weathering (Beven and Germann, 2013; Hendrickx and Flury, 2001; Uhlenbrook, 2006)

Preferential infiltration is typically triggered when thresholds in rain intensity and soil moisture levels are exceeded (Rahman and Rosolem, 2017; Tritz et al., 2011). When activated, preferential infiltration pathways may enhance groundwater recharge while limiting surface runoff (e.g. Bargués Tobella et al., 2014). In karst systems, permeability is often so high that surface runoff is negligible, and virtually all precipitation infiltrates (Contreras et al., 2008; Fleury et al., 2007; Hartmann et al., 2014). Furthermore, preferential infiltration pathways can affect the temporal dynamics of recharge. For instance, Cuthbert et al. (2013) showed that macro-pores in the soil can generate quick responses in the water table, and Arbel et al. (2010) observed that dripping rates in a karst cave can fluctuate following inter-seasonal and intra-seasonal precipitation variations.

Changes in weather patterns, specifically in the intensity and frequency of precipitation, may alter the activation of preferential flow pathways. In karst areas, several modelling studies have shown that groundwater recharge (Hartmann et al., 2012a; Loáiciga et al., 2000), spring discharge (Hao et al., 2006) and streamflow (Samuels et al., 2010) respond to changes in climate. However, to our knowledge, only the study by Hartmann et al. (2017) quantified the sensitivity of simulated karst groundwater recharge to specific precipitation characteristics, namely the mean precipitation and the intensity of daily heavy precipitation events. This modelling study was conducted across carbonate rock areas in Europe, the Middle East and northern Africa for different climate change projections. The study results suggested that, due to the presence of preferential flow pathways, recharge in karst systems tends to exhibit a higher sensitivity to changes in mean precipitation and in the intensity of heavy precipitation events in dry climates, in addition to a lower sensitivity in wet climates, compared to non-karst systems. Hartmann et al. (2017) also showed that the intensity of heavy precipitation can have both a positive and a negative effect on recharge for both karst and non-karst systems. However, other observational studies in non-karst areas associate increases in extreme rainfall with a higher recharge amount, e.g. in a semi-arid tropical region (Taylor et al., 2013) and in a seasonally humid tropical region (Owor et al., 2009). The discrepancy between these results might be explained by the fact that Hartmann et al. (2017) only tested recharge sensitivity against precipitation properties, while ignoring their interactions with other meteorological variables such as temperature or humidity.

In addition to climate change, land cover/use change is also expected to have a major impact on hydrological processes in the future (DeFries and Eshleman, 2004; Vörösmarty, 2002). Changes in land cover/use can impact the partitioning between green water (water that can be lost through evapotranspiration) and blue water (water potentially available for human activities, namely groundwater recharge and runoff) (Falkenmark and Rockström, 2006). Green water tends to be higher for forested areas than for shorter vegetation (e.g. Brown et al., 2005), which has also been confirmed for karst areas (Ford and Williams, 2007; Williams, 1993). Significant land cover/use changes are expected to occur in the future, including in European and Mediterranean karst areas. These changes will partly be caused by modifications in socio-economic and technological factors, such as changes to food and wood demand or changes in agricultural management practices that could enhance agricultural yields (see e.g. Holman et al., 2017, for a European assessment; Hurtt et al., 2011, for a global assessment). Future vegetation will also be impacted by other changing environmental conditions such as increases in atmospheric CO2 leading to differences in plant behaviour and leaf area index (LAI), as well as to differences in climate and weather patterns, which in turn change the frequencies of natural disturbances such as wildfires, storms or bark beetle infestations (Seidl et al., 2014; Zhu et al., 2016).

The above review of the literature reveals that changes in climate characteristics (e.g. precipitation intensity and frequency) and in land cover properties can be expected to have significant impacts on the hydrology of karst regions. However, the joint effect of changes in these boundary conditions has not been studied systematically, and such assessment is needed at large scales to inform water resources management plans (Archfield et al., 2015). In this study we introduce a novel large-scale hydrological model that includes an explicit representation of both karst and vegetation properties; furthermore, we systematically explore the sensitivity of this model's groundwater recharge predictions to meteorology and vegetation properties. Our model extends an existing karst hydrology model, called VarKarst, which was recently developed for large-scale applications and was tested over European and Mediterranean carbonate rock areas (Hartmann et al., 2015). However, VarKarst only contains a simplistic and empirical representation of evapotranspiration processes that does not explicitly include land cover properties, which prevented its application in land cover change impact studies.

Our study has two objectives. First, we add an explicit representation of land cover properties into VarKarst by improving its evapotranspiration (ET) estimation. While we seek to keep the model structure parsimonious, we want the new version of the model (V2Karst V1.1) to be appropriate for assessing the hydrological impacts of combined land cover and climate change at a large scale. The model should ultimately simulate groundwater recharge, which is the amount of renewable groundwater available for human consumption and ecosystems, at both seasonal and annual timescales (e.g. Döll and Fiedler, 2008; Scanlon et al., 2006; Wada et al., 2012). Second, we test the plausibility of the V2Karst model behaviour by comparing its predictions against observations available at carbonate rock FLUXNET sites, and by analysing the sensitivity of simulated recharge using both measured and synthetic data to force the model. In particular, the use of synthetic data in virtual experiments allows us to control variations in climate and vegetation inputs, so that we can isolate their individual and combined effects on simulated recharge and overcome some of the issues found by Hartmann et al. (2017) as discussed above.

2 New version of VarKarst with explicit representation of land cover properties (V2Karst)

2.1 Rationale for our approach to land cover representation in VarKarst

The new version of the VarKarst model should be appropriate to assess the impact of climate and land cover change on karst groundwater recharge. It should also consider the range of challenges related to modelling ET at large scales, namely a lack of ET observations to compare with model predictions and a lack of observations of vegetation properties (e.g. rooting depth, stomatal resistance and canopy interception storage capacity) to constrain the model parameters (see Sect. S1 in our Supplement). We define the three following criteria to develop an enhanced version of the VarKarst model with explicit representation of land cover properties:

  1. The new model should assess the three main ET components (bare soil evaporation in presence of sparse canopy, transpiration and evaporation from canopy interception) separately and explicitly. In fact, these fluxes exhibit different dynamics and sensitivity to environmental conditions; therefore, they are likely to respond differently to climate and land cover changes (Gerrits, 2010; Maxwell and Condon, 2016; Savenije, 2004; Wang and Dickinson, 2012).

  2. The model should use a Penman–Monteith formulation to estimate the potential evapotranspiration (PET; Monteith, 1965), so that it can separate the effects of climate and land cover on each of the ET components. In fact, empirical PET formulations such as the Priestley–Taylor equation (Priestley and Taylor, 1972) do not allow for such separations as they do not explicitly include land cover properties.

  3. All processes should be represented parsimoniously in accordance with the modelling philosophy underpinning the first version of VarKarst (Hartmann et al., 2015). This criterion aims to avoid over-parameterisation given the limited amount of available information to constrain and test model simulations – specifically at large scales (Abramowitz et al., 2008; Beven and Cloke, 2012; Haughton et al., 2018; Hong et al., 2017; IPCC, 2013, chap. 9, 790–791; Young et al., 1996). In particular, parameters that account for the physical properties of the system (e.g. soil and vegetation properties) are commonly taken from look-up tables; however, the physical meaning of these parameters has come into question, and they may actually not be commensurate with field measurements as discussed in Hogue et al. (2006) and Rosero et al. (2010). Therefore, it has been suggested that even these “physical” parameters should be calibrated in order to optimise the model performance (Chaney et al., 2016; Rosolem et al., 2013). Importantly, parsimony also limits the computational time for model simulations and allows for the assessment of the impact of modelling choices and the uncertainty and sensitivity of model output using Monte Carlo simulation (Hong et al., 2017; Young et al., 1996).

With respect to previous modelling studies of karst systems, to our knowledge, only four have used models that explicitly include land cover properties, all of which were applied at the local scale with detailed on-site information. Three of these studies (Canora et al., 2008; Doummar et al., 2012; and Zhang et al., 2011) used generic hydrological models that were not specifically developed for karst areas but included enough flexibility in their spatially distributed parameters to represent the variability in soil and bedrock properties. The large number of parameters in these models hampers their application at large scales and does not comply with criterion 3 (parsimony). The fourth study (Sauter, 1992) used a lumped model that is much more parsimonious but does not represent soil evaporation and uses empirical PET equations, which does not allow for the separation of the effect of climate and land cover (criteria 1 and 2).

As for large-scale models, we can identify two main types: hydrological models, which focus on the assessment of hydrological fluxes, and land surface models, which also evaluate many other fluxes such as sensible heat, latent heat, ground heat flux, radiation and carbon fluxes (for a review, see Bierkens, 2015). Land surface models do not usually comply with criterion 3 because they have many parameters, including a number of empirical parameters that are difficult to constrain (Mendoza et al., 2015). Moreover, it has been shown that land surface models could be simplified when the objective is to assess hydrological fluxes only. For instance, Cuntz et al. (2016) demonstrated that a large number of parameters of the Noah land surface model are non-influential or have very little influence on simulated runoff. In contrast, hydrological models focus on the representation of hydrological processes and include far fewer parameters. However, our literature review (summarised in Tables A1–A3) showed that we cannot directly adopt any of their ET representations into VarKarst. In fact, as shown in Tables A1–A3, the most parsimonious models (WBM, WaterGAP and mHM) neglect some ET components and/or use empirical PET equations, which contradicts criteria 1 and 2, while models that comply with criteria 1 and 2 (PCR-GLOBWB, VIC and the model of Kergoat, 1998) use heavily parameterised schemes, such as a Jarvis type parameterisation of surface resistance (Jarvis, 1976; Stewart, 1988); therefore, these models do not satisfy criterion 3 (parsimony). Moreover, we found that large-scale models include empirical schemes with no clear origin, such as the reference crop formulation used in the PCR-GLOBWB model for PET calculation or the interception model used in the LPJ dynamic global vegetation model and in the model of Kergoat (1998). Importantly, our review revealed the tremendous variability of approaches used in large-scale models to represent ET processes. A detailed list of all parameters involved in the representation of ET in the models in Tables A1–A3 can be found in Sect. S2 of our Supplement. Consequently, no clear indication has emerged regarding a “best way” to parameterise the different ET processes at large scales, which leaves us with a large range of different formulations to choose from to implement an explicit representation of land cover processes into VarKarst.

Figure 1Schematic representation of (a) the VarKarst model (Hartmann et al., 2015) and (b) the new version of the model V2Karst using six vertical compartments. Model parameters are in green (see their definition in Table 1), inputs are in blue (P – precipitation, Rn – net radiation, T – temperature, RH – relative humidity, WS – wind speed), model fluxes are in black (Epot – potential total evapotranspiration, Tpot – potential transpiration, Ecpot – potential evaporation from canopy interception, Espot – potential soil evaporation, Eact – total actual ET, Tact – actual transpiration, Ecact – actual evaporation from canopy interception, Esact – actual bare soil evaporation, Tf – throughfall, Qlat,23 – lateral flow from the second to the third compartment, Qsurf – surface runoff and Qepi – recharge) and state variables are in red.


Table 1Description of V2Karst parameters, unconstrained ranges used in the application at the four FLUXNET sites to capture the variability across soil, epikarst and vegetation types, and the category of the parameters (which indicates whether the parameters depend on soil, epikarst or vegetation properties). Parameters a, Vsoil, Vepi and Kepi were already present in the previous version of the model (VarKarst). More information on how the ranges were determined is provided in Sect. S3 of our Supplement.

Download Print Version | Download XLSX

2.2 Previous representation of evapotranspiration (ET) processes in VarKarst

VarKarst (Hartmann et al., 2015) is currently the only karst recharge model developed for large-scale applications. It is a conceptual semi-distributed model that simulates karst potential recharge (Fig. 1a). The representation of karst processes in VarKarst is based on a previous karst model developed for applications at the local scale introduced in Hartmann et al. (2012b). VarKarst includes two horizontal subsurface layers, a top layer called “soil” and a deeper layer called “epikarst”. The soil layer corresponds to the layer from which ET can occur. The epikarst layer corresponds to the uppermost layer of weathered carbonate rocks where it is assumed that water cannot be lost through ET. VarKarst represents karst processes because for each model grid cell, the water balance is evaluated separately over a number of vertical compartments with varying soil and epikarst properties. Recharge is the flux leaving the bottom of the epikarst and is produced in saturated or unsaturated compartments where there is moisture stored in the epikarst. Recharge can occur at faster rates in deeper compared to shallower compartments. Additionally, lateral flow concentrates the infiltrating water from saturated to unsaturated compartments. Conceptually, the direct contribution of precipitation to infiltration and recharge can be associated with the diffuse fraction, while the contribution of lateral flow can be associated with the concentrated fraction. The model has been shown to be more appropriate for applications over karst areas than other large-scale hydrological models that do not represent karst processes (Hartmann et al., 2017). The ET component of the VarKarst model is very simple and does not include explicit representation of land cover properties. ET is lumped in the soil layer, is estimated from PET and is reduced by a water stress factor, which is estimated as a linear function of soil moisture. The PET rate is calculated with the empirical Priestley–Taylor equation (Priestley and Taylor, 1972) using a spatially and temporally uniform value of the empirical coefficient. This approach does not allow for the separation of the effects of climate and land cover, as the empirical coefficient reflects both climate and vegetation characteristics simultaneously. Therefore, the ET component of VarKarst needs to be modified if the model is to be used for large-scale land cover change impact assessment.

2.3 V2Karst: the new version of VarKarst for integrated vegetation–recharge simulations over karst areas

In this section, we propose a new version of the VarKarst model, called V2Karst (Sarrazin et al., 2018; Fig. 1b). In accordance with criteria 1 and 2 defined in Sect. 2.1, the new V2Karst model includes a physically based PET equation, separates the evapotranspiration flux into three components (transpiration, bare soil evaporation and evaporation from canopy interception) and comprises three soil layers. In agreement with criteria 3 of Sect. 2.1 (parsimony), we sought to parsimoniously represent the different ET processes in VarKarst. In fact, V2Karst uses 12 parameters to represent ET and vegetation seasonality, namely 11 newly introduced parameters that replace the Priestley–Taylor empirical coefficient α used in Varkarst and the soil water capacity parameter Vsoi already present in VarKarst (model parameters are described in Table 1 and Fig. 1). This is less than other existing large-scale models that use the Penman–Monteith equation and separate the three ET components, as these models have over 15 parameters in their ET component (PCR-GLOBWB, VIC and the model of Kergoat, 1998, shown in Tables A1–A3). We assumed homogeneous above ground vegetation properties across model compartments.

The new model is forced by time series of precipitation (P), air temperature (T) and net radiation (Rn) in the same manner as VarKarst. Additionally, time series of relative humidity (RH) and wind speed (WS) are now needed for PET calculation. The V2Karst model can be run at both daily and sub-daily time steps. In this study, we present simulation results obtained using a daily time step, while results from hourly simulations are reported in Sect. S8 of our Supplement. We did not find significant differences between daily and hourly simulations for assessing recharge at monthly and annual timescales, which is the focus of our study. Therefore, it is reasonable to apply the model at a daily time step, which significantly reduces the computational requirements.

2.3.1 Definition of soil and epikarst properties in V2Karst

The computation of the water storage capacity of the entire soil column VS,i (mm) and of the epikarst VE,i (mm), and the epikarst outflow coefficient KE,i (days) for the ith model compartment is carried out in the same fashion as in VarKarst:


where Vmax,S (mm) is the maximum soil storage capacity over all model compartments; Vmax,E (mm) is the maximum epikarst storage capacity; Kmax,E (days) is the maximum outflow coefficient; nc (–) is the number of model compartments, which is set to 15 following (Hartmann et al., 2013, 2015); and a (–) is the spatial variability coefficient. A previous study showed that VS,i, VE,i and KE,i can be determined using the same distribution coefficient a (Hartmann et al., 2013). In V2Karst, Vmax,S, Vmax,E and KE,i are computed as a function of the average properties of the cell using the following formulas:


where Vsoi (mm) is the mean soil storage capacity, Vepi (mm) is the mean epikarst storage capacity and Kepi (mm) is the mean epikarst outflow coefficient. We note that the definition of the three parameters Vsoi, Vepi and Kepi is revised compared to VarKarst.

As in VarKarst, we neglect ET from the epikarst. Several studies have shown that in the presence of shallow soil and a dry climate, plants can take up water in the weathered bedrock where soil pockets can sustain root development (Schwinning, 2010). However, given the uncertainty in soil depth for large-scale applications, V2Kast does not allow ET from the epikarst to avoid over-parameterisation. Therefore, the V2Karst soil layer must be interpreted as a conceptual layer that does not exactly correspond to the physical soil layer (layer of loose material) but is defined as the portion of the subsurface where ET losses can occur.

In V2Karst, the soil layer is further divided into a shallow top layer from which water can be lost from both evaporation and transpiration, a second middle layer where only transpiration can occur and a third deeper layer below the root zone where transpiration can only take place when the first two layers are depleted. The maximum storage capacity of the first layer is noted as Ve (mm), and the maximum storage capacity of first and second layers combined is noted as Vr (mm), which corresponds to the maximum storage capacity of the root zone. The model assumes that Ve is smaller than Vr, which is in turn smaller than the storage capacity of the deeper model compartment VS,n.

2.3.2 Soil water balance

The soil water storage Vsoi,ij(t) (mm) in the ith compartment and the jth soil layer j=1,2,3 is updated at the end of each time step t as follows:


where Tf(t) (mm) is the throughfall i.e. the fraction of precipitation that reaches the ground (Eq. 8), Qlat,i-1i(t) (mm) is the lateral flow from the (i−1)th to the ith model compartment (Sect. 2.3.4), Esact,i(t) (mm) is the actual soil evaporation (Eq. 9), Tact,ij(t)(mm) is the actual transpiration in the jth soil layer (Eqs. 11–12), R12,i(t) (mm) is the downward flow from the first to the second soil layer, R23,i(t) (mm) is the downward flow from the second to the third soil layer and Repi,i(t) (mm) is the downward flow from the soil to the epikarst.

It is assumed that percolation from the unsaturated soil to the epikarst is negligible due to the low permeability of the soil. This assumption seems reasonable as karst soils usually have a high clay content (Blume et al., 2010; Clapp and Hornberger, 1978). However, clayey soil typically presents cracks (Lu et al., 2016); therefore, when the soil reaches saturation, preferential flow starts to occur in the soil cracks, which causes all saturation excess to quickly infiltrate to the epikarst. Just as in VarKarst, such preferential vertical flow is represented by the variable Repi,i(t) (used in Eq. 3) and is set equal to the saturation excess in the (lowest) soil layer. In V2Karst, a similar approach is also used to assess the other vertical flows from one soil layer to another (R12,i(t) and R23,i(t) in Eq. 3).

2.3.3 Evapotranspiration

We adopt the representation of sparse vegetation proposed by Bohn and Vivoni (2016) for the VIC model which is referred to as a “clumped” vegetation scheme. Each model compartment is divided into a vegetated and a non-vegetated fraction using a canopy cover fraction coefficient fc(t) (–). The uptake of soil moisture for transpiration and soil evaporation is coupled in a way that, for each model compartment, we evaluate an overall water balance over the two fractions. Using a coupled approach such as this facilitates the representation of the seasonal variations in vegetated and non-vegetated fractions compared to an uncoupled “tile” approach, in which a separate soil moisture state is represented for vegetated and bare soil fractions. Consistent with other existing large-scale models, aerodynamic interactions between both fractions are neglected to keep the number of parameters to a minimum (Table A3).

The canopy coefficient fc(t) is estimated in V2Karst using the Beer–Lambert's law as in Van Dijk and Bruijnzeel (2001) and Ruiz et al. (2010). This law was originally used to separate the fraction of incident radiation (and by extension of net radiation) absorbed by the canopy from the fraction penetrating the canopy (Kergoat, 1998; Ross, 1975; Shuttleworth and Wallace, 1985). The canopy cover fraction at time t is expressed as a function of the cell average leaf area index LAI(t) (m2 m−2) and an extinction coefficient k (–), which is understood to vary across vegetation types as it accounts for leaf architecture (Ross, 1975):

(4) f c ( t ) = 1 - e - k LAI ( t ) .

Notice that Eq. (4) allows for the description of the seasonal variations in the canopy cover fraction without introducing additional parameters in the model, given that it will simply follow the seasonal variations in the LAI.

Canopy interception

The interception storage capacity over the vegetated fraction Vcan,max(t) (mm) depends (1) on the leaf area index over the vegetated fraction, which is estimated by rescaling the cell average leaf area index LAI(t) by the vegetation cover fraction fc(t) (as in Bohn and Vivoni, 2016) and (2) on the canopy storage capacity per unit of leaf area index, denoted by Vcan, which is understood to depend on the vegetation type since it accounts for leaf architecture (Gerrits, 2010). Specifically, it is expressed as

(5) V can,max ( t ) = V can LAI ( t ) f c ( t ) .

The interception storage over the vegetated fraction Ic(t) (mm) is then updated at each time step as follows:

(6) I c ( t ) = min P ( t ) + I c ( t - 1 ) - E c act ( t ) f c ( t ) , V can,max ( t ) .

The actual evapotranspiration from canopy interception Ecact(t) (mm) is computed as

(7) Ec act ( t ) = f c ( t ) min Ec pot ( t ) , P ( t ) + I c ( t - 1 ) , V can,max ( t ) ,

where Ecpot(t) (mm) is the potential evaporation from canopy interception (Eq. 14) and P(t) (mm) is the precipitation. The factor fc(t) in Eq. (7) accounts for the fact that evaporation from canopy occurs over the vegetated fraction only. Finally, the throughfall is calculated as

(8) T f ( t ) = max P ( t ) - Ec act ( t ) - f c ( t ) I c ( t ) - I c ( t - 1 ) , 0 .

Previous studies suggest that daily simulations of interceptions can provide reasonable results (Gerrits, 2010; De Groen, 2002; Savenije, 1997). For daily simulations, the model does not account for the carry-over of the interception storage from one day to the next, which means that Ic(t) is set to zero at the end of each day and that all precipitation which is not evaporated from the interception store becomes throughfall as in Gerrits (2010), De Groen (2002) and Savenije (1997). This assumption can be justified by the fact that the interception process is highly dynamic at a sub-daily timescale, because the canopy can go through several wetting–drying cycles within a day (Gerrits, 2010). Therefore, at a daily time step, the canopy layer must be interpreted as a conceptual layer, where the storage capacity does not exactly correspond to the physical storage capacity of the canopy (i.e. the amount of water that can be held at a given time) but to the cumulative amount of water that can be held by the canopy over one day (Gerrits, 2010).

Bare soil evaporation

It is assumed that soil evaporation is a faster process than transpiration, which is consistent with general knowledge on ET processes (Wang and Dickinson, 2012). Therefore, soil moisture can be first evaporated and then transpired if some available moisture remains for plant water uptake. Soil evaporation is withdrawn for the first soil layer as a function of the potential rate and soil moisture, similar to the previous version of VarKarst:


where Espot(t) is the potential soil evaporation (Eq. 14). The factor (1−fc(t)) in Eq. (9) accounts for the fact that soil evaporation occurs from the non-vegetated fraction only; therefore, the potential rate has to be weighted by the bare soil cover fraction. The right term of the equation (Vsoi,i1(t-1)+Tf(t)) is not weighted because we assume that the soil moisture is uniform over the fractions of each model compartment (we compute a unique water balance); this means that the total moisture present in the first soil layer is available to soil evaporation because the vegetated fraction can supply moisture to the bare soil fraction.

Transpiration over the vegetated fraction

Transpiration mainly occurs in the first and second soil layers, and it switches to the third soil layer when the first two layers are depleted. The extraction of water by the roots below the root zone has been documented in studies such as Penman (1950). We account for this process by representing a soil layer below the root zone, which can provide water to the root zone through capillary rise as in the ISBA model (Boone et al., 1999). In V2Karst, the rate at which transpiration occurs in the two first soil layers Trate,i12(t) (mm) and in the third soil layer Trate,i3(t) (mm) are assessed as follows:


where Tpot(t) is the potential transpiration (Eq. 14), twet(t) (–) is the fraction of the time step with a wet canopy (Eq. 13) and fred (–) is a reduction factor which accounts for the fact that moisture below the root zone is less easily accessible to the roots than moisture in the root zone (Penman, 1950); fred is expected to vary across soil type since it is linked to the soil capability to supply water to the root zone. It is assumed that transpiration occurs in the two first soil layers when Trate,i12(t) is higher than Trate,i3(t), and that transpiration is drawn from the third soil layer otherwise. The actual transpiration in the two first soil layers Tact,i12(t) (mm) and in the third soil layer Tact,i3(t) (mm) are therefore calculated as follows:


Actual transpiration in the upper two layers Tact,i12(t) is partitioned between the two soil layers within the root zone as is used in the PCR-GLOBWB model (Van Beek, 2008). In V2Karst, the transpiration is attributed to the two first soil layers proportional to their storage content. This simple representation assumes that the roots can equally access the moisture stored in the first and second layer. Actual transpiration from the first layer Tact,i1(t) (mm) and the second layer Tact,i2(t) (mm) are computed as follows:


The fraction of the day with a wet canopy twet(t) (–) is estimated as the fraction of available energy that was used to evaporate water from the interception store similar to Kergoat  (1998):

(13) t wet ( t ) = Ec act ( t ) f c ( t ) Ec pot ( t ) .

Potential evapotranspiration

We replace the Priestley–Taylor potential evaporation equation originally used in VarKarst with the Penman–Monteith equation (Monteith, 1965), which has been shown to be applicable at both a daily and sub-daily time step (e.g. Allen et al., 2006; Pereira et al., 2015). The potential transpiration rate over the vegetated fraction of the cell Tpot(t) (mm) is estimated from the canopy aerodynamic resistance ra,can(t) (s m−1) and surface resistance rs,can(t) (s m−1), the potential evaporation from interception over the vegetated fraction of the cell Ecpot(t) (mm) is assessed assuming that the surface resistance is equal to 0 following studies such as Shuttleworth (1993), and the potential bare soil evaporation rate over the bare soil fraction of the cell Espot(t) (mm) is calculated from the soil aerodynamic resistance ra,soi(t) (s m−1) and surface resistance rs,soi (s m−1), using the following equations:


where Rn(t) (MJ m−2Δt−1) is the net radiation (Δt is the simulation time step), G(t) (MJ m−2Δt−1) is the ground heat flux, λ(t) (MJ kg−1) is the latent heat of the vaporisation of water, Δ(t) (kPa C−1) is the gradient of the saturated vapour pressure–temperature function, γ(t) (kPa C−1) is the psychrometric constant, ρa(t) (kg m−3) is the air density, cp (MJ kg−1C−1) is the specific heat of the air and is equal to 1.013×10-3 MJ kg−1C−1, es(t) (kPa) is the saturation vapour pressure, ea(t) (kPa) is the actual vapour pressure and Kt (s Δt−1) is a time conversion factor which corresponds to the number of seconds per simulation time step (equal to 86 400 s d−1 for daily simulations). In this study, we neglect ground heat flux for the daily time step, which seems to be a reasonable assumption (see e.g. Allen et al., 1998; Shuttleworth, 2012).

The aerodynamic resistances of the canopy (ra,can(t)) and the soil (ra,soi(t)), which depend on the properties of the land cover and the soil, respectively, are computed using the formulation of Allen et al. (1998). To assess ra,can(t), roughness lengths and the zero displacement plane for the canopy are estimated from the vegetation height hveg (m) (Allen et al., 1998). To calculate ra,soi(t)), the zero plane displacement height is equal to zero (d=0) and the roughness length for momentum and for heat and water vapour transfer are assumed to be equal, as in Šimůnek et al. (2009) and denoted as z0 (m).

Finally, the canopy surface resistance is computed by scaling the stomatal resistance rst (s m−1) to the canopy level using the leaf area index over the vegetated fraction (as in Eq. 5 to assess canopy interception capacity); therefore, a homogeneous response across all stomata in the canopy is assumed (Allen et al., 1998; Liang et al., 1994):

(15) r s,can ( t ) = r st LAI ( t ) f c ( t ) .

In other large-scale models, rs,can is also often expressed as a function of the LAI, which allows for the direct representation of its seasonality following the variations in the LAI.

Seasonality of vegetation

We represent the seasonality of vegetation by describing the seasonal variation of the cell average LAI. We use two parameters, the maximum LAImax (m2 m−2), which is the annual maximum value of the LAI during the growing season (assumed to be from June to August in this study) and LAImin (%), which is the percentage reduction in the LAI during the dormant season (assumed to be from December to February in this study). The monthly value of the leaf area index LAIm (m2 m−2) for the mth month is computed using a continuous, piecewise linear function of LAImax and LAImin, which allows for a smooth transition between dormant and growing seasons and is similar to the function proposed by Allen et al. (1998) to assess the seasonality in crop factors:

(16) LAI m = LAI min 100 LAI max when m = 1 , 2 , 12 LAI m = LAI min 100 LAI max 4 6 - m + LAI max 4 ( m - 2 ) when m = 3 , 4 , 5 LAI m = LAI max when m = 6 , 7 , 8 LAI m = LAI min 100 LAI max 4 ( m - 8 ) + LAI max 4 ( 12 - m ) when m = 9 , 10 , 11 .

The advantage of using this simple parameterisation is that it permits one to easily analyse the effect of vegetation seasonality by studying the sensitivity of the model predictions to the LAImin parameter, which captures the strength of the seasonal variation in LAI. Timings of the four phases of the seasonality model reported in Eq. (16) are adapted for the application at the sites used in the present study, which are all located in Europe (Sect. 3.1).

2.3.4 Water storage in the epikarst

Epikarst water storage Vepi,i(t) (mm) for the ith compartment is updated at the end of each time step t as follows:


where Qepi,i(t) (mm) is the potential recharge to the groundwater (Eq. 18), Qlat,ii+1(t) (mm) is the lateral flow from the ith to the (i+1)th model compartment and Qsurf,nc(t) (mm) is the surface runoff generated by the ncth compartment. When soil and epikarst layers are saturated, the concentration flow component of the model is activated. The ith model compartment generates lateral flow towards the (i+1)th compartment Qlat,ii+1(t) (mm) equal to its saturation excess. Lateral flow from the ncth compartment is lost from the cell as surface runoff, while the other model compartments do not produce any surface runoff. The epikarst is simulated as a linear reservoir (Rimmer and Hartmann, 2012) with the outflow coefficient KE,i (days):


Figure 2The four carbonate rock FLUXNET sites selected for the analyses. Mean annual precipitation (P) and mean annual temperature (T) were estimated over the period from 1 January 2001 to 17 December 2009 for the German site, 1 January 2006 to 31 December 2008 for the Spanish site (dry years), 1 January 2009 to 30 December 2011 for the Spanish site (wet years), 1 January 2010 to 30 December 2011 for the French 1 site and 1 April 2003 to 31 March 2009 for the French 2 site. Sources of the photos: Pinty et al. (2011) for the German site, Alcalá et al. (2011) for the Spanish site, (last access: 30 November 2018) for the French 1 site, (last access: 30 November 2018) for the French 2 site. Source of the carbonate rock and country map: Williams and Ford (2006) (country map obtained from Terraspace, Russian space agency).


3 Sites and data for model testing

3.1 Description of study sites

We test the model with plot scale measurements from sites of the FLUXNET network (Baldocchi et al., 2001). We identified four FLUXNET sites across European and Mediterranean carbonate rock areas for which sufficient data were available to force V2Karst and to test the model (see Sect. 3.2). A short summary of each site's characteristics is provided in Fig. 2, and more detailed information can be found in Table B1.

The sites have different climate and land cover properties. The first site (Hainich site, referred to as “German site”) is located in the protected Hainich National Park, Thuringia, central Germany, and is characterised by a suboceanic–submountain climate and a tall and dense deciduous broadleaf forest. The second site (Llano de los Juanes site, referred to as “Spanish site”) is located on a plateau of the Sierra de Gádor mountains, south-eastern Spain, has a semi-arid mountain Mediterranean climate and is an open shrubland. The third site (Font-Blanche site, referred to as “French 1 site”) is located in south-eastern France, has a Mediterranean climate and its land cover is medium-height mixed evergreen forest. The fourth site (Puéchabon site, referred to as “French 2 site”) is located in southern France and is characterised by a Mediterranean climate with a short evergreen broadleaf forest. Overground vegetation properties are well characterised at all sites, but subsurface properties are more uncertain. In particular, the rooting depth water capacity was only well investigated at the French 2 site. The four sites are appropriate for testing V2Karst, as they satisfy the model assumptions: a karstified or fissured and fractured bedrock, overall high infiltration capacity with limited surface runoff and high clay content in the soil (Table B1).

3.2 Data description and preprocessing

Data available at the four FLUXNET sites include measurements of precipitation, temperature, net radiation, relative humidity and wind speed to force the model. Eddy-covariance measurements of latent heat flux (density) and at the German and Spanish sites measurements of soil moisture are also available to estimate the model parameters (Sect. 4.1). Specifically, at the German site, soil moisture was measured in one vertical soil profile at three different depths (5, 15 and 30 cm) with Theta probes (Knohl et al., 2003). We selected the measurement at 30 cm depth; we deemed this depth to be most representative of the entire soil column which has a depth between 50 and 60 cm. At the Spanish site, soil moisture was assessed at a depth of 15 cm using a water content reflectometer (Pérez-Priego et al., 2013).

Table 2Simulation period at the four FLUXNET sites, and number of months for which latent heat measurements and soil moisture measurements are available to calibrate the model. Soil moisture measurements are not provided at the two French sites.

Download Print Version | Download XLSX

Data to force the model were gap-filled and aggregated from a 30 min to a daily timescale. V2Karst output observations, namely latent heat flux and soil moisture measurements, were aggregated from a 30 min to a monthly timescale and we discarded the months when more than 20 % of the 30 min data were missing. We discarded the monthly observations of latent heat flux and soil moisture for months in which the forcing data contained many gaps; therefore, the impact of the gap-filling of the data on the simulation results is likely to be too significant to sensibly compare simulated and observed soil moisture and latent heat flux. Additionally, we removed monthly aggregated latent heat flux measurements when the mismatch in the energy balance closure was higher than 50 %, similar to Miralles et al. (2011). We corrected latent heat flux measurements to force the closure of the energy balance assuming that measured latent heat flux (LE; MJ m−2 month−1) and sensible heat flux (H; MJ m−2 month−1) have similar errors. The corrected evapotranspiration estimate (Eact,cor; mm month−1) is calculated as proposed in Foken et al. (2012) and Twine et al. (2000):

(19) E act,cor = R n λ 1 + H LE .

Further details on the data processing and on the analysis of the uncertainty in latent heat flux measurements are reported in Sect. S4 of our Supplement. In particular, we showed that, while the actual value of the latent heat flux can have large uncertainties, we have a much higher confidence in its temporal variations.

Table 2 reports the simulation period and the number of monthly latent heat flux and soil moisture observations that were used to estimate the model parameters at the four FLUXNET sites. We extracted a continuous time series of forcing data covering about 10 years at the German site, 7 years at the Spanish site, 3 years at the French 1 site and 8 years at the French 2 site, while latent heat flux and soil moisture measurements were not available over the entire simulation time series. All model simulations were performed using a 1-year warm-up period, which we found to be sufficient to remove the impact of the initial conditions on the simulation results (see Sect. S5 of our Supplement).

4 Methods

To test the plausibility of V2Karst predictions at the FLUXNET sites, we run three sets of analyses: we first estimate the model parameters using actual ET and soil moisture observations available at the four sites (methods in Sect. 4.1); we then analyse the sensitivity of simulated annual recharge to the model parameters using measured forcing data to understand the controlling processes (Sect. 4.2); finally, we investigate the sensitivity of simulated recharge to precipitation properties and land cover type in virtual experiments, in order to understand the controlling processes in recharge when the forcing data is varied more widely than observed at the study sites (Sect. 4.3). All of the analyses were performed using the SAFE toolbox for global sensitivity analysis (Pianosi et al., 2015).

Table 3Site-specific constrained parameter ranges at the four FLUXNET sites for the vegetation parameters (hveg, rst, LAImin, LAImax and Vr) and the soil storage capacity (Vsoi). More information on how the ranges were determined is provided in Sect. S3 of our Supplement. Parameters are defined in Table 1.

Download Print Version | Download XLSX

4.1 Parameter constraints at the FLUXNET sites using soft rules

We investigate whether it is possible to estimate parameter values that produce plausible simulations based on information available at each FLUXNET site. To this end, and similarly to Hartmann et al. (2015), we use “soft rules” to accept or reject parameter combinations based on the consistency between monthly model simulations on one side, and monthly observations and a priori information on model fluxes on the other side. Using soft rules instead of “hard rules” (i.e. the minimisation of the mismatch between observations and simulations) allows for the identification of a set of plausible parameter sets and accounts for the fact that (1) the observed soil moisture is not strictly commensurate with simulated soil moisture, (2) observations are affected by uncertainties (see Sect. 3.2) and (3) it is not expected that V2Karst simulations closely match site-specific data, as the model structure is based on general understanding of karst systems for large-scale applications and may not account for some site specificities. We define five soft rules to identify acceptable (“behavioural”) parameter combinations:

  1. The bias between observed and simulated actual ET is below 20 %:

    (20) Bias = t M ET ( E act,sim ( t ) - E act,cor ( t ) ) t M ET E act,cor ( t ) < 20 % ,

    where Eact,sim(t) (mm) is the simulated actual ET for month t (sum of transpiration, soil evaporation and evaporation from canopy interception), Eact,cor(t) (mm) is the corrected observed actual ET (Eq. 19) and MET is the set of months for which latent heat measurements are available. This rule allows one to constrain the simulated water balance.

  2. The correlation coefficient (ρET) between observed monthly actual ET (Eact,bw) and simulated total actual ET (Eact,sim) is above 0.6. This rule ensures that the temporal pattern of simulated ET follows the observed pattern.

  3. The correlation coefficient (ρSM) between observed monthly soil moisture (SMobs – % soil saturation) and simulated monthly soil moisture (SMsim – m3 m−3 soil volume) is above 0.6. Simulated soil moisture SMsim for month t is calculated as the average soil moisture within the root zone over all model compartments. This rule guarantees that soil moisture variations are consistent with observations.

  4. Total simulated surface runoff (Qsurf) is less than 10 % of precipitation, in accordance with a priori information on the carbonate rock sites, which attests that runoff is negligible (see Sect. 3.1).

  5. Soil and vegetation parameter values are consistent with a priori information, i.e. they fall within constrained (site-specific) ranges. This rule applies to the parameters for which a priori information is available at the FLUXNET sites, namely hveg, rst, LAImin, LAImax, Vr and Vsoi and the constrained ranges are reported in Table 3. This rule ensures that acceptable model outputs are produced using plausible parameter values.

For each site, we derived a sample size of 100 000 for the 15 parameters of V2Karst using Latin hypercube sampling and unconstrained (wide) ranges for the model parameters to explore a large range of soil and vegetation types. We applied the above rules in sequence to either reject or accept the sampled parameter combinations. We sampled the constrained parameter ranges used in rule 5 more densely so that a sufficiently large number of parameterisations remain after applying rule 5. Similarly to Hartmann et al. (2015), a priori information on parameter ranges (rule 5) is applied last so that we can first assess the constraint of the parameter space based on information on model output only (rules 1–4); we then assess the consistency of the constraint using a priori information (rule 5).

We also note that the thresholds used in rules 1 to 3 are stricter than in the study by Hartmann et al. (2015), in which the threshold for the bias rule 1 was set to 75 % and for the correlation rules 2 and 3 was set to 0. The reason for this is that in Hartmann et al. (2015) behavioural parameter sets had to be consistent with observations at all sites within each climate zone defined in the study, while here we perform the parameter estimation for each site separately and therefore we expect better model performances.

4.2 Global sensitivity analysis

We use a global sensitivity analysis called the elementary effect test (Saltelli et al., 2008), also referred to as the method of Morris (1991), to analyse sensitivity across the entire space of parameter variability (Saltelli et al., 2008). We chose this method in particular because it is well suited for identifying uninfluential parameters (Campolongo et al., 2007; Saltelli et al., 2008), which is one of the objectives of our analysis, and because it can be applied to dependent parameters (in V2Karst it is assumed that VeVrVS,n as explained in Sect. 2.3.1).

The method requires the computation of the elementary effects (EEs) of each parameter in n different baseline points in the parameter space. The EE of the ith parameter xi at a given baseline point x1j,x2j,,xi-1j,xij,,xMj and for a predefined perturbation Δ is assessed as follows:

(21) EE i j = y x 1 j , x 2 j , , x i - 1 j , x i j + Δ , x M j - y x 1 j , x 2 j , , x i - 1 j , x i j , x M j Δ ,

where M is the number of parameters and y is the model output (simulated annual recharge in our case). We analyse the mean of the absolute values of the EEs (denoted by μi*) introduced in Campolongo et al. (2007), which is a measure of the total effect of the ith parameter; we also analyse the standard deviation of the EEs (σi) proposed in Morris (1991), which is an aggregate measure of the intensity of the interactions of the ith parameter with the other parameters and of the degree of non-linearity in the model response to changes in the ith parameter.

The total number of model executions required to compute these two sensitivity indices is n(M+1), where n is the number of baseline points chosen by the user. The baseline points and the perturbation Δ of Eq. (21) were determined following the radial design proposed by Campolongo et al. (2011). The baseline points and the perturbation were randomly selected using Latin hypercube sampling for the 15 parameters of V2Karst, and dropping the parameter sets that did not meet the condition VeVrVS,n. In our application, we used n=500 points, which means that we needed 8000 model executions for each sensitivity analysis for each of the four FLUXNET sites. We derived confidence intervals on the sensitivity indices using 1000 bootstrap resamples, and checked the convergence of the results at the chosen sample size, as suggested in Sarrazin et al. (2016).

4.3 Virtual experiments to analyse sensitivity to climate and land cover change

Our last analysis consists of a set virtual experiments to investigate the sensitivity of recharge simulated by V2Karst to changes in (1) the precipitation properties, specifically monthly total precipitation and daily temporal distribution (i.e. frequency and intensity), which are likely to change (IPCC, 2013) and (2) land cover, specifically from forest to shrub or vice versa.

Virtual experiments permit full control of experimental conditions, thus, changes in model outputs can be unequivocally attributed to changes in model inputs (see e.g. Pechlivanidis et al., 2016; Weiler and McDonnell, 2004). Different studies have used virtual experiments to analyse the impact of precipitation spatial and temporal variability on hydrologic model outputs. In fact, using historical precipitation time series or future projections only allows for the exploration of a limited range of possible realisations, which makes it difficult to disentangle the effects of different precipitation properties on model outputs. Instead, synthetic precipitation time series can be tailored to analyse the impact of specific precipitation characteristics, for instance precipitation spatial distribution (Pechlivanidis et al., 2016; Van Werkhoven et al., 2008) and precipitation temporal distribution, namely frequency and intensity (Jothityangkoon and Sivapalan, 2009; Porporato et al., 2004), storminess (Jothityangkoon and Sivapalan, 2009), and seasonality (Botter et al., 2009; Jothityangkoon and Sivapalan, 2009; Laio et al., 2002; Yin et al., 2014).

In this study, we create a synthetic precipitation time series where the same precipitation event is periodically repeated. The precipitation time series is characterized by the intensity of precipitation events Ip (mm d−1) and the interval between two wet days Hp (days). The duration of each precipitation event is set here to 1 day. Therefore, the average monthly precipitation Pm (mm month−1) for an average month with 30 days is equal to:

(22) P m = 30 I p 1 + H p .

Table 4Values of V2Karst parameters and weather variables used in the virtual experiments. Values for the virtual forest site and the virtual shrub site are based on the characteristics of the German FLUXNET site and Spanish FLUXNET site, respectively. Values of the model parameters (parameters are defined in Table 1) correspond to behavioural parameterisations obtained when calibrating the model at FLUXNET sites. Values of the weather variables (Rn – net radiation, T – temperature, RH – relative humidity and WS – wind speed) correspond to the average values calculated at FLUXNET sites.

* At the virtual shrub site, WS was recalculated at a height of 43.5 m because the original measurement provided at a height of 2.5 m at the Spanish site was too low to simulate a change of land cover to tall vegetation (forest). More details on this are reported in Sect. S4 of our Supplement.

Download Print Version | Download XLSX

To determine the possible range of variation of the three variables, Pm, Ip and Hp, we analysed their distributions at the four FLUXNET sites and over all European and Mediterranean carbonate rock areas using Global Land Data Assimilation System (GLAS) data (Rodell et al., 2004; distributions are reported in Sect. S6 of our Supplement). We found that wide but plausible ranges are as follows: Pm varies between 0 and 500 mm month−1, Ip varies between 0 and 200 m d−1 and Hp varies between 0 and 89 days (note that Hp=0 means that it rains every day). We then derived a set of 2266 precipitation time series by deterministically sampling Pm and Hp within those ranges (and consequently deriving a sampled value of Ip from Eq. 22). We sampled more densely closer to the lower bound of the ranges as this is where lower values of Pm and Hp are more likely to occur. For each of the precipitation time series obtained in this manner, we ran the V2Karst model until a steady-state was reached (i.e. periodic oscillations of all state and flux variables) and we analysed the steady-state monthly average of recharge.

The experiments are conducted at two “virtual sites” that are designed based on the characteristics of the FLUXNET sites. Specifically, we use a virtual “forest site” that has the characteristics of the German site (i.e. its behavioural parameterisations for the soil, epikarst and vegetation parameters, and its climate characteristics) as the baseline and a virtual “shrub site” that has the characteristics of the Spanish site. We do not investigate the effects of varying temperature, net radiation, relative humidity or wind speed characteristics as we did for precipitation, because these variables are correlated (see e.g. Ivanov et al., 2007) and, therefore, they cannot be varied independently. We account for their overall combined effect and we vary them jointly to reproduce two conditions, i.e. winter (low energy for ET) and summer (high energy for ET). For winter (summer) conditions we forced the model using constant values of temperature, net radiation, relative humidity and wind speed that were taken as the average values measured at the FLUXNET sites during winter (summer). To investigate the impact of a change in land cover, we swapped the vegetation parameters between the two virtual sites. Table 4 reports the values of the parameters and the values of the weather variables used at the two virtual sites.

Figure 3Reduction in the number of behavioural parameterisations of the V2Karst model at the four FLUXNET sites, when sequentially applying the five soft rules defined in Sect. 4.1 (no rule – initial sample; rule 1 – ET bias; rule 2 – ET correlation; rule 3 – soil moisture correlation; rule 4 – runoff; and rule 5 – a priori information). Rule 3 could not be applied to the French sites where soil moisture observations were not available.


5 Results

5.1 Parameter constraints at FLUXNET sites

We first present the results of applying the soft rules defined in Sect. 4.1 at the four FLUXNET sites. Figure 3 shows that behavioural parameterisations consistent with all rules can be identified at all sites, although the number of parameterisations is very different from one site to another. Specifically, out of the initial 100 000 randomly generated parameter samples, we found 36 838 behavioural parameterisations at the German site, 147 at the Spanish site, 6354 at the French 1 site and 4077 at the French 2 site. From Fig. 3, we also see that the application of each rule reduces the number of behavioural parameterisations, except for rule 4 (value of total surface runoff <10 % of precipitation), as all model simulations produce less than 7 % surface runoff at all sites. This can be explained by the fact that V2Karst gives priority to recharge production over surface runoff. Therefore, the latter only occurs under extremely wet conditions when all model compartments are saturated.

Figure 4Parallel coordinate plots representing V2Karst behavioural parameterisations, and their corresponding simulated output values, identified when sequentially applying the five soft rules defined in Sect. 4.1 at (a) the German site, (b) the Spanish site, (c) the French 1 site and (d) the French 2 site. Parameters are defined in Table 1. The bias absolute mean error between observed and simulated total actual ET (rule 1), the ρET correlation coefficient between observed and the simulated total actual ET (rule 2), the ρSM correlation coefficient between observed and simulated soil moisture (rule 3) and the Qsurf surface runoff (rule 4) are displayed on the right of the figure. Rule 5 corresponds to application of a priori information on parameter ranges (black vertical bars).


Figure 4 reports a parallel coordinate plot of the behavioural parameter sets and associated values of the output metrics after sequential application of the soft rules. The application of rules 1 to 4 does not significantly reduce the parameter ranges, but it allows for low values of parameters Vr and Vsoi to be discarded at all sites (dark blue lines in Fig. 4). In contrast, the application of rule 5 (a priori parameter ranges, black vertical lines in Fig. 4) permits for a significant reduction in parameter ranges, not only for the parameters that are directly constrained by this rule (hveg, rst, LAImin, LAImax, Vr and Vsoi) but also for the spatial variability coefficient a. Specifically, behavioural values of the a parameter are found to be between 0 and 3.2 at the French 1 site and between 0 and 2.8 at the French 2 site. At the Spanish site, we also observe that the behavioural simulations (red lines) cover some portions of the ranges more densely, specifically higher values of parameters rs,soi and a, and lower values of z0 and Ve. This means that the value for these parameters is more likely to be within these sub-ranges.

Figure 5(a) Simulated recharge (Qepi) and actual ET (Eact) expressed as a percentage of total precipitation, and (b) simulated actual transpiration (Tact), actual soil evaporation (Esact) and actual evaporation from interception (Ecact) expressed as a percentage of Eact. The figure reports the ensemble mean and 95 % confidence intervals calculated over the behavioural simulation ensemble of the V2Karst model at the four FLUXNET sites. Simulated fluxes were evaluated over the period from 1 January 2001 to 17 December 2009 for the German site, 1 January 2006 to 31 December 2008 for the Spanish site (dry years), 1 January 2009 to 30 December 2011 for the Spanish site (wet years), 1 January 2010 to 30 December 2011 for the French 1 site and 1 April 2003 to 31 March 2009 for the French 2 site.


Figure 6Monthly time series of observed variables, namely precipitation input P, actual ET Eact,cor (in blue) and soil moisture SMobs (in green), and monthly time series of simulated variables including recharge Qepi, actual ET Eact,sim and soil moisture in the root zone SMsim (non-behavioural simulations are in grey, and behavioural simulations are in black) at (a) the German site, (b) the Spanish site, (c) the French 1 site and (d) the French 2 site. Blue and green shaded areas correspond to the periods in which observations of ET and soil moisture, respectively, were selected to apply the soft rules of Sect. 4.1 (further details on data processing are provided in Sect. 3.2).


We also analyse the repartition of the simulated water fluxes when using the behavioural parameterisations. The top panel (Fig. 5a) compares the total simulated recharge and the total actual ET, expressed as a percentage of the total precipitation at the four FLUXNET sites (mean and 95 % confidence interval across the behavioural parameterisations). At the Spanish site, we present the results over two different time periods that have very different precipitation values, namely a drier period from 1 January 2006 to 31 December 2008 and a wetter period from 1 January 2009 to 30 December 2011 (see Fig. 2). Figure 5 shows that, apart from extremely wet periods at the Spanish site, in all other cases the fraction of recharge (Qepi) is significantly lower than the fraction of actual ET (ETact). Figure 5b shows the partitioning of ET among its different components (transpiration, soil evaporation and interception). We observe that transpiration (Tact) is the largest component at all sites, while the relative importance of evaporation from canopy interception (Ecact) and soil evaporation (Esact) varies across sites. In particular, at the German site Ecact is remarkably high on average compared to the other sites.

Finally, Fig. 6 reports monthly time series of observed variables, namely precipitation input P, actual ET Eact,cor (in blue) and soil moisture SMobs (in green). It also reports time series of simulated variables, including recharge Qepi, actual ET Eact,sim and soil moisture in the root zone SMsim (non-behavioural simulations are in grey and behavioural simulations in black). We see that the soft rules allow for a significant reduction of the uncertainty in model outputs at all sites. In fact, the width of the behavioural ensemble, i.e. the ensemble of simulations obtained by the application of the rules (black lines), is much narrower than the non-behavioural ensemble (grey lines). Simulated actual ET (Eact,sim) is also closer to the observations (blue line) in the behavioural ensemble compared to the non-behavioural ensemble. This means that the application of the soft rules and a priori information on parameter ranges allows not only for an improvement of the precision of the simulated states and fluxes (reduced uncertainty ranges of the simulations), but also for an improvement of the accuracy of simulated actual ET (simulations close to observations).

From Fig. 6, we also observe that the seasonal variations in model predictions are consistent with our understanding of the sites over the entire simulation horizon and not only over the months for which ET and soil moisture observations are used to estimate the parameters (blue and green areas in the plot). Specifically, at the German site we find a marked seasonality of simulated Eact,sim and SMsim: low Eact,sim and high SMsim in winter and high Eact,sim and low SMsim in spring and summer. In fact, in winter, the energy available for ET is low and the deciduous vegetation is not able to transpire or intercept large amounts of precipitation. On the contrary, in spring and summer more energy is available for ET and the vegetation has a higher LAI value; therefore, losses due to ET can occur during this time and deplete the soil moisture. At the other sites we observe a similar pattern for SMsim, while Eact tends to peak in spring and be lower in summer when the ET fluxes are more water-limited than at the German site.

Figure 7Sensitivity indices of the V2Karst parameters (μ* is the mean of the absolute elementary effects and σ is the standard deviation of the elementary effects) for total simulated recharge (expressed as a percentage of total precipitation) at the four FLUXNET sites, when constrained (site-specific) parameter ranges are used (ranges of Table 3) and when unconstrained ranges are used (ranges of Table 1). Sensitivity indices were computed over the period from 1 January 2001 to 17 December 2009 for the German site, 1 January 2006 to 31 December 2008 for the Spanish site (dry years), 1 January 2009 to 30 December 2011 for the Spanish site (wet years), 1 January 2010 to 30 December 2011 for the French 1 site and 1 April 2003 to 31 March 2009 for the French 2 site.


5.2 Global sensitivity analysis of the model parameters

The sensitivity analysis results are reported in Fig. 7 and refer to the sensitivity of the total simulated recharge (expressed as a percentage of total precipitation) to the 15 parameters of the V2Karst model. For each parameter, the plots in Fig. 7 report the absolute mean of the elementary effects (μ*, total effect of the parameters) on the horizontal axis and their standard deviation (σ, degree of linearity of the effect of the parameters) on the vertical axis. In all plots, we observe that the bootstrap confidence intervals of the sensitivity indices are narrow and show little overlap, which gives confidence that the sensitivity results are robust. Similar to the analysis of the simulated fluxes in Sect. 5.1 (Fig. 5), at the Spanish site we present the results for two different time periods with different precipitation amounts.

We first examine the left panels in Fig. 7, which show the sensitivity results when hveg, rst, LAImin, LAImax, Vr and Vsoi are sampled within the constrained ranges (Table 3), while the ranges for other parameters are taken from Table 1. The objective is to inform model calibration in future model applications, as such parameter ranges capture the uncertainty in parameter values left after considering site-specific information. We first note that μ* and σ take a non-zero value for all parameters at all sites, which means that all parameters are influential and have a non-linear effect on recharge, possibly through interactions with other parameters.

We observe that the spatial variability coefficient a has the largest influence by far, followed by the Vsoi and Vr parameters. In fact, their value of μ* is significantly higher than the other parameters at all sites. The implication for model calibration in future applications of V2Karst is that efforts should primarily seek to reduce the uncertainty in a, Vsoi and Vr parameters. These three parameters also have a significantly large value of σ, which indicates non-linearities in the model response to variations in these parameters. Interestingly, Vcan, which controls evaporation from interception, has an impact on recharge at most sites, and rs,soi, which controls soil evaporation, has an impact on recharge at the Spanish site during wet years. This shows that the processes of evaporation from interception and soil evaporation can be important for recharge simulations.

Moreover, we observe that the LAImin, z0, k and Ve parameters have a very small impact on total recharge at all sites (μ*<3 %). However, Sect. S7 of our Supplement reports additional sensitivity analysis results for other model outputs and shows that the most influential parameters which should be the focus of the calibration strategy vary depending on the output of interest.

We then examine the right panels of Fig. 7, which show the sensitivity indices when sampling parameters within unconstrained ranges (Table 1). This analysis allows one to test the plausibility of the model structure through the assessment of the model sensitivity across a large spectrum of soil and vegetation conditions.

The most apparent difference with respect to the previous results is that vegetation parameters (hveg, rst, LAImin, LAImax and Vr) now have a much higher value of the sensitivity indices (both μ* and σ). More specifically, LAImax has a very high sensitivity index at all sites (μ*>10.5 %), which can be attributed to the fact that this parameter is used to calculate different model components. Interestingly, the seasonality of the LAI appears to play an important role in V2Karst as μ* for LAImin, although always lower than μ* for LAImax, stands out at all sites.

The considerable variations in parameter sensitivities observed across sites can be interpreted by considering their climatic differences. We would expect transpiration to be mainly energy-limited at the German site, given that it has a suboceanic–submountain climate; in contrast, we would expect it to be mainly water-limited at the French sites, which have a Mediterranean climate, and at the Spanish site, which has a semi-arid Mediterranean climate. In this regard, the most influential parameter at the Spanish site is by far the a parameter (high μ*), which has an impact on the water storage in the soil and therefore on the amount of water available to sustain ET between rain events, while at the German site, parameters rst and hveg, which control PET, are more influential compared to the other sites.

Finally, we observe that, the parameters that specifically control the volume of transpiration (rst and Vr) have a significantly higher value of μ* than the parameters that specifically control soil evaporation (z0, rs,soi and Ve) and evaporation from interception (Vcan). Moreover, z0, rs,soi and Ve have a very small impact (μ*<3 %), while parameter Vcan can have an important effect at the German site (μ*=5.7 %). Additionally, we see that parameter fred, which controls transpiration from the third soil layer, has an impact on recharge simulated at the Spanish site.

Figure 8Average monthly recharge (Qepi) simulated with V2Karst for different values of average monthly precipitation Pm (mm month−1) and the interval between wet days Hp (days) of the synthetic periodic precipitation input used to force the model, at the virtual forest and shrub sites and under winter and summer conditions.


5.3 Virtual experiments

After showing that the V2Karst model behaves reasonably at the four FLUXNET sites, in this section we use virtual experiments to further learn about the sensitivity of simulated recharge to precipitation characteristics and land cover using virtual sites (see Sect. 4.3). Figure 8 shows the monthly average value of simulated recharge Qepi, for a range of synthetic precipitation inputs with different values of the precipitation monthly amount Pm (x axis) and of the interval between rainy days Hp (y axis) at the virtual forest and shrub sites. We do not report Qepi values in the top right of the plots because this region corresponds to very intense precipitation events (higher than 200 mm d−1) that have a very low probability of occurrence (see Sect. 4.3).

From the top left panel of Fig. 8, we see that winter Qepi is mostly sensitive to Pm, in fact simulated recharge increases along the x axis from left to right but shows little variations along the y axis (when Hp is varied). This result is due to the fact that the actual ET is very limited in winter because of the low energy available. We indeed estimated that the maximum value of total ET across the different precipitation inputs is 13 mm month−1 at the forest site and 35 mm month−1 at the shrub site. Therefore, a large part of precipitation becomes recharge rather independently of its temporal distribution.

Figure 9Change in monthly recharge (ΔQepi=Qepishrub-Qepiforest) simulated with V2Karst when the land cover is set to shrub compared to forest for different values of average monthly precipitation Pm (mm month−1) and the interval between wet days Hp (days) of the synthetic periodic precipitation input used to force the model, at the virtual forest and shrub sites and under winter and summer conditions.


From the right panel of Fig. 8, we observe a systematic reduction in summer Qepi compared to winter at both virtual sites. Moreover, summer recharge is generally highly sensitive not only to Pm but also to Hp, as it increases when moving along the y axis from bottom to top, i.e. when the same amount of monthly precipitation falls in less frequent but more intense events. This result can be explained by the fact that in summer potential ET is larger and therefore, if events are less intense, a larger part of the precipitation is lost via ET. In contrast, if events are more intense, the canopy and soil stores reach saturation. In this case, a saturation excess flow to the epikarst is generated – hence, there is more recharge and less ET. Moreover, in summer, Qepi shows a limited sensitivity to Pm and Hp when these quantities take low values (brown and red dots on the left of the plots). This is due to the fact that only a few soil compartments reach saturation under drier conditions and therefore little recharge can be generated. We also note that Qepi is a significant flux (Qepi>5 mm) for smaller values of Pm and Hp at the shrub site compared to the forest site. This may be due to the fact that the soil water capacity (Vsoi) is much smaller at the shrub site, and the soil compartments can therefore reach saturation under drier conditions.

Finally, Fig. 9 reports the results of another virtual experiment which is focused on the impact of land cover change. Specifically, the panels in Fig. 9 show the variation in simulated recharge when land cover is changed from forest to shrub at the virtual forest site (and vice versa at the virtual shrub site), and more specifically, Fig. 9 reports ΔQepi=Qepishrub-Qepiforest. We see that in all plots ΔQepi is positive, which means that recharge is larger and actual ET is lower under shrub land cover compared with forest land cover for both sites. From the left panels of Fig. 9, we observe that ΔQepi is very limited in winter, which is expected as ET fluxes are small in winter as explained above.

Conversely, the right panels of Fig. 9 show that summer ΔQepi is much higher than in winter conditions and is sensitive to both Pm and Hp. The sensitivity of ΔQepi is highly variable across the different precipitation inputs, and more specifically an increase in Hp can have a different effect on ΔQepi depending on the value of Pm (no variation, increase or decrease in ΔQepi). In fact, when Pm is low, ΔQepi is always low and does not vary sensibly when Pm and Hp are varied (brown area in the left end of the plot), since recharge is always low under these precipitation conditions as shown in Fig. 8. For intermediate values of Pm, ΔQepi has a similar pattern at both sites and increases when either Hp or Pm increases. Instead, for high values of Pm, we see that for both sites ΔQepi decreases when Hp increases. ΔQepi is largest when the monthly precipitation Pm is high and the interval between wet days Hp is low (green dots at the virtual forest site and dark blue dots at the virtual shrub site), because under these precipitation conditions the amount of moisture available for ET is maximum.

Importantly, our results also show that the impact of a change in land cover can vary greatly across sites, as at the virtual shrub site summer ΔQepi reaches much higher values and is sensitive to Pm and Hp over a larger range of values of Pm and Hp compared to the virtual forest site.

6 Discussion

6.1 Plausibility of V2Karst simulations against site specific information

We tested the model by evaluating its ability to reproduce observations at four carbonate rock FLUXNET sites, which is a standard approach to model testing, used for instance to test the previous version of the model VarKarst (Hartmann et al., 2015) and large-scale ET products (Martens et al., 2017; McCabe et al., 2016; Miralles et al., 2011). We demonstrated that V2Karst is able to produce behavioural simulations consistent with observations and a priori information at FLUXNET sites. In addition, the time series of simulated water balance components are coherent with our understanding of the sites. A different number of behavioural parameterisations was identified at the different sites, with the highest number found at the more humid German site and the lowest number at the semi-arid Spanish site, which is coherent with previous findings that a better fit to observations can be obtained at wetter locations (Atkinson et al., 2002; Bai et al., 2015).

Interestingly, the behavioural parameter sets for the French 1 site corroborate the hypothesis that root water uptake is likely to extend below the physical soil layer, as communicated by Guillaume Simioni (investigator of the site). In fact, we found behavioural values of the Vr parameter above 59 mm, while site-specific information indicates that the physical soil layer has a storage capacity of 49 mm (Table B1). This result further attests to the realism of V2Karst structure.

The results of the global sensitivity analysis of simulated recharge to the model parameters are easily interpreted in light of the different climatic conditions at the four FLUXNET sites. Parameters that control PET and the soil water storage capacity generally have a large impact on simulated recharge. The interception capacity is also very important at the German site. These findings are in line with the few sensitivity analysis studies of large-scale hydrological models (Güntner et al., 2007; Werth et al., 2009), which have examined the sensitivity of continental water storage and river discharge simulated by the WaterGAP model. The importance of the maximum leaf area index and, to a lesser extent, of the seasonality of vegetation in the V2Karst model is also consistent with previous studies. For example, Tesemma et al. (2015) found that assimilating year-to-year monthly LAI in the VIC model can significantly improve runoff simulations compared with using long-term average LAI; furthermore, Rosero et al. (2010) determined that LAI has a large influence on simulated latent heat in the Noah land surface model. Therefore, our sensitivity analysis results generally suggest that the model behaves sensibly and consistently with our understanding of the key vegetation–recharge processes that we aim to reproduce.

6.2 Sensitivity of simulated groundwater recharge to changes in climate and vegetation characteristics

Our results also provide valuable insights about possible vulnerabilities of karst systems. For example, the fact that the LAI is a highly influential parameter at all FLUXNET sites suggests that the increasing trend in global LAI documented by Zhu et al. (2016) could be reflected as a significant impact on groundwater recharge and on the partitioning between green and blue water in karst areas.

Variations in simulated recharge were also systematically examined under combined changing land cover type and climate conditions (precipitation overall amount and temporal distribution), as the model should be appropriate for climate and land cover change impact studies. We used virtual experiments with synthetic data because they allow one to unequivocally attribute changes in recharge to changes in climate and land cover. Importantly, we showed that precipitation characteristics and land cover have an interacting effect on simulated recharge; these findings are consistent with the global and observational analysis by Kim and Jackson (2012) of non-karst areas, which revealed interactions among vegetation and climate (in particular the overall amount of precipitation) in producing recharge. We found that simulated recharge is lower (and ET is higher) for forest than for the shorter vegetation type considered, i.e. shrubs, which is in line with expectations and with previous studies undertaken in karst areas (Ford and Williams, 2007; Williams, 1993). We demonstrated that the overall amount and daily intensity of precipitation have a positive effect on simulated recharge, i.e. the higher these factors are the higher the recharge. This is consistent with previous observational studies, which found that heavy precipitation events enhance groundwater recharge in non-karst areas (e.g. Taylor et al., 2013, in a semi-arid tropical region; Owor et al., 2009, in a seasonally humid tropical region). Conversely, the modelling study by Hartmann et al. (2017) showed that heavy precipitation events can have a negative effect on recharge simulated by the VarKarst model in some carbonate rock areas in Europe, the Middle East and northern Africa. However, Hartmann et al. (2017) analysed the sensitivity of simulated recharge using forcing data from global circulation models (GCMs). Given that climate variables from GCMs show complex patterns, it might be difficult to isolate the effect of a specific climate property (e.g. heavy precipitation events) on the simulated fluxes, which could explain the differences between their findings and ours.

Precipitation patterns are more complex than the simple periodic variations we used here, and the steady-state conditions may never be reached in practice. Nevertheless, we believe that performing virtual experiments similar to those proposed in the present study is a complementary approach to the application of climate projections provided by GCMs and future land cover change scenarios (e.g. Holman et al., 2017; Hurtt et al., 2011), to systematically assess the sensitivity of a model to changes in input characteristics, to test the soundness of the model sensitivities and to determine which aspects of a model input would be worth further investigating.

6.3 Applying V2Karst over larger domains

The results of our global sensitivity analyses suggest that all newly introduced processes in V2Karst are relevant for applications over large domains because the parameters that control these processes can affect simulated recharge, depending on the climatic, soil and land cover conditions. Specifically, transpiration and vegetation seasonality are important processes under the climate of the four FLUXNET sites, evaporation from canopy interception is important at all forested sites (German site and two French sites), the contribution of water stored below the root zone is important under the climate of the Spanish site and soil evaporation is important at the semi-arid site with sparse and short vegetation (Spanish site). The significance of representing canopy interception, in particular for forested land covers, has already been mentioned in previous studies (Gerrits, 2010; Savenije, 2004), whilst the importance of separating transpiration and soil evaporation was reported in Maxwell and Condon (2016) and Wang and Dickinson (2012).

Regarding parameter estimation, the use of soft rules similar to this study and to the large-scale study of Hartmann et al. (2015) can be envisaged for applications over large domains. We showed that a priori information on parameter ranges is needed to constrain the simulations. At large scales, a priori information on vegetation parameters can be derived from large-scale databases of vegetation properties (more details in Sect. S1 of our Supplement), and a priori information on parameters that control the subsurface properties (Vsoi) can be derived following Hartmann et al. (2015). It is also particularly important to assess the sensitivity of model parameters across the modelling domain to test the suitability of fixing model parameters, as done in this study at FLUXNET sites. We indeed found that parameters Vcan and rs,soi, which are typically fixed in the other large-scale models detailed in Table A1, do have an impact on simulated recharge at least one FLUXNET site. Moreover, as mentioned in Sect. 2.3, Vcan and rs,soi are understood to vary across land cover types and soil types, respectively, even if no clear ranges of these parameters have been established across either land cover or soil types. Therefore, fixing these two parameters could potentially introduce large uncertainties in V2Karst simulations. Other studies have reported on this issue; Cuntz et al. (2016), in particular, showed that some constant parameters of the Noah-MP land surface model can be highly influential for some model outputs. Finally, parameter estimation for applications over large domains will need to account for the large variability in hydraulic properties and recharge patterns observed across karst systems reported, for example, in Hartmann et al. (2014) and Klimchouk and Ford (2000), for instance using a simplified classification of karst landscapes based on climate and topography as in Hartmann et al. (2015).

We note that although soils may be absent in some karst areas (e.g. karren field in high mountain areas as reported in Hartmann et al., 2014), V2Karst always includes a soil layer. However, this soil layer can be very thin and has a limited impact on recharge. Additionally, when the simulation units have a large extent, such as 0.25×0.25 in Hartmann et al. (2015), it is reasonable to assume that soil layers can always be found.

The V2Karst model can account for the sub-grid heterogeneity in the vegetation type using a “tile” approach as in Mac-PDM (Gosling and Arnell, 2011) and VIC (Bohn and Vivoni, 2016), which consists of subdividing each model grid cell into a number of independent units (tiles), each of which has a specific land cover (e.g. short or tall vegetation). The model is then evaluated separately over each tile and the overall simulated fluxes are computed as the area weighted average of the fluxes calculated over the tiles. In V2Karst, each “tile” includes a vegetated and non-vegetated fraction as explained in Sect. 2.3.3.

The V2Karst model implements a simplified representation of ET and karst processes, given the limited amount of data available at large scales as discussed in Sect. 2.1. In particular, further information on typical karstic features over large domains would be valuable to refine the representation of karst processes in V2Karst in future developments of the model, to consider additional processes explicitly (the formation of a perched aquifer in the epikarst, which can occur in highly weathered karst systems as described e.g. in Williams (1983, 2008); this process is represented e.g. in a previous karst model developed for applications at the local scale introduced in Hartmann et al., 2012b).

7 Conclusions

The objectives of the present study were (1) to develop and test an ET component with the explicit representation of land cover processes for the karst recharge model VarKarst, so that the model can be used for combined climate and land cover change impact assessment at large scales, and (2) to evaluate the mechanisms of recharge production in the model as well as the model sensitivity to temporal precipitation patterns and land cover using observations and synthetic data to force the model.

Many different approaches are used to represent ET in large-scale hydrologic models, and the lack of in situ ET observations makes it difficult to assess and compare the performance of these different formulations. Moreover, some models use a large number of parameters that can only be poorly constrained by the few available observations. High model complexity also makes Monte Carlo simulation computationally expensive and hampers uncertainty and sensitivity analysis. The new version of the VarKarst model developed here, V2Karst (V1.1), is the first large-scale hydrological model to include explicit representations of both karst and land cover processes. We sought to include parsimonious process descriptions that are understood to be relevant for climate and land cover impact assessment, namely, (1) a representation of the three ET components (transpiration, soil evaporation in presence of sparse canopy and evaporation from canopy interception) and (2) a physically based PET equation (Penman–Monteith). The model also comprises a parsimonious representation of vegetation seasonality.

We demonstrated that V2Karst produces credible simulations at four carbonate rock FLUXNET sites, by showing that it reproduces observations of latent heat and soil moisture and that the parameters which dominate the model sensitivity are in accordance with our perception of expected controls on recharge. We also established that the model has plausible sensitivities to climate and land cover changes when using virtual experiments with synthetic precipitation and land cover scenarios. Additionally, we showed that the newly introduced processes in V2Karst can have varying impacts on simulated recharge quantities depending on the climate, the soil properties and the land cover.

Overall, our study demonstrates the value of a model development and evaluation process that considers both how well a model reproduces historical observations as well as how this performance is achieved. The latter is examined using global sensitivity analysis to understand which processes and inputs dominate the model output. Through the proposed process we have a higher chance of obtaining the right results for the right reasons than we would have using performance analysis alone (Kirchner, 2006). We further believe – given the lack of observations of many hydrological fluxes at large scales – that an integrated analysis using historical data, our expectations derived from previous studies and synthetic data (to understand model behaviour more generally) can provide a more holistic view of a model than using observations alone.

Code availability

The code of the V2Karst model is open source and freely available under the terms of the GNU General Public License version 3.0. The model code is written in Matlab and is provided through a Github repository (, last access: 30 November 2018); it is assigned a DOI (; Sarrazin et al., 2018).

Appendix A: Review of ET component in large-scale hydrological models

Table A1Characteristics of selected large-scale models: simulation time step (Δt), representation of sub-grid variability of soil moisture, solving of the energy balance, ET processes represented (overstorey transpiration Tactover, understory transpiration Tactunder, soil evaporation Esact, evaporation from canopy interception Ecact and carbon cycle i.e. vegetation dynamic model), and number of parameters for ET estimation. The models were selected based on the following criteria: (1) explicit representation of land cover properties, (2) calculation of ET and soil water balance at a daily or sub-daily time step, and (3) applications in previous studies over a wide range of climate and land cover types. Tables A2 and A3 present the parameterisations used in these models, while a detailed list of all parameters involved in the representation of ET in these models can be found in Sect. S2 of our Supplement.

a None of these models account for karst processes as done by the VarKarst model (Hartmann et al., 2015). b Number of parameters for a given land cover type, excluding parameters used in the representation of vegetation seasonality, carbon cycle (vegetation dynamic), sublimation from snowpack and snowmelt evaporation to make models more comparable. c Number of parameters considering three soil layers. d This number includes the parameters used for the computation of both understory and overstorey (grass) transpiration. e Number of parameters assuming tall vegetation (interception is considered for tall vegetation only).

Download Print Version | Download XLSX

Table A2Representation of potential evapotranspiration (PET) and stress model for actual evapotranspiration (ET) calculation from PET in the large-scale models of Table A1.

a T: Thornthwaite (1948); HS: Hargreaves and Samani (1985); PT: Priestley–Taylor (1972); PM: Penman–Monteith (Monteith, 1965); SW: Shuttleworth and Wallace (1985). b This approach consists of calculating a value of PET for a reference grass surface with known properties and to adjust this potential rate using land cover specific empirical crop factors. This formulation avoids the specification of the stomatal resistance whose value is largely uncertain (see Sect. S1 of our Supplement). Tabulated values of the crop factors for agricultural crops are provided in Allen et al. (1998). However, the origin of the crop factor formulation for non-agricultural crops is not clear. c This number includes the parameters used for the computation of PET for both understory and overstorey (grass). d Number of parameters assuming tall vegetation (interception is considered for tall vegetation only).

Download Print Version | Download XLSX

Table A3Representation of sparse vegetation, soil layers, evaporation from canopy interception and seasonality of vegetation in the large-scale models of Table A1.

a Uncoupled approaches consist of separately assessing the water balance for the vegetated and bare soil fractions (overstorey and understory fractions for Mac-PDM). Therefore, this approach is based on the simplifying assumption that the vegetation roots do not extent beyond the surface area covered by the vegetation canopy. Instead, coupled approaches evaluate the overall water balance over both fractions, thus allowing for interactions for soil moisture uptake between vegetated and bare soil fractions. All models neglect aerodynamic interactions between vegetation and bare soil. This can be accounted for using for instance the Shuttleworth–Wallace PET equation (Shuttleworth and Wallace, 1985), which requires the specification of further resistance parameters compared to the Penman–Monteith equation. The Shuttleworth–Wallace equation was used anecdotally in the WBM model for a few applications. b Esact: actual soil evaporation; Tact: actual vegetation transpiration. c This number includes the parameters used for the computation of PET for both understory and overstorey (grass). d Number of parameters assuming tall vegetation (interception is considered for tall vegetation only).

Download Print Version | Download XLSX

Appendix B: Additional information on the four carbonate rock FLUXNET sites

Table B1Description of the four carbonate rock FLUXNET sites.

* Between wilting point and field capacity.

Download Print Version | Download XLSX


The supplement related to this article is available online at:

Author contributions

FS, AH, RR and TW designed the model equations of V2Karst. FS, AH, FP and TW contributed to the methodology to test the model. FS led the analyses, performed the numerical experiments and processed the data, and developed the code of the V2Karst model; the V2Karst V1.1 code was built upon the code of the previous version of the model (VarKarst), which was developed by AH. AH contributed more particularly to the karst aspect of the work, RR to the land surface modelling aspect of the work and to processing the data, and FP and TW to the sensitivity analysis aspect of the work. The paper was prepared by FS and edited by all of the authors.

Competing interests

The authors declare that they have no conflict of interest.


Fanny Sarrazin was supported by a PhD scholarship programme from the University of Bristol Development and Alumni Relations Office. Support to Andreas Hartmann was provided by the Emmy Noether Programme of the German Research Foundation (DFG; grant no. HA 8113/1-1; project “Global Assessment of Water Stress in Karst Regions in a Changing World”). Francesca Pianosi was partially supported by a UK Engineering and Physical Sciences Research Council fellowship (EPSRC; grant no. EP/R007330/1). Support to Rafael Rosolem was provided by the NERC “A MUlti-scale Soil moisture Evapotranspiration Dynamics study” (AMUSED project; grant no. NE/M003086/1) and the “Brazilian Experimental datasets for MUlti-Scale interactions in the critical zone under Extreme Drought” (BEMUSED project, grant no. NE/R004897/1). Thorsten Wagener was partially supported by a Royal Society Wolfson Research Merit Award. We thank Timothy Foster for advice on ET modelling. We thank Yoshihide Wada for providing information on ET representation in the PCR-GLOBWB model. We thank the investigators of the four FLUXNET sites (Hainich, Llano de los Juanes, Font-Blanche and Puéchabon) for allowing us to use their data, and in particular Martina Mund and Manfred Fink (De-Hai), Penelope Serrano Ortiz, Francisco Domingo Poveda and Andrew Kowalski (ES-Lju), Richard Joffre and Serge Rambal (FR-Pue), and Guillaume Simioni (FR-FBn) for sharing information on the characteristics of the FLUXNET sites. We appreciate the support from Eleonora Canfora and Dario Papale from the FLUXNET European Fluxes Database Cluster. We are grateful to two anonymous referees for the detailed comments that helped us to significantly improve the paper.

Edited by: Didier Roche
Reviewed by: two anonymous referees


Abramowitz, G., Leuning, R., Clark, M., and Pitman, A.: Evaluating the performance of land surface Models, J. Clim., 21, 5468–5481,, 2008. 

Alcalá, F. J., Cantón, Y., Contreras, S., Were, A., Serrano-Ortiz, P., Puigdefábregas, J., Solé-Benet, A., Custodio, E., and Domingo, F.: Diffuse and concentrated recharge evaluation using physical and tracer techniques: results from a semiarid carbonate massif aquifer in southeastern Spain, Environ. Earth Sci., 62, 541–557,, 2011. 

Allen, R. G., Pereira, L. S., Raes, D. and Smith, M.: Crop evapotranspiration: Guidelines for computing crop requirements, FAO Irrigation and Drainage Paper 56, Food and Agriculture Organization (FAO), Rome, Italy, 1998. 

Allen, R. G., Pruitt, W. O., Wright, J. L., Howell, T. A., Ventura, F., Snyder, R., Itenfisu, D., Steduto, P., Berengena, J., and Yrisarry, J. B.: A recommendation on standardized surface resistance for hourly calculation of reference ETo by the FAO56 Penman-Monteith method, Agric. Water Manag., 81, 1–22,, 2006. 

Arbel, Y., Greenbaum, N., Lange, J., and Inbar, M.: Infiltration processes and flow rates in developed karst vadose zone using tracers in cave drips, Earth Surf. Process. Landforms, 35, 1682–1693,, 2010. 

Archfield, S. A., Clark, M., Arheimer, B., Hay, L. E., McMillan, H., Kiang, J. E., Seibert, J., Hakala, K., Bock, A., Wagener, T., Farmer, W. H., Andréassian, V., Attinger, S., Viglione, A., Knight, R., Markstrom, S., and Over, T.: Accelerating advances in continental domain hydrologic modeling, Water Resour. Res., 51, 10078–10091,, 2015. 

Arnell, N. W.: A simple water balance model for the simulation of streamflow over a large geographic domain, J. Hydrol., 217, 314–335,, 1999. 

Atkinson, S. E., Woods, R. A., and Sivapalan, M.: Climate and landscape controls on water balance model complexity over changing timescales, Water Resour. Res., 38, 1314,, 2002. 

Bai, P., Liu, X., Liang, K., and Liu, C.: Comparison of performance of twelve monthly water balance models in different climatic catchments of China, J. Hydrol., 529, 1030–1040,, 2015. 

Baldocchi, D., Falge, E., Gu, L., Olson, R., Hollinger, D., Running, S., Anthoni, P., Bernhofer, C., Davis, K., Evans, R., Fuentes, J., Goldstein, A., Katul, G., Law, B., Lee, X., Malhi, Y., Meyers, T., Munger, W., Oechel, W., Paw, K. T., Pilegaard, K., Schmid, H. P., Valentini, R., Verma, S., Vesala, T., Wilson, K., and Wofsy, S.: FLUXNET: A New Tool to Study the Temporal and Spatial Variability of Ecosystem-Scale Carbon Dioxide, Water Vapor, and Energy Flux Densities, B. Am. Meteorol. Soc., 82, 2415–2434,<2415:FANTTS>2.3.CO;2, 2001. 

Bargués Tobella, A., Reese, H., Almaw, A., Bayala, J., Malmer, A., Laudon, H., and Ilstedt, U.: The effect of trees on preferential flow and soil infiltrability in an agroforestry parkland in semiarid Burkina Faso, Water Resour. Res., 50, 3342–3354,, 2014. 

Beven, K. and Germann, P.: Macropores and water flow in soils revisited, Water Resour. Res., 49, 3071–3092,, 2013. 

Beven, K. J. and Cloke, H. L.: Comment on “hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth's terrestrial water” by Eric F. Wood et al., Water Resour. Res., 48, W01801,, 2012. 

Bierkens, M. F. P.: Global hydrology 2015: State, trends, and directions, Water Resour. Res., 51, 4923–4947,, 2015. 

Blume, H.-P., Brümmer, G. W., Horn, R., Kandeler, E., Kögel-Knabner, I., Kretzschmar, R., Stahr, K., and Wilke, B.-M.: Lehrbuch der Bodenkunde, Springer-Verlag, Berlin Heidelberg,, 2010. 

Bohn, T. J. and Vivoni, E. R.: Process-based characterization of evapotranspiration sources over the North American monsoon region, Water Resour. Res., 52, 358–384,, 2016. 

Boone, A., Calvet, J.-C., and Noilhan, J.: Inclusion of a third soil layer in a land surface scheme using the force–restore method, J. Appl. Meteorol., 38, 1611–1630,<1611:IOATSL>2.0.CO;2, 1999. 

Botter, G., Porporato, A., Rodriguez-Iturbe, I., and Rinaldo, A.: Nonlinear storage-discharge relations and catchment streamflow regimes, Water Resour. Res., 45, W10427,, 2009. 

Brown, A. E., Zhang, L., Mcmahon, T. A., Western, A. W., and Vertessy, R. A.: A review of paired catchment studies for determining changes in water yield resulting from alterations in vegetation, J. Hydrol., 310, 28–61,, 2005. 

Calder, I. R. (Ed.): Evaporation in the Uplands, John Wiley & Sons Ltd., Chichester, UK, 1990. 

Campolongo, F., Cariboni, J., and Saltelli, A.: An effective screening design for sensitivity analysis of large models, Environ. Model. Softw., 22, 1509–1518,, 2007. 

Campolongo, F., Saltelli, A., and Cariboni, J.: From screening to quantitative sensitivity analysis. A unified approach, Comput. Phys. Commun., 182, 978–988,, 2011. 

Canora, F., Fidelibus, M. D., Sciortino, A., and Spilotro, G.: Variation of infiltration rate through karstic surfaces due to land use changes: A case study in Murgia (SE-Italy), Eng. Geol., 99, 210–227,, 2008. 

Cantón, Y., Villagarcía, L., José Moro, M., Serrano-Ortíz, P., Were, A., Javier Alcalá, F., Kowalski, A. S., Solé-Benet, A., Lázaro, R., and Domingo, F.: Temporal dynamics of soil water balance components in a karst range in southeastern Spain: estimation of potential recharge, Hydrol. Sci. J., 55, 737–753,, 2010. 

Chaney, N. W., Herman, J. D., Ek, M. B., and Wood, E. F.: Deriving global parameter estimates for the Noah land surface model using FLUXNET and machine learning, J. Geophys. Res., 121, 13218–13235,, 2016. 

Chen, Z., Auler, A. S., Bakalowicz, M., Drew, D., Griger, F., Hartmann, J., Jiang, G., Moosdorf, N., Richts, A., Stevanovic, Z., Veni, G., and Goldscheider, N.: The World Karst Aquifer Mapping project: concept, mapping procedure and map of Europe, Hydrogeol. J., 25, 771–785,, 2017. 

Clapp, R. B. and Hornberger, G. M.: Empirical equations for some soil hydraulic properties, Water Resour. Res., 14, 601–604,, 1978. 

Contreras, S., Boer, M. M., Alcala, F. J., Domingo, F., Garcia, M., Pulido-Bosch, A., and Puigdefabregas, J.: An ecohydrological modelling approach for assessing long-term recharge rates in semiarid karstic landscapes, J. Hydrol., 351, 42–57,, 2008. 

COST: Cost action 65 – Hydrogeological aspects of groundwater protection in karstic areas, Report EUR 16547, European Commission, Directorate-General XII Science, Research Development, Luxembourg, 1995. 

Coxon, C.: Agriculture and Karst, in Karst Management, edited by van Beynen, P. E., Springer Netherlands, Dordrecht, 103–138, 2011. 

Cuntz, M., Mai, J., Samaniego, L., Clark, M., Wulfmeyer, V., Branch, O., Attinger, S., and Thober, S.: The impact of standard and hard-coded parameters on the hydrologic fluxes in the Noah-MP land surface model, J. Geophys. Res., 121, 10676–10700,, 2016. 

Cuthbert, M. O., Mackay, R., and Nimmo, J. R.: Linking soil moisture balance and source-responsive models to estimate diffuse and preferential components of groundwater recharge, Hydrol. Earth Syst. Sci., 17, 1003–1019,, 2013. 

DeFries, R. and Eshleman, K. N.: Land-use change and hydrologic processes?: a major focus for the future, Hydrol. Process., 18, 2183–2186,, 2004. 

De Groen, M. M.: Modelling interception and transpiration at monthly time steps?: introducing daily variability through Markov chains, PhD thesis, Delft University of Technology, Delft, the Netherlands, ISBN: 9058093786, 2002. 

Döll, P. and Fiedler, K.: Global-scale modeling of groundwater recharge, Hydrol. Earth Syst. Sci., 12, 863–885,, 2008. 

Döll, P., Kaspar, F., and Lehner, B.: A global hydrological model for deriving water availability indicators: Model tuning and validation, J. Hydrol., 270, 105–134,, 2003. 

Doummar, J., Sauter, M., and Geyer, T.: Simulation of flow processes in a large scale karst system with an integrated catchment model (Mike She) – Identification of relevant parameters influencing spring discharge, J. Hydrol., 426–427, 112–123,, 2012. 

Ecofor: Site atelier de Font Blanche, available at:, last access: 13 December 2017. 

Falkenmark, M. and Rockström, J.: The new blue and green water paradigm?: breaking new ground for water resources planning and management, J. Water Resour. Plan. Manag., 132, 129–132,, 2006. 

Federer, C. A.: Transpirational Supply and Demand: plant, soil, and atmospheric effects evaluated by simulation, Water Resour. Res., 18, 355–362,, 1982. 

Federer, C. A., Vörösmarty, C., and Fekete, B.: Sensitivity of Annual Evaporation to Soil and Root Properties in Two Models of Contrasting Complexity, J. Hydrometeorol., 4, 1276–1290,<1276:SOAETS>2.0.CO;2, 2003. 

Fleury, P., Plagnes, V., and Bakalowicz, M.: Modelling of the functioning of karst aquifers with a reservoir model?: Application to Fontaine de Vaucluse (South of France), J. Hydrol., 345, 38–49,, 2007. 

Foken, T., Leuning, R., Oncley, S. R., Mauder, M., and Aubinet, M.: Corrections and Data Quality Control, in Eddy Covariance: A Practical Guide to Measurement and Data Analysis, edited by: Aubinet, M., Vesala, T., and Papale D., Springer Netherlands, Dordrecht, 85–131, 2012. 

Ford, D. and Williams, P. (Eds.): Karst Hydrogeology and Geomorphology, John Wiley & Sons Ltd., Chichester, UK, 2007. 

Gash, J. H. C.: An analytical model of rainfall interception by forests, Q. J. Roy. Meteor. Soc., 105, 43–55,, 1979. 

Gea-Izquierdo, G., Guibal, F., Joffre, R., Ourcival, J. M., Simioni, G., and Guiot, J.: Modelling the climatic drivers determining photosynthesis and carbon allocation in evergreen Mediterranean forests using multiproxy long time series, Biogeosciences, 12, 3695–3712,, 2015. 

Gerrits, M.: The role of interception in the hydrological cycle, PhD thesis, Delft University of Technology, Delft, the Netherlands, available at: (last access: 30 November 2018), 2010. 

Gerten, D., Schaphoff, S., Haberlandt, U., Lucht, W., and Sitch, S.: Terrestrial vegetation and water balance – Hydrological evaluation of a dynamic global vegetation model, J. Hydrol., 286, 249–270,, 2004. 

Gosling, S. N. and Arnell, N. W.: Simulating current global river runoff with a global hydrological model: Model revisions, validation, and sensitivity analysis, Hydrol. Process., 25, 1129–1145,, 2011. 

Güntner, A., Stuck, J., Werth, S., Döll, P., Verzano, K., and Merz, B.: A global analysis of temporal and spatial variations in continental water storage, Water Resour. Res., 43, W05416,, 2007. 

Hao, Y., Yeh, T. C. J., Gao, Z., Wang, Y., and Zhao, Y.: A gray system model for studying the response to climatic change: The Liulin karst springs, China, J. Hydrol., 328, 668–676,, 2006. 

Hargreaves, G. H. and Samani, Z. A.: Reference crop evapotranspiration from temperature, Appl. Eng. Agric., 1, 96–99,, 1985. 

Hartmann, A. and Baker, A.: Modelling karst vadose zone hydrology and its relevance for paleoclimate reconstruction, Earth-Sci. Rev., 172, 178–192,, 2017. 

Hartmann, A., Lange, J., Vivó Aguado, À., Mizyed, N., Smiatek, G., and Kunstmann, H.: A multi-model approach for improved simulations of future water availability at a large Eastern Mediterranean karst spring, J. Hydrol., 468–469, 130–138,, 2012a. 

Hartmann, A., Lange, J., Weiler, M., Arbel, Y., and Greenbaum, N.: A new approach to model the spatial and temporal variability of recharge to karst aquifers, Hydrol. Earth Syst. Sci., 16, 2219–2231,, 2012b. 

Hartmann, A., Barberá, J. A., Lange, J., Andreo, B., and Weiler, M.: Progress in the hydrologic simulation of time variant recharge areas of karst systems – Exemplified at a karst spring in Southern Spain, Adv. Water Resour., 54, 149–160,, 2013. 

Hartmann, A., Goldscheider, N., Wagener, T., Lange, J., and Weiler, M.: Karst water resources in a changing world: Review of hydrological modeling approaches, Rev. Geophys., 52, 218–242,, 2014. 

Hartmann, A., Gleeson, T., Rosolem, R., Pianosi, F., Wada, Y., and Wagener, T.: A large-scale simulation model to assess karstic groundwater recharge over Europe and the Mediterranean, Geosci. Model Dev., 8, 1729–1746,, 2015. 

Hartmann, A., Gleeson, T., Wada, Y., and Wagener, T.: Enhanced groundwater recharge rates and altered recharge sensitivity to climate variability through subsurface heterogeneity, P. Natl. Acad. Sci. USA, 114, 2842–2847,, 2017. 

Haughton, N., Abramowitz, G., and Pitman, A. J.: On the predictability of land surface fluxes from meteorological variables, Geosci. Model Dev., 11, 195–212,, 2018. 

Hendrickx, J. M. H. and Flury, M.: Uniform and Preferential Flow Mechanisms in the Vadose Zone, in: Conceptual Models of Flow and Transport in the Fractured Vadose Zone, edited by: National Research Council, The National Academies Press, Washington, DC, 149–188, 2001. 

Hogue, T. S., Bastidas, L. A., Gupta, H. V., and Sorooshian, S.: Evaluating model performance and parameter behavior for varying levels of land surface model complexity, Water Resour. Res., 42, W08430,, 2006. 

Holman, I. P., Brown, C., Janes, V., and Sandars, D.: Can we be certain about future land use change in Europe? A multi-scenario, integrated-assessment analysis, Agric. Syst., 151, 126–135,, 2017. 

Hong, E.-M., Pachepsky, Y. A., Whelan, G., and Nicholson, T.: Simpler models in environmental studies and predictions, Crit. Rev. Environ. Sci. Technol., 47, 1669–1712,, 2017. 

Hurtt, G. C., Chini, L. P., Frolking, S., Betts, R, A., Feddema, J., Fischer, G., Fisk, J. P., Hibbard, K., Houghton, R. A., Janetos, A., Jones, C. D., Kindermann, G., Kinoshita, T., Klein Goldewijk, K., Riahi, K., Shevliakova, E., Smith, S., Stehfest, E., Thomson, A., Thornton, P., van Vuuren, D. P., and Wang, Y. P.: Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands, Clim. Change, 109, 117–161,, 2011. 

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M. M. B., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M., Cambridge University Press, 2013. 

Ivanov, V. Y., Bras, R. L., and Curtis, D. C.: A weather generator for hydrological, ecological, and agricultural applications, Water Resour. Res., 43, W10406,, 2007. 

Jarvis, P. G.: The Interpretation of the Variations in Leaf Water Potential and Stomatal Conductance Found in Canopies in the Field, Phil. Trans. R. Soc. Lond. B., 273, 593–610,, 1976. 

Jothityangkoon, C. and Sivapalan, M.: Framework for exploration of climatic and landscape controls on catchment water balance, with emphasis on inter-annual variability, J. Hydrol., 371, 154–168,, 2009. 

Kergoat, L.: A model for hydrological equilibrium of leaf area index on a global scale, J. Hydrol., 212–213, 268–286,, 1998. 

Kim, J. H. and Jackson, R. B.: A Global Analysis of Groundwater Recharge for Vegetation, Climate, and Soils, Vadose Z. J., 11,, 2012. 

Kirchner, J. W.: Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology, Water Resour. Res., 42, W03S04,, 2006. 

Klimchouk, A. B. and Ford, D. C.: Types of karst and evolution of hydrogeologic setting, in Speleogenesis, Evolution of Karst Aquifers, edited by: Klimchouk, A. B., Ford, D. C., Palmer, A., and Dreybrodt, W., National Speleological Society, Huntsville, Alabama, USA, 45–53, 2000. 

Knohl, A., Schulze, E. D., Kolle, O., and Buchmann, N.: Large carbon uptake by an unmanaged 250-year-old deciduous forest in Central Germany, Agric. Forest Meteorol., 118, 151–167,, 2003. 

Kumar, R., Samaniego, L., and Attinger, S.: Implications of distributed hydrologic model parameterization on water fluxes at multiple scales and locations, Water Resour. Res., 49, 360–379,, 2013. 

Laio, F., Porporato, A., Ridolfi, L., and Rodriguez-Iturbe, I.: On the seasonal dynamics of mean soil moisture, J. Geophys. Res., 107, 4272,, 2002. 

Li, X. Y., Contreras, S., and Solé-Benet, A.: Spatial distribution of rock fragments in dolines: A case study in a semiarid Mediterranean mountain-range (Sierra de Gádor, SE Spain), Catena, 70, 366–374,, 2007. 

Li, X. Y., Contreras, S., Solé-Benet, A., Cantón, Y., Domingo, F., Lázaro, R., Lin, H., Van Wesemael, B., and Puigdefábregas, J.: Controls of infiltration-runoff processes in Mediterranean karst rangelands in SE Spain, Catena, 86, 98–109,, 2011. 

Liang, X., Lettenmaier, D. P., Wood, E. F., and Burges, S. J.: A simple hydrologically based model of land surface water and energy fluxes for general circulation models, J. Geophys. Res., 99, 14415–14428,, 1994. 

Loáiciga, H. A., Maidment, D. R., and Valdes, J. B.: Climate-change impacts in a regional karst aquifer, Texas, USA, J. Hydrol., 227, 173–194,, 2000. 

Lu, Y., Liu, S., Weng, L., Wang, L., Li, Z., and Xu, L.: Fractal analysis of cracking in a clayey soil under freeze-thaw cycles, Eng. Geol., 208, 93–99,, 2016. 

Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., Fernández-Prieto, D., Beck, H. E., Dorigo, W. A., and Verhoest, N. E. C.: GLEAM v3: satellite-based land evaporation and root-zone soil moisture, Geosci. Model Dev., 10, 1903–1925,, 2017. 

Maxwell, R. M. and Condon, L. E.: Connections between groundwater flow and transpiration partitioning, Science, 353, 377–380,, 2016. 

McCabe, M. F., Ershadi, A., Jimenez, C., Miralles, D. G., Michel, D., and Wood, E. F.: The GEWEX LandFlux project: evaluation of model evaporation using tower-based and globally gridded forcing data, Geosci. Model Dev., 9, 283–305,, 2016. 

Mendoza, P. A., Clark, M. P., Barlage, M., Rajagopalan, B., Samaniego, L., Abramowitz, G., and Gupta, H.: Are we unnecessarily constraining the agility of complex process-based models?, Water Resour. Res., 51, 716–728,, 2015. 

Miralles, D. G., Gash, J. H., Holmes, T. R. H., De Jeu, R. A. M., and Dolman, A. J.: Global canopy interception from satellite observations, J. Geophys. Res., 115, D16122,, 2010. 

Miralles, D. G., Holmes, T. R. H., De Jeu, R. A. M., Gash, J. H., Meesters, A. G. C. A., and Dolman, A. J.: Global land-surface evaporation estimated from satellite-based observations, Hydrol. Earth Syst. Sci., 15, 453–469,, 2011. 

Monteith, J. L.: Evaporation and environment, Symp. Soc. Exp. Biol., 19, 205–234, 1965. 

Morris, M. D.: Factorial sampling plans for preliminary computational experiments, Technometrics, 33, 161–174,, 1991. 

Müller Schmied, H., Eisner, S., Franz, D., Wattenbach, M., Portmann, F. T., Flörke, M., and Döll, P.: Sensitivity of simulated global-scale freshwater fluxes and storages to input data, hydrological model structure, human water use and calibration, Hydrol. Earth Syst. Sci., 18, 3511–3538,, 2014. 

Mund, M., Kutsch, W. L., Wirth, C., Kahl, T., Knohl, A., Skomarkova, M. V., and Schulze, E. D.: The influence of climate and fructification on the inter-annual variability of stem growth and net primary productivity in an old-growth, mixed beech forest, Tree Physiol., 30, 689–704,, 2010. 

Owor, M., Taylor, R. G., Tindimugaya, C., and Mwesigwa, D.: Rainfall intensity and groundwater recharge: empirical evidence from the Upper Nile Basin, Environ. Res. Lett., 4, 35009,, 2009. 

Pechlivanidis, I. G., McIntyre, N., and Wheater, H. S.: The significance of spatial variability of rainfall on simulated runoff: an evaluation based on the Upper Lee catchment, UK, Hydrol. Res., 48, nh2016038,, 2016. 

Penman, H. L.: The dependance of transpiration on weather and soil conditions, J. Soil Sci., 1, 74–89,, 1950. 

Pereira, L. S., Allen, R. G., Smith, M., and Raes, D.: Crop evapotranspiration estimation with FAO56: Past and future, Agric. Water Manag., 147, 4–20,, 2015. 

Pérez-Priego, O., Serrano-Ortiz, P., Sánchez-Cañete, E. P., Domingo, F., and Kowalski, A. S.: Isolating the effect of subterranean ventilation on CO2 emissions from drylands to the atmosphere, Agric. Forest Meteorol., 180, 194–202,, 2013. 

Pianosi, F., Sarrazin, F., and Wagener, T.: A Matlab toolbox for Global Sensitivity Analysis, Environ. Model. Softw., 70, 80–85,, 2015. 

Pinty, B., Jung, M., Kaminski, T., Lavergne, T., Mund, M., Plummer, S., Thomas, E., and Widlowski, J. L.: Evaluation of the JRC-TIP 0.01 products over a mid-latitude deciduous forest site, Remote Sens. Environ., 115, 3567–3581,, 2011. 

Porporato, A., Daly, E., and Rodríguez-Iturbe, I.: Soil water balance and ecosystem response to climate change, Am. Nat., 164, 625–632,, 2004. 

Priestley, C. H. B. and Taylor, R. J.: On the Assessment of Surface Heat Flux and Evaporation Using Large-Scale Parameters, Mon. Weather Rev., 100, 81–92,<0081:OTAOSH>2.3.CO;2, 1972. 

Rahman, M. and Rosolem, R.: Towards a simple representation of chalk hydrology in land surface modelling, Hydrol. Earth Syst. Sci., 21, 459–471,, 2017. 

Rambal, S.: Quercus ilex facing water stress: a functional equilibrium hypothesis, in Quercus ilex L. ecosystems: function, dynamics and management, Advances in vegetation science, edited by: Romane, F. and Terradas, F., Springer, Dordrecht, the Netherlands, AIVS, 13, 147–153, 1992. 

Rambal, S.: Le Paradoxe hydrologique des écosystèmes méditerranéens sur des sols karstiques, in: Numéro spécial des Annales de la Société d'Horticulture et d'Histoire Naturelle de l'Hérault, 61–67, 2011. 

Rambal, S., Ourcival, J. M., Joffre, R., Mouillot, F., Nouvellon, Y., Reichstein, M., and Rocheteau, A.: Drought controls over conductance and assimilation of a Mediterranean evergreen ecosystem: Scaling from leaf to canopy, Glob. Chang. Biol., 9, 1813–1824,, 2003. 

Reichstein, M., Tenhunen, J. D., Roupsard, O., Ourcival, J. M., Rambal, S., Miglietta, F., Peressotti, A., Pecchiari, M., Tirone, G., and Valentini, R.: Severe drought effects on ecosystem CO2 and H2O fluxes at three Mediterranean evergreen sites: Revision of current hypotheses?, Glob. Change Biol., 8, 999–1017,, 2002. 

Rimmer, A. and Hartmann, A.: Simplified conceptual structures and analytical solutions for groundwater discharge using reservoir equations, in Water resources management and modeling, edited by: Nayak, P., InTech, Kakinada, 217–238, 2012. 

Rodell, M., Houser, P. R., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C.-J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., Entin, J. K., Walker, J. P., Lohmann, D., and Toll, D.: The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381–394,, 2004. 

Rosero, E., Yang, Z. L., Wagener, T., Gulden, L. E., Yatheendradas, S., and Niu, G.-Y.: Quantifying parameter sensitivity, interaction, and transferability in hydrologically enhanced versions of the Noah land surface model over transition zones during the warm season, J. Geophys. Res.-Atmos., 115, D03106,, 2010. 

Rosolem, R., Gupta, H. V., Shuttleworth, W. J., Gonçalves de Gonçalves, L. G., and Zeng, X.: Towards a comprehensive approach to parameter estimation in land surface parameterization schemes, Hydrol. Process., 27, 2075–2097,, 2013. 

Ross, J.: Radiative transfer in plant communities, in Vegetation and the Atmosphere, volume I Principles, edited by: Monteith, J., Academic Press, London, 13–55, 1975. 

Ruiz, L., Varma, M. R. R., Kumar, M. M. S., Sekhar, M., Maréchal, J.-C., Descloitres, M., Riotte, J., Kumar, S., Kumar, C., and Braun, J.-J.: Water balance modelling in a tropical watershed under deciduous forest (Mule Hole, India): Regolith matric storage buffers the groundwater recharge process, J. Hydrol., 380, 460–472,, 2010. 

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S. (Eds.): Global Sensitivity Analysis, The Primer, John Wiley & Sons Ltd., Chichester, UK, 2008. 

Samaniego, L., Kumar, R., and Attinger, S.: Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale, Water Resour. Res., 46, W05523,, 2010. 

Samaniego, L., Brenner, J., Demirel, C. M., Jing, M., Kaluza, M., Kumar, R., Langenberg, B., Rakovec, O., Schäfer, D., Schrön, M., Schweppe, R., and Thober, S.: The mesoscale Hydrologic Model mHM, Documentation for version 5.9, Helmoltz Centre for Environmental Research (UFZ), Leipzig, Germany, 2018. 

Samuels, R., Rimmer, A., Hartmann, A., Krichak, S., and Alpert, P.: Climate Change Impacts on Jordan River Flow: Downscaling Application from a Regional Climate Model, J. Hydrometeorol., 11, 860–879,, 2010. 

Sarrazin, F., Pianosi, F., and Wagener, T.: Global Sensitivity Analysis of environmental models: Convergence and validation, Environ. Model. Softw., 79, 135–152,, 2016. 

Sarrazin, F., Hartmann, A., Pianosi, P., Rosolem, R., and Wagener, T.: V2Karst version v1.1,, 2018. 

Sauter, M.: Quantification and Forecasting of Regional Groundwater Flow and Transport in a Karst Aquifer (Gallusquelle, Malm, SW. Germany), PhD thesis, Tübinger Universität, Tübinger, Germany, 1992. 

Savenije, H. H. G.: Determination of evaporation from a catchment water balance at a monthly time scale, Hydrol. Earth Syst. Sci., 1, 93–100,, 1997. 

Savenije, H. H. G.: The importance of interception and why we should delete the term evapotranspiration from our vocabulary, Hydrol. Process., 18, 1507–1511,, 2004. 

Scanlon, B. R., Keese, K. E., Flint, A. L., Flint, L. E., Gaye, C. B., Edmunds, W. M., and Simmers, I.: Global synthesis of groundwater recharge in semiarid and arid regions, Hydrol. Process., 20, 3335–3370,, 2006. 

Schwinning, S.: The ecohydrology of roots in rocks, Ecohydrology Bearings – Invited Commentary, Ecohydrology, 3, 238–245,, 2010. 

Seidl, R., Schelhaas, M.-J., Rammer, W., and Verkerk, P. J.: Increasing forest disturbances in Europe and their impact on carbon storage, Nat. Clim. Change, 4, 806–810,, 2014. 

Serrano-Ortiz, P., Kowalski, A. S., Domingo, F., Rey, A., Pegoraro, E., Villagarcía, L., and Alados-Arboledas, L.: Variations in daytime net carbon and water exchange in a montane shrubland ecosystem in southeast Spain, Photosynthetica, 45, 30–35,, 2007. 

Shuttleworth, W. J.: Evapotranspiration, in: Handbook of Hydrology, edited by: Maidment, D. R., McGraw-Hill Inc., New York, 4.1–4.53, 1993. 

Shuttleworth, W. J. (Eds.): Terrestrial Hydrometeorology, John Wiley & Sons Ltd., Chichester, UK, 2012. 

Shuttleworth, W. J. and Wallace, J. S.: Evaporation From Spare Crops – An Energy Combination Theory, Q. J. R. Meteorol. Soc., 111, 839–855,, 1985. 

Simioni, G., Durand-Gillmann, M., and Huc, R.: Asymmetric competition increases leaf inclination effect on light absorption in mixed canopies, Ann. Forest Sci., 70, 123–131,, 2013. 

Šimůnek, J., Šejna, M., Saito, H., Sakai, M., and van Genuchten, M. T.: The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, Version 4.08, University of California Riverside, Riverside, USA, 2009. 

Sitch, S., Smith, B., Prentice, I. C., Arneth, a., Bondeau, a., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., Thonicke, K., and Venevsky, S.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Change Biol., 9, 161–185,, 2003. 

Smith, K. A.: Investigating Uncertainty in Global Hydrology Modelling, PhD thesis, University of Nottingham, Nottingham, UK, 2016. 

Sperna Weiland, F. C., Vrugt, J. A., Van Beek, R. L. P. H., Weerts, A. H., and Bierkens, M. F. P.: Significant uncertainty in global scale hydrological modeling from precipitation data errors, J. Hydrol., 529, 1095–1115,, 2015. 

Stewart, J. B.: Modelling surface conductance of pine forest, Agric. Forest Meteorol., 43, 19–35,, 1988. 

Sutanudjaja, E. H., van Beek, L. P. H., de Jong, S. M., van Geer, F. C., and Bierkens, M. F. P.: Large-scale groundwater modeling using global datasets: a test case for the Rhine-Meuse basin, Hydrol. Earth Syst. Sci., 15, 2913–2935,, 2011. 

Taylor, R. G., Todd, M. C., Kongola, L., Maurice, L., Nahozya, E., Sanga, H., and MacDonald, A. M.: Evidence of the dependence of groundwater resources on extreme rainfall in East Africa, Nat. Clime Chang., 3, 374–378,, 2013. 

Tesemma, Z. K., Wei, Y., Peel, M. C., and Western, A. W.: The effect of year-to-year variability of leaf area index on Variable Infiltration Capacity model performance and simulation of runoff, Adv. Water Resour., 83, 310–322,, 2015. 

Thornthwaite, C. W.: An Approach toward a Rational Classification of Climate, Geogr. Rev., 38, 55–94,, 1948. 

Tritz, S., Guinot, V., and Jourde, H.: Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model, J. Hydrol., 397, 250–262,, 2011. 

Twine, T. E., Kustas, W. P., Norman, J. M., Cook, D. R., Houser, P. R., Meyers, T. P., Prueger, J. H., Starks, P. J., and Wesley, M. L.: Correcting eddy covariance flux underestimates over grassland, Agric. Forest Meteorol., 103, 279–300,, 2000. 

Uhlenbrook, S.: Catchment hydrology – a science in which all processes are preferential, Hydrol. Process., 20, 3581–3585,, 2006. 

Valente, F., David, J. S., and Gash, J. H. C.: Modelling interception loss for two sparse eucalypt and pine forests in central Portugal using reformulated Rutter and Gash analytical models, J. Hydrol., 190, 141–162,, 1997. 

Van Beek, R.: Forcing PCR-GLOBWB with CRU data, Utrecht University, the Netherlands, available at: (last access: 30 November 2018), 2008. 

Van Beek, L. P. H. and Bierkens, M. F. P.: The Global Hydrological Model PCR-GLOBWB: Conceptualization, Parameterization and Verification, Report Department of Physical Geography, Utrecht University, Netherlands, available at: (last access: 30 November 2018), 2008. 

Van Dijk, A. I. J. M. and Bruijnzeel, L. A.: Modelling rainfall interception by vegetation of variable density using an adapted analytical model – Part 1: Model description, J. Hydrol., 247, 230–238,, 2001. 

Van Werkhoven, K., Wagener, T., Reed, P., and Tang, Y.: Rainfall characteristics define the value of streamflow observations for distributed watershed model identification, Geophys. Res. Lett., 35, L11403,, 2008. 

Vörösmarty, C. J.: Global change, the water cycle, and our search for Mauna Loa, Hydrol. Process., 16, 135–139,, 2002. 

Vörösmarty, C. J., Moore, B., Grace, A. L., Gildea, M. P., Melillo, J. M., Peterson, B. J., Rastetter, E. B., and Steudler, P. A.: Continental scale models of water balance and fluvial transport: An application to South America, Global Biogeochem. Cy., 3, 241–265,, 1989. 

Vörösmarty, C. J., Federer, C. A., and Schloss, A. L.: Potential evapotranspiration functions compared on US watersheds: implications for global-scale water balance and terrestrial ecosystem modeling, J. Hydrol., 207, 147–169,, 1998. 

Wada, Y., Van Beek, L. P. H., and Bierkens, M. F. P.: Nonsustainable groundwater sustaining irrigation: A global assessment, Water Resour. Res., 48, W00L06,, 2012. 

Wang, K. and Dickinson, R. E.: A review of global terrestrial evapotranspiration: observation, modelling, climatology, and climatic variability, Rev. Geophys., 50, 1–54,, 2012. 

Weiler, M. and McDonnell, J.: Virtual experiments: A new approach for improving process conceptualization in hillslope hydrology, J. Hydrol., 285, 3–18,, 2004. 

Werth, S., Güntner, A., Petrovic, S., and Schmidt, R.: Integration of GRACE mass variations into a global hydrological model, Earth Planet. Sci. Lett., 277, 166–173,, 2009. 

Williams, P. W.: The role of the subcutaneous zone in karst hydrology, J. Hydrol., 61, 45–67,, 1983. 

Williams, P. W. (Ed.): Environmental change and human impact on karst terrains: an introduction, in: Karst Terrains – Environmental changes and human impact, Catena Verlag, Cremlingen-Destedt, Germany, 1–19, 1993. 

Williams, P. W.: The role of the epikarst in karst and cave hydrogeology?: a review, Int. J. Speleol., 37, 1–10,, 2008. 

Yin, J., Porporato, A., and Albertson, J.: Interplay of climate seasonality and soil moisture-rainfall feedback, Water Resour. Res., 50, 6053–6066,, 2014.  

Young, P., Parkinson, S., and Lees, M.: Simplicity out of complexity in environmental modelling: Occam's razor revisited, J. Appl. Stat., 23, 165–210,, 1996. 

Zhang, Z., Chen, X., Ghadouani, A., and Shi, P.: Modelling hydrological processes influenced by soil, rock and vegetation in a small karst basin of southwest China, Hydrol. Process., 25, 2456–2470,, 2011. 

Zhu, Z., Piao, S., Myneni, R. B., Huang, M., Zeng, Z., Canadell, J. G., Ciais, P., Sitch, S., Friedlingstein, P., Arneth, A., Cao, C., Cheng, L., Kato, E., Koven, C., Li, Y., Lian, X., Liu, Y., Liu, R., Mao, J., Pan, Y., Peng, S., Peñuelas, J., Poulter, B., Pugh, T. A. M., Stocker, B. D., Viovy, N., Wang, X., Wang, Y., Xiao, Z., Yang, H., Zaehle, S., and Zeng, N.: Greening of the Earth and its drivers, Nat. Clim. ChangE, 6, 791–795,, 2016. 

Short summary
We propose the first large-scale vegetation–recharge model for karst regions (V2Karst), which enables the analysis of the impact of changes in climate and land cover on karst groundwater recharge. We demonstrate the plausibility of V2Karst simulations against observations at FLUXNET sites and of controlling modelled processes using sensitivity analysis. We perform virtual experiments to further test the model and gain insight into its sensitivity to precipitation pattern and vegetation cover.