Keywords: Petrophysics, Inversion, Time-lapse, Seismics, Groundwater

The hydrodynamic characterization of the epikarst, the shallow part of the unsaturated zone in karstic systems, has always been challenging for geophysical methods. This work investigates the feasibility of coupling time-lapse refraction seismic data with petrophysical and hydrologic models for the quantitative determination of water storage and residence time at shallow depth in carbonate rocks. The Biot -Gassmann fluid substitution model describing the seismic velocity variations with water saturation at low frequencies needs to be modified for this lithology. I propose to include a saturation dependent rock-frame weakening to take into account water-rock interactions. A Bayesian inversion workflow is presented to estimate the water content from seismic velocities measured at variable saturations. The procedure is tested first with already published laboratory measurements on core samples, and the results show that it is possible to estimate the water content and its uncertainty. The validated procedure is then applied to a time-lapse seismic study to locate and quantify seasonal water storage at shallow depth along a seismic profile. The residence time of the water in the shallow layers is estimated by coupling the time-lapse seismic measurements with rainfall chronicles, simple flow equations and the petrophysical model. The daily water input computed from the chronicles is used to constraint the inversion of seismic velocities for the daily saturation state and the hydrodynamic parameters of the flow model. The workflow is applied to a real monitoring case and the results show that the average residence time of the water in the epikarst is generally around three months, but it is only eighteen days near an infiltration pathway.

During the winter season the residence times are three times shorter in response to the increase of the effective rainfall.

Introduction

Understanding the flow and transport properties through the unsaturated zone of karst systems has always been a great challenge in hydrology. Traditional problems such as borehole positioning or karstification intensity and infiltration potential mapping [START_REF] Kavouri | PaPRIKa: a method for estimating karst resource and source vulnerability-application to the Ouysse karst system (southwest France)[END_REF] require the identification of targets like perched aquifers or preferred infiltration pathways. More recently the development of unconventional hydrocarbon resources is drawing major concerns about groundwater management and protection. When the exploration is targeting formations overlain by karstified surface carbonates, the concerns are even more serious because of the intrinsic vulnerability of the karstic aquifers to contaminant. It is therefore all the more important to find or develop the techniques which can monitor the saturation state of the unsaturated layer. The vadose zone of the karst plays a major role in water recharge and distribution to the deep saturated zone through the low-permeability volume (LPV). It is the focus of many studies using geochemical or isotopic tracers sampled in cave drips [START_REF] Arbel | Infiltration processes and flow rates in developed karst vadose zone using tracers in cave drips[END_REF][START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the relationship[END_REF], but geophysical methods have been seldom used for this purpose. Time-lapse gravity measurements [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF] or magnetic resonance soundings [START_REF] Mazzilli | Applicability of MRS Soundings for the Characterisation of the Unsaturated Zone of Karst Systems[END_REF] have been attempted, but the former deliver only non-localised water balances, and the weak response of the later is difficult to interpret quantitatively. [START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF] used time-lapsed electrical and seismic tomographies (respectively ERT and SRT) and showed that electrical and acoustical properties of the epikarst are changing after rainfall. Variations in electrical resistivity and seismic velocity were easily observed in field data and on the inverted models, but penetration depth, resolution and repeatability were much better for SRT than ERT. In the present paper I extend this work and show that quantitative estimates of water content can be inverted from the seismic time-lapse data in karstic hydrogeophysics. Quantitative time-lapse analysis is already used for hydrocarbons recovery [START_REF] Landrø | Uncertainties in quantitative time-lapse seismic analysis[END_REF][START_REF] Veire | Stochastic inversion of pressure and saturation changes from time-lapse multi component data[END_REF] or for CO 2 injection monitoring [START_REF] Ivanova | Monitoring and volumetric estimation of injected CO 2 using 4D seismic, petrophysical data, core measurements and well logging: a case study at Ketzin, Germany[END_REF], and I investigate its applicability in karstic hydrology. I review first the saturation-related properties of rocks and introduce a model of saturation -velocities relationship taking into account the specificities introduced by the interaction of water with carbonates at shallow depth. This model is used to invert time-lapse seismic velocities for saturation parameters within a Bayesian framework. The effectiveness of the procedure is tested first on existing laboratory data available in the literature, and then a real case application for temporary water storage mapping is presented. Going one step further, I present an inversion workflow coupling time-lapse seismic measurements with simple hydrologic modelling, and show that the residence time of the water inside the epikarst can be estimated locally from seismic velocities.

Modelling water saturation effects in carbonates

I review in this section the dependencies of the seismic compression and shear waves velocities V p , V s upon the variable saturation S w . The seismic response of a medium is characterised by the bulk modulus, K, the shear modulus, , and the bulk density, , which are related to seismic velocities through the fundamental equations

𝜌𝑉 𝑃 2 = 𝐾 + 4 3 𝜇, (1)
𝜌𝑉 𝑆 2 = 𝜇.

(2)

When replacing air by water in the pore space, seismic velocity is normally computed using effective medium theories [START_REF] Mavko | The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media[END_REF] and the quasi static equation for fluid saturation effects [START_REF] Gassmann | Über die elastizität poröser medien[END_REF]Biot, 1956a). The bulk modulus K is replaced at any saturation by K sat defined by

𝐾 𝑠𝑎𝑡 𝐾 0 -𝐾 𝑠𝑎𝑡 = 𝐾 𝑑𝑟𝑦 𝐾 0 -𝐾 𝑑𝑟𝑦 + 𝐾 𝑓𝑙 ∅(𝐾 0 -𝐾 𝑓𝑙 ) , (3)
where K dry is the dry-rock modulus, K fl the bulk modulus of the water-air mixture, K 0 is the rock frame modulus, and  is the total porosity. The shear modulus  is invariant during fluid substitution. When air and water are mixed at fine scale (this scale issue is discussed later in this section), the fluid bulk modulus is described by the Reuss average

1 𝐾 𝑓𝑙 = 𝑆 𝑤 𝐾 𝑤 + 1 -𝑆 𝑤 𝐾 𝑎 , (4)
with K w and K a the moduli of water and air respectively.

The model predicts that Vp velocity decreases with saturation (trend A in Figure 1) before a sharp increase near full saturation, while Vs decreases monotonically. The linear decrease is mainly a density effect (increasing water content) and does not exceed 15% in general. The Biot-Gassmann model provides a robust fluid prediction tool as long as the main assumptions behind the theory are satisfied. I discuss now two of them, which are the most critical for the lithology considered in this work.

The first requirement is that there is no interaction between the rock frame and the fluid, an assumption, which is obviously false during water saturation of carbonates. Several authors [START_REF] Clark | Effect of volatiles on attenuation (Q-1) and velocity in sedimentary rocks[END_REF][START_REF] Cadoret | Effet de la saturation eau/gaz sur les propriétés acoustiques des roches: étude aux fréquences sonores et ultrasonores[END_REF] pointed out that during the initial introduction of moisture the velocity drops sharply by a few percent (B in Figure 1) for different types of samples.

This apparent softening of the rock has been attributed to various mechanisms such as the dissolution of cement coating the grains, or the disruption of cohesive forces between very dry surfaces. From a chemical point of view, [START_REF] Vanorio | The effect of chemical and physical processes on the acoustic properties of carbonate rocks[END_REF] showed on carbonates samples that after 10 days of immersion in cold (10°C) carbonated water at room pressure, the saturated Vp and Vs ultrasonic velocities are smaller than the values predicted by the Gassmann model. Moreover, after drying the samples, measurable changes in the microstructure of the samples occurred; porosity and permeability increased by 3%, and up to 400mD respectively. A circulation of an aqueous CO 2 -rich solution leads with time to a decrease in elastic moduli by up to 60% at room pressure [START_REF] Vanorio | Rock physics analysis and time-lapse rock imaging of geochemical effects due to the injection of CO 2 into reservoir rocks[END_REF], while the decrease is smaller (5%) at higher pressure [START_REF] Grombacher | Time-lapse acoustic, transport, and NMR measurements to characterize microstructural changes of carbonate rocks during injection of -rich water[END_REF]. Most of the rock weakening is taking place during the early stages of water circulation and the effect is levelling off with time. The conclusion of all these experiments is that it is expected to observe large velocity variations (at least 30 %?) when the saturation of the epikarst changes, because the confining pressure is small and the water reactive surface is large in this highly fractured shallow layer. It is therefore necessary to include these chemo-mechanical mechanisms in the description of saturation effects.

An equivalent form of equation ( 3) is [START_REF] Mavko | The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media[END_REF])

1 𝐾 𝑠𝑎𝑡 = 1 𝐾 0 + 𝜙 𝐾 𝜙 + 𝐾 0 𝐾 𝑓𝑙 𝐾 0 -𝐾 𝑓𝑙 , (5)
where K Φ is the pore space compressibility. If the rock weakening is originating from a change in the rock frame properties, I can add an extra compliance

1 𝐾 0 * = 1 𝐾 0 + 1 𝐾 𝑐ℎ𝑒𝑚 , (6)
and because K fl is much smaller than K 0 and K 0 * for air-water mixtures,

1 𝐾 𝑠𝑎𝑡 * = 1 𝐾 0 * + 𝜙 𝐾 𝜙 + 𝐾 0 * 𝐾 𝑓𝑙 𝐾 0 * -𝐾 𝑓𝑙 ≅ 1 𝐾 𝑠𝑎𝑡 + 1 𝐾 𝑐ℎ𝑒𝑚 . (7)
In the same way, if the porosity increases (dissolution), writing *=+d and substituting in (5) yield the same kind of relation. I therefore do not attempt to describe exactly the rock softening with a model using some microstructural parameters (eg some kind of reverse diagenetic trend using the contact-cement model of [START_REF] Dvorkin | Elasticity of high-porosity sandstones: Theory for two North Sea data sets[END_REF], but introduce instead a single extra saturationdependent compliance 1/K chem to take into account all possible rock weakening chemo-mechanical effects. For the sake of simplicity I define K chem by

𝐾 𝑐ℎ𝑒𝑚 = 𝐹 𝑐ℎ𝑒𝑚 (𝑐ℎ𝑒𝑚, 𝑆 𝑤 ) 𝐾 0 , (8)
where F chem is a simple mathematical function increasing from 0 to a maximum value describing the levelling off above a threshold saturation arbitrarily set to 0.2 (see Figure 2). I proceed in the same way with the shear modulus, using the same function

𝜇 𝑐ℎ𝑒𝑚 = 𝐹 𝑐ℎ𝑒𝑚 (𝑐ℎ𝑒𝑚, 𝑆 𝑤 ) 𝜇 0 , (9)
1 𝜇 0 * = 1 𝜇 0 + 1 𝜇 𝑐ℎ𝑒𝑚 . ( 10
)
Note that according to Gassmann theory we will have  sat = µ 0 *. The model uses a single parameter chem to describe the chemo-mechanical effects on the elastic moduli K and µ. This is certainly a crude approximation but I want to keep the model simple and tractable for inversion purposes.

Figure 2 shows that smaller value of chem yield larger chemical compliance; the parameter will describe at least qualitatively the reactivity of the rock to water.

The second main requirement is that air and water should be mixed at fine scale to compute the bulk modulus of the mixture. Complex saturation patterns are common in carbonate rock, owing to their fabric heterogeneities. [START_REF] Sharma | Sensitivity of flow and elastic properties to fabric heterogeneity in carbonates[END_REF] studied several examples of natural and forced saturation, and showed that up to 44% of the pore volume was not invaded during natural (benchtop) imbibition, while up to 20% was still unsaturated after the forced (pressure) saturation stage. Some samples of the [START_REF] Cadoret | Effet de la saturation eau/gaz sur les propriétés acoustiques des roches: étude aux fréquences sonores et ultrasonores[END_REF] study used in this paper show some departures from the Gassmann trend and differences in velocities between saturation by depressurization and desaturation by simple drying. When heterogeneities (patches) occur in rocks at mesoscopic scale (size larger than grain size but smaller than the seismic wavelength, i.e. a few metres inside the epikarst), and if the seismic frequency is high enough, pressure gradients build up between patches and wave-induced fluid flow (WIFF) takes place in the pore space to equilibrate the pressure of the passing seismic wave (Biot, 1956a(Biot, , 1956b;;Dutta andOdé, 1979a, 1979b;[START_REF] Mavko | Estimating grain-scale fluid effects on velocity dispersion in rocks[END_REF].

WIFF is now considered a major cause of seismic attenuation and velocity dispersion in porous media [START_REF] Pride | Seismic attenuation due to wave-induced flow[END_REF][START_REF] Müller | Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks -A review[END_REF], even for the low (10 Hz -occur because of the complexity of porosity in carbonates rocks [START_REF] Cadoret | Fluid distribution effect on sonic attenuation in partially saturated limestones[END_REF][START_REF] Adam | Seismic wave attenuation in carbonates[END_REF], the non-uniform (patchy) partial saturation [START_REF] Toms | Seismic attenuation in porous rocks with random patchy saturation[END_REF], or the presence of water-saturated cracks in a porous rock frame. It is beyond the scope of this paper to report on the saturation-related seismic attenuation effects; this topic has been discussed for the data used in this work in [START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF]. We focus here on the velocity dispersion effect because I plan to invert the velocities for the water content. Several theories are available in the literature to model the effects of WIFF (see [START_REF] Müller | Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks -A review[END_REF] for a review). As an example, Figure 1 displays the velocity curve obtained using the dual porosity model of [START_REF] Chapman | Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity[END_REF] modified to take in account the chemo-mechanical effect at low saturation. We can see that the velocity starts to increase at partial saturation (C in Figure 1), when the squirt flow between pores and cracks takes place. Such kind of models requires the parameterization in some way of the microstructure of the rock, which is not easy to achieve. In this work I simplify the problem by considering that patchy saturation occurs even with a small amount of water, and I just replace the effective fluid modulus of Eq. 4 by the Voigt average [START_REF] Mavko | The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media[END_REF])

𝐾 𝑓𝑙 = 𝑆 𝑤 𝐾 𝑤 + (1 -𝑆 𝑤 )𝐾 𝑎 . ( 11
)
In such a case, the velocity will increase monotonically with saturation (dashed line in Figure 1). I have described a model providing a proxy between seismic velocities of carbonates and their water storage capabilities described by porosity and saturation. In the following sections the accuracy of the model is checked against experimental data available in the literature, and I present two realcase applications in the study of the vadose zone in karstic systems.

Validation of the model

I validate first the modified Biot -Gassmann model using the experimental data provided by the reference work of [START_REF] Cadoret | Effet de la saturation eau/gaz sur les propriétés acoustiques des roches: étude aux fréquences sonores et ultrasonores[END_REF]. This dataset consists of benchtop velocity measurements at sonic frequencies of 1 kHz or 2 kHz with resonant bar (1 m long cylindrical samples) and controlled partial water saturation, for several carbonates lithologies. The samples were chosen homogenous enough to limit high frequency rock stiffening complications; porosity values were computed from X-rays tomography maps on small-size plugs. Cadoret's work used sonic frequencies of 1 kHz or 2 kHz, but I anticipate that the conclusions of the inversion tests will be valid in the lower frequency band of the field seismic records. The modified Biot -Gassmann model m and probabilistic Bayesian inversion [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] are used to invert the experimental data d (measured seismic velocities Vp and Vs in the dry and partially saturated states) for the vector of rock properties m=(,S w ,, chem). Writing Bayes theorem,

𝑃(𝐦|𝐝) = 𝑃𝑃𝐷 = 𝑃(𝐝|𝐦)𝑃(𝐦) 𝑃(𝐝) , (12)
I compute P(m|d) the posterior probability density (PPD) over the model space, from P(d|m) the likelihood (misfit, including data uncertainty) and P(m) the prior knowledge about the solution. In practice the quantity P(d) on the right-hand side of equation ( 12) is regarded as a normalization factor [START_REF] Bachrach | Joint estimation of porosity and saturation using stochastic rock-physics modeling[END_REF] and dropped in the calculations. The results can be displayed using m*, the a posteriori most probable set of parameters (MAP); however the benefit of Bayesian framework is that uncertainties for a given parameter m i can be studied by computing P(m i |d), the marginal probability density distribution of parameter m i knowing data d. The uncertainties can be described by a credibility interval computed from such marginal distributions. In the same way, joint probability density maps can be used to study any possible correlations for a couple of parameters. Within this framework, it is possible to handle this nonlinear problem and calculate the most likely values of porosity and saturation given the seismic velocities [START_REF] Bachrach | Joint estimation of porosity and saturation using stochastic rock-physics modeling[END_REF].

Figure 3 displays the prediction results obtained with the Estaillades limestone of [START_REF] Cadoret | Effet de la saturation eau/gaz sur les propriétés acoustiques des roches: étude aux fréquences sonores et ultrasonores[END_REF].

The solutions space is sampled with a Markov Chain Monte Carlo (MCMC) approach using 40000 random draws of the parameter vector m for homogenous saturation and pure calcite mineralogy.

Figure 3a displays the joint a posteriori probability distribution between S w and ; the black dot represents the measured solution with a horizontal error bar indicating the porosity heterogeneity.

The distribution shows that it is impossible to accurately determine each parameter; the highest probability trend (smaller misfit, darkest area in the figures) follows a constant S w (water content) trend (dotted lines). This is an expected result, for the model is based on Gassmann equation, which

shows that at partial saturation the velocity decreases mostly because the bulk density increases while adding water to the pore space. Consequently, I will invert from now only for water content S w and chemical compliance Chem.

I proceed in the following way to introduce some a priori knowledge in the inversion and to sample the model space. Firstly, I determine from all measurements for all samples of [START_REF] Cadoret | Effet de la saturation eau/gaz sur les propriétés acoustiques des roches: étude aux fréquences sonores et ultrasonores[END_REF] two functions S1(Vs) and S2(Vs), where Vs is the measured shear waves velocity at any partial saturation; these functions provide upper and lower bounds for the S w solutions. This step is performed only once. For real field data, such bounds could be computed from samples taken in the area in order to better fit with local lithology. Nevertheless in this paper I used the same bounds for the field examples discussed in the following sections of the paper. In the next step, for any set of velocities (Vp, Vs) (dry, wet) to be inverted, I determine S min and S max from S1, S2, and Vs in the wet state. For each value of S drawn from the MCMC sampling of this interval I generate several (,S)

values and combine with the random sampling of the Chem values to get the 40000 solutions used to build the probability distributions.

By this way I attempt to solve the ambiguity near full saturation (two different saturations for the same Vp value, see Figure 1). When computing the probability distribution, the misfit term may include or not the Vs value in the wet state. In the latter case, it is possible to use only some rough estimate for this Vs value to set the bounds; this is particularly useful when dealing with field data because Vs is seldom available for both dry and wet states. Figure 3b displays the joint probability distribution between S w and chem for the same saturation as Figure 3a. We can see now that S is estimated separately from the chemical compliance value; there is a prominent probability peak close to the experimental value (dashed line). The 95% credibility interval (dotted lines) describes the uncertainty associated with this estimation. Note that the chemical compliance is displayed as 1/chem to show increasing chemical effect (see Figure 2).

I compare in Figure 3c the estimated water content with the value computed from the measured porosity and partial saturations. We can see that it is possible to estimate correctly the water content with its uncertainty provided by the error bars; for saturation values close to 100%, it is nevertheless necessary to include Vs wet in the misfit function to get correct values. We can draw the same conclusion from the results for all samples available in Cadoret's work; therefore I claim that the two parameters model is describing correctly the water saturation effects for homogenous saturation.

Figure 3c suggests that a small bias toward higher water content may exist; one possible explanation for this bias is that the experimental porosity value is too low. Porosity has been computed from small plugs, which are usually cut away from macro porosity heterogeneities for obvious reason of plugs quality. I used in the example shown an average porosity of 0.3 computed on several plugs, but values of 0.34 or higher are found; it is therefore possible that the plugs average may not be representative of the 1m long sample used in velocity measurements.

I conclude this section by discussing the critical issue of the dry sample requirement. It is necessary to have Vp and Vs values for the dry state in order to perform the calculation. This requirement is somewhat easy to meet for a laboratory experiment where samples can be oven-dried, but for velocities measured in the field, it may be difficult to record; it is always expected to have at least some moisture (bound water) in the ground. Figure 3d displays the same comparison as Figure 3c, but the velocities measured at 7% saturation have been used for the inversion instead of the correct dry velocities. If I replace in Figure 3d S by S, the water content variation, the result is still acceptable. In other words, for field velocity measurements, I consider that the inversion is robust enough (at least up to 10% residual water) to estimate water content variations if it is possible to identify the driest state between the two measurements. The most secure way to find this state is to look at the Vs data and choose the state with the highest Vs velocity. In the following sections, I show two real cases using this approach applied to time-lapse seismic surveys shot for hydrologic purposes in southern France.

Water content estimation in the epikarst

The first example deals with the mapping of the temporary water storage locations inside the epikarst. [START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF] The water content S w in April is estimated for each cell of the velocity model using the October measurements as dry reference (if some residual water was present in October, it is rather the variation in water content, see previous discussion). The input data are Vp and Vs velocities provided by independent traveltimes inversions. The stochastic inversion is carried out independently in each cell using 10000 random draws to build the probability function (considering here only homogenous saturation). The estimated variations in water content are displayed at the bottom of Figure 4. The results have been slightly smoothed to attenuate the blurry aspect generated by the stochastic inversion process. The picture shows that at the end of the wet season (April) the water is stored mainly inside a near-surface layer bound at depth by an irregular limit (dashed line in Figure 4 bottom). This layer which is less than 10m-thick is in good agreement with the seven metres of highly fractured epikarst found in boreholes located a few metres off-line from the middle of the profile.

The water storage area is thicker at the foot of the slope, suggesting that some water is flowing laterally downslope inside the epikarst instead of infiltrating directly the LPV. The presence of an unusual local accumulation of water at shallow depth after the dry season (black line of Figure 4) could be explained in the following way. The October seismic survey took place one month after the drilling of the two wells located some 15 m away from this location. The driller reported that the holes remained dry until the down-the-hole hammer reached the targeted karstic conduit 94 m below ground surface; at that time a large amount of water was raised by the compressed air from the 38 metre deep static level and flushed out around the holes for cleaning and testing purposes.

This artificial recharge of water was probably still present in the epikarst during the October survey.

Estimating water residence time

The previous example showed that it is possible to map the saturation in the epikarst. I show now that the geophysical quantities (the velocities) can be converted to hydrodynamic parameters using a simple hydrologic model of the epikarst and the actual volume of water (the rainfall) entering the system. I invert now for p=(,Ch, epk), i.e. respectively total porosity, chemical interaction, and the epikarst parameters defined hereafter, by fitting the velocities predicted by the petrophysical model to the time-lapse experimental values, knowing the water quantity S w . I use in this work data collected during the study of the Durzon karst system; details about the geological and hydrogeological setting are found in [START_REF] Bruxelles | Dépôts et altérites des plateaux du Larzac central : causses de l'Hospitalet et de Campestre[END_REF] and [START_REF] Jacob | Absolute gravity monitoring of water storage variation in a karst aquifer on the larzac plateau (Southern France)[END_REF]. The vadose zone, which is at least 100m-thick, has been periodically monitored by geophysical methods since the year 2006 [START_REF] Jacob | Absolute gravity monitoring of water storage variation in a karst aquifer on the larzac plateau (Southern France)[END_REF][START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF][START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF]. This paper uses the results of four seismic surveys (both 2D and 3D) shot between 2009 and 2012. The 150m x 150m area studied in this work is sitting on outcropping middle Jurassic dolostones; three boreholes indicate that the thickness of the epikarst is 7m, calcite is commonly sealing fractures in the upper LPV, and a sinkhole with a cave entrance has developed inside the dolostones.

The main outcome of the time-lapse study of the epikarst by [START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF] was that this unit must be divided into two sub-units when considering the response to rainfall events. There is an upper layer where seismic velocity markedly decreases with rainfall, and a deeper one showing some increase during the same elapsed time. Hydrogeologists [START_REF] Bakalowicz | Contribution de la géochimie des eaux à la connaissance de l'aquifère karstique de la karutification[END_REF][START_REF] Perrin | Epikarst storage in a karst aquifer: a conceptual model based on isotopic data, Milandre test site, Switzerland[END_REF] have already described this duality and they conceptually divide the unsaturated zone in two layers; the epikarst, the uppermost weathered layer of the karstic system with substantially enhanced porosity and permeability [START_REF] Klimchouk | Towards defining, delimiting and classifying epikarst: Its origin, processes and variants of geomorphic evolution[END_REF], and the underneath infiltration zone corresponding to the unsaturated part of the LPV. It is the result of the interactions between water, carbonates and carbon dioxide (CO 2 ); at shallow depth in the epikarst, CO 2 is continuously fed by the atmosphere or by the roots of the vegetation, leading to rock dissolution, while some calcite precipitation occurs at the top of the LPV. I postulate that the observations of [START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF] are the geophysical response to these complex interactions, and investigate the possibility to use velocity measurements as a proxy to the transit time of water through the upper unsaturated zone. In karst hydrology, reservoirs models are commonly used at the catchment scale to model spring hydrographs [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] Bailly-Comte | Water exchange and pressure transfer between conduits and matrix and their influence on hydrodynamics of two karst aquifers with sinking streams[END_REF][START_REF] Tritz | Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model[END_REF][START_REF] Ke | Modeling the groundwater recharge in karst aquifers by using a reservoir model[END_REF]; it is a convenient way to hide the complexity of the flow through the karst system while still providing useful insight with simple equations. Following this approach, the epikarst is conceptualized using a two reservoirs model to reflect the duality underlined by [START_REF] Valois | Caractérisation structurale de morphologies karstiques superficielles et suivi temporel de l'infiltration à l'aide des méthodes électriques et sismiques[END_REF]; the goal is to take advantage of the local information provided by the seismic to estimate locally the parameters describing the transfer function of the epikarst. I assume that I know from the seismic refraction results the geometry of the two layers (Figure 5) involved in the model and I build upon this assumption. The top surface layer is generally thin (at most a few metres, thickness E1), its velocity V1 is very low and seems to be negatively correlated to the shortterm cumulative rain. The second layer is thicker (thickness E2) with intermediate velocity V2

apparently positively correlated to the rain. The third layer is much faster and velocity variations with time are smaller and complex to interpret. The base of the epikarst (higher porosity layer) corresponds to the limit between layer 2 and 3 in the studied example; we call the first layer "Soil" because it is including the pedologic cover, and the second layer will be "Epikarst", albeit it is actually the lower epikarst. The associated hydrologic conceptual model is made of two reservoirs (Figure 5); input to the upper (soil) reservoir is the effective rain, i.e. the balance between rainfall R and predicted evapo-transpiration (PET), and the output is Q1. The volume of water H1 residing in the reservoir is therefore given by the balance equation

𝑑𝐻1 𝑑𝑡 = 𝑅 -𝑃𝐸𝑇 -𝑄1 , (13)
with the condition dH1/dt = 0 when H1=0. In practice the equation is discretized with an explicit Euler method using a daily time step because R and PET data are provided by a nearby weather station on a daily basis (in mm or l/m3, as usual in hydrology, hence the name water "height"). The output (discharge) Q1 from the soil reservoir is modelled according to a classical linear law of discharge through a porous outlet [START_REF] Maillet | La vidange des systèmes de réservoirs[END_REF]. In this case the discrete form of the daily volume of water leaving the reservoir is finally given by

𝑄1(𝑛) = 1 𝑇1 (𝐻1(𝑛) -ℎ1) , ( 14
)
where T1 is the recession constant of the flow exponential decay with time , usually given in days, H1

is the water height in the reservoir expressed in mm and the constant h1 has been introduced to describe some non-linearity; it is a threshold, the minimal water height necessary to trigger the discharge. The water height H1 for day n is computed from the effective rain and the previous height, 𝐻1(𝑛) = 𝐻1(𝑛 -1) + 𝑅(𝑛) -𝑃𝐸𝑇(𝑛) -𝑄1(𝑛) .

(15)

The discharge Q1 is the input to the deeper epikarst reservoir, which is described by the same law,

𝑄2(𝑛) = 1 𝑇2
(𝐻2(𝑛) -ℎ2) , ( 16)

𝐻2(𝑛) = 𝐻2(𝑛 -1) + 𝑄1(𝑛) -𝑄2(𝑛) . ( 17
)
The decay constants T1 and T2 are describing the flow through respectively the epikarst reservoir and the underneath low-permeability volume (LPV), while the thresholds h1 and h2 are describing flow delays for respectively the soil and the epikarst reservoir. [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France)[END_REF] split the discharge into fast and slow discharges. This conceptual partition can be related to the division of the pore space into well-connected ("open fractures") and less-connected ("vuggy") porosity commonly described in carbonates. I choose not to model this partition because the time delays between the geophysical measurements largely exceed the short-term delay involved in fast discharge.

Nevertheless at least one measurement was made very shortly after a flood event and some insight into fast transmission could be gained.

Finally we can compute for each layer the water saturation Sw(n) from the daily water height H(n)

and the thickness e. I assume that total porosity  is homogenous over the considered rock volume and we have for a surface area of 1 m 2 ,

𝑆𝑤(𝑛) = 𝐻(𝑛)/(1000𝜙 𝑒) , ( 18
)
when H is given in mm (l/m 2 ). It is necessary to define the initial state (water height) of each layer to run the model; I assign full saturation (Sw=1) to the layers in January 1996, when temporary lakes and streams appeared in the vicinity of the survey location [START_REF] Bruxelles | Lacs temporaires et circulation de surface sur le Causse de l'Hospitalet du Larzac (Aveyron) en 1996 : fonctionnement et implications géomorphologiques[END_REF]. I use first a single layered velocity model computed from all the seismic first break picks gathered over a distance of approximately 150m; this coarse spatial scale is in the range of the sampling which could be used for building vulnerability maps for instance. Figure 6 displays an example of inversion result: the horizontal axis is calendar time, and from bottom to top, we can see the effective rain (bar graph, in mm), the water saturation, computed from the water height of Maillet law and the thickness provided by the geophysics, then an example of the estimate of the residence time R (in days) inside the soil reservoir, which is time varying to reflect the variable water input to the system.

The quantity R(n) for a single layer or the two layers is computed from the daily water volume H(n)

and the discharge Q(n) (expressed in mm/day) by finding the necessary time (in days) to drain completely the upper layers, namely for the two reservoirs,

𝐻1(𝑛) + 𝐻2(𝑛) = ∑ ∆𝑡 𝑄2(𝑖) 𝑛+𝑅(𝑛) 𝑖=𝑛+1
, with ∆t= 1day . (

Figure 6 displays also the computed seismic velocity time serie from the MAP solution and the experimental data for both the soil and the lower reservoir. The columns M1, M2 and H1 of Table 1 display the results for three different locations chosen according to surface morphology features and of seismic velocity variations.

The results for the soil reservoir indicates that porosity is significant, between 0.24 and 0.40, and the mixing of water and air is at a fine scale because permeability is large enough to yield homogenous saturation; consequently seismic velocity decreases with saturation and this decrease is large because there is a very significant chemical weakening of the rock frame when introducing water in the pores. The situation at M2 is completely different; porosity is much lower, the saturation is patchy, an indication of lower permeability and velocity increases with saturation. This specific behaviour is logically related to the presence of clay in the soil reservoir at this location. For the deeper reservoir, the saturation model fitting the data is always the patchy model, in response to the permeability drop occurring at the base of the epikarst, and chemical weakening is almost insignificant. Despite the differences in thickness and saturation type, the average residence time at locations M1, M2 and H1 inside the model is always larger than three months, with a decrease during high flow events and an increase during dry periods. It is possible to make the inversion at a finer spatial scale; namely for location X108 in Table 1 a local layered model is built from velocities and thicknesses extracted locally from time-lapse seismic tomography velocity sections. In this case the horizontal analysis scale is about one seismic wavelength, i.e. five to ten metres. Location X108 was chosen because it is located just outside the rim of a sinkhole, there is a cave entrance nearby, and a low velocity anomaly is present and extends inside the LPV [START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF]. The result at this location is peculiar, in the sense that saturation is always homogenous throughout the model, even at depth, an indication that permeability is higher than in previous locations, and consequently the decay constants and the residence time are much smaller. Location X108 is likely a preferred infiltration pathway through the LPV.

Discussion

I review now some of the main issues involved in the uncertainties affecting the results of the inversions. The first cause of uncertainty is the precision in the input velocity data, which is itself the result of an inversion procedure. Confidence intervals can be defined when using least-square fitting for standard refraction, but it is not straightforward for tomographic inversion. In the latter case I simply ran a few inversions with randomly perturbed travel times and computed separately an average standard deviation over all cells of the low velocity layer and its substratum. From these tests I set the velocity uncertainty to 5% and this value was incorporated in the Bayesian inversion.

Another source of uncertainty is the petrophysical model itself, because I made several simplifying assumptions in order to get a tractable model. The most critical issue is perhaps the dry-rock properties needed to perform Gassmann fluid substitution. The robustness of the model was tested against residual saturation, but the model is likely to fail for residual water in excess of 10%. I used coarse determinations of Vs velocities from surface waves to find the driest measurement and choose a single representative value of dry Vp for each layer. Lastly, the model does not handle the most complicated porosity patterns, when the saturation switches from homogenous to patchy during imbibition (strong dual porosity). The model handles homogenous or patchy saturations, not a combination of both. I postulate this case is quite unlikely to happen at seismic frequencies at least in the highly permeable epikarst.

Regarding the hydrologic conceptual model used in the previous section, the water input to the epikarst is not accurately known. The actual evapotranspiration may be quite different from the predicted PET, and fast infiltration is not considered, therefore the daily volume of water entering the system could be overestimated. The epikarst model is simpler than that of [START_REF] Fleury | Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France)[END_REF] or [START_REF] Tritz | Modelling the behaviour of a karst system catchment using non-linear hysteretic conceptual model[END_REF] because it is desirable to limit the number of parameters to fit, given the scarcity of the data; most of the complexity introduced by hydrogeologists is generally motivated by the need to allow non-linearity to model the hydrograph at the outlet, which is not the goal in the present work.

One other point is that I assume implicitly a pure vertical flow through the epikarst, but the actual flow may be quite different with some lateral component toward the nearest vertical pathway through the LPV.

Despite all the shortcomings that may exist, the model is able to provide useful insight to the water transfer at shallow depth. It is not within the scope of this work to discuss in details the hydrologic significance of the results for each location; this paper only shows that they are consistent with local surface geologic features and the general functioning of epikarst and LPV. The results show that the velocities decrease after rainfall in the uppermost high permeability layer of the epikarst, in agreement with the theoretical framework [START_REF] Gassmann | Über die elastizität poröser medien[END_REF]Biot, 1956a) of partial saturation in porous rocks at low frequencies. Likewise lower permeability clayey sands or the layer overlying the LPV exhibits the patchy saturation behaviour (the P-waves velocity increases after the rain) because air pockets are trapped and their size is greater than the maximum size allowed at seismic frequencies to consider iso-stress conditions, a prerequisite for the fine-scale mixing effective fluid model in Biot-Gassmann theory. In any cases the shallow layers are not fully saturated, in agreement with geochemical data showing that the unsaturated zone is an atmospheric-open system [START_REF] Aquilina | Water storage and transfer in the epikarst of karstic systems during high flow periods[END_REF][START_REF] Peyraube | Geochemical evolution of groundwater in the unsaturated zone of a karstic massif, using the relationship[END_REF]. The results show too that the main epikarstic water buffer resides near the base of the low velocity layer where the permeability drops [START_REF] Perrin | Epikarst storage in a karst aquifer: a conceptual model based on isotopic data, Milandre test site, Switzerland[END_REF], and it is possible to provide quantitative estimates of water volumes. The computed values of residence times in the epikarst (several months) are well within the estimations commonly reported by geochemical tracers [START_REF] Aquilina | Water storage and transfer in the epikarst of karstic systems during high flow periods[END_REF], but coupling a simple hydrologic model with the geophysical measurements allow to spatially locate these values. The provided example of short residence time at X108 corresponds in the field to a nearby sinkhole and vertical cave, two surface morphologic features indicating fast infiltration. The seasonal variability of the transit time through the epikarst found in this work is in agreement with classical hydrologic investigations using the analysis of piezometric levels [START_REF] Delbart | Temporal variability of karst aquifer response time established by the sliding-windows crosscorrelation method[END_REF]. Nevertheless these studies consider the unsaturated zone as a whole, while the present results indicate that this variability originates from the epikarst.

It is not the first time that simple hydrological modelling is coupled to geophysical measurements in karst studies. [START_REF] Deville | On the impact of topography and building mask on time varying gravity due to local hydrology[END_REF] attempted to invert for the parameters of a simple two reservoirs model from three years gravimetric time series, but reported that it was difficult to correlate the results with surface geological features. [START_REF] Mazzilli | On the inclusion of ground-based gravity measurements to the calibration process of a global rainfall-discharge reservoir model: the case of the Durzon karst sysem (Larzac, southern France)[END_REF] showed with the same data that the assimilation of gravity information was actually degrading the prediction of the terminal spring discharge using a rather complex two reservoirs model with temporary secondary outlets. This failure is probably related among other reasons to the fact that gravimetric measurements provide a global water budget and it is difficult to infer the vertical localization of the water storage. Using seismic data overcomes this limitation because the results are intrinsically spatially localized with a theoretical resolution related to the wavelength, and this work shows that the results are strongly correlated to surface morphology. It seems therefore possible to use seismic velocities in the shallow layers of the karst as quantitative proxies for temporary water storage volumes and transfer rates, two hydrologic variables which in turn could be used for larger scale hydrologic modelling or the mapping of karstification intensity and infiltration potential during environmental studies [START_REF] Kavouri | PaPRIKa: a method for estimating karst resource and source vulnerability-application to the Ouysse karst system (southwest France)[END_REF].

I conclude the discussion with a few words about the underlying LPV. The present work uses only first breaks seismic data as input (actually the velocities inverted from first breaks), which is suitable for the study of the shallow low velocity layer and its fast basement. Reliable velocity models for the deeper karst are difficult to build with traveltimes only [START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF], and it is necessary to use the amplitudes of the seismic wavefield to achieve this goal. Seismic attenuation is theoretically a better proxy for the study of the deep partially saturated LPV, because WIFF-related attenuation is enhanced by the low permeability of the medium and its partial saturation. Consistent increase of attenuation with saturation has been observed for the seismic data used in this paper [START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF] and the result was confirmed with borehole measurements. However, extracting reliable and quantitative attenuation information is still challenging, and models are not yet very tractable and robust to relate attenuation to permeability and partial saturation.

Conclusions

Estimating the seasonal variations of water saturation in the epikarst is possible when combining time-lapse seismic refraction data with a suitable petrophysical model. It is necessary to introduce a chemo-mechanical effect in the Biot-Gassmann fluid substitution model to account for the large velocity variations observed on field data; this rock frame weakening is probably the consequence of the chemical interaction between water and carbonates. Inversion of the time-lapse seismic velocities for water content and rock properties is performed using Bayesian inference and the modified Biot-Gassmann model with a chemical compliance parameter. Testing the procedure with laboratory measurements shows that the water content estimated from dry and partially saturated velocities is consistent with the experimental data. A real case application along a seismic profile shows a detailed picture of the water buffering in the epikarst; the storage locations agree with boreholes information and the conceptual model of a shallow temporary karstic reservoir.

The residence times of the water inside the epikarstic reservoir can be estimated if hydrologic modelling is further introduced and coupled with the time-lapse seismic data and the petrophysical model. The epikarst is jointly described by the seismic layering and a conceptual hydrologic model using simple flow equations, and the daily water input in the layers is computed from the rainfall chronicles. This water input provides a constraint for the inversion of four sets of seismic velocities corresponding to four hydrodynamic states during the year. The inverted parameters are in this case porosity, saturation type and water residence time for each layer. Numerical results show that the epikarst must be divided into two layers with different seismic properties: a high permeability upper layer where homogeneous air-water mixing occurs (homogeneous for the low frequency scale of surface seismic), and a second deeper low permeability layer where patchy saturation applies.

Calculations have been performed at different locations chosen for their characteristic karstic morphologies, and the values of residence times agree with these features and the results reported in the literature using geochemical tracers.

Finally, this work has shown that time-lapse seismic velocities are proxies for hydrologic parameters used in applications like quantitative vulnerability mapping of aquifers or distributed flow modelling. with the upper and lower bounds of the 95% credibility interval. The threshold constant h is poorly constrained. Last row displays the average residence time of the water in the model (in days), and the higher and lower values during the four years period.

  found during his time-lapse study of several karst systems that velocity was consistently decreasing by a large amount in the upper epikarst during the wet season, while the lower epikarst was showing a velocity increase. I use the velocity models from a seismic line shot during April 2009 (end of wet season) and repeated in October six months later at the end of the dry summer. Details about the survey, its location and the tomographic inversions are reported in[START_REF] Galibert | Seismic study of the lowpermeability volume in southern France karst systems[END_REF]; the line is laid out on the lower slopes of a dry valley carved in a limestone plateau. The local geology consists of a shallow layer of weathered brecciated limestone (epikarst) overlying micritic limestones with intercalations of thin marly beds (LPV). Both surveys included compression and shear waves records, therefore we can apply the inversion procedure using Vp and Vs data. The wet state Vp model from the baseline tomographic inversion and the velocity variations during summer are displayed in Figure4(top and middle). We can see that Vp is increasing during the dry summer almost everywhere at shallow depth, in response to the drying-up of the sub-surface when rainfall is scarce and the stored water percolates downward to the LPV. There is however an anomaly near abscissa 85m showing a decrease in Vp; when using the Vs velocity variations in each cell of the model to identify locally without ambiguity the driest situation, this anomaly corresponds to a drier state (higher Vs) during April.

  I use daily Rand PET values from a Meteo-France weather station located 15 km to the south of the site. Once the daily saturation is computed, the final step of the forward model is to compute the seismic velocities from Sw,  and the elastic properties of the rock frame using the modified Gassmann model.I invert this time for the parameters p=(, Ch, T 1 ,h 1 , T 2 , h2) of the epikarst model, i.e. respectively total porosity, chemical interaction, and Maillet's decay times and saturation thresholds for the two epikarst layers, by fitting the modelled velocities to the experimental velocities d, the daily rain providing an external constraint driving the model. The experimental data d consists of the four seismic velocities Vp and Vs computed for each layer in four different hydrodynamic conditions for the same rock volume. The Vs information is actually provided by simple surface-waves inversion to identify the driest state; it is not included in the cost function, and I assume a constant Vp/Vs ratio based on calcite properties for the forward computations. The MCMC sampling building the probability distribution uses 50000 random draws of the parameter vector p for each type of saturation models, i.e. the model space includes homogenous and patchy saturations, and therefore the saturation type is another inverted parameter.

Figure 1 Figure 2

 12 Figure 1 Example of acoustic velocity variations with partial saturation. Experimental data is from

Figure 3

 3 Figure 3 Bayesian inference of water content (see text). The experimental data is for Estaillades

Figure 4

 4 Figure 4 Example of water content estimation in the epikarst from seismic time-lapse results. Top: P-

Figure 5

 5 Figure 5 Conceptual model of epikarst (left) and the dual layer characterisation from refraction

Figure 6

 6 Figure 6 Determination of the parameters for the conceptual model of the epikarst from seismic

Table 1

 1 This conclusion is likely valid for different geological settings if the petrophysical model is adapted. Estimated parameters of the conceptual epikarst model. The decay constant T is provided

		M1	M2	H1	X108
	Surface lithology	Sand	Clayey sand	Dolostones	Infiltration
		Dolostones			point
	Soil reservoir	8m	2.5m	3.5m	2m
	Saturation	homogenous	patchy	homogenous	homogenous
	Porosity	0.40	0.10	0.35	0.28
	T days (lower, upper)	25 (20, 26)	2 (1, 3)	5 (4, 6)	6 (5, 7)
	h (mm)	-	0	1	1
	Lower reservoir	20m	11m	11m	13m
	Saturation	patchy	patchy	patchy	homogenous
	Porosity	0.09	0.10	0.05	0.10
	T days (lower, upper)	51 (17, 121)	84 (62,91)	86 (62, 240)	6 (5, 9)
	h (mm)	10	-	9	3
	R days (lower, upper)	97 (27,189)	90 (15,154)	101 (23,174)	18 (5, 73)

Acknowledgments

This work was partially supported by the program Hydrokarst G2 from the Agence Nationale de la Recherche, France. I am grateful to Thierry Cadoret for providing the numerical results of his laboratory sonic measurements. I thank the anonymous reviewers for their comments and suggestions that helped improve this paper.