Chloé Ollivier

Naomi Mazzilli

Albert Olioso

Konstantinos Chalikakis

Simon D Carrière

Charles Danquigny

Christophe Emblanch

Karst recharge-discharge semi distributed model to assess spatial variability of flows

Keywords: Karst, Recharge, Rainfall-discharge modelling

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Groundwater resources is affected by global change. Climate change can greatly affect groundwater recharge and hydrodynamics [START_REF] Taylor | Ground water and climate change[END_REF]. Groundwater depletion may be due to a decrease in recharge or exacerbated by human activities that overexploit groundwater [START_REF] Konikow | Groundwater depletion: a global problem[END_REF][START_REF] Leduc | Impacts of hydrological changes in the mediterranean zone: environmental modifications and rural development in the merguellil catchment, central tunisia[END_REF][START_REF] Xanke | Numerical long-term assessment of managed aquifer recharge from a reservoir into a karst aquifer in jordan[END_REF][START_REF] Meixner | Implications of projected climate change for groundwater recharge in the western united states[END_REF]. Groundwater use planning and management may benefit greatly from quantitative information on groundwater recharge rates. Evaluation of the water flux from the atmosphere to the underlying aquifer is difficult because it is due to multiple factors [START_REF] Lerner | Groundwater recharge: a guide to understanding and estimating natural recharge[END_REF][START_REF] De Vries | Groundwater recharge: an overview of processes and challenges[END_REF]. Aquifer recharge depends on climatic conditions (precipitation, temperature), land use (vegetation, cropland), and hydraulic properties of rocks.

When considering a karst aquifer, the evaluation is even more complex because karst systems are highly heterogeneous and anisotropic [START_REF] Bakalowicz | Karst groundwater: a challenge for new resources[END_REF]. Still, it is important to investigate them because they supply water to a quarter of the world's population [START_REF] Ford | Karst Hydrogeology[END_REF]. Karst hydrosystems are a product of morphological evolution due to rock dissolution and the evacuation of dissolved solids by underground flows [START_REF] Bakalowicz | Karst groundwater: a challenge for new resources[END_REF]. Karst hydrodynamic is very heterogeneous. The conduits network ensures rapid transfer through the system. Slower transfer of water also occurs within carbonate rock.

Karst aquifer permeability ranges from extremely large values in karst conduits, between 10 -1 m.s -1 to several m.s -1 [START_REF] Jeannin | Modeling flow in phreatic and epiphreatic karst conduits in the hölloch cave (muotatal, switzerland)[END_REF] to small values of between 10 -6 m.s -1 and 10 -12 m.s -1 in micropores [START_REF] Borgomano | Petrophysics of lower cretaceous platform carbonate outcrops in provence (se france): implications for carbonate reservoir characterisation[END_REF]. Karstogenesis also affects landscape and land surface properties. Atop carbonate rock, soils are usually poorly developed, but karst features such as dolines may promote soil devellopement, thus enabling agriculture.

The transfer of water from the surface to groundwater ranges in velocity depending on variations in soil thickness, lithology, and structure of the vadose zone. Geophysical field investigations have highlighted wide variability of recharge processes in the vadose zone of karst systems [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF][START_REF] Mazzilli | Contribution of magnetic resonance soundings for characterizing water storage in the unsaturated zone of karst aquifers[END_REF][START_REF] Watlet | Imaging groundwater infiltration dynamics in karst vadose zone with long-term ert monitoring[END_REF].

Karst aquifer heterogeneity and its influence on underground flow can be modelled

Preprint submitted to STOTEN and published DOI : 10.1016/j.scitotenv.2019.134368 October 5, 2019 using different approaches based on water transfer equations that range in complexity.

These approaches may be lumped, semi-distributed, or distributed. Lumped models use simple linear or nonlinear equations that directly relate system discharge to precipitation [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF][START_REF] Tritz | Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model[END_REF][START_REF] Caballero | Caractérisation de la recharge des aquifères et évolution future en contexte de changement climatique[END_REF]. Lumped models are usually applied at aquifer scale and they require a limited amount of data. However, by nature, they are unable to assess spatial variability of recharge and flows. Distributed models have high data requirements, which are quite difficult to satisfy for karst systems [START_REF] Scanlon | Can we simulate regional groundwater flow in a karst system using equivalent porous media models? case study, barton springs edwards aquifer, usa[END_REF]. The uncertainty of permeability distribution is a major limitation [START_REF] Scanlon | Impact of land use and land cover change on groundwater recharge and quality in the southwestern US[END_REF][START_REF] Hartmann | Karst water resources in a changing world: Review of hydrological modeling approaches[END_REF]. In particular, knowledge relative to the conduit network is rarely available at watershed scale. Therefore, applications of distributed models for entire karst aquifers remain limited [START_REF] Jukić | Groundwater balance estimation in karst by using a conceptual rainfall-runoff model[END_REF][START_REF] Robineau | A double medium approach to simulate groundwater level variations in a fissured karst aquifer[END_REF]. Whereas the heterogeneity of karst aquifers cannot be described in fully distributed models, semi-distributed models seem capable of shedding light on flow distribution with reduced knowledge requirements (Pardo-Igúzquiza et al., 2018a).

Semi-distributed karst models combine the advantages of lumped models, in particular parsimony, with an assessment of the spatial distribution of principal aquifer properties and principal forcing variables. Two principal sources of heterogeneity can be considered in semi-distributed karst models: i) surface / sub-surface heterogeneity that controls aquifer recharge [START_REF] Andreo | Methodology for groundwater recharge assessment in carbonate aquifers: application to pilot sites in southern Spain[END_REF][START_REF] Hughes | Evaluation of distributed recharge in an upland semi-arid karst system: the west bank mountain aquifer, middle east[END_REF][START_REF] Bailly-Comte | A conceptual semidistributed model of the coulazou river as a tool for assessing surface water-karst groundwater interactions during flood in mediterranean ephemeral rivers[END_REF][START_REF] Malard | A novel approach for estimating karst groundwater recharge in mountainous regions and its application in switzerland[END_REF]Bittner et al., 2018a;Pardo-Igúzquiza et al., 2018b) and ii) the underground karst network that controls flow paths [START_REF] Ladouche | Semi-distributed lumped model of a karst system under active management[END_REF].

The number of parameters and the possibility to determine their values is a common issue in environmental models [START_REF] Ebel | Physics-based hydrologic-response simulation: Seeing through the fog of equifinality[END_REF][START_REF] Bittner | On the relation between parameters and discharge data for a lumped karst aquifer model[END_REF]. This concern is made exacerbated for karst systems, because difficulties in characterizing karst heterogeneity widen the gap between available information and the number of model parameters required for distributed or semi-distributed modelling [START_REF] Hartmann | Karst water resources in a changing world: Review of hydrological modeling approaches[END_REF][START_REF] Parente | Bayesian Calibration and Sensitivity Analysis for a Karst Aquifer Model Using Active Subspaces[END_REF]. Integration of measurable hydraulic parameters into semi-distributed models is also challenging. Semi-distributed models are based on physically sound structures and equations, but their parameters are not direct equivalents of measurable flow properties. However, linking flow properties to model parameters appears to be a promising way to facilitate the identification of parameter values and to improve the predictive capacity of models. Such an approach was proposed by [START_REF] Hartmann | A new approach to model the spatial and temporal variability of recharge to karst aquifers[END_REF], which shows that considering the statistical distribution of surface and epikarst properties (porosity, and hydraulic conductivity) may improve recharge estimates [START_REF] Hartmann | Enhanced groundwater recharge rates and altered recharge sensitivity to climate variability through subsurface heterogeneity[END_REF].

In this study, we address the following questions:

• Can we develop a semi-distributed model of a karst aquifer with parameters dependent on environmental factors (either measurable karst aquifer properties or qualitative indexes)?

• If so, does spatial distribution of environmental factors improve aquifer recharge and discharge modelling and provide a better understanding of associated aquifer behavior?

• Is aquifer discharge depending on recharge location?

We present a new semi-distributed model for karst aquifers, the KaRaMel model (Karst Recharge and discharge Model). Key features of the model are :

• Spatial parallelization of a lumped model whose structure is based on a widely accepted karst conceptual model and that is commonly used to model karst aquifers (e.g. [START_REF] Pinault | Inverse modeling of the hydrological and the hydrochemical behavior of hydrosystems: Characterization of Karst System Functioning[END_REF][START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF][START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systemes karstiques[END_REF][START_REF] Mazzilli | Karstmod: a modelling platform for rainfall-discharge analysis and modelling dedicated to karst systems[END_REF][START_REF] Baudement | Groundwater management of a highly dynamic karst by assessing baseflow and quickflow with a rainfall-discharge model (Dardennes springs, SE France)[END_REF];

• Direct integration of Soil Available Water Capacity (SAWC, Vörösmarty et al., 1989) estimates as a model parameter;

• Indirect integration of intrinsic vulnerability indices. These indices are related to the proportion of preferential underground flows, on the basis of information regarding aquifer properties such as lithology, karst network evidences, tracer tests [START_REF] Zwahlen | Vulnerability and risk mapping for the protection of carbonate (karst) aquifers[END_REF][START_REF] Kavouri | Paprika: a method for estimating karst resource and source vulnerability-application to the ouysse karst system (southwest france)[END_REF].

The test site for KaRaMel is the Fontaine de Vaucluse aquifer. The large watershed area (about 1160 km 2 ) of this aquifer is associated with broad heterogeneity of climate variables and aquifer properties, which makes it suitable for characterizing the impact of different levels of heterogeneity on aquifer flows and discharge. After calibration and validation of the model on this test case, two assumptions on possible modifications of future precipitation distribution are applied to a past precipitation dataset. Such a test makes it possible to characterize aquifer sensitivity to spatial and temporal precipitation distribution.

Model rationale and organisation of the study

KaRaMel, the proposed Karst Recharge and discharge Model, is designed to provide a distributed simulation of karst aquifer recharge and underground flow. Here, we define karst recharge as water flows that penetrate the unsaturated zone and that are not released to the atmosphere by evapotranspiration (Figure 1). Due to the water that transit througth the unstaurated zone, the recharge of unstarurated zone may be different of recharge of saturated zone [START_REF] Caballero | Caractérisation de la recharge des aquifères et évolution future en contexte de changement climatique[END_REF]. The direct observation of recharge is usually lacking in karst systems, while in many cases discharge through karst springs is recorded. Aquifer flow processes are generally considered to account for spring discharge, which is thus valuable information for the evaluation of karst models. KaRaMel simulates aquifer discharge by distributing recharge and partitioning it into rapid and slow flow pathways. The recharge distribution is evaluated against discharge observations. Recharge depends primarily on precipitation [START_REF] Fiorillo | A model to simulate recharge processes of karst massifs[END_REF][START_REF] Malard | A novel approach for estimating karst groundwater recharge in mountainous regions and its application in switzerland[END_REF], evapotranspiration [START_REF] Jasechko | The pronounced seasonality of global groundwater recharge[END_REF][START_REF] Doble | Review: Current and emerging methods for catchmentscale modelling of recharge and evapotranspiration from shallow groundwater[END_REF], to surface properties that promote runoff or infiltration and to the soil ability to hold water [START_REF] Kim | A global analysis of groundwater recharge for vegetation, climate, and soils[END_REF][START_REF] Ries | Recharge estimation and soil moisture dynamics in a mediterranean karst aquifer[END_REF] (Figure 1). All these recharge determinants vary in space. Therefore, the first layer of KaRaMel is spatially distributed and accounts for distributed aquifer-atmosphere exchange. The soil reservoir is filled by precipitation and enables evapotranspiration. Excess precipitation percolates toward the underground aquifer. Reservoir capacity is usually calibrated [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF][START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systemes karstiques[END_REF]Chang et al., 2017, e.g.). In KaRaMel, no calibration was done and reservoir capacity was defined according to the Soil Available Water Capacity (SAWC, Vörösmarty et al., 1989). SAWC represents the maximum amount of water that can be lost to evapotranspiration. We determined SAWC for each model grid cell from the information contained in available soil surveys using digital soil mapping methodology [START_REF] Mcbratney | On digital soil mapping[END_REF].

The methodology for mapping SAWC is presented in section 3 and the resulting map in section 4.

Underground flows and stocks are modelled on the basis of a semi-distributed model.

Flow parameters for rapid and slow pathways are calibrated but we constrain their distribution through the distribution of intrinsic vulnerability indices. Indeed, intrinsic vulnerability summarises aquifer properties that promote or delay flow paths through the unsaturated zone [START_REF] Zwahlen | Vulnerability and risk mapping for the protection of carbonate (karst) aquifers[END_REF]. The used intrinsic vulnerability method is PaPRIKa [START_REF] Dörfliger | Cartographie de la vulnerabilité des aquifères karstiques. guide méthodologique de la méthode paprika[END_REF]. Calibration of model parameters is based on the comparison between observed and simulated discharges (section 5). Consideration of these two parameter distribution levels makes it possible to simulate the spatial variability of hydrosystem stocks and the discharge time series. Daily recharge distribution over the last ten years is simulated with KaRaMel. Daily rainfall locations over the last ten years are partially modified to determine the sensitivity of the aquifer response to the recharge location.

The KaRaMel model

Model description

KaRaMel specifications are the following.

• The model should simulate spatial distribution of recharge at a short time scale (e.g.: daily) in order to represent the fast response of the karst system to precipitation [START_REF] Ollivier | Contribution of Hydrogeological Time Series Statistical Analysis to the Study of Karst Unsaturated Zone (Rustrel, France)[END_REF].

• The model should use easily accessible data to facilitate its application. model units have the same three-reservoir structure (Figure 2-a and Table 1), which is commonly used to model karst aquifers [START_REF] Pinault | Inverse modeling of the hydrological and the hydrochemical behavior of hydrosystems: Characterization of Karst System Functioning[END_REF][START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systemes karstiques[END_REF][START_REF] Schmidt | Characterisation and modelling of conduit restricted karst aquifers -Example of the Auja spring, Jordan Valley[END_REF][START_REF] Baudement | Groundwater management of a highly dynamic karst by assessing baseflow and quickflow with a rainfall-discharge model (Dardennes springs, SE France)[END_REF][START_REF] Schmidt | Characterisation and modelling of conduit restricted karst aquifers -Example of the Auja spring, Jordan Valley[END_REF]. The upper reservoir (S) represents soil behaviour. Upper reservoir collects precipitation, hosts evapotranspiration, and governs infiltration. These flows are controlled by the reservoir holding capacity (SAWC). Infiltration (I) occurs when the amount of water in the soil exceeds SAWC.

Evapotranspiration (ET) is equal to climatic demand (reference evapotranspiration ETo) as long as there is enough water in the soil reservoir. Water that percolates from the soil reservoir is divided into two flows. One flow goes to a slow discharge reservoir (M) and the second goes to a rapid discharge reservoir (C). The rapid discharge reservoir may represent a karst conduit network. The slow discharge reservoir represents the less altered part of the aquifer that mostly corresponds to the limestone matrix. Reduction of the natural complexity of flows into two end-members (slow, matrix flows and fast, conduit flows) is a common and successful approximation in karst hydrology (e.g. [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF][START_REF] Mazzilli | Karstmod: a modelling platform for rainfall-discharge analysis and modelling dedicated to karst systems[END_REF] Aquifer discharge (Qsim) is the sum of discharge from the slow and rapid reservoirs of all model units.

The model grid should be adapted to the studied aquifer and respect the following assumptions.

• The grid should be coarse enough to respect the assumption of the absence of water exchange between model grids. Runoff may thus occur, but its characteristic length should be less than unit length.

• The grid should be fine enough to assume homogeneous landscape properties within the grid.

Governing equations

Flow within each model unit i (i=1,...n, where n is the total number of model units) is governed by the following budget equations:

ET i =    ET o i , if Hs i + P i ≥ ET o i Hs i + P i , if Hs i + P i < ET o i (1)
I i =    P i -ET i , if HS i + P i -ET i ≥ SAW C i 0, if Hs i + P i -ET i < SAW C i (2) dHs i dt =    P i -ET i -I i , if Hs i ≥ SAW C i max(P i -ET i , 0), if Hs i < SAW C i (3) dHm i dt = x i × I i -qm i (4) dHc i dt = (1 -x i ) × I i -qc i (5)
Where P i is incoming rainfall (in L.T -1 ), ET i is evapotranspiration, and ETo i the reference evapotranspiration (in L.T -1 , eq.1). As long as the water level of the soil reservoir Hs i (in L) is lower than the threshold value of Soil Available Water Capacity (SAWC i in L), infiltration is null. Infiltration occurs when the water level of the soil reservoir exceeds the SAWC i (eq.2). I i is the infiltration from the soil reservoir to the lower reservoirs (in L.T -1 ). The coefficient x i (dimensionless, between 0 and 1) represents the infiltration portion that moves towards the reservoir M i (eq.4). Thus, 1-x i represents the infiltration part that moves toward the reservoir C i . Hs i , Hm i and Hc i (in L) are the water reservoir levels of S i , M i and C i , respectively. qm i and qc i are the discharges per unit of area of reservoirs M i and C i (both in L.T -1 ), defined as:

qm i = km i × Hm i ( 6
)
qc i = kc i × Hc i (7)
where km i and kc i are reservoir-specific discharge coefficients (in T -1 ). The total unit outflow Q i is the sum of slow reservoir and rapid reservoir flows (eq.8). Aquifer discharge Q sim (in L.T -1 , eq.9) is the sum of all model unit outflows scaled by the unit area (A i in L 2 ) relative to the watershed area (A, in L 2 ):

Q i = qm i + qc i (8) Q sim = n i (Q i * Ai A ) (9)
To prevent long-term error accumulation, model equations are solved using the analytical solution instead of Euler solution used by [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF].

Model inputs, parameters, output variables, and calibration strategy

Model inputs are distributed precipitation and reference evapotranspiration. This information may be derived from re-analyses of atmospheric conditions computed by major meteorological centers such as METEO-FRANCE (SAFRAN, Quintana-Seguí et al., 2008), and European Centre for Medium-Range Weather Forecasts (ERA-interim, ERA-5, WEB1). Reference evapotranspiration can be easily estimated from gridded meteorological variables in the above-mentioned re-analyses by using various equations such as that of Thornthwaite (e.g. [START_REF] Quesada-Montano | Can climate variability information constrain a hydrological model for an ungauged Costa Rican catchment?[END_REF] or Penman-Monteith (e.g. [START_REF] Guo | Impact of evapotranspiration process representation on runoff projections from conceptual rainfall-runoff models[END_REF]. Several gridded precipitation products are also available based on meteorological re-analysis and/or remote sensing data (e.g. MSWEP in [START_REF] Beck | MSWEP: 3-hourly 0.25° global gridded precipitation (1979–2015) by merging gauge, satellite, and reanalysis data[END_REF]), CHIRPS in Funk et al. (2015)). A gridded dataset offers the advantage of being available over large areas with spatial resolution between 0.05 • and 0.20 • .

The KaRaMel semi-distributed model is composed of n units (Figure 2). Each unit has four parameters (SAWC i , x i , km i and kc i ) that must be distributed over the entire propose a parsimonious model and to limit the number of parameters to be calibrated, the parameter values are spatially constrained according to the commonly available aquifer characteristics. Indeed, expansion of digital technology, soil data collection, spatial data management, and methods to describe soil patterns provide opportunities for predicting soil properties using readily available products such as soil property maps produced by the European Soil Data Centre [START_REF] Ballabio | Mapping topsoil physical properties at european scale using the lucas database[END_REF]WEB2) or GlobalSoilMap [START_REF] Sanchez | Digital soil map of the world[END_REF]. Here we propose to use a similar methodology for distributing SAWC i (detailed in section 4.4). We also use vulnerability map patterns to constrain x i and kc i parameters, based on the assumption that vulnerability maps of karst aquifers are qualitative measurements of the importance of rapid water transfers [START_REF] Kavouri | A coupled groundwater-flow-modelling and vulnerability-mapping methodology for karstic terrain management[END_REF]. Uniform values are assumed for km. Calibration is performed against spring discharge at the aquifer outlet. However, spatially variable inner flows and saturation states may also be considered as valuable model output.

Application and evaluation of KaRaMel

Test site

The Fontaine de Vaucluse karst system is located in the south-eastern France (Fig-

ure 3). The area has a Mediterranean climate: a hot summer, mild winter, and most precipitation falls in spring and autumn (Figure 4). The majority of precipitation is produced by storms resulting from air streams that arrive from the Mediterranean Sea. The average monthly temperature is 2 • C for the coldest month (February) and 20 • C for the warmest month (August). From 2004 to 2015, annual precipitation ranges from about 660 to 1260 mm.y -1 , with an average of 960 mm.y -1 , daily precipitation ranges from 0 to 80 mm.d -1 .

The area is mainly covered by forest and bush (85%), 14% of the watershed is used for farming and 1% is built environment (after the land use map of Provence-Alpes-Côted'Azur area of 2014, derived from remote sensing data, WEB3). The aquifer is made up

of Cretaceous limestone with a potential thickness of 1500 m [START_REF] Blavoux | Bilan, fonctionnement du système karstique de la fontaine de vaucluse (sud-est de la france)[END_REF]. The Cretaceous limestone constitutes an important topographic unit delimited at the north by a mountain chain that contains Mont Ventoux (1910 m a.s.l.) as the summit and bordered on the south by the Apt syncline (Figure 3). Eastern and western boundaries correspond to the Rhône and Durance river valleys. For further insights the regional geological map with the fault network is presented in Appendix A.14. The impluvium is large, with an area of 1162 km 2 (after Puig, 1989, updated with GIS). The area has a mean elevation of 880 m a.s.l. There is no permanent river over the watershed area. The main outlet of the aquifer is the Fontaine de Vaucluse spring (84 m a.s.l.). The mean thickness of the unsaturated zone is about 800 m. Due to thickness and lithology, the unsaturated zone plays an important hydrological role as a buffer stock of water [START_REF] Puig | Le système karstique de la Fontaine de Vaucluse[END_REF][START_REF] Emblanch | Carbon 13 of tdic to quantify the role of the unsaturated zone: the example of the vaucluse karst systems (southeastern france)[END_REF]. The unsaturated zone supplies on average 55% of the discharge over the hy-drological cycle, and up to 80% during flooding [START_REF] Garry | Etude des processus d'écoulements de la zone non saturée pour la modélisation des aquifères karstiques[END_REF]. The Fontaine de Vaucluse aquifer has a mature karst network, with well-developed conduits and a large component of drowned karst downstream [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] El-Hakim | Significance and origin of very large regulating power of some karst aquifers in the Middle East. Implication on karst aquifer classification[END_REF]. The development of a major karst network several hundred metres below the actual groundwater table is attributed to the Messinian crisis [START_REF] Audra | The effect of the Messinian Deep Stage on karst development around the Mediterranean Sea. Examples from Southern France[END_REF][START_REF] Gilli | Les lithophages pliocènes de la fontaine de Vaucluse (Vaucluse, France). Un argument pour une phase messinienne dans la genèse du plus grand karst noyé de France[END_REF].

Tracer tests evaluating mass transfer through the saturated zone from different locations show a range of maximum flow velocity between 12 and 208 m.h -1 [START_REF] Couturaud | Hydrogéologie de la partie occidentale du système karstique de Vaucluse (karstification et aquifère sous couverture[END_REF].

Rapid pressure transfer through the saturated zone is demonstrated by the 6 hour delay that occurs between spring discharge changes and piezometric head changes observed at a well located 30 km east of the Fontaine de Vaucluse spring (LSBB site in Figure A.14).

Rapid pressure transfer occurring through the aquifer is demonstrated by a maximum precipitation-discharge cross-correlation for delays ranging from 1 to 6 days depending on rainfall event characteristics and the aquifer's saturation state [START_REF] Ollivier | Contribution of Hydrogeological Time Series Statistical Analysis to the Study of Karst Unsaturated Zone (Rustrel, France)[END_REF].

Fontaine de Vaucluse discharge records

Fontaine de Vaucluse discharge has been accurately monitored since November 2003 by a gauging station installed in the river, 450 m downstream of the spring [START_REF] Ollivier | Caractérisation et spatialisation de la recharge des systèmes karstiques. Application au système de la Fontaine de Vaucluse[END_REF]. Water levels are recorded every hour and converted into hourly discharge using the station's calibration curve. Hourly data is used for computing the spring's daily discharge. Records of Fontaine de Vaucluse discharge prior to 2003 are mainly instantaneous measurements of the spring's water levels. These records are specific measurements collected by an observer on a nearly daily basis [START_REF] Cognard-Plancq | Historical monthly rainfallrunoff database on fontaine de vaucluse karst system: review and lessons, in: Karst, cambio climatico y aguas submediterraneas[END_REF]. Analysis of hourly gauging station records since 2003 shows that discharge has varied from day to day over 80% of this period [START_REF] Ollivier | Caractérisation et spatialisation de la recharge des systèmes karstiques. Application au système de la Fontaine de Vaucluse[END_REF]. Thus, older records may be used for qualitative characterisation but not for quantitative studies.

In this study we use daily discharge data from September 2004 to August 2015 (Figure 4). Over this period, discharge ranged from 2.8 to 63 m 3 .s -1 , with an average of 15 m 3 .s -1

and a median of 10 m 3 .s -1 .

Precipitation and reference evapotranspiration dataset

Meteorological data were extracted from the SAFRAN meteorological dataset provided by METEO-FRANCE, the French national meteorological service. SAFRAN is a meteorological re-analysis system combining meteorological observations with estimates provided by numerical weather models [START_REF] Quintana-Seguí | Analysis of Near-Surface Atmospheric Variables: Validation of the SAFRAN Analysis over France[END_REF][START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the Safran system[END_REF].

The main objective of SAFRAN is to produce an accurate estimate of meteorological variables for climate and hydrologic analysis. Daily data of the primary relevant atmospheric parameters (i.e., air temperature, wind speed, humidity, precipitation, and incident radiation) have been available for France since 1958, with a mesh of 8 by 8 km (also available over Switzerland and Spain). The SAFRAN product includes reference evapotranspiration calculated using the Penman-Monteith equation.

The Fontaine de Vaucluse watershed area is covered by 34 SAFRAN cells (Figure 3).

Daily precipitation varied from 0 to 90 mm.d -1 (Figure 4) and reference evapotranspiration from 0.1 to 9 mm.d -1 .

Soil Available Water Capacity (SAWC)

SAWC is generally defined as the water amount intermediate between field capacity and wilting point that is held in the soil. This information is difficult to obtain over large areas and in particular in karst environments. Karst surfaces (epikarst or soil) are complex environments where soil matrix is mixed with rock elements and there are notable spatial variations; in this medium, root systems develop differently from place to place [START_REF] Williams | The role of the epikarst in karst and cave hydrogeology: a review[END_REF][START_REF] Bakalowicz | Epikarst, in: Encyclopedia of Caves[END_REF][START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF]. Because no detailed SAWC map was available for hydrological modelling for our area, we set up a specific procedure to build a new map based on local estimates and principles of digital soil mapping. The proposed methodology consists of: (1) creating a set of SAWC values from soil pit observations, (2) identifying SAWC spatial covariates, (3) establishing the best mathematical equation to predict SAWC.

Local SAWC estimates

Local estimations of SAWC were obtained from soil pit information by considering the description of soil horizons in terms of depth, texture, and rock content and by using pedotransfer functions. We computed SAWC from soil layer thickness (d j in m, j being the soil layer number), fraction of stony coarse elements in each layer (CC j ), and an estimate of the water holding capacity of the soil layer matrix (WHC j in m 3 .m -3 ) in each layer (eq. 10). We derived the latter from the pedotransfer table provided by [START_REF] Baize | Guide pour la description des sols. 2ème édition revue et corrigée[END_REF] which uses Jamagne's soil textural triangle [START_REF] Jamagne | Quelques données sur la variabilité dans le milieu naturel de la rèserve en eau des sols[END_REF][START_REF] Bruand | Use of class pedotransfer functions based on texture and bulk density of clods to generate water retention curves[END_REF][START_REF] Bruand | Estimation des propriétés de rétention en eau des sols à partir de la base de données solhydro: une première proposition combianant le type d'horizon, sa texture et sa densité apparente[END_REF]) with 14 soil classes (Table 2). A water holding capacity was established according to the texture of fine materials in each soil layer. Total SAWC was computed by summing the values for all n soil layers:

SAW C = n j=1 W HC j * d j * (1 -CC j ) (10)
Information from 323 pits over the impluvium were extracted from two soil databases (Figure 5). The first database was provided by the French National Forest Inventory (IFN). This database mainly covers natural lands and forests [START_REF] Morneau | Les forêts mélangées en France métropolitaine[END_REF]. It contains 261 experimental points that were acquired between 2005 and 2013 over the study area following a mostly regular grid. The second database included information on 62 soil pits spread over agricultural areas. These data were available from the French National Soil Inventory [START_REF] Grolleau | Le système d'information national sur les sols : DONESOL et les outils associés[END_REF].

SAWC mapping function

Local SAWC estimations have to be extrapolated to build a soil map of the area.

Because surface properties are heterogeneous over the impluvium, a linear extrapolation of soil observations was not suitable for representing the spatial distribution of SAWC.

Digital soil mapping approaches assume that it is possible to identify factors that produce soil patterns [START_REF] Pachepsky | Artificial neural networks to estimate soil water retention from easily measurable data[END_REF][START_REF] Grunwald | Multi-criteria characterization of recent digital soil mapping and modeling approaches[END_REF], so that local estimates over soil pits can be related to descriptors spatially available over the watershed (e.g. climate, vegetation, morphology). To predict soil properties the SCORPAN approach (eq. 11,

McBratney
The review of soil modelling methods led by [START_REF] Grunwald | Multi-criteria characterization of recent digital soil mapping and modeling approaches[END_REF] showed that depending on the study site not all factors were suitable for soil property mapping. With its plateaus, mountains and canyons, the Fontaine de Vaucluse impluvium presents a wide diversity of landscapes (Figure 3) over a limited geographical extent; as a result we retained only The SAWC model. Various quantitative methods ranging from statistical, geostatistical, and hybrid methods to mechanistic simulation models can be used to predict soil properties. The most popular methods are based on regressions, classification/discrimination methods, kriging, and tree-based models. Other methods such as GIS-based explicit modelling, neural networks, fuzzy logic-based models, and stochastic simulations are also used but to a lesser extent [START_REF] Grunwald | Multi-criteria characterization of recent digital soil mapping and modeling approaches[END_REF].

In this study, we tested two types of mathematical models to relate SAWC to R and O factors: (1) a multilinear regression model (LM) and (2) an artificial neural network (ANN). ANN models are more flexible than linear models because they are able to catch non-linear behaviours without a priori knowledge of the shape of the relation. The robustness of the ANN models was verified by using a split sampling test [START_REF] Klemeš | Operational testing of hydrological simulation models[END_REF]. This [START_REF] Venables | Modern Applied Statistics with S. Fourth ed[END_REF]. In both cases, two different models were developed, one for natural land areas and another for agricultural areas. Modelling natural and agricultural areas together provided lower performances for both types of models.

Infiltration distribution coefficients

Vulnerability is an expression of aquifer geological and hydrogeological properties that define its susceptibility to surface pollution, independently of the nature of a contaminant [START_REF] Margat | Vulnérabilité des nappes d'eau souterraine à la pollution[END_REF][START_REF] Vrba | Continental scale models of water balance and fluvial transport: An application to South America[END_REF][START_REF] Foster | Groundwater protection: the science and practice of land surface zoning[END_REF]. Because the vulnerability map identified aquifer properties that may promote rapid infiltration, we used it to constrain the distribution of model parameters that control rapid flows: partition coefficients x i and discharge coefficients kc i . Discharge coefficients of slow reservoirs km i could not be related to vulnerability or any other qualitative indeces and a uniform value was used over the aquifer. The vulnerability indices of the Fontaine de Vaucluse watershed were mapped using the PaPRIKa method [START_REF] Dörfliger | Cartographie de la vulnerabilité des aquifères karstiques. guide méthodologique de la méthode paprika[END_REF] by Ollivier et al. (2019a,b) at 50 m spatial resolution and upscaled to the model resolution (1 km 2 ). The resulting vulnerability index distribution is presented in Figure 6 : 1. "Low" over 13% of the area, mainly in the steepest zones such as mountainsides.

2. "Moderate" over 48% of the area, mainly gentle slopes, surface karst features, and secondary fault networks.

3. "High" over 36% of the area, mainly steeper slope zones with thin soil.

4. "Very high" over 3% of the area due to major fault networks and geological discontinuities.

Because the "low" and "very high" classes represent a small part of the watershed they were combined with the "moderate" and "high" counterparts, ending up with two vulnerability classes : a) a less vulnerable part composed of areas with "low" and "moderate" indices and b) a more vulnerable part made up of areas with "high" and "very high" indices. Parameters x i and kc i are assumed to be homogeneous within each vulnerability zone. We denote by x L , kc L and x H , kc H the calibrated x i , kc i values for the less vulnerable and more vulnerable parts of the system, respectively.

Calibration and evaluation

Calibration was performed using the Particle Swarm Optimization algorithm [START_REF] Clerc | Particle swarm optimization[END_REF]; Zambrano-Bigiarini and Rojas, 2013, PSO). Parameter sampling followed a Latin hypercube procedure. We used a single objective function, the Kling Gupta efficiency (KGE) after [START_REF] Gupta | Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling[END_REF]. KGE is a multi-objective function with a limited high discharge bias compared to the usual Nash and Sutcliffe Efficiency (NSE, [START_REF] Nash | River flow forecasting through conceptual models part i -a discussion of principles[END_REF]. The KGE function images the Euclidian distance from the ideal point in the scaled space of the three components: r, β and α. r is the correlation coefficient, β is the ratio between mean simulated flow and mean observed flow and αis the ratio between standard deviation of simulated and observed values.

KGE = 1 -(r -1) 2 + (α -1) 2 + (β -1) 2 (12)
The set of parameters (x L , x H , kc L , kc H , km) providing the best performance with respect to the KGE objective function was retained. Model performances were evaluated with KGE, Root Mean Square Error (RMSE), and NSE. A 2-year warm-up period, from September 2004 to August 2006 (Figure 4), was used to avoid calibration bias due to initial conditions, which may be significant for lumped models [START_REF] Mazzilli | Sensitivity and uncertainty in karst hydrosystem modelling[END_REF] Moreover, discharge simulation with KaRaMel was compared to discharge simulation with the KDM model which is a lumped model of the Fontaine de Vaucluse aquifer proposed by [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF]. KDM has the same structure as the KaRaMel unit.

N SE = (Qo -Qs) 2 (Qo -Qo) 2 (13) RM SE = 1 n (Qo -Qs) 2 (

Scenarios of modification of rainfall distribution

The major goal of modelling the distribution of aquifer recharge is the ability to forecast the response of the aquifer to location of recharge.To identify the sensitivity of recharge distribution to temporal and spatial precipitation distribution, and its impacts on aquifer discharge, we modified the precipitation pattern over the last ten years. We tested two modifications of the precipitation pattern:

• rainfall events of the western part are shifted toward the eastern part of the watershed,

• rainfall events are shifted from lower elevation areas to higher (> 1000 m. a.s.l.).

These two scenarios are gradually applied; first modification involved only the location of the high intensity rainfall event. Then we decreased the threshold value of modified rainfall events. Finally, the scenario involved all rainfall events of more than 5mm.d -1 .

Results

SAWC mapping

The analysis of soil available water capacity (SAWC) and covariates variance led to the identification of principal explicative covariates. The covariates used for agricultural areas are the curvature and roughness of the surface, elevation, western exposure and the vegetation index. The natural areas covariates are latitude, longitude, elevation, north or south mountain side, mountain amplitude, curvature and vegetation index.

The SAWC distribution was calculated using two types of regression models (multilinear or neural network) that were trained on soil pit information and spatial covariates following a local application of the SCORPAN method. The SAWC was distributed over the watershed on a 1 km 2 grid. A regression model was established separately over natural and agricultural areas (Figure 7). They were evaluated over a part of the dataset that was not used for training. Multilinear regression models did not correctly simulate SAWC observations (r 2 of 0.1), whereas ANN models were much more efficient (r 2 from 0.6 to 0.9). The ANN model of natural areas had nine neurons with one blind layer, which enabled a good simulation with an r 2 of 0.58 and an RMSE of 17 mm (Figure 7). The ANN model of agricultural areas had seven neurons and one blind layer; it explained the full range of SAWC variation (about 0 to 250 mm) with an r 2 of 0.98 and an RMSE of 13 mm. The two ANN models were used for mapping SAWC over the watershed with a 1km 2 spatial resolution (Figure 8). The highest SAWC values were located over the central plateau, in particular in the Sault graben where soil developed on Oligocene fill and over flat areas with a high density of dolines filled by accumulations of clay, where agricultural activities are well developed. Mountainsides with slope gradients and natural vegetation had lower SAWC values. There, Cretaceous limestones were subject to high levels of denudation and large parts of the surface were covered by rocks [START_REF] Thomas | Limited influence of climatic gradients on the denudation of a Mediterranean carbonate landscape[END_REF].

Model calibration results

Calibrated values are given in Table 3 and calibration performances in Table 4. The The calibrated values of coefficient x i are consistent with the vulnerability indices;

the higher the vulnerability of the area, the higher the fast flow ratio. Rapid flows are 24% higher in high vulnerability areas, as compared to moderate vulnerability areas.

Lower kc H than kc L means that rapid transfer is more rapid for areas with moderate vulnerability. However, kc H and kc L have the same order of magnitude. The calibrated value of the slow reservoir discharge coefficient is two orders of magnitude lower than kc i .

This value is strongly constrained by low flow periods because analysis of Fontaine de

Vaucluse recession curves provided recession coefficients between 0.005 d -1 and 0.006 d -1 , which are on the same order of magnitude [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF].

KaRaMel performances obtained over calibration and validation periods are given in

Sensitivity of the model to the distribution of its properties

The step-by-step procedure highlights the impact of the distribution of the three heterogeneity sources on hydrosystem discharge. Figure 9-a shows various model setups according to the introduction of distributed information into the model. Model setup number 1 corresponds to the semi-distributed model with homogeneous values only for all parameters (SAWC i , x i , kc i ) and all input variables (P i and ETo i ). This is equivalent to the KDM [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF]. Model set number 8 corresponds to the KaRaMel model with all variables and parameters spatially distributed. Figure 9-b presents model performances; the best performances are located near the outside of the radar chart.

The best performances were those of the completely distributed KaRaMel inputs (see also Table 4). Less optimal performances were obtained when no parameters are distributed (model number 1). This model exhibits in the same performance as the KDM model calibrated as presented here. The overall performance of KaRaMel (expressed by KGE) is strongly dependent on the distribution of rapid flow parameters. The lowest KGE values were obtained when x i and kc i were not distributed. In detail, performances of high and low flow periods improved only when the rapid flow parameter distribution was combined with the distribution of another heterogeneity source: either climate variables (model number 6), or SAWC (model number 4). The use of spatialized meteorological variables and the distribution of SAWC had a greater impact on low flow than on high flow periods. Improvements were significant when two sources of variability were accounted for, in particular if one of them was SAWC.

Aquifer responses to recharge distribution

We tested two scenarios for modifying the distribution of precipitation over the watershed. Figure 10 Differences between the simulated discharge with KaRaMel and observed discharges (in grey). Trend lines of differences between simulated and observed discharges are computed with a moving average.

Differences between simulated discharge with the KDM and observed discharge is represented by the blue trend line.

Discharge simulations

The discharge simulation obtained using the parameters retrieved in the previous section is shown in the upper part of Figure 11. The statistical evaluation in Table 4 highlights the excellent quality of simulated discharge, particularly for high flow periods (RMSE of 6.7 m 3 .s -1 ). However, comparison of the discharge time series shows that the first high flow peaks that occurred in autumn after the long summer recession period is overestimated. The difference between simulated discharge and observation, presented in the lower part of Figure 11, confirms a discharge overestimation trend at the beginning of the wet season (autumn). It is followed by an underestimation of high flow peaks during the remaining part of the wet season. The cycle of overestimation and underestimation of discharge was already visible for simulations by the 1D model (blue line), which probably indicates that this trend is due to the lumped structure of the model. Underestimates and overestimates tend to balance each other over one hydrologic year, so that model performances were good at the annual scale. The absolute difference between simulated and observed mean annual discharge ranges only from 6 to 22 mm.y -1 , while annual discharges range from 288 to 571 mm.y -1 .

Temporal distribution of recharge

Daily reference evapotranspiration and precipitation from September 1958 to August 2015 were used to simulate annual recharge variations. Precipitation, simulated evapotranspiration, and recharge are presented in Figure 12 at annual scale (upper part). Note that annual quantities were computed for the hydrologic year, from September to August.

Annual recharge and evapotranspiration are strongly related to annual precipitation.

Annual precipitation varied greatly from one year to the next; evapotranspiration is considerably more stable than both precipitation and recharge (Figure 12). This implies that recharge-precipitation and evapotranspiration-precipitation ratios changed over hydrological years (Figure 12, lower part). Recharge-precipitation ratio ranged from about 20% to 60%. Figure 12 shows that the higher the annual precipitation, the higher the recharge and the ratio of precipitation/recharge. This ratio depends on temporal precipitation distribution. When winter precipitation represents more than 60% of annual precipitation, annual recharge is higher. Conversely, if winter precipitation represents less than 40% of annual precipitation, the ratio is much lower.

The evapotranspiration-precipitation ratio also follows the annual precipitation trend: the higher the precipitation, the lower the ratio. This ratio also depends on the temporal precipitation distribution. When most annual precipitation occurs during winter, the ratio is lower than when most precipitation occurs during spring and summer.

The evapotranspiration-precipitation ratio and the recharge-precipitation ratio change yearly, depending first on the annual precipitation amount and then on the temporal precipitation distribution. This highlights the important role of evapotranspiration. When precipitation occurs during spring or summer, climatic demand is higher than during winter, resulting in decrease of recharge.

Discussions and perspectives

Model performances

In this article, we presented a new model for the description of water flow in karst systems. This model is a semi-distributed model which consist in three main part:

• the upper reservoir corresponding to soil and epikarst, which is the place for evapotranspiration and recharge,

• a rapid flow pathway which represents the main pathway related to the conduit network,

• a slow flow pathway which represent the transfer through the calcareous matrix.

All these compartments are spatially distributed by including as much as possible distributed information. This includes soil characteristics, elevation, remote sensing of vegetation activity, climate, geology and hydrogeology properties.

The semi distributed model was tested over the Fontaine de Vaucluse watershed showing satisfactory results in the simulation of spring outflow and improvement compared to the KDM model (lumped model previously developed by [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF]. In a previous study, [START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systemes karstiques[END_REF] evaluated the capacity of four lumped models to simulate Fontaine de Vaucluse spring flow: GR4J (4 parameters, [START_REF] Perrin | Improvement of a parsimonious model for streamflow simulation[END_REF], KDM (4 parameters, [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF], HBV-6p (6 parameters, Bergström and Forsman, 1973) and TOPMO (modified version of TOPMODEL of Beven and Freer, 2001, , 6 It is difficult to compare these evaluations to the results of KaRaMel. Indeed, the objective functions, calibration periods, discharge observations, and watershed surfaces used by [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF] and [START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systemes karstiques[END_REF] were different than those used in our work.

However, the former modelling exhibited that the precipitation-discharge relationship is already well described by conceptual model with one dimension. This conceptual model make it possible to simulate discharge satisfactorily with performance values around 0.8.

Considering that these performance criteria have an ideal value of 1, which can never be achieved [START_REF] Gupta | Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling[END_REF], the possibility of significantly improving discharge modeling is small. The comparison of KaRaMel simulation with KDM shows that KaRaMel bring improvements of discharge modelling (KGE increased from 0.8 to 0.9). Thus the consideration of the hydrosystem heterogeneity enables significant improvements. Perhaps modelling performances around 0.9 have reached an upper limit. Spring discharge modelling may not be more improved until the trend of overestimation-underestimation is solved.

When KaRaMel simulations of spring discharge were analysed in details (Figure 11) a cycle of overestimation -underestimation appeared. This cycle may be due to either the mathematical model or to a hydrologic process within the aquifer that we have not conceptualized and represented. A review of literature focusing on karst aquifer discharge modelling with distributed [START_REF] Robineau | A double medium approach to simulate groundwater level variations in a fissured karst aquifer[END_REF] and lumped models [START_REF] Hosseini | Spring Hydrograph Simulation of Karstic Aquifers: Impacts of Variable Recharge Area, Intermediate Storage and Memory Effects[END_REF] indicates that this trend has been present in numerous cases. We therefore believe that this cycle stems from an internal process of the karstic aquifer. Based on a study performed at the local scale on the same aquifer, [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF] hypothesize that flows that connect different water bodies within the unsaturated zone of the aquifer may change during the year. During the dry period, weak flows drain only a few water masses, but during wet periods, increases in flows promote hydraulic connections between different water stocks that feed the flows more abundantly [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF]. We tried to account for this trend and modified the KaRaMel model to store more water at the beginning of humid season and then release water by, for instance, adding an exchange flow between reservoir M and C governed by saturation in reservoir M, or by modifying the x i in function of the filling of M. These attempts did not remove the trend.

Heterogeneity of recharge

The annual recharge (calendar year) of the Fontaine de Vaucluse aquifer was estimated with seven usual methods by [START_REF] Lanini | Groundwater recharge and associated uncertainty estimation combining multi-method and multi-scale approaches[END_REF]. This study used the time series The difference between annual recharge estimation may be due to the distribution of climate variables and SAWC. In this study, climate variables exert a major control over the spatial distribution of recharge. However, the resolution of climate data from SAFRAN is coarse compared to the resolution of the SAWC map (8 x 8 km vs. 1 x 1 km). Moreover SAFRAN data results from the combination of information with even lower resolutions (30 km x 30 km for the atmospheric model and some thousands of km 2 for homogeneous climatic area classification) and meteorological network measurements [START_REF] Vidal | A 50-year high-resolution atmospheric reanalysis over France with the Safran system[END_REF]. The impact of the lower resolutions is visible on rain maps in Figure 13, particularly for 2011-2012. Thus, the resolution of climatic data masks a part of recharge heterogeneity (in particular in relation to rainfall). The representation of recharge heterogeneity could be improved if a finer resolution of climatic data was used. This is particularly true as our study site presented large range of elevation (even within one single SAFRAN mesh).

After climate variables, the SAWC is the second most important factor controlling the spatial distribution of recharge. We chose to constrain this value by using observations made in soil pits. Soil pit data collected in the natural environment and in agricultural areas provide a relevant assessment of the spatial variability of SAWC. However, in karst environments, soils are often quite rich in coarse elements which limits the depth of soil pits. Therefore, measurements in the natural area are not always representative of total SAWC. For example, in forested areas, trees can extend their roots rock fractures until more than 5 or 10 m [START_REF] Canadell | Maximum rooting depth of vegetation types at the global scale[END_REF]Carrière et al., 2019b,a). This deeper water resource is not accounted for by observations made in a soil pit. SAWC measurements including shallow and deep water ressource exist but their are limited to experimental site [START_REF] Rambal | Variabilité des propriétés hydrodynamiques du sol à l'échelle d'un versant karstique[END_REF]. To overcome this problem, the SAWC is usually multiplied by a factor of two for ecophysiological modelling [START_REF] Cailleret | A synthesis of radial growth patterns preceding tree mortality[END_REF]. In the future it will be relevant to determine the amount of water that trees can take up at depth to improve the water balance, with inverse ecophysiological modelling, inverse hydrological modelling [START_REF] Jukić | Groundwater balance estimation in karst by using a conceptual rainfall-runoff model[END_REF] or isotopic tracing [START_REF] Barbeta | The combined effects of a long-term experimental drought and an extreme drought on the use of plant-water sources in a Mediterranean forest[END_REF]. Soil information could also benefit from the world soil map at finer resolution announced by [START_REF] Mulder | Globalsoilmap france: High-resolution spatial modelling the soils of france up to two meter depth[END_REF].

We have shown that recharge distribution is highly variable inside the watershed. In fact, the results show wide variations of contributory zones to recharge from one year to the next. Previous studies demonstrated also the variation of recharge area of karst aquifers by using hydro-chemical monitoring [START_REF] Ravbar | The study of hydrodynamic behaviour of a complex karst system under low-flow conditions using natural and artificial tracers (the catchment of the Unica River, SW Slovenia)[END_REF] or hydrological modelling [START_REF] Hartmann | Progress in the hydrologic simulation of time variant recharge areas of karst systems -Exemplified at a karst spring in Southern Spain[END_REF]. However, they were not locating spatially this variation.

Our work has endorsed this concept and seeks to locate and quantify the participation level of each pixel of watershed in overall recharge. Our work goes further; we show that recharge distribution has a significant impact on the hydrodynamic response of karst.

Understanding global change impacts on karst hydrosystem response

The Mediterranean region has been identified as a "hotspot" of current and future climate change suggesting that changes in precipitation patterns may be significant [START_REF] García-Ruiz | Mediterranean water resources in a global change scenario[END_REF][START_REF] Pachauri | Climate change 2014: synthesis report. Intergovernmental Panel on Climate Change[END_REF]. Our study shows that annual precipitation varies considerably from one year to another while evapotranspiration remains more stable.

This highlights that mediterranean vegetation is acclimated to drought and preserve a relatively stable transpiration as illustrated in Figure 12. Increased ETo due to climate change is expected and will cause increased hydraulic failure of plants and may lead to widespread forest die-off. Therefore, KaRaMel implemented with evapotranspiration model taking into account plant activity may be appropriate for predicting the impacts of climate change evapotranspiration [START_REF] Ollivier | Estimation of recharge in karst aquifer using improved evapotranspiration monitoring thanks to remote sensing data[END_REF]. A sustainable prediction of climate changes implication for groundwater resources may not be led with 1D lumped model because it can not supported the land cover changes and even precipitation patterns, but it can be done with KaRaMel. A change in spatial distribution of annual precipitation would have an immediate effect on the hydrodynamics of a karst system.

For equivalent cumulative rainfall, we have shown that a shift of storm location from the western to the eastern part of the watershed could have a drastic impact on aquifer discharge (Figure 10). Likewise, an accentuation of the precipitation gradient will also lead to a modified hydrological response of the hydrosystem.

Recharge sensitivity to temporal precipitation distribution has been characterized by various approaches including water balance modelling [START_REF] Hughes | Evaluation of distributed recharge in an upland semi-arid karst system: the west bank mountain aquifer, middle east[END_REF], chloride mass balance [START_REF] Marei | Estimating groundwater recharge using the chloride mass-balance method in the West Bank, Palestine[END_REF] and hydrological modelling [START_REF] Pulido-Velazquez | Assessment of future groundwater recharge in semi-arid regions under climate change scenarios (Serral-Salinas aquifer, SE Spain). could increased rainfall variability increase the recharge rate?[END_REF].

The hydrological modelling in this study contributes to the characterisation of temporal recharge dynamics. Our study highlights that in a Mediterranean climate, heavy rainfall during autumn or winter promotes greater recharge than precipitation that occurs in the spring or summer. This effect is partly due to lower evapotranspiration during the colder autumn and winter. Thus, recharge is closely linked to land use and plant transpiration [START_REF] Scanlon | Global models underestimate large decadal declining and rising water storage trends relative to grace satellite data[END_REF]. Climate change may impact recharge in two ways, changes of

Figure 1 :

 1 Figure 1: Scheme of main elements of recharge.

Figure 2 :

 2 Figure 2: (a) Structure of one model unit. Unit parameters are indicated in red. Index "i" relates to the grid cell number. (b) Semi-distributed model architecture enables the distribution of flows and reservoirwater levels over the aquifer; For the sake of clarity only 3 units are drawn (i=1, 13, and 20)

Figure 3 :

 3 Figure 3: Overview of the study area. The Fontaine de Vaucluse aquifer is located in southeastern France, 30 km east of the city of Avignon. The watershed area of about 1160 km 2 , delimited by the red line, is covered by 34 cells of the meteorological product SAFRAN (green grid). The altitude variation is illustrated by the grey scale, the lowest altitude is the spring of Fontaine de Vaucluse at 89 m and the highest is the summit of Mont Ventoux at 1912 m. The orange grid represents the limestones of the Cretaceous period.

Figure 4 :

 4 Figure 4: Daily time series of precipitation (P) and reference evapotranspiration (ETo) of the studied area. Fontaine de Vaucluse discharge time series with modelling periods.

  et al. (2003)) identifies seven types of soil genesis factors: soil properties (S), climate and environment (C), organisms and vegetation (O), topography and landscape (R), parent material and lithology (P), age (A), spatial position (N). The soil attribute

R

  and O factors for mapping SAWC. Selection of appropriate covariates of SAWC for agricultural areas and natural areas was based on the analysis of variance. R factor. Topography and landscape factors were derived from a principal component analysis (PCA) of a Digital Elevation Model (DEM, Figure3). We used the DEM BD Alti from the French National Institute of Geographic and Forest Information with a spatial resolution of 25 m. PCA components were for example North-South exposure, western exposure, and mountain amplitude.

Figure 5 :

 5 Figure 5: Soil pit locations over the Fontaine de Vaucluse impluvium.

  approach evaluated the model dependence to the set of training data. Repeated training of the ANN model with various datasets provided similar results, thus highlighting ANN model constancy and limited dependence on the training data set. ANN and LM were established using the nnet and lm functions in the R programming language (Venables

Figure 6 :

 6 Figure 6: Intrinsic vulnerability map of Fontaine de Vaucluse with resolution of 1 km 2 , after Ollivier et al. (2019a,b).

  14) where Q O is observed discharge, ( QO ) is the mean observed discharge, Q S is the simulated discharge, n is the number of time steps. For more details, the performances were tiered for low flow periods (LF) and high flow periods (HF). Base flow periods match discharges mainly supported by slow flow. High flow periods match discharges mainly supported by rapid flows. Altogether, five indices KGE, NSE, and RMSE of low flow periods, and NSE and RMSE of high flow periods, were used to evaluate the simulations.

  . Model calibration was performed over three hydrologic years from September 2006 to August 2009. The 2006-2007 hydrologic year presented dry conditions with only 750 mm of precipitation and a mean daily discharge of 8 m 3 .s -1 . The 2008-2009 year was humid with 950 mm of precipitation and a mean daily discharge of 20 m 3 .s -1 . The 2007-2008 year presented an intermediate stage with about 850 mm of precipitation and a mean daily discharge of 12 m 3 .s -1 . Model evaluation was performed over the period of September 2009 to August 2015.4.7. Evaluation of parameter spatial variabilityKaRaMel has three sources of spatial variability: distribution of climatic variables (P and ETo), distribution of soil property (SAWC), and distribution of underground flow parameters (x i and kc i ). Due to the lack of direct measurements, which is the most common case in karst systems, direct validation of these distributed outputs is not possible. We thus consider a step by step procedure, with gradual introduction of the variability sources. The model was run with different combinations of variability sources : overall, eight combinations were tested. Model calibration was performed for each model setup. Comparison of model performances over the validation period, depending on the eight model setups, made it possible to hierarchize the effect of variability sources.

Figure 7 :

 7 Figure 7: Simulated versus observed SAWC values for forested and natural areas (green triangles) and agricultural areas (brown dots); (a) with multilinear regression models and (b) artificial neural network models.

Figure 8 :

 8 Figure 8: Soil available water capacity map obtained with an artificial neural network.

  soil water reservoir capacity of the semi-distributed model was directly constrained by the SAWC map. Parameters corresponding to rapid flows (x and kc) were distributed depending on the intrinsic vulnerability index. The slow reservoir discharge coefficient (km) is considered uniform over the watershed.

Figure 9 :

 9 Figure 9: Influence of model parameter distribution on model performance during the validation period: precipitation (P) and reference evapotranspiration (ETo), soil available water capacity (AWC), infiltration distribution coefficients (x i ) and specific discharge coefficients of rapid reservoirs (kc i ). a) Eight model setup were tested and b) their performances. The closer the points are to the outside of the radar chart, the better the model performs.

  shows the effects of these modifications on aquifer discharge and a comparison of results with observations over the period from September 2009 to August 2015. The wide range of discharge response is due to the percentage of precipitation affected by the scenario.The two proposed scenarios imply highest rainfall over selected parts of the aquifer: highest part or eastern part of the watershed. These scenarios increased recharge heterogeneity resulting in less recharge where precipitation was lower and higher recharge where precipitation was concentrated. The overall effect of these changes affects a larger volume of rapid flow through a part of the karst. Overall evapotranspiration decreases, recharge increases, thus discharge increases. Modification of the precipitation distribution has a direct impact on aquifer discharge. This result highlights the rapid hydrological response of the karst system to climate modifications.

Figure 10 :

 10 Figure 10: The orange areas correspond to the range of results obtained with various percentages of storms that move from the western to the eastern part of the watershed. Red areas correspond to the range of results obtained with various percentages of precipitation that moves from highlands to mountains in the watershed. The blue line represents observed discharge. At left, the effects of climate modification on daily discharge from September 2009 to August 2015. At right, cumulate discharge over this period.

Figure 11 :

 11 Figure 11: Upper part: observed (grey) and simulated (red) discharges with KaRaMel. Lower part:

Figure 12 :

 12 Figure 12: Annual recharge (Ry), evapotranspiration (ETy), and precipitation (Py) of hydrological years from September 1958 to August 2015. Winter precipitation is defined as precipitation that occurs between October and February.

Figure 13 :

 13 Figure 13: Spatial distribution of annual precipitation (P), evapotranspiration (ET), recharge (R) over the watershed area of the Fontaine de Vaucluse aquifer (blue scale bar) and the R/P and ET/P ratios(red scale bar).

  parameters). Models performances were evaluated with the objective function NSE calculated on the square root of discharge. The performances ranged between 0.88 to 0.92 over the 2000-2004 period. Fleury et al. (2007) indicates a NSE of 0.92 for KDM simulation over 1996-2005 period. Note that our application of KDM over 2009-2015 led to lower NSE of 0.79 when using the original calibration by Fleury et al. (2007) and 0.83 with our calibration.

  of precipitation and evapotranspiration calculated for the watershed from SAFRAN data, the same variables as those used in this work. However, they used a slightly different watershed surface and an unique value of SAWC of 20 mm when the mean SAWC used in KaRaMel is 70 mm for the whole watershed. Over the period 1996-2011, they estimated the annual recharge from 89 mm.y -1 (in 2007) to 623 mm.y -1 (in 2008) with an average value of 335 mm.y -1 . The annual recharges simulated with KaRaMel from 1996 to 2011 (calendar year) range from 78 mm.y -1 (2007) to 695 mm.y -1 (2008), with an average value of 412 mm.y -1 . The average difference between the two modeling experiments is 100 mm.y -1 , and the annual recharges simulated with KaRaMel are higher. Surprisingly, with a weaker SAWC,[START_REF] Lanini | Groundwater recharge and associated uncertainty estimation combining multi-method and multi-scale approaches[END_REF] would simulate a higher recharge.The difference between the annual recharge estimated by the two modelling experiments may be due to evapotranspiration models and the distribution of variables. Empirical or conceptual recharge models based solely on precipitation records can simulate a wide range of recharge, but it is impossible to verify the consistency of the estimate. The use of discharge records enables the validation of simulated recharge.

  spatial and temporal precipitation distribution may modify and alter the nature of land cover. Future model development will have to implement climate change, land use, and land cover change in order to anticipate the impact of global change on groundwater resources.7. ConclusionsEnvironmental factors (i.e. climate, soil, and geology) explain the diversity of recharge paterns. Constraining a semi-distributed karst model with readily available surface description data was both necessary and challenging. The Karst Recharge and discharge Model (KaRaMel) was designed to embed the primary knowledge available about karst hydrosystems. Our innovation is that parameters are based on aquifer characteristics that have been quantified or constrained using soil available water capacity and intrinsic vulnerability mapping. The use of measurable parameters in the field reinforces the realism of the model and minimizes minimizes the number of parameters to be calibrated. The proposed link between rapid flow pathways and the vulnerability index makes it possible to easily integrate a synthesis of voluminous data based on geologist and hydrogeologist expertise. The proposed soil mapping method is robust for large scale studies. It makes it possible to consider the hydrological role of soil heterogeneity. Rainfall-discharge records are used for model calibration. This method may help resource managers and researchers incorporate a large amount of knowledge of the studied aquifer with the great advantage of needing few calibrated parameters.The major advantage of KaRaMel is that it provides assumptions on aquifer storage distribution that the 1D model does not offer. If the purpose of future studies is to simulate discharge of a karstic system similar to that of the Fontaine de Vaucluse, lumped models are sufficient. If the purpose of the study is to locate recharge areas, KaRaMel is an efficient tool that has limited data demands and requires minimal calculation time (< 1 min. for simulation of daily flows over ten years).The application of KaRaMel to the Fontaine de Vaucluse karst system shows that a large part of discharge variation at the outlet may be explained by recharge distribution. The daily distribution of precipitation, reference evapotranspiration, and soil available water capacityt are the primary sources of the overall improvement of discharge modelling in term of quantities. The evapotranspiration demand calculated with the Penman-Monteith model[START_REF] Allen | Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[END_REF] based on SAFRAN data assumes that evapotranspiration is mainly driven by weather conditions. Future development might include a more refined representation of the evaporation and transpiration processes. Over last decade, annual recharge has changed radically from one year to the next over the watershed area. Half of the years studied show a complex situation with the coexistence of areas of high recharge and areas of almost no recharge. Therefore, a homogeneous average recharge value for the entire watershed area is an imprecise approximation of reality.The average overlooks the complex and heterogeneous mechanisms at stake. In 2011-2012 despite an intermediate cumulative rainfall (630 mm) on the Fontaine de Vaucluse impluvium, the eastern part had a strong water deficit in terms of recharge. An understanding of spatial and temporal recharge distribution is crucial for sustainable groundwater management. KaRaMel may be a useful tool to help water resource managers face future issues in the context of global change. Our simulations show that a slight modification of rainfall distribution on the Fontaine de Vaucluse impluvium rapidly increased or decreased discharge at its outlet. The steep precipitation gradient and the heterogeneity of soil and hydraulic properties, and karst characteristics dominate the hydrological setting in our study area but similar characteristics might be found for a majority of other aquifers. The coexistence of flows with contrasting (slow and fast) dynamics is also observed in fissured and porous aquifers. It could be relevant to apply KaRaMel to other aquifer types where hydraulic parameter are poorly documented if the modelling objective is to understand the spatial distribution of aquifer stocks.

Table 1 :

 1 Overview of the model parameters, their description and physical units.

	Variables Description	Unit	Layer
	ETo	Reference evapotranspiration	L	Atmosphere
	ET	Evapotranspiration	L
	P	Precipitation	L
	S	Soil reservoir		Soil
	SAWC	Soil Available Water Capacity	L
	Hs	Filling rate of soil reservoir	L
	I	Infiltration from soil to the lower reservoirs	L
		(corresponding to the recharge of the aquifer)
	x	Infiltration part that feed reservoir M	(-)
	M	Slow discharge reservoir		Underground
	Hm	Filling rate of reservoir M	L
	km	Specific discharge coefficient of reservoir M	T -1
	qm	Outflow of reservoir M	L.T -1
	C	Rapid discharge reservoir
	Hc	Filling rate of reservoir C	L
	kc	Specific discharge coefficient of reservoir C	T -1
	qc	Outflow of reservoir C	L.T -1
	Q	Sum of outflows of reservoir C and M	L.T -1
	Qsim	Aquifer discharge	L.T -1 Atmosphere
		(corresponding to the discharge of the aquifer)
	aquifer. The number of parameters to be calibrated increases in proportion to the number
	of units of the model, theoretically there are 4*n parameters to calibrate. In order to

Table 2 :

 2 Texture of fine earth constituents (particles ranging from 0 to 2mm) of soil according to the Jamagne classification and corresponding water holding capacity (WHC).

	Soil class	Soil texture	WHC
	S	Sand	0.07
	SL	Loamy sand	0.1
	SA	Clay Sand	0.135
	LlS	Silty sandy loam	0.12
	LS	Sandy loam	0.145
	LMS	Silt loam (medium) 0.16
	LSA	Clay sandy loam	0.165
	LAS	Sandy clay loam	0.175
	Ll	Silty loam	0.13
	LM	Silt loam (fine)	0.175
	LA	Silty clay loam	0.195
	AS	Sandy clay	0.17
	A	Clay	0.175
	AL	Loamy clay	0.18

(Sa) model can be written as follows (x, y and t: spatial and time coordinates): S a = f (S x,y,t ; C x,y,t ; O x,y,t ; R x,y,t ; P x,y,t ; A x,y,t ; N x,y,t )

Table 3 :

 3 Distribution of infiltration distribution coefficients and reservoir specific discharge coefficients according to intrinsic vulnerability class.

Table 4 .

 4 Evaluation statistics were similar regardless of the period (KGE of 0.91 and 0.92), showing that the calibration period was long enough to be representative of a wide range of aquifer behaviours. KaRaMel performances over the validation period were higher than performances of the 1D model, showing that spatial distribution of aquifer properties and climate variables improved the modelling. The KGE increased from 0.87 to 0.92; most of the improvement was obtained for high flow periods (NSE increased from

0.66 to 0.75). The spatial distribution of aquifer flows resulted in larger improvements for high flow than for low flow periods.

Table 4 :

 4 Performances of KaRaMel over calibration and validation periods. KDM performances using the original calibration by[START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to fontaine de vaucluse (south of france)[END_REF].

			Low flows	High flows
			KGE
		Period	NSE	RMSE	NSE	RMSE
				(m 3 .s -1 )		(m 3 .s -1 )
		Calibration 0.91 0.75	2.68	0.75	6.78
	KaRaMel
		Validation 0.92 0.83	2.56	0.75	6.68
	KDM	Validaion	0.87 0.79	2.97	0.66	7.48

Acknowledgments

This research was supported by the Conseil Régional Provence-Alpes-Côte d'Azur (no grant number) and the Syndicat Mixte du Bassin des Sorgues (www.laSorgue.com). This work benfited from fruitful discussion within the CERISE LEFE-INSU project and within the KARST observatory network (SNO KARST) initiative from the INSU/CNRS. SNO KARST aims to strengthen knowledge sharing and to promote cross-disciplinary research on karst systems at the national scale. SAFRAN data were provided by METEO FRANCE, the authors apprciate the support from Frédéric Huard (INRA). The authors thank the three anonymous reviewers for their constructive and helpful comments.

Appendix A. Geological setting of Fontaine de Vaucluse aquifer water eventually infiltrates into the limestone, [START_REF] Ollivier | Caractérisation et spatialisation de la recharge des systèmes karstiques. Application au système de la Fontaine de Vaucluse[END_REF]