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.

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)

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:

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).

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.

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.

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.

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:

VS,i=Vmax,Sinca(1)VE,i=Vmax,EincaKE,i=Kmax,Enc-i+1nca,

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-1→i(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)fc(t)=1-e-kLAI(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)Vcan,max(t)=VcanLAI(t)fc(t).

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

(6)Ic(t)=min⁡P(t)+Ic(t-1)-Ecact(t)fc(t),Vcan,max(t).

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

(7)Ecact(t)=fc(t)min⁡Ecpot(t),P(t)+Ic(t-1),Vcan,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)Tf(t)=max⁡P(t)-Ecact(t)-fc(t)Ic(t)-Ic(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)twet(t)=Ecact(t)fc(t)Ecpot(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−1∘C−1) is the specific heat of the air and is
equal to 1.013×10-3 MJ kg−1∘C−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)rs,can(t)=rstLAI(t)fc(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:

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,i→i+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,i→i+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):

(18)Qepi,i(t)=min⁡Vepi,i(t-1)KE,i,Vepi,i(t-1)+Repi,i(t).

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,
http://www.gip-ecofor.org/f-ore-t/fontBlanche.php (last access: 30 November 2018) for the French 1 site,
http://puechabon.cefe.cnrs.fr/ (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.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.

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)Eact,cor=Rnλ1+HLE.

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).

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.

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:

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

(20)Bias=∑t∈MET(Eact,sim(t)-Eact,cor(t))∑t∈METEact,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.

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.

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).

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 Ve≤Vr≤VS,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)EEij=yx1j,x2j,…,xi-1j,xij+Δ,…xMj-yx1j,x2j,…,xi-1j,xij,…xMjΔ,

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 Ve≤Vr≤VS,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)Pm=30Ip1+Hp.

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.

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.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.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).

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.

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).

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).

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).

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.

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.

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, https://doi.org/10.1098/rstb.1976.0035, 1976.

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, https://doi.org/10.1029/2005WR004362, 2006.

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.

We propose the first large-scale vegetation–recharge model for karst regions (V2Karst), which...