# Data-model comparisons of isotopic and hydrological variabilities of the karstic vadose zone above Villars Cave, SW-France based on 20 years' monitoring record

Jian Zhang, Dominique Genty, Colette Sirieix, Cécile Verdet, Sylvain Mateo, Ludovic Devaux, Giovanni Sgubin, Stéphane Bujan, Bénédicte Minster, Edouard Régnier, et al.

► To cite this version:

Jian Zhang, Dominique Genty, Colette Sirieix, Cécile Verdet, Sylvain Mateo, et al. Data-model comparisons of isotopic and hydrological variabilities of the karstic vadose zone above Villars Cave, SW-France based on 20 years' monitoring record. Chemical Geology, 2023, pp.121802. 10.1016/j.chemgeo.2023.121802. hal-04264851

HAL Id: hal-04264851

https://hal.science/hal-04264851v1

Submitted on 28 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Data-model comparisons of isotopic and hydrological variabilities of**
2 **the karstic vadose zone above Villars Cave, SW-France based on 20**
3 **years' monitoring record**

4 Jian Zhang<sup>a, b, c*</sup>, Dominique Genty<sup>a*</sup>, Colette Sirieix<sup>b, c</sup>, Cécile Verdet<sup>b, c</sup>, Sylvain

5 Mateo<sup>b, c</sup>, Ludovic Devaux<sup>a</sup>, Giovanni Sgubin<sup>a, d</sup>, Stéphane Bujan<sup>a</sup>, Bénédicte

6 Minster<sup>e</sup>, Edouard Régnier<sup>e</sup>, Ting-Yong Li<sup>f</sup>, François Bourges<sup>g</sup>, Ming-Qiang

7 Liang<sup>h</sup>, Chaojun Chen<sup>f</sup>, Vincent Marieu<sup>a</sup>

8 <sup>a</sup> *Environnements et Paléoenvironnements Océaniques et Continentaux (EPOC)*, UMR CNRS,

9 5805, Université de Bordeaux, 33615 Pessac Cedex, France

10 <sup>b</sup> Université de Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France

11 <sup>c</sup> Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, INRAE, I2M, UMR 5295, F-33400

12 Talence, France

13 <sup>d</sup> Dipartimento di Scienze della Terra e del Mare, Università degli Studi di Palermo, via Archirafi

14 22, 90123 Palermo, Italy

15 <sup>e</sup> Laboratoire de Sciences du Climat et de l'Environnement, Université Paris-Saclay, Orme des

16 Merisiers, 91191, Gif-sur-Yvette Cedex, France

17 <sup>f</sup> Yunnan Key Laboratory of Plateau Geographical Processes & Environmental Changes, Faculty

18 of Geography, Yunnan Normal University, 650500, Kunming, China

19 <sup>g</sup> Géologie Environnement Conseil, 09200, Saint-Girons, France

20 <sup>h</sup> Georges Lemaître Center for Earth and Climate Research, Earth and Life Institute, Université

21 catholique de Louvain, 1348, Louvain-la-Neuve, Belgium

22

23

24

25

26

27

28 * The corresponding authors: Jian Zhang and Dominique Genty

29 E-mail: jian.zhang@u-bordeaux.fr; dominique.genty@u-bordeaux.fr

30 Postal address : Allée Geoffroy Saint-Hilaire Bât. B18N, Université de Bordeaux,

31 33615 Pessac Cedex, France
32 **Abstract**

Because of the great heterogeneity of the limestone fissures, conduits, and reservoirs in karst areas, the understanding of the infiltration in caves and karst terrains is complicated. For example, the time delays observed between rainfall and the various flows into caves, such as shaft flow and seepage flow, exhibit significant variability. To better understand the infiltration above a typical shallow cave, Villars Cave, South-west France, we used both a very long monitoring (>20 years) isotopic data set of rainfall and dripping rates at two different cave gallery levels (14 m and 36 m below the surface, respectively) and a conceptual model (KarstFor) that mix water flows and isotopes. In order to determine the most appropriate reservoir and flow settings, we used one time-lapse electrical resistivity tomography (ERT) image made above the cave, which were incorporated into the KarstFor model. The modeled variations in storage water and infiltration flows are in a good agreement with the observed dripping discharge rates, indicating seasonal variability and various time delays between water excess ((Rainfall-Evapotranspiration)>0) and drip rates in the different cave levels. The fact that the model strikingly simulates the small, but constant, isotopic difference of drip water $\delta^{18}$O ($\delta^{18}$O$_d$) values between the upper and the lower galleries of Villars Cave, not only reinforces reliability in the model itself but also reveals the impact of multiple infiltration routes and karst reservoir dynamics on the $\delta^{18}$O$_d$. On an inter-annual scale, the observed $\delta^{18}$O$_d$ values follow the variations of the water excess and of the cave temperature, which, in the period analyzed, show a marked trend, possibly induced by the climate variability patterns above the cave. Overall, these results are not only important to better understand the infiltration in karstic zones, but also to better interpret speleothem isotopic records widely used in paleo-climate studies.

57

**Keywords:** Karst reservoirs; Infiltration; Unsaturated zones; Stable isotopes; Electrical resistivity tomography; Drip rate; Speleothems

60 **Main manuscript content**

61 **1. Introduction**

The process of rainfall infiltration in karst zones is challenging to describe and model, despite the fact that this knowledge is of utmost importance for preserving freshwater quality and protecting cave heritage (Ford and Williams, 1989; Hartmann et al., 2012; Hartmann and Baker, 2017; Bourges et al., 2014; 2020). Modeling the infiltration zone is essential for understanding the transmission of the climatic and geochemical/isotopic signals from the outside to cave deposits (speleothems, cave sediments) which provide remarkable archives for paleo-climatological reconstructions (Baker et al., 2013).
Cave monitoring studies have been conducted to better understand the relationship between rainfall and seepage water in the karst zone in several cave sites (Kaufman et al., 2003; Spötl et al., 2005; Cruz et al., 2005; Fairchild et al., 2006; Pape et al., 2010; Mattey et al., 2008; 2010; Kluge et al., 2008, 2010; Genty et al., 2014; Rau et al., 2015; Markowska et al., 2015; Duan et al., 2016; Baker et al., 2019; 2020). It is recognized that the infiltration duration and its routing, controlled by the reservoirs and infiltration types, influence cave drip water oxygen isotopes (δ¹⁸O_d) value. Generally, due to mixing in the soil zone and epikarst, δ¹⁸O_d variability is damped compared to rainfall δ¹⁸O (δ¹⁸O_r), and it is usually considered to be close to the annual weighted (rainfall amount) average δ¹⁸O_r (Yonge et al., 1985; McDermott, 2004; Li et al., 2011; Riechelmann et al., 2011; Zhang and Li, 2019). In some cases, where there is less storage and faster flow, the seasonality is observed in δ¹⁸O_d values in the Mawmluh Cave, India (Breitenbach et al., 2015) and the Black Chasm Cavern, California (Oster et al., 2012). Recently, it has been suggested that δ¹⁸O_d variations may be considered an averaged value of weighted δ¹⁸O_r for a much longer period up to several years, like in Villars Cave (SW-France) (Genty et al., 2014; Jean-Baptiste et al., 2019). It is worth pointing out that more recently, Ellis et al. 2020 used a mixing model to make δ¹⁸O_d variations closely match δ¹⁸O_r variations to calculate the residence time in epikarst, ranging from 3 to 18 months in Mulu Caves, Northern Borneo. Consequently, within-cave responses to outside climate changes may be significantly variable, depending on the geometry and the dynamics of reservoirs above the cave. Rainfall infiltration in karst terrains is heterogeneous and combines slow infiltrations and fast ones that may be a mixture of “old and young” water in reservoirs, resulting in different responses of the δ¹⁸O_d. The δ¹⁸O_d may represents various outside-cave climate dynamics: 1) seasonality (Oster et al., 2012); 2) atmospheric circulation patterns, such as El Niño-Southern Oscillation (ENSO) (McDonald et al., 2004; Tadros et al., 2016; Chen and Li, 2018; Zhang and Li, 2019; Ellis et al., 2020; Zhang et al., 2022; Xue et al., 2023); 3) heavy rain events in the wet season (Pape et al., 2010), or 4) the degree of enrichment of drip water evaporation between replenishment events in arid regions (Cuthbert et al., 2014).

Modeling infiltration through the reservoir in karst terrains has been a challenging topic. For example, based on a lumped parameter approach to model the complexities of karst hydrogeology, the KarstFor model was developed by Bradley et al. (2010) and enhanced with multiple adaptations in the different climate systems, such as monsoon, temperate, and Mediterranean climates. The initial schematic illustration of the KarstFor model consisted of climatic inputs and a monthly time-step for three reservoirs: the soil, the epikarst, and the vadose zone (Bradley et al., 2010). Cuthbert et al. (2014) previously applied a single-reservoir (soil moisture balance – karst) model to the drip recharge, and this method was then used to compare modeled and measured δ¹⁸O values for two modern speleothems (Markowska et al., 2020). Meanwhile, Dee et al. (2015) proposed a simple convolution of the δ¹⁸O_r time series
to simulate $\delta^{18}$O$_d$ values using drip water transit time distributions (pulse response function) in their PRoxY System Model (PRYSM) (Dee et al., 2015; Hu et al., 2021).

Regarding cave drip discharge, Healy (2010) reported the use of the high frequency and continuous discharge time series of drip water infiltrating into caves as a physical method for quantifying rainfall thresholds in the karstic unsaturated zone (Healy, 2010). Automated drip rate loggers monitor hydrological variability in karst zones, capturing continuous data in infiltration zones and enabling recharge event comparisons (Baker et al., 2020). Generally, in the time frame of several months to several years, the fluctuations of drip flow are usually driven by seasonal changes in rainfall availability (Hu et al., 2008; Genty and Deflandre 1998). However, on a daily to weekly time scale, the drip rate may respond to a single rainfall event (Baldini et al., 2012) and cave (and outside) air pressure (Genty and Deflandre, 1998; Tremaine and Froelich, 2013). Flow rates may also influence speleothem climate proxies, such as $\delta^{18}$O values and trace elements, through the dilution and mixing of percolation waters before reaching the cave (Treble et al., 2013). Despite the very few studies that compare the measurement of $\delta^{18}$O$_d$ values with changes in drip rate, we observed that in most cases there is no correlation between both measurements (Genty et al., 2014; Riechelmann et al., 2017).

The exceptionally long monitoring record of $\delta^{18}$O$_r$ and $\delta^{18}$O$_d$ since 1997 (over 20 years) in Villars Cave, SW-France, clearly provides insights into the study of the relationship between rainfall, drip water, and climate parameters. This research is a continuation of the earlier work presented by Genty et al. (2014) with extended $\delta^{18}$O$_r$ and $\delta^{18}$O$_d$ data to 2019 coupled with an infiltration conceptual-model analysis (KarstFor). Although prior works in Villars Cave have provided insight into infiltration seasonality, key questions remain about reservoir dynamics and geometry. By utilizing 2D time-lapse electrical resistivity tomography (ERT), a non-destructive and minimally invasive method (Xu et al., 2015; Leopold et al., 2021), we were able to add a clear setting of the geometry of the reservoirs along a profile above the cave.

Thanks to the information provided by the ERT regarding the numbers and locations of the reservoirs, we were able to significantly improve the settings of KarstFor. This better-informed version of KarstFor was able to reasonably reproduce the observed data, leading to a better identification of the processes driving the long-term trends. Consequently, several questions were addressed in this study: (1) How do the modeled drip rates and $\delta^{18}$O values align with actual cave observations, (2) What are the implications in term of reservoir dynamics (water time transfer, mixing, overflows) ?

## 2. Study area and monitoring sites

The Villars Cave (45°26′18″N, 0°47′2″E), SW-France, belongs to the maritime climate region with mild winters and warm summers (Figure 1A). During the monitoring period (1997–2019 AD) the annual precipitation (P) above the Villars
Cave is 1014 mm (Standard Deviation: SD =197 mm), and the mean annual temperature is 12 °C (SD = 5.6 °C) (Genty, 2008; Genty et al., 2014). Local temperature in the region of Villars Cave is characterized by a strong seasonal-to-interannual variability (Genty et al. 2014), which is primarily influenced by large-scale circulation patterns, such as the North Atlantic Oscillation (NAO) and the Atlantic Multidecadal Oscillation (AMO) (Zhang et al., 2020). The NAO, an atmospheric mode of variability of flow patterns over the North Atlantic, exerts a significant influence on weather conditions in the North and Western Europe, notably during the winter season (Hurrel et al. 2013). In contrast, the AMO, an oceanic mode of variability affecting the Atlantic Sea Surface Temperature, significantly impacts Western Europe during the summer. These large-scale patterns are therefore essential drivers of seasonal-to-interannual variability in the Villars Cave region, influencing local temperature and moisture levels (Zhang et al., 2020).

Villars Cave is developed within Middle Jurassic limestone, specifically the Bajocian and Bathonian formations (Genty, 2008). It is situated less than 3 kilometers south of the crystalline rocks of the limestones. The surrounding landscape features hills with a height of approximately 50 meters, reaching an elevation of around 200 meters above sea level. The intricate galleries of Villars Cave are characterized by their small size and muddy nature. They form a complex network that spans approximately 10 kilometers in length. The direction of the galleries is influenced by St. Pardoux fault, which acts as the boundary between limestone and shale rock formations located a few kilometers north of Villars Cave. This fault follows a southeast-northwest direction, influencing the development of the cave system (Genty, 2008; Genty et al., 2014). The vegetation above the Villars Cave is mainly composed by oak and hornbeam tree. The soil thickness above the cave is unevenly distributed ranging from 20 to 40 centimeters (Genty, 2008; Genty et al., 2014).

## 3. Methods and data

### 3.1 Meteorological data

During the monitoring period, spanning from 1997 to 2019, we collected external air temperature and precipitation data from a meteorological station located just above the cave. The potential evapotranspiration (ETP) and water excess ((Rainfall-Evapotranspiration) >0, also called efficient rainfall) have been calculated using meteorological measurements provided by the nearby Meteo-France station (Saint-Martin de Freissengeas, 4.8 km East) during the period 2006-2019. The Penman-Monteith method (Allen et al., 1998) for estimating ETP incorporates variables such as minimum and maximum temperature, average wind speed at 10 m, average vapor pressure, and global radiation, which are crucial factors in determining the water loss from the surface through evaporation and transpiration.
## 3.2 Drip rate

Inside the cave, we monitor four stations, namely Vil-#10A, Vil-#10B in the upper galleries (14.38 ±0.23 m depth) and Vil-#1A, and Vil-#1B in the lower galleries (36.31 m±0.48 m depth) (Zhang et al., 2023). Vil-#10A and Vil-#10B are situated approximately 4 meters apart while Vil-#1B and Vil-#1A are separated by 0.5 meters (Figure 1C; Genty, 2008). The depth of the different stations was determined by using magnetic probes which gave very accurate results (Appendix Figure 1). On a horizontal plane, the Vil-#10A station is positioned 4.3 meters west away from the ERT profile (Figure 1C; Appendix Figure 1). The locations of Vil-#1B and Vil-#3 are situated approximately 100 meters east away from the ERT profile presented in Figure 1C and Appendix Figure 1.

Drip rates, presented in this work, were recorded between October 2013 and April 2019, which coincides with a continuous period without any stop due to technical problem. Drip loggers (name: STALAGMATE counter) were used to record the frequency of drip rate onto the surface of a sealed box containing an acoustic sensor.

## 3.3 Stable isotopic compositions of rainwater and drip water

Rainfall and drip water were collected for δ¹⁸O-δD isotopes measurements from 1997 to 2019 (Genty et al., 2014; Zhang et al., 2020). The δ¹⁸O time series from 2006 to 2019 provides the most continuous data, with sampling intervals occurring on average every 1-2 months. There were intermittent gaps in the data collection between 2003 and 2006. However, before 2003, the sampling intervals were conducted approximately every 1 to 2 months. Water samples were sent to the LSCE laboratory (Laboratoire de Sciences du Climat et de l'Environnement, Gif-sur-Yvette, France) and analyzed following two methods. Rainfall and drip water δD and δ¹⁸O were measured using a PICARRO laser spectrometer analyzer, with the same standards, leading to an accuracy of 0.5 ‰ for δD and δ¹⁸O values. However, due to insufficient precision in the δ¹⁸O measurements for PICARRO laser spectrometer, we also used a MAT 252 mass spectrometer for δ¹⁸O measurements, which provided a 2-sigma error of ±0.05‰, significantly improving the accuracy compared to the PICARRO measurements.

## 3.4 ERT measurement

To visualize the unsaturated zone above Villars Cave, we made electrical resistivity tomography (ERT) images made in January 2021 that display more or less humid zones in the unsaturated zone. In this experiment, 72 electrodes, with an electrode spacing of 1m, resulting in a total length of 71m, have been planted in the soil above the Villars Cave. Then, two electrical arrays were used Gradient Array, and Pole-Dipole Array (forward and reverse) (Xu et al., 2017; Figure 1C). We take into
account the topography of the profile at each electrode before the inversion of the data.
A correction has been applied to minimize the effects of the remote electrode (Xu et
al., 2015). The inverted resistivity models were calculated using the Res2Dinv®
software, employing the L1 norm (robust inversion) method with the refinement of
the mesh model, as described by Loke (2004). Additionally, we applied topography
and temperature corrections to the ERT images, following the methodology outlined
by Xu et al. (2017), as these corrections are essential to ensure accurate interpretation
of the data by considering the influence of surface features and temperature
fluctuations, respectively.

The variation of the resistivity is essentially linked to the variation of the water
content of the rock (Archie law) (Xu et al., 2017). Therefore, a variation of resistivity
is linked to the hydrological change in the unsaturated zone (Xu et al., 2015; Leopold
et al., 2021; Figure 2A). In this article, we show one ERT image whose purpose is
only to identify the position and the number of reservoirs that can be linked to the
conceptual model (KarstFor) we used (Zhang, 2023 and Zhang et al., 2023; Figure
2A).

3.5 Modeling infiltration process and lag time calculation

Drip water hydrological flux and δ¹⁸O value were modeled using a version of
a karst hydrology model (KarstFor) first developed by Bradley et al., 2010. This
model has been applied to multiple climatic regions so that we can determine the
theoretical range of drip δ¹⁸O in some locations: northern temperate (Tartair Cave,
Scotland), monsoonal (Ethiopia), and Mediterranean (St. Michaels Cave, Gibraltar,
South Spain Peninsula), sub-humid to semi-arid (central Texas caves, United States of
America) (Bradley et al., 2010; Baker et al., 2013; Treble et al., 2013; 2019; Hu et al.,
2021).

Here, our final configuration parameter settings are based on Bradley et al., (2010), with some adjustments made in the context of the temperate oceanic climate (Table 1). Four storage reservoirs were set in the conceptual model: Soil store, Epikarst store, Karst store 1, and Karst store 2 (Figure 2B) based on identified reservoirs from the ERT image (Figure 2A).

All the main hydrological parameters of each reservoir are summarized in Table 1. KarstFor model is based on a local physical understanding of each parameter and then adjusting the parameters one by one. For example, we set the different drainage fluxes (F1, F4, F5) and overflow (Ovicap and Epicap, Table 1) of karst reservoirs to control the variance of δ¹⁸O in the Epikarst and Karst stores 1 and 2. The term “Ovicap” is used to refer to the setting parameter that represents the overflow capacity threshold of the karst store. Conversely, “Epicap” pertains to the setting parameter that represents the overflow capacity threshold of the epikarst (Treble et al., 2013). The high/low proportion of karst water available for evapotranspiration (k_eevap) constrains the mean value of δ¹⁸O<sub>d</sub> (Figure 2B; Bradley et al., 2010; Treble et al., 2013; 2019; Hu et al., 2021).
The former work in modelling infiltration suggested that (1) the rainfall effectively infiltrates (water excess) over a longer period than expected (i.e., most of the year) and (2) the average age of the Villars drip water was first determined by Genty et al., 2014 who used a simple isotopic mixing model using rainfall amount and rainfall $\delta^{18}$O value for the input. Later, Jean-Baptiste et al., 2019 used tritium measurements in drip water combined with $\delta^{18}$O$_r$ value to produce a more precise model. Thus, the average transition time, between rainfall and drip stations, has been evaluated to be 7.0±0.5 years for the upper gallery and 11.1±0.5 years for the lower gallery. (3) the previous simulation model of the infiltration was a simple 1-D model of the downward water percolation without considering the complex fact that there may be multiple reservoirs in the unsaturated zone (Jean-Baptiste et al., 2019). (4) previous studies have observed a long-term increase in $\delta^{18}$O$_d$ variation, but it is unclear why this occurs (Genty et al., 2014).

Our new modeling approach incorporates multiple reservoirs to represent the soil, epikarst, and unsaturated zone. Additionally, we utilized monthly climate data as inputs to simulate the oxygen isotopes in drip water. KarstFor simulations include $\delta^{18}$O fractionation during evaporation, mixing water in matrix flow with conduit flow in the reservoirs, and overflows. The model simulates rainfall infiltration through the soil reservoir, epikarst reservoir, and two underlying reservoirs before reaching the top of the cave. As a result, this new model incorporates several additional parameters and reservoirs, leading to improved results compared to the previous models proposed by Genty et al. (2014) and Jean-Baptiste et al. (2019).

In this study, the variable infiltration routing of percolating waters is conceived between two sub-surface water stores (Soil Store and Epikarst Store, Figure 2B) and two lower-level karst water reservoirs (Karst Store 1 and 2). Each of them potentially represents the source of drip water in the cave (Bradley et al., 2010; Treble et al., 2013; 2019). The flow route or hydrology dynamic of the conceptual model is as follows (Figure 2B):

(1) In the Soil Store, soil thickness is ranging from 20 cm to 40 cm based on field observations. We set an initial storage volume of only 50 mm (while it can go up to 150 mm for more developed soil (Treble et al., 2013)). The local precipitation amount is considered an input parameter into the “Soil Store” receiving precipitation directly. Outputs are evapotranspiration and drainage to Epikarst (F1).

(2) In the Epikarst, the hydrological flux receives drainage from Soil Store (F1), and outflows to Karst Store 2 (F4) or Karst Store 1 (F3) according to the proportional drainage function: Epikarst × (Epikarst/Epikarst capacity) × 0.25, i.e., flow is 25% of volume when Epikarst is at capacity.

(3) In Karst Store 2, the hydrological flux receives drainage from Soil Store and Epikarst. “Karst Store 2” receives F4 from “Epikarst” and F8 from a “bypass” (Table 1) when the drainage is over its capacity.

(4) Karst Store 1 receives hydrological fluxes provided by Epikarst (F3) and Karst
Store 2 (F6 and F7).

The δ¹⁸O value in each store is updated with the time step based on the water balance calculations and isotopic evaporation coefficient in terms of rainfall, evapotranspiration, drainage, and overflow. Three outputs, representing multiple karst configurations, are modeled similarly to Bradley et al., 2010. The karst component of the model follows the KarstFor version presented in Bradley et al., 2010, with water levels in the reservoirs (Soil Store, Epikarst Store, Karst Store 1, and Karst Store 2) given by an initial storage value (mm), and then recalculated (starting from the top store) at a monthly time step via:

$$ \Delta V = (\Sigma In - \Sigma Out) \Delta t \quad (1) $$

where ΔV = net change in monthly storage volume (mm), In = input flow, out = output flow (mm/month), Δt = 1 month. Besides the model analysis, we also studied the time lags between drip rate and water excess at the different drip sites.

To determine the lag time between drip rate and water excess at different drip sites, we employed the Pearson correlation coefficient, which is commonly used to assess the strength of a linear relationship between two variables (x and y). By calculating the maximum correlation coefficient (r) for specific days, we were able to identify the time lag in days between the drip rate and water excess. We utilized the following equation:

$$ r = \frac{\Sigma_i (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\Sigma_i (x_i - \bar{x})^2 \Sigma_i (y_i - \bar{y})^2}} \quad (2) $$

where, r represents the Pearson correlation coefficient, x represents water excess, y represents drip rate, and the number of lag days corresponds to the days with the highest correlation coefficient.

# 4. Results

## 4.1 Rainfall and drip water δ¹⁸O time series

The δ¹⁸O<sub>r</sub> values show clear seasonal changes, with higher values in summer and lower values in winter (Figure 3A; Genty et al. 2014). Based on the δ¹⁸O and δD values of the rainwater samples above the Villars Cave, we established an updated equation of the local meteoric water line (LMWL): δD= (7.01±0.12) δ¹⁸O+ (4.7±0.74) (n = 201), with a correlation coefficient R = 0.98 (p <0.01) (Figure 4; Zhang et al., 2020).

The δ¹⁸O<sub>d</sub> values for all drip water monitoring sites exhibit a range of -5.37‰ to -6.76‰, as shown in Table 2. Interestingly, Vil-#10A from the upper cave gallery exhibits a high variability in δ¹⁸O<sub>d</sub> values (Figure 3B), and shows a similar distribution with rainfall isotopic composition variations (Figure 4A), while the other three cave δ¹⁸O<sub>d</sub> time series display almost no fluctuations during years without any visible seasonal variability (Figures 3C and 4B). An important observation is that there is a mean difference of approximately 0.2‰ between the lower and the upper
348 galleries. This difference falls in the analytical uncertainty of ±0.05‰ for the δ¹⁸O
349 values (Figure 3C). While a difference of 0.2‰ may appear negligible, it can be
350 actually significant as shown in Figure 3C. Specifically, the lower gallery stations
351 (Vil-#1A and Vil-#1B) have the lowest mean values of δ¹⁸O<sub>d</sub>, which are -6.35‰
352 (SD=0.08, n=106) and -6.36‰ (SD=0.08, n=112), respectively. Simultaneously, the
353 upper gallery stations (Vil-#10A and Vil-#10B) have mean values of -6.13‰
354 (SD=1.33, n=129) and -6.19‰ (SD=0.21‰, n=96), respectively (Table 2). All these
355 values clearly show a large damping of the rainfall isotopic signal that is the
356 consequence of the mixing of the recharged water for several months to several years
357 (Genty et al., 2014; Jean-Baptiste et al., 2019).

## 4.2 Rainfall amount, evapotranspiration, and drip rate relationships

359 The climate around the Villars Cave is temperate and maritime. The local
360 monthly precipitation shows that there is no significant seasonality above the Villars
361 Cave (Appendix Figure 2). The mean annual value is 1014 mm in the 1997-2019
362 period (SD =192 mm). But there is a strong seasonal signal in the quantity of water
363 that infiltrates the karst, which can be estimated by the water excess (Figure 3E). The
364 main cause is the seasonality of the outside air temperature (Figure 3D), leading to
365 seasonality in the ETP (potential evapotranspiration) whose mean annual value is 837
366 mm (SD = 58 mm).

367 Between October 2013 and April 2019, drip loggers recorded a large variability
368 at different time-scales, but they are characterized by a well-marked variability
369 (Figure 5). The averaged drip rates over the entire hydrological years at the three drip
370 water sites (Vil-#10A, Vil-#10B, Vil-#1B) are 21, 3.5, and 16.5 drip/mins,
371 respectively, with a variability characterized by their standard deviation values (SD)
372 of 28.4, 3.7 and 13.9 drip/min, respectively (Figure 5). More specifically, during the
373 first cycle of 2014-2015, it was observed that there was a water excess of more than
374 zero in mid-September 2014 in Vil-#10A (Appendix Figure 3). It takes approximately
375 two to three months for the reservoirs' water to be filled up and reach the percolation
376 threshold, which then enables water to start moving down into the cave, forming drip
377 water. Consequently, there was a rise in drip water rate in December 2014, which
378 persisted until around May 2015 in Vil-#10A.

379 The time delay between water infiltration, once it reaches the percolation
380 threshold (quantified here as water excess), and the drip rates has been estimated
381 through cross-correlations (Figure 6). Results show that the time lag is in the interval
382 of 59-180 days for Vil-#1B station located in the lower gallery, 3-36 days for Vil-
383 #10B, and 4-19 days for Vil-#10A located in the upper gallery (Figure 6), showing
384 that time lags range between a few days and several months depending on the depth
385 of the drip stations and of its inertia linked to reservoir behaviors.

386 We also perform a manual one-to-one comparison by comparing peaks of the
387 different time series of water excess and drip rate in different years (Appendix Figures
388 3 and 4) to verify the cross-correlation results. Results, clearly visible on the time
series, show that the lag time between the water excess and Vil-10A# (upper gallery) drip rate is about 3-18 days (Appendix Figures 3 and 4), which is consistent with our cross-correlation results (4-19 days, Figure 6A-6E). The water excess peak, typically observed in January to February each year at the Vil-#1B station, leads to a time delay of generally more than three months before the maximum drip rate occurs around April to May (Appendix Figures 3 and 4; Figures 6F-6J).

## 4.3 ERT imaging of the karst reservoirs above the cave

The ERT resistivity model, generated on January 7, 2021 (Figure 2A), provides a comprehensive depiction of the various karst zones above Villars Cave. This model offers a broad visualization of the cave galleries and the saturation of the karst reservoirs (Figure 2A). The following information is found in this image: (1) The resistivity higher than 1000 Ω·m can be related to the unsaturated limestone. This value increases over 3240 Ω·m in three red high resistivity zones (altitude: 165-175 m), which is related to the upper galleries of Villars Cave (Figure 2A); (2) The low green resistivity ranges from 48 Ω·m to 550 Ω·m, indicating areas that can be filled with more water in winter, likely soil, epikarst, and some reservoir zones (Figure 2A).

When comparing the cave topographical map with the ERT image (Figures 1C & 2A), we observe a very good agreement between cave galleries (map) and high resistivity zones (ERT). In particular, Figure 2A shows three red zones (Zones A-C) with high resistivity values, which correspond to the three cave galleries depicted in Figure 1C. However, the ERT image does not provide precise information about the exact location of the lower gallery stations (Vil-#1A and Vil-#1B). These stations are situated in a lower vertical section, potentially located near the side of Karst Store 1 and at a greater depth (Figure 2B).

## 4.4 Simulations of the flowing rate of drip water

The links between different flow paths and water reservoirs can be described as follows: (1) F1 flow is from Soil Store, and outflows to Epikarst. (2) Karst store 2 can receive water flow from F4 and F8 and outflows in F7. (3) Karst store 1 can receive water flow from F7 and drain in F5. Modeled variations of the quantity of water stored in each reservoir during the 2006-2019 period is plotted in Figures 7C, 7D, and 7E. Results show, firstly, that there are considerable differences in water storage dynamics between individual water reservoirs. Secondly, the seasonal flowing successfully simulates higher/lower hydrological response discharge, agreeing with observed seasonal drip rates. Both reach a maximum drip rate, following winter water excess high values (Figures 5 & 7).

Compared to the behavior of the different reservoirs (Appendix Figure 5), the modeled output with the same vertical axis shows differences in the relative magnitude of individual water fluxes and obvious seasonal variability in all reservoirs. The Soil Store exhibits significant seasonality, ranging from 0 to 9 mm in magnitude,
and serves as the primary source for the seasonal variations in the following reservoirs (Appendix Figure 5). Karst Store 1 shows a seasonal variation from 0 to 60 mm, and the magnitude variability of the Epikarst Store (from 150 to 180 mm) and of the Karst Store 2 (from 180 to 200 mm) is between the initial and maximum setting storage capacity (Appendix Figure 5). Once the maximum storage capacity is exceeded, water infiltration overflows to the next reservoirs, explaining the flat shape of maximum flow periods (Appendix Figure 5).

To assess the reliability of the model output, we conducted a correlation analysis between the simulated flux and the observed drip rate. This correlation is based on monthly values for all the studied periods from 2013-2019 (Figure 8). Results show that the correlations between Vil-#10A and F4 flow (from soil and epikarst stores), between Vil-#10B and F7 flow (from Karst store 2), and between Vil-#1B and F5 flow (from Karst store 1), are significant (p < 0.01) with a high coefficient of correlation: R = 0.8, R = 0.75 and R=0.61, respectively (Figures 8A& 8B& 8C).

## 4.5 Simulations of isotopes of drip water

The Epikarst δ¹⁸O values display different isotope variability amplitudes when we set different F1 and Epicap (Epikarst capacity) values (Figure 9A). The amplitude is reduced when we set F1=0.1 (mm/month) and Epicap value=160mm, leading to isotopic values close to the one observed at Vil-#10A station (about 1.5 ‰ maximum amplitude) (Figure 9A). When we test another parameter setting in F1 = 0.95/0.5 (mm/month) and Epicap values = 100/50 mm, the modeled Epikarst has a much larger variance (Figure 9A). We evaluated the impact of the evapotranspiration (ETP) on Epikarst δ¹⁸O by examining the proportion of karst water (k_eevap) across varying parameters (0.1, 0.04, 0.01 mm/month). By applying a correction factor (K_d18O_Epi= 0.04) to account for the effects of ETP, we were able to generate a mean value of Epikarst δ¹⁸O (Appendix Figure 6).

Karst Store 2 simulation suggests that water feeding the Vil-#10B station (mean value=-6.14‰, SD=0.09 ‰) has a low variability amplitude (Figure 9B). Modeled data are in agreement with observed δ¹⁸O_d when the low F4 (F4=0.1) and high Ovicap value (Ovicap=180) were set in the model (Figure 9B). Here too, when we also tested another parameter for F4 = 0.95/0.5 (mm/month) and Ovicap values = 100/50 mm, the modeled Karst Store 2 has a much larger variance (Figure 9B). At Karst Store 1, the simulated δ¹⁸O_d values give a similar trend with lower gallery dripping points (Vil-#1A and Vil-#1B stations, mean value=-6.35‰, SD=0.08‰), when the low F5 value (F5=0.1, mm/month)) and high Ovicap value (Ovicap=180 mm) were considered (Figure 9C).
465 **5. Discussion**

466 **5.1. Linking observed drip stations to ERT image and KarstFor**
467 **model**

468 One challenging question is to attribute each dripping station in Villars Cave to a
469 specific reservoir exit or flux conduit in the KarstFor model pattern. We utilized the
470 time series of drip rates, where the variability indicates its more or less connectivity to
471 water excess changes. Results show that Vil-#10A is the only dripping station that
472 exhibits a highly variable drip rate, which can be attributed to its direct connection to
473 the Soil or Epikarst store (Figures 2B & 7). Therefore, we propose to locate the Vil-
474 #10A station at the drainage exit of the Epikarst reservoir (F4 flow; Figures 2B &
475 7C).

476 The upper ERT image, characterized by low resistivity values, suggests
477 reservoirs that can be assimilated to Karst Store 2 (Figures 2A; 2B). According to the
478 drip rate change of Vil-#10B, more damped flows are observed compared to Vil-
479 #10A. For this reason, we suggest putting it at the exit of the Karst Store 2 reservoir
480 (F7; Figures 2B &7C).

481 The deepest part of the ERT profile (below 162m NGF) has low resistivity
482 values (Figure 2B), which may be explained by the presence of a deeper reservoir
483 storage area. In the lower gallery, we observed that Vil-#1B dripping station never
484 stops, even during dry seasons. Therefore, we proposed that this site should be put at
485 the drainage exit of the Karst Store 1 reservoir (F5 flow; Figures 2B and Figure 7). In
486 fact, the lower cave gallery (Vil-#1B) is situated beneath Karst Store 1, receiving its
487 water supply from Karst Store 1 (Figure 2B). Through these attributions, we not only
488 validate the accuracy of the modeled flowing processes, but also show that, even
489 when situated at the same depth and close to each other within the upper gallery, two
490 drip stations may be supplied by different flow paths and reservoir systems (Vil-#10A
491 and Vil-#10B for example).

492 **5.2. Long-term changes of Villars Cave drip water $\delta^{18}$O**

493 At an inter-annual time-scale, we compared modeled and observed data, showing
494 that the KarstFor model's simulations are in coherence with the observations,
495 revealing the $\delta^{18}$O<sub>d</sub> increasing trend between 2006 and 2019 (Figure 10A). The
496 enriched trend of $\delta^{18}$O<sub>d</sub> is parallel to the interannual decline of the water excess
497 (Figure 10C), which can be attributed to the rising trend of potential
498 evapotranspiration (ETP) (Figure 10D). Our observations reveal a consistent increase
499 in temperature within the cave (Figure 10B). Specifically, in the upper gallery, we
500 observed a temperature gradient of 0.1°C/10 years, while in the lower gallery, the
501 temperature gradient was 0.3°C/10 years (Figure 10B). This finding suggests that the
502 interannual trend in $\delta^{18}$O<sub>d</sub> may be attributed to the temperature increase, associated
503 with global warming, which led to the rising of the evapotranspiration (ETP)
(Figures 10B and 10D).

The way by which the ETP controlled the dripping δ¹⁸O is still not clear. It could
theoretically produce a rise in soil isotope fractionation which could affect the mean
seepage water δ¹⁸O. However, it is important to note that this phenomenon will have
a negligible impact on the Villars site because it is continuously covered by vegetation
which prevent evaporation (Figure 10D). But an ETP increase means that less rain
water infiltrates which also potentially can modify the mean annual infiltrated δ¹⁸O by
a more or less complex way. The study of other sites with similar observed δ¹⁸O
trends in the drip water may give more precise explanations.

5.3 Infiltration “time delay” versus isotopic mixing process

The karst hydrological model results suggest that drip water at different
monitoring sites comes from different reservoirs. Vil-#10A, which is connected to
sub-surface reservoirs and is in close proximity to the soil and epikarst (Figures 2A
and 2B), exhibits the highest variability in δ¹⁸O_d values compared to other drip sites
(Figure 3B). Vil-#10B from the upper gallery and Vil-#1B from the lower gallery
(Figure 3C) provide support for the initial hypothesis that these reservoirs are distinct
(Genty et al., 2014; Jean-Baptiste et al., 2019). Because of their low δ¹⁸O_d variability
and their very stable values, it is likely that they are fed by reservoirs with a long
infiltration mixing isotopic process. This hypothesis is further supported by the drip
rate characteristics, including a long-time delay in response to excess water and low
amplitude.

As mentioned before, one of the characteristics of the Villars Cave isotopic
survey is the small (0.2‰), but significant, difference between the upper and lower
gallery δ¹⁸O_d value (Vil-#10B and Vil-#1A/1B; Figure 3C). This difference might be
controlled by reservoir overflows in the infiltration zone as suggested by Genty et al.,
(2014). Lower values of δ¹⁸O_r occur during the high recharge period in winter. When
the reservoir storage exceeds the water level, it triggers an overflow mechanism that
supplies water to the deeper reservoir. This overflow subsequently mixes with the
Karst Store 1 water and provides supply to the drip stations (Vil-#1A and Vil-#1B) in
the lower gallery. Conversely, during the summer, the δ¹⁸O_r is much higher, which
coincides with a low recharge period, leading to the absence of low δ¹⁸O_d overflow

The time delays between the water excess and drip rates are a few days to several
months (Section 4.2). This “time delay” is not the actual time of the water infiltration,
but rather the time necessary to “fill” the karst reservoirs and “push” the water
downwards towards the dripping stations in a “piston flow” system. Actually, the
mean age of the drip water is much older and may reach several years (Genty et al.,
2014; Jean-Baptiste et al., 2019).

To resolve this apparent inconsistency between infiltration “time delay” and
isotopic mixing time, we have to imagine a piston flow system in which the change in
drip rate is due to changes in hydraulic charge related to the filling of the reservoir
that feeds the drip stations. Then, hydraulic pressure, linked to water excess changes,
determines the dripping rate variations. However, the mixing process and overflow
545
determine the isotope values of the water, which is always mixed with "old water"
546
inside the reservoirs. Consequently, we can distinguish the "time delay" that applies
547
to the water excess/drip rate variations, while the "isotopic mixing transit time"
548
concerns the mean age of the water that infiltrates through the reservoirs (Genty et al.,
549
2014).
550

551
5.4. Implications for speleothem records

The karst hydrological model that we used here allows us to simulate and
552
understand the sensitivity of the dripping stalactites feeding stalagmites to the rainfall
553
input signal. Consequently, it could be used to assess whether speleothems from
554
different sampling points in the cave are suitable for reconstructing paleoclimate
555
change. For example, if a stalagmite is fed by a very variable dripping point with a
556
highly seasonal isotopic signal, we should be able to reconstruct past seasonal
557
climatic signals. On the opposite, if the feeding of a stalagmite is well-damped, with
558
no seasonality in the δ¹⁸O_d values, then it will be more adapted to study longer
559
climatic trends.
560

Model simulations show that deeper drip stations fed by reservoir overflows are
561
subject to isotopically biased lower values, which is important for
562
palaeoclimatological studies. Our data from Villars Cave upper and lower galleries
563
confirm this result. The isotopic differences in the drip water, like the observed and
564
modeled one in Villars Cave between lower and upper galleries (−0.2‰), will
565
produce differences in the isotopic composition of contemporaneous speleothems,
566
adding more variabilities to other factors like kinetics fractionations. Indeed,
567
conducting additional tests on actively growing calcite from various levels within
568
Villars Cave will provide fresh insights into this matter, not only within this cave but
569
also across other cave systems. On another note, comprehending the processes and
570
variations in the isotopic composition of dripping water is pivotal for the
571
interpretation of speleothem fluid inclusions. In the scenario of a well-mixed seepage
572
reservoir supplying a dripping station, the absence of seasonal δ¹⁸O variation suggests
573
that fluid inclusions should reflect the average annual (or multi-year) composition of
574
rainfall water.
575

The KarstFor model is a real tool for the understanding of the infiltration
576
processes and for formulating hypotheses about the presence and connectivity of
577
reservoirs in the karst system. For instance, we can investigate the effects of different
578
infiltration patterns into the reservoir, varying seasonality, and frequency of heavy
579
rainfall events. Such simulations allow us to understand the reasons behind the
580
discrepancies in dripping rates and isotopic composition. Furthermore, forecasting
581
with models can extend to predicting stalagmite isotopic composition using an
582
internal-process model like the ISOLUTION model proposed by Deininger and
583
Scholz, 2019.
584
The drip water δ¹⁸O enrichment can be attributed to the global warming over the past 20 years. However, it is also crucial to consider the impact of cave air temperature on the δ¹⁸O value of speleothems. Variations in cave air temperature have the potential to affect the fractionation between drip water and speleothem calcite, thereby influencing the isotopic composition of speleothems. In the future, by incorporating future global warming scenarios, the KarstFor model indicates that oxygen isotopes in drips and stalagmites may exhibit a slight enrichment.

## 6. Conclusions

Based on extensive hydrological and isotopic monitoring of rainfall and drip water within Villars Cave, we conducted a comprehensive study using a hydrological-isotopic conceptual model to analyze rainfall infiltration behavior in a karst environment. Villars Cave served as an ideal case study due to its multiple dripping stations, which have been monitored for over 20 years. These stations are located at varying depths within the cave and exhibit distinct hydrological characteristics driven by diverse pathways and reservoirs. The geometrical information obtained from the first direct images of the reservoirs, acquired through the ERT measurement above the cave, allowed us to associate each dripping station with a specific location in the conceptual model (KarstFor). After having fixed the parameters used in the karst hydrological model (Bradley et al., 2010; Baker et al., 2013; Treble et al., 2013; 2019), simulation outputs lead to the following conclusions:

1) the seasonality of rainfall infiltration, as characterized by water excess (or efficient rainfall), is effectively simulated in the epikarst zones and within various reservoir overflows. The closest dripping station, Vil-#10A, located in the upper gallery, displays the highest isotopic variability, underscoring its connection to these subsurface reservoirs.

2) the “apparent” time lag between water excess and various dripping stations, calculated using cross-correlation methods, varies from a few days to several months. The extent of these time lags depends on the station's depth and its inertia linked to reservoir behaviors. The KarstFor model accurately simulates this “apparent” time delay as a consequence of recharge that push infiltrated water into a “piston flow” system. However, it's worth noting that the actual water transit time is considerably longer, spanning up to several years, as demonstrated by prior studies.

3) the small yet consistent isotopic difference between the upper and lower gallery drip waters suggests that the -0.2 ‰ lower isotopic values in the deeper gallery are likely attributable to overflows, which are more likely to happen during the winter season when the δ¹⁸O values are lower. In this regard, model simulations and
625 observations align coherently.

626 4) all the drip stations from the Villars Cave display a δ¹⁸O₈ long-term
627 (decadal) increase, possibly due to higher global external
628 temperature leading to more ETP and less water excess, but more
629 studies are necessary.

630 The combination of a long drip rate and isotope monitoring, and modeling linked
631 to an actual pattern of reservoirs seen by ERT measurements, is a novel and
632 promising approach to studying the origin of differences observed in the flows and the
633 isotopic compositions of dripping sites in caves. Indeed, this approach appears useful
634 to estimate the consequences of any modifications of the climatic conditions above a
635 cave, especially changes in rainfall intensities, seasonality, and temperatures, due to
636 climatic changes.

637 Acknowledgments

638 We thank the Versaveaud family, owner of the cave site, who are always
639 interested in our research and whose constant help makes the long-term monitoring
640 work in Villars Cave possible. We would like to express our sincere gratitude to the
641 editor, Christian France-Lanord and three anonymous reviewers for the invaluable
642 assistance in enhancing the quality of this article. We want to express sincere thanks
643 to Olivier Mestre for the help with meteorological data, Thierry Baritaud for his help
644 in the cave exploration and mapping We also would like to express sincere thanks to
645 Manon Sabourdy for the help of ArcGIS software. We thank Jean Pierre CASSOU for
646 his helps in constructing the new map of the Villars Cave (Figure 1B and 1C), and
647 Daniel Chailloux for the measurement of the different dripping station depth with a
648 magnetic probe. This research was supported by the China Scholarship Council (CSC)
649 to Jian. Zhang (201906990014). This long monitoring work (DG) has been funded,
650 since 1996, by several programs from the CNRS (INSU), and three laboratories
651 (GEOPS, Orsay, Paris-Saclay University; Laboratoire des Sciences du Climat et de
652 l'Environnement, LSCE, CEA/CNRS/UVSQ and Environnements et
653 Paléoenvironnements Océaniques et Continentaux, EPOC, University of Bordeaux,
654 France). Le Mas and Villars stations are now part of the new French rainfall isotopic
655 network RENOIR SOERE Observatory of the INSU/CNRS. I2M programs on
656 decorated caves for geophysical surveys also supported this work. Part of this work is
657 linked to the new ANR DECACLIM research program (coord. B. Lartiges, GET,
658 Toulouse) and constitutes a contribution to the CNRS scientific group (GIS) GEMS
659 (Groupe d'Etude des milieux Souterrains).

660 CRediT authorship contribution statement

661 **Jian Zhang**: Conceptualization, Methodology, Formal analysis, Writing - review & editing. **Dominique Genty**: Writing - review & editing, Supervision, Funding acquisition, Monitoring data acquisitions, Water sampling. **Colette Sirieix**: Revising, Funding acquisition and Supervision. **Bénédicte Minster**:
water isotopic analysis. **Edouard Régnier:** Monitoring data acquisitions and Water
sampling. **Giovani Sgubin:** Data and revising. **Ting-Yong Li :** Revising and
discussion, **François Bourges:** Methodology and discussion, **Ming-Qiang Liang:**
Software, **Chaojun Chen:** Software. **Cécile Verdet:** Software and cave mapping,
**Sylvain Mateo, Ludovic Devaux, Stéphane Bujan, Vincent Marieu:** Field works.

670 Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported in this
paper.
671
672
673

674 References

675 Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. 1998. Crop evapotranspiration-Guidelines for
676 computing crop water requirements-FAO Irrigation and drainage paper 56. Fao, Rome, 300(9),
677 D05109.

678 Almalki, S.J. and Nadarajah, S., 2014. Modifications of the Weibull distribution: A review.
679 Reliab. Eng. Syst. Saf., 124, pp.32-55.

680
681
682
Ayalon, A., M. Bar-Matthews, and E. Sass. 1998., Rainfall-recharge relationships within a karstic terrain in the Eastern Mediterranean semi-arid region, Israel: δ¹⁸O and δD characteristics, J.Hydrol, 207(1), 18-31.

683
684
685
Baker, A., C. Bradley, S. Phipps, M. Fischer, I. J. Fairchild, L. Fuller, C. Spötl, and C. Azcurra
.2012., Millennial-scale forward models and pseudoproxies of stalagmite δ¹⁸O: An example
from NW Scotland, Clim. Past, 8, 1153–1167, doi:10.5194/cp-8-1153-2012.

686
687
688
Baker, A. Berthelin, R., Cuthbert, M.O., Treble, P.C., Hartmann, A. and the KSS Cave Studies
Team. Rainfall recharge thresholds in a subtropical climate determined using a regional cave
drip water monitoring network. J. Hydrol, 2020: 125001.

689
690
691
Baker, A., Scheller, M., Oriani, F., Mariethoz, G., Hartmann, A., Wang, Z., & Cuthbert, M. O.
2021. Quantifying temporal variability and spatial heterogeneity in rainfall recharge thresholds
in a montane karst environment. J. Hydrol, 594, 125965.

692
693
694
Baker, A., A. Asrat, I. J. Fairchild, M. J. Leng, L. Thomas, M. Widmann, C. N. Jex, Buwen Dong,
P. van Calsteren, and C. Bryant .2010., Decadal-scale rainfall variability in Ethiopia recorded
in an annually laminated, Holocene-age, stalagmite, The Holocene, 20(6), 827-836.

695
696
Baker, A., C. Bradley, and S. J. Phipps.2013., Hydrological modeling of stalagmite δ¹⁸O response
to glacial-interglacial transitions, Geophys. Res. Lett., 40(12), 3207-3212.

697
698
699
Baker, A., Hartmann, A., Duan, W., Hankin, S., Comas-Bru, L., Cuthbert, M.O., Treble, P.C.,
Banner, J., Genty, D., Baldini, L., Bartolomé, M., Moreno, A., and Pérez-Mejías, C., 2019.
Global distribution and controls on cave drip water oxygen isotope composition. Nat.
Commun., 10, Article number: 2984.
701

702

703

Baker, A., Mockler, N.J., Barnes, W.L., 1999. Fluorescence intensity variations of speleothem-forming groundwaters: implications for paleoclimate reconstruction. Water Resour. Res. 35, 407-413.

704

705

706

707

Baldini, J. U. L., F. McDermott, L. M. Baldini, C. J. Ottley, K. L. Linge, N. Clipson, and K. E. Jarvis. 2012.. Identifying short-term and seasonal trends in cave dripwater trace element concentrations based on a daily-scale automatically collected dripwater dataset, Chem. Geol., 330-331, 1-16.

708

709

710

711

Bar-Matthews, M., Ayalon, A., Matthews, A., Sass, E., & Halicz, L. 1996. Carbon and oxygen isotope study of the active water-carbonate system in a karstic Mediterranean cave: Implications for paleoclimate research in semiarid regions. Geochim. Cosmochim. Acta, 60(2), 337-347.

712

713

714

Bourges, F., Genthon, P., Genty, D., Lorblanchet, M., Mauduit, E., d'Hulst, D. 2014. Conservation of prehistoric caves and stability of their inner climate: lessons from Chauvet and other French caves. Sci. Total Environ., 493, 79-91.

715

716

717

Bourges, F., Genty, D., Perrier, F., Lartiges, B., Régnier, É., François, A., et al. 2020. Hydrogeological control on carbon dioxide input into the atmosphere of the Chauvet-Pont d'Arc cave. Sci. Total Environ., 716, 136844.

718

719

720

Bradley, C., A. Baker, C. N. Jex, and M. J. Leng. 2010.. Hydrological uncertainties in the modelling of cave drip-water δ¹⁸O and the implications for stalagmite palaeoclimate reconstructions, Quat. Sci. Rev., 29(17-18), 2201-2214.

721

722

723

Breitenbach, S. F., Lechleitner, F. A., Meyer, H., Diengdoh, G., Mattey, D., & Marwan, N., 2015. Cave ventilation and rainfall signals in dripwater in a monsoonal setting-a monitoring study from NE India. Chem. Geol., 402, 111-124.

724

725

726

Crawford, J., Hughes, C. E., Lykoudis, S., 2014. Alternative least squares methods for determining the meteoric water line, demonstrated using GNIP data. J.Hydrol., 519, 2331-2340.

727

728

729

Cobb, K.M., Adkins, J.F., Partin, J.W., Clark, B., 2007. Regional-scale climate influences on temporal variation of rainwater and cave dripwater oxygen isotopes in northern Borneo. Earth Planet. Sci. Lett. 263, 207-220

730

731

732

Coleborn, K., Rau, G.C, Cuthbert, M.O., Baker, A and Navarre, O., 2016. Solar forced diurnal regulation of cave drip rates via phreatophyte evapotranspiration. Hydrol. Earth. Syst. Sci., 20, 4439-4455

733

734

735

Cruz, F.W., Karmann, I., Viana, O., Burns, S.J., Ferrari, J.A., Vuille, M., Sial, A.N., Moreira, M.Z., 2005. Stable isotope study of cave percolation waters in subtropical Brazil: implications for paleoclimate inferences from speleothems. Chem. Geol. 220, 245-262

736

737

Cuthbert, M.O., Baker, A., Jex, C.N., Graham, P.W., Treble, P., Andersen, M.S and Acworth, R.I., 2014. Drip water isotopes in semi-arid karst: implications for speleothem paleoclimatology.
774 37, 173–191.

775 Genty, D., Labuhn, I., Hoffmann, G., Danis, P.A., Mestre, O., Bourges, F., Wainer, K., Massault, M., Van Exter, S., Régnier, E., Orengo, Ph., Falourd, S., Minster, B., 2014. Rainfall and cave water isotopic relationships in two South-France sites. Geochim. Cosmochim. Acta. 131 (5), 323–343.

779 Hartmann A., Lange J., Weiler M., Arbel Y., & Greenbaum N., 2012. A new approach to model the spatial and temporal variability of recharge to karst aquifers. Hydrol. Earth. Syst Sci., 16(7), 2219–2231. 10.5194/hess-16-221.

782 Hartmann, A., Baker, A., 2017. Modelling karst vadose zone hydrology and its relevance for paleoclimate reconstruction. Earth-Sci. Rev. 172, 178–192. https://doi.org/10.1016/j.earscirev.2017.08.001.

785 Healy, Richard W. Estimating groundwater recharge. Cambridge University Press, 2010.

786 Hendy, C. H. 1971, The isotopic geochemistry of speleothems—I. The calculation of the effects of different modes of formation on the isotopic composition of speleothems and their applicability as palaeoclimatic indicators, Geochim. Cosmochim. Acta, 35(8), 801–824.

789 Hu, C., Henderson, G. M., Huang, J., Chen, Z., and Johnson, K. R.: Report of a three-year monitoring programme at Heshang Cave, Central China, Int. J. Speleol., 37, 143–151, doi:10.5038/1827-806X.37.3.1, 2008.

792 Hu, J., Dee, S. G., Wong, C. I., Harman, C. J., Banner, J. L., & Bunnell, K. E. 2021., Assessing proxy system models of cave dripwater $\delta^{18}$O variability. Quat. Sci. Rev, 254, 106799.

794 Jean-Baptiste, P., Genty, D., Fourre, E., Regnier, E., 2019. Tritium dating of dripwater from Villars cave (SW-France). Appl. Geochem. 107, 152–158.

796 Kaufman A., Bar-Matthews M., Ayalon A. and Carmi I. 2003. The vadose flow above Soreq Cave, Israel: a tritium study of the cave waters. J. Hydrol. 273, 155–163.

798 Kluge T., Riechelmann D. F. C., Wieser M., Spotl C., Sultenfuss J., Schroder-Ritzrau A., Niggemann S. and Aeschbach-Hertig W. 2010., Dating cave drip water by tritium. J. Hydrol. 394, 396–406.

801 Kluge, T., T. Marx, D. Scholz, S. Niggemann, A. Mangini, and W. Aeschbach-Hertig. 2008, A new tool for palaeoclimate reconstruction: Noble gas temperatures from fluid inclusions in speleothems, Earth Planet. Sci. Lett., 269(3–4), 408–415.

804 Lachniet, M. S. 2009. Climatic and environmental controls on speleothem oxygen-isotope values. Quatern. Sci. Rev, 28(5-6), 412-432.

806 Leopold, M., Gupanis-Broadway, C., Baker, A., Hankin, S., Treble, P. 2021. Time lapses electric resistivity tomography to portray infiltration and hydrologic flow paths from surface to cave. J. Hydrol, 593, 125810.
809

810

811

812

Li, T.-Y., Shen, C.-C., Li, H.-C., Li, J.-Y., Chiang, H.-W., Song, S.-R., Yuan, D.-X., Lin, C.D.J., Gao, P., Zhou, L.-P., Wang, J.-L., Ye, M.-Y., Tang, L.-L., Xie, S.-Y., 2011. Oxygen and carbon isotopic systematics of aragonite speleothems and water in Furong Cave, Chongqing, China. Geochim. Cosmochim. Acta. 75(15), 4140-4156.

813

814

Loke, M., 2004. Tutorial: 2-D and 3-D Electrical Imaging Surveys Instructions for Use of RS2DINV.

815

816

Loke, M.H., Dahlin, T., Rucker, D.F., 2014. Smoothness-constrained time-lapse inversion of data from 3D resistivity surveys. Near Surf. Geophys. 12 (1):5-24.

817

818

819

Mahmud, K., G. Mariethoz, A. Baker, P. C. Treble, M. Markowska, and E. McGuire 2016, Estimation of deep infiltration in unsaturated limestone environments using cave lidar and drip count data., Hydrol. Earth Syst. Sci., 20(1), 359-373.

820

821

822

Markowska, M., Andersen, M.S., Treble, P.C., Baker, A., Tadros, C., Hankin, S., and Jex, C.N., 2015. Unsaturated zone hydrology and cave drip discharge water response: Implications for speleothem palaeoclimate record variability. J. Hydrol, 529, 662-675.

823

824

825

Markowska, M., Cuthbert, M.O., Baker, A., Treble, P.C., Andersen, M.S., Adler, L., Griffiths, A., and Frisia, S., 2020. Modern speleothem oxygen isotope hydroclimate records in water-limited SE Australia. Geochim. Cosmochim. Acta, 270, 431-448.

826

827

828

Mattey, D., Lowry, D., Duffet, J., Fisher, R., Hodge, E., Frisia, S., 2008. A 53 year seasonally resolved oxygen and carbon isotope record from a modern Gibraltar speleothem: reconstructed drip water and relationship to local precipitation. Earth Planet. Sci. Lett. 269, 80-95

829

830

831

832

833

Mattey, D.P., Fairchild, I.J., Atkinson, T.C., Latin, J.-P., Ainsworth, M., Durell, R., 2010. Seasonal microclimate control of calcite fabrics, stable isotopes and trace elements in modern speleothem from St Michaels Cave, Gibraltar. In: Pedley, H.M., Rogerson, M. (Eds.), Tufas and Speleothems: Unravelling the Microbial and Physical Controls. The Geological Society of London, London, pp. 323-344.

834

835

836

McDonald, J. 2004.: The 2002-2003 El Niño recorded in Australian cave drip waters: Implications for reconstructing rainfall histories using stalagmites, Geophys. Res. Lett., 31, L22202.

837

838

Michelle, S. 2002. Numerical Analysis: A Mathematical Introduction, Chapters 4 and 6. Clarendon Press, Oxford.

839

840

Oster, Jessica L., Isabel P. Montañez, and Neil P. Kelley. 2012. Response of a modern cave system to large seasonal precipitation variability. Geochim. Cosmochim. Acta 91 92-108.

841

842

843

Pape, J. R., J. L. Banner, L. E. Mack, M. Musgrove, and A. Guilfoyle. 2010. Controls on oxygen isotope variability in precipitation and cave drip-waters, central Texas, USA, J. Hydrol, 385(1-4), 203-215.

844

Proctor, C. J., Baker, A., Barnes, W. L., and Gilmour, M. A. 2000. A thousand-year speleothem
proxy record of North Atlantic climate from Scotland, Clim. Dyna., 16, 815-820.

Rau, G.C., Cuthbert, M.O., Andersen, M.S., Baker, A., Rutlidge, H., Markowska, M., Roshan, H., Marjo, C.E., Graham, P.W. and Acworth, R.I. 2015. Controls on cave drip water temperature and implications for speleothem-based paleoclimate reconstructions. Quat. Sci. Rev., 127, pp.19-36.

Riechelmann, D. F. C., Schröder-Ritzrau, A., Scholz, D., Fohlmeister, J., Spötl, C., Richter, D. K., & Mangini, A. 2011. Monitoring Bunker Cave (NW Germany): A prerequisite to interpret geochemical proxy data of speleothems from this site. J. Hydrol, 409(3-4), 682-695.

Scholz, D., C. Mühlinghaus, and A. Mangini. 2009., Modelling δ¹³C and δ¹⁸O in the solution layer on stalagmite surfaces, Geochim. Cosmochim. Acta, 73(9), 2592-2602.

Soil Guidelines-soil water storage capacity and available soil moisturehttps://www.droughtmanagement.info/literature/BC_MA_Soil_Water_Storage_Capacity_2005.pdf

Spötl, C., I. J. Fairchild, and A. F. Tooth.2005., Cave air control on dripwater geochemistry, Obir Caves (Austria): Implications for speleothem deposition in dynamically ventilated caves, Geochim. Cosmochim. Acta, 69(10), 2451-2468.

Stoll, H. M., A. Mendez-Vicente, S. Gonzalez-Lemos, A. Moreno, I. Cacho, H. Cheng, and R. L. Edwards.2015., Interpretation of orbital scale variability in mid-latitude speleothem δ¹⁸O: Significance of growth rate controlled kinetic fractionation effects. Quat. Sci. Rev. 127, 215-228.

Tadros, Carol V., Pauline C. Treble, Andy Baker, Ian Fairchild, Stuart Hankin, Regina Roach, Monika Markowska, and Janece McDonald. 2016. ENSO-cave drip water hydrochemical relationship: a 7-year dataset from south-eastern Australia. Hydro. Earth. Syst Sci., 20, no. 11: 4625-4640.

Tooth, A. F., & Fairchild, I. J. 2003. Soil and karst aquifer hydrological controls on the geochemical evolution of speleothem-forming drip waters, Crag Cave, southwest Ireland. J. Hydrol, 273(1-4), 51-68.

Treble, P. C., C. Bradley, A. Wood, A. Baker, C. N. Jex, I. J. Fairchild, M. K. Gagan, J. Cowley, and C. Azcurra. 2013., An isotopic and modelling study of flow paths and storage in Quaternary calcarenite, SW Australia: implications for speleothem paleoclimate records, Quat. Sci. Rev, 64, 90-103.

Treble, P., Mah, M., Griffiths, A., Baker, A., Deininger, M., Kelly, B., Scholz, D. and Hankin, S., 2019. Separating isotopic impacts of karst and in-cave processes from climate variability using an integrated speleothem isotope-enabled forward model. 2019. EarthArXiv, 10.31223/osf.io/j4kn6c.

Treble, P.C., Baker, A., Abram, N.J., Hellstrom, J.C., Crawford, J., Gagan, M.K., Borsato, A., Griffiths, A.D., Bajo, P., Markowska, M., Priestley, S.C., Hankin, S., Paterson, D., 2022.
Ubiquitous karst hydrological control on speleothem oxygen isotope variability in a global study. Commun Earth Environ 3 (1), 29.

Tremaine, D. M., & Froelich, P. N. 2013. Speleothem trace element signatures: A hydrologic geochemical study of modern cave dripwaters and farmed calcite Geochim. Cosmochim. Acta, 121, 522-545.

Wong, C. I., and D. O. Breecker.2015, Advancements in the use of speleothems as climate archives, Quatern. Sci. Rev, 127, 1-18.494

Yonge, C.J., Ford, D.C., Gray, J., Schwarcz, H.P., 1985. Stable isotope studies of cave seepage water. Chem. Geol. 58, 97-105

Xu, S., Sirieix, C., Ferrer, C., Lacanette-Puyo, D., Riss, J., Malaurent, P., 2015. A geophysical tool for the conservation of a decorated cave. - a case study for the Lascaux cave. Archaeol. Prospect. http://dx.doi.org/10.1002/arp.1513 (online in Wiley OnlineLibrary)

Xu, S., Sirieix, C., Riss, J., Malaurent, P., 2017. A clustering approach applied to time-lapse ERT interpretation case study of Lascaux cave. J. Appl. Geophys. 144, 115-124

Xue, Y. X., Zhang, J., Su, Z., Wu, Y., Liang, Q. S., Liang, M. Q., Liu, X-Z., Chen, C-J, Qiu, H-Y Li, J-Y., Li, T. Y., 2023. Quantifying source effects based on rainwater $\delta^{18}$O from 10-year monitoring records in Southwest China. Appl. Geochem. 105706.

Zhang, J., Li. Ting-Yong., 2019. Seasonal and inter-annual variation of hydro-chemical characteristics and stable isotopic compositions of drip waters in Furong Cave, Southwest China based on 12 years' monitoring. J. Hydrol. 572:40-50.

Zhang, J., Genty, D., Sirieix, C., Michel, S., Minster, B., & Régnier, E. 2020. Quantitative assessments of moisture sources and temperature governing rainfall $\delta^{18}$O from 20 years' monitoring records in SW-France: Importance for isotopic-based climate reconstructions. J. Hydrol., 591, 125327.

Zhang, J., Liang, M., Li, T., Chen, C. & Li, J. 2022. Asian-Australian monsoon evolution over the last millennium linked to ENSO in composite stalagmite $\delta^{18}$O records. Quat. Sci. Rev. 281, 107420.

Zhang, J., Sirieix, C., Genty, D., Verdet, C., Salmon, F., Xu, S., Mateo, S., Larcanché, M., Bujan, S., Devaux, L. 2023 Imaging hydrological dynamics in karst unsaturated zones by time-lapse electrical resistivity tomography. Sci. Total Environ. (accepted)

Zhang. J., Characterization of the infiltration in the cave of Villars (Dordogne) by geophysical and isotopic methods: Importance for paleo-climatic studies based on speleothems as well as for the conservation of heritage sites. 26th. June. 2023., Ph.D.

Caption
**Figure 1 (A)** Location of Villars Cave. **(B)** (B) Three-dimensional diagram depicting the cave structure and the distribution of surveillance sites (marked as red circles) within Villars Cave (Vil-#1A, Vil-#1B, Vil-#10A, and Vil-#10B) is shown (J.P. Cassou, 2020, unpublished). For detailed site descriptions and maps, please refer to (Genty, 2008; Genty et al., 2014). **(C)** Profile of ERT (Red line) above the Villars Cave map (made by Jean-Pierre Cassou, 2020) and from SCP (Speleo Club de Périgueux, Th. Baritaud and P. Vidal). The red line represents the ERT profile that was monitored. In Figure 2C, the dotted circles labeled as Zone A, Zone B, and Zone C correspond to the ERT image identifying three cave galleries in Figure 2A. The different colors in Panel C indicate the cave galleries.

**Figure 2.** ERT image revealing karst reservoirs. This image helps to set four storage reservoirs (Soil store, Epikarst store, Karst store 1, and Karst store 2) which entered the isotope model. **(A)** Resistivity model of the 07th of January 2021 with an electrode spacing of 1 m. The inversion was performed based on the concatenation of the data of all the two arrays (Gradient Array and Pole-Dipole Array), using Res2Dinv software to process data. **(B)** Conceptual model of karst hydrology in individual karst water reservoirs. The boxes represent different reservoirs, blue lines represent water levels, blue arrows represent fracture flow, and the red arrow represents diffuse flow (Bradley et al., 2010, Baker et al., 2013). A brief overview of the overall model is given in the main text, and the details of the KarstFor model are given in (Bradley et al., 2010).

**Figure 3.** Combined time series of meteoric conditions and isotope hydrology of waters in Villars Cave. **(A)** Comparison with rainfall and drip water δ¹⁸O values of sites (Vil#1A, Vil#1B, Vil#10A and Vil#10B) in Villars Cave. **(B)** Detailed figure of the δ¹⁸O values for Vil-#10A. **(C)** Detailed figure of the δ¹⁸O values for Vil-#1A, Vil-#1B, and Vil-#10B. **(D)** Mnthly average air temperature outside of Villars Cave. **(E)** Water excess outside Villars Cave. The pink shading represents the winter months.

**Figure 4.** Local meteoric water line (LMWL) (red solid line) based on the δ¹⁸O and δD values of rainfall outside Villars Cave. **(A)** δ¹⁸O and δD values of drip water Vil-#10A station; note its distribution around the LMWL. **(B)** δ¹⁸O and δD values of drip water Vil-#1B, Vil-#10A and Vil-#10B stations in Villars Cave. Note their damped amplitude signal compared to that of precipitation. Global Meteoric Water Line (GMWL) in black (Craig, 1961) is plotted for comparison.

**Figure 5.** Comparison of the variability of daily R-ETP and drip rate during 2013-2019. R-ETP>0 is called water excess. A smoothed (10 points) average of the R-ETP appears in blue.
**Figure 6.** Time delay between drip rates and water excess during the 2014-2018 period. Pearson's linear correlation coefficients are from MATLAB coherence package. (A-E) Vil-#1A & Water-excess, (F-J) Vil-#1B & Water-excess, (L-M) Vil-#10B & Water-excess. For the R-ETP parameter, the date of minimum values was chosen. The delay time (X is the highest correlation value) is shown in red. We performed a five-day smoothing on water excess and drip rates to remove small peak effects.

**Figure 7.** Hydrological model outputs for the 2006-2019 period. (A) Comparison of local precipitation and ETP in Villars Cave. (B) Water excess (Precipitation minus potential evapotranspiration (R-ETP)>0). (C) Modelled discharge (F1, F4, F7) in the unsaturated zone of Villars Cave. (D) Modelled discharge (F5) in the unsaturated zone of Villars Cave. (E) Modelled discharge (F3, F6) in the unsaturated zone of Villars Cave. Hydrological output is implemented as a python module, see examples at : https://github.com/agriff86/karstolution. We set the parameters in F1=0.95, F4=0.6, F5=0.5, F6=0.001, F7=0.6.

**Figure 8.** Correlation and comparison between the simulated fluxes and the observed drip rates during the 2013-2019 period. (A) Vil-#10A& F4 flow, (B) Vil-#1B & F5 flow, (C) Vil-#10B& F7 flow. Correlation is calculated on a monthly time scale.

**Figure 9.** Optimization and sensitivity tests using different parameters in path processes (F) and overflow (Epicap and Ovicap) based on KarstFor model for Villars Cave. (A) Comparisons of modeled epikarst δ¹⁸O values and observed Vil-#10A (upper gallery) in Villars Cave based on adjusting F1 and overflow (Epicap) in Epikarst. (B) Comparisons of modeled Karst store 2 δ¹⁸O value and observed Vil-#10B (upper gallery) by adjusting F4 and overflow (Ovicap) in karst store 2 in Villars Cave. (C) Comparisons of modeled Karst store 1 and observed Vil-#1A and #1B δ¹⁸O values by adjusting F5 and overflow (Ovicap) in karst store 2 in Villars Cave (lower gallery). There is a 1 °C temperature difference between the upper gallery (12.5°C) and the lower gallery (11.5°C). To account for this, we assigned different evaporation parameters (k_d18O_epi) in the upper galleries (Vil-#1B) with a higher value of 0.06‰, and in the lower galleries (Vil-#1A and #1B) with a lower value of 0.04‰. We did not consider the influence of humidity as it remained above 90% saturation throughout.

**Figure 10.** (A) Modeled δ¹⁸O values in the Karst Store 1 & 2 and observed Vil-#10B (Upper gallery) and Vil-#1A/#1B δ¹⁸O (lower gallery) in Villars Cave. (B) Inside cave temperature (Vil-#1B and Vil-#10B), (C) Mean monthly and annual outside water excess calculated on a monthly time scale. (D) Mean annual outside ETP calculated on a monthly time scale. The red-fitted linear regression indicates its inter-annual change trend.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig 1</figcaption></figure>

1001

1002

1003
<figure><img src="image_4.png" /><figcaption>Fig 2</figcaption></figure>

1004

1005

1006

Date: January 07th, 2021

(A)

(B)
<figure><img src="image_5.png" /><figcaption>Fig 3</figcaption></figure>

1007

1008

1009
<figure><img src="image_6.png" /><figcaption>Fig 4</figcaption></figure>

1010

1011

1012
<figure><img src="image_7.png" /><figcaption>Fig 5</figcaption></figure>

1013

1014

1015

1016
<figure><img src="image_8.png" /><figcaption>Fig 6</figcaption></figure>

1017

1018

1019

1020
<figure><img src="image_9.png" /><figcaption>Fig 7</figcaption></figure>

1021

1022

1023

1024
<figure><img src="image_10.png" /><figcaption>Fig 8</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig 9</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG 10</figcaption></figure>
<table><caption>**Table 1.** Compilation of selected model parameter values and notes about its implementation and configuration. All parameter names are from conceptual Figure 2B.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Final configuration*</td></tr></thead><tbody><tr><td>Soil store (mm)</td><td>Set to a small configuration, no soil evaporation. To determine the storage volume, the available water storage capacity was also taken into account. Soil thickness range from approximately 20 to 40 cm. According to the soil handbook, the clay component of the soil has a water storage capacity of about 21% of its volume (Soil guidelines, 2005). From these values, a selection of 50 mm was made as the initial soil storage volume.</td><td>Init: 50<br>Max: 150</td></tr><tr><td>Epikarst (mm)</td><td>Set to a small configuration, no epikarst evaporation</td><td>Init: 100<br>Max: 180<br>Epicap: 50/100/160</td></tr><tr><td>Karst store 1 (mm)</td><td>The main karst storage that receives flux from epikarsts by k_diffuse and/or F3; may also receive the overflow from Karst store 2.</td><td>Init: 200<br>Max: 200</td></tr><tr><td>Karst store 2 (mm)</td><td>Secondary karst storage, filled through F4, through F6 or allows modeling the switchable overflow of KS1 via F7. The amount of overflow is proportional to the ratio of KS2 / KS1.</td><td>Init: 200<br>Max: 200<br>Ovicap: 50/100/180</td></tr><tr><td>F1 (mm/month)</td><td>Determines steady state values in soil store.</td><td>0.1/0.5/0.95</td></tr><tr><td>F3 (mm/month)</td><td>The flux represents the fracture flow from Epikarst to KS1. The setting is based on field observation, and the flow is zero diffusion flow.</td><td>0.0008</td></tr><tr><td>F4 (mm/month)</td><td>Flux from Epikarst to KS2 is activated when threshold 'Epicap' is reached.</td><td>0.1/0.5/0.6/0.9</td></tr><tr><td>F5 (mm/month)</td><td>Drainage flux of KS1.</td><td>0.01/0.1/0.2/0.5</td></tr><tr><td>F7 (mm/month)</td><td>Overflow from KS2 back into KS1 once 'ovicap' is exceeded.</td><td>0.4</td></tr><tr><td>F6 (mm/month)</td><td>Drainage flux of KS2.</td><td>0.001</td></tr></tbody></table>
<table><caption>**Table 2.** Average and ranges for the δ¹⁸O/δD values for rainfall and drip water (Vil-#10A, Vil-#10B, Vil-#1A, Vil-#1B, V-SMOW) in Villars Cave.</caption><tr><td>F8 (mm/month)</td><td>Bypass from the surface to KS2. Used to test the configuration used by Treble et al. (2013), rainfall events> 7mm/month give priority to charging KS2, but it is set to zero in the final configuration.</td><td>0.0001</td></tr><tr><td>k_diffuse (mm/month)</td><td>The flux is to simulate the diffusion flow through the pdf function.</td><td>0.001</td></tr><tr><td>k_eevap</td><td>Fraction of water remaining in epikarst reservoir available to evaporate</td><td>0.016</td></tr><tr><td>k_d18O_soil</td><td>Isotopic evaporation coefficient for soil reservoir</td><td>0.06 ‰ month<sup>-1</sup> mm<sup>-1</sup></td></tr><tr><td>k_d18O_epi</td><td>Isotope evaporation coefficient of epikarst reservoir</td><td>0.01/0.04/0.06/0.1 ‰ month<sup>-1</sup> mm<sup>-1</sup></td></tr></table>

<table><thead><tr><th></th><th></th><th>Rainfall</th><th>Vil-#10A</th><th>Vil-#10B</th><th>Vil-#1A</th><th>Vil-#1B</th></tr></thead><tbody><tr><td rowspan="4">δ¹⁸O</td><td>Average</td><td>-6.13</td><td>-6.19</td><td>-6.14</td><td>-6.35</td><td>-6.36</td></tr><tr><td>SD</td><td>1.33</td><td>0.21</td><td>0.09</td><td>0.08</td><td>0.08</td></tr><tr><td>Max</td><td>-3.30</td><td>-5.37</td><td>-5.97</td><td>-6.00</td><td>-6.06</td></tr><tr><td>min</td><td>-9.93</td><td>-6.76</td><td>-6.48</td><td>-6.73</td><td>-6.62</td></tr><tr><td rowspan="4">δD</td><td>Average</td><td>-37.99</td><td>-37.43</td><td>-37.03</td><td>-39.24</td><td>-39.44</td></tr><tr><td>SD</td><td>9.37</td><td>1.68</td><td>0.97</td><td>1.15</td><td>0.76</td></tr><tr><td>Max</td><td>-19.60</td><td>-30.81</td><td>-34.50</td><td>-33.70</td><td>-36.20</td></tr><tr><td>min</td><td>-67.40</td><td>-42.50</td><td>-39.90</td><td>-42.80</td><td>-41.70</td></tr></tbody></table>