# Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity

Thomas Jacob, Jean Chery, Roger Bayer, Nicolas Lemoigne, J.-P. Boy, Philippe Vernant, Federick Boudin

► To cite this version:

Thomas Jacob, Jean Chery, Roger Bayer, Nicolas Lemoigne, J.-P. Boy, et al. Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity. Geophysical Journal International, 2009, 177 (2), pp.347-360. 10.1111/j.1365-246X.2009.04118.x. hal-00413528

HAL Id: hal-00413528

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

Submitted on 10 Jun 2021

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 CC BY 4.0 - Attribution - International License
# Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity

Thomas Jacob,<fnref n="1" /> Jean Chery,<fnref n="1" /> Roger Bayer,<fnref n="1" /> Nicolas Le Moigne,<fnref n="1" /> Jean-Paul Boy,<fnref n="2" /> Philippe Vernant<fnref n="1" /> and Frédéric Boudin<fnref n="1" />

Accepted 2009 January 13. Received 2008 November 26; in original form 2008 August 22

## SUMMARY

In this study we attempt to understand the water storage variations in a karst aquifer on the Larzac Plateau (South of France) using ground-based gravimetry. Surface to 60 m depth gravity measurements are performed three times a year since 2006 down a pothole, in complement to monthly absolute gravity (AG) measurements at three sites. The time variations of the surface to depth (STD) gravity differences are compared with the AG variations. Using a simple Bouguer plate model, we find that the STD gravity differences are very similar to the AG variations. The STD gravity differences are then used to determine apparent density values. These integrative density values are compared with measured grain densities from core samples to obtain an apparent porosity and saturation change representative of the investigated depth. The apparent porosity ranges from 4.8 to 7.3 per cent. We then discuss on the repartition of the apparent physical properties with respect to the epikarst and infiltration zone karst structures. We argue that AG and STD differences monitor epikarst water storage variations. Within this scope, we discuss the fact that seasonal scale water storage variation occurs predominantly in the epikarst.

**Key words:** Time variable gravity; Hydrogeophysics; Permeability and porosity.

## 1 INTRODUCTION

The near-surface water content and its variation on continents in- duce both surface deformation and gravity changes. Therefore, hy- drological effects, both natural and anthropogenic, can be detected by precise geodetic techniques such as space radar interferometry (Amelung et al. 1999), GPS (Burbey 2003; King et al. 2007), tilt- meters (Dal Moro & Zadro 1998; Kümpel et al. 2001) and gravime- try (Naujoks et al. 2007). Because modern gravimetric tools reach µGal (10⁻⁸ m s⁻<fnref n="2" />) repeatability, they allow detecting water content changes as small as a few cm of equivalent water height. Sev- eral gravimetric systems have been used to study hydrological pro- cesses. At large scale, space borne gravimetry such as the Gravity Recovery and Climate Experiment (GRACE) allows quantifying, with unprecedented accuracy, water storage variation of continental scale basins, thus bringing new constraints on global hydrological cycle (Seo et al. 2006). However, µGal accuracy of GRACE gravity field is limited to wavelengths larger than 1000 km, therefore pre- cluding local hydrological studies. On the ground, superconducting gravimeters offer continuous gravity monitoring (Goodkind 1999). Because of their high resolution (0.01 µGal) and nearly continuous sampling, they are well suited to study the local hydrology-induced gravimetric effects such as those due to rainfall, soil moisture, evap-

otranspiration and water table change on measurement site vicin- ity (Bower & Courtier 1998; Kroner 2001; Takemoto et al. 2002; Abe et al. 2006; Harnisch & Harnisch 2006; Imanishi et al. 2006; Van Camp et al. 2006). Unfortunately, superconducting gravimeters have not specifically been set up for hydrological studies, rather for studying geodynamics and Earth's structure (Hinderer & Crossley 2000).

The FG5 absolute gravity (AG) meter manufactured by Micro-g LaCoste provides 1 µGal (10⁻⁸ m s⁻<fnref n="2" />) accuracy (Niebauer et al. 1995) and can detect an equivalent of 5 cm water table variation. Despite its lower accuracy compared to superconducting gravime- ters, it has the advantage to be portable and therefore allows measurements at different sites of a hydrological basin. This leads to successful detection of water storage variations with a single in- strument, especially in complex hydrological situations, such as on karst aquifers (Jacob et al. 2008).

In spite of the insights brought by gravimetry to the understanding of the hydrological cycle, all the instruments above cannot help solving the depth source. Indeed, the inverse problem of inferring mass distribution at depth is ill conditioned if only surface gravity measurements are performed. This aspect severely limits the use of gravimetry for hydrological purposes, as it is often desirable to know not only horizontal but also vertical distribution of water storage

<footnote n="1">¹Géosciences Montpellier, UMR CNRS/UM2 5243, Université Montpellier II, Montpellier, France. E-mail: Thomas.Jacob@gm.univ-montp2.fr</footnote>

<footnote n="2">²Institut de Physique du Globe de Strasbourg, UMR CNRS/ULP 7516, Université L. Pasteur, Strasbourg, France</footnote>
variation. To overcome this intrinsic limitation of the gravimetric method, both surface and underground measurements are needed to measure the attraction of time-varying perturbing masses (here the water) from above and below.

In this paper, we investigate the time evolution of surface to depth (STD) gravity differences on a karst system using gravity measurements at the surface and at 60 m depth at the base of a pothole, using a portable Scintrex CG-5 Autograv relative gravimeter. On the studied karst system, AG measurements are performed on a monthly basis at three sites since January 2006. Observed gravity variations have been linked to water storage variations within the karst through mass balance modelling (Jacob *et al.* 2008). One of the aims of this paper is to determine at what depth within the karst water storage variations occur. Indeed, previous hydrochemistry studies (Plagnes 1997) invoke the saturated zone as the main storage entity in this karst system. It has, however, been argued that water storage variations driving AG variations must occur in the unsaturated zone (Jacob *et al.* 2008). Within the unsaturated zone, the epikarst zone, the uppermost weathered zone of the karst, is of particular interest because it can potentially be a storage entity (Mangin 1975; Klimchouk 2004; Williams 2008). The STD gravity experiment monitors gravity variations related to water storage change within a known depth, which comprises the epikarst. Using a Bouguer plate approximation, we compare the STD differences with AG variations and find that they are of the same order.

We then use the STD gravity measurements to determine apparent physical properties for the investigated depth. We first determine apparent densities derived from STD gravity measurements and compare the latter with grain densities measured from rock samples. In this manner, an apparent porosity range and apparent saturation changes representative of the investigated depth are determined. We then discuss the implications of an epikarst overlying an infiltration zone on porosity, depth and saturation change. Finally, we discuss the fact that AG monitors most water storage variations occurring in the epikarst on the studied karst. A conceptual model of karst functioning is given.

## 2 THE DURZON KARST SYSTEM

The Durzon karst system is located in the Grands Causses area, southern French Massif Central. This aquifer is embedded in a 400-m-thick formation of middle-to-upper Jurassic limestones and dolomites, deposited on top of a 200 m thick upper-Liassic marl formation (Bruxelles 2001a; Fig. 1). This latter formation acts as an impermeable barrier that defines the lower limit of the saturated zone of the karst system. Middle Jurassic formations are predominant on the recharge area (see Fig. 1), particularly a 200 m thick dolomite formation of Bathonian age, which largely outcrops.

On the Durzon karst aquifer, the main recharge comes from rainfall, which infiltrates at the surface of the spring catchment (750 m elevation on average), and the discharge occurs at the Durzon Spring (533 m elevation). Thanks to its monitoring by the Parc National des Grands Causses, this perennial spring is known to have a mean daily discharge of 1.4 m³ s⁻¹ (calculated over the 2002–2007 period), with maximum daily discharges reaching 18 m³ s⁻¹ during high-flow events.

During exceptionally long rainfall periods, temporary lakes (shown in blue in Fig. 1) appear at the south of the recharge area (Plagnes 1997; Bruxelles 2001a), their last appearance being in 1996. These lakes reflect that a large amount of precipitation may sometimes exceed the infiltration capacity of the epikarst and infiltration zone towards the saturated zone (Ricard & Bakalowicz 1996; Plagnes 1997; Bruxelles 2001a). The vadose zone, including the epikarst and infiltration zone, has therefore a different behaviour in the North and South of the recharge area. It is well drained in the North and inefficiently drained in the South (Fig. 1; Ricard & Bakalowicz 1996; Bruxelles 2001b). The STD gravimetric site study is located in a well-drained zone, as attested by some deep potholes and caves reaching depth of 60–170 m below the surface.

In agreement with an accepted karst model (Mangin 1975), the karst system may be horizontally layered into three zones (see Fig. 2): (1) the epikarst zone including both soil and weathered rock has 5–30 m depth. This zone has a high secondary porosity and is expected to be an important water reservoir (Williams 1983,

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Hydrogeological location map of the studied area, modified after Jacob *et al.* (2008). Absolute gravity sites CANA, BLAQ, SALV and surface to depth gravity site BEAU are indicated by black dots.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic functional diagram of the Durzon karst, showing the main karst structures, see text for explanations. Not to scale for horizontal distances.</figcaption></figure>

2008). (2) the infiltration zone below is mostly composed of massive rock with penetrative fissures and conduits, therefore making a fast vertical water transfer possible. (3) The saturated or phreatic zone formed by large conduits insures a mostly horizontal water flow to the outlet.

Karst spring hydrographs are characterized by high flow events consecutive to important precipitations and a base flow component. High flow events necessitate fast water infiltration to the phreatic zone through enlarged drainage shafts (see Fig. 2) and are generally short-lived. Base flow discharge is sustained by groundwater storage within the karst. This water storage's location differs from karst to karst and lies within the epikarst and/or in the phreatic zone (Bakalowicz 2005).

STD gravity measurements appear well suited to probe water storage variation in the epikarst.

## 3 SURFACE TO DEPTH GRAVITY EXPERIMENT

The STD gravity experiment consists of measuring the time-lapse gravity difference between the surface and depth at a given site.

### 3.1 Experimental setup

The experiment site is called La Beaumelle (BEAU) and is located to the northwest of BLAQ site (see Fig. 1). It corresponds to a vertical shaft or pothole 20 m long, 5 m wide and more than 50 m deep and allows an easy descent (see Fig. 3) at this depth. The surface lithology consists of Bathonian-aged homoclinal dolomites. We performed six STD measurements between September 2006 and April 2008 (see Table 1) using a Scintrex CG-5 relative spring gravimeter. We used CG-5#167 for all measurements except $t_3$, for which we operated #323. Scintrex relative gravimeters have been used in the past for precise microgravimetry surveys (Jousset et al. 2000; Ballu et al. 2003; Bonvalot et al. 2008; Merlet et al. 2008).

The CG-5 instrument has a reading resolution of 1 µGal and a repeatability of less than 10 µGals (Scintrex 2007). Its gravity sensor is based on a capacitive displacement transducer electrostatic feedback system to detect movements of the fused quartz proof mass and to force the mass back to a null position. The drift of the CG-5 sensor is caused by an unavoidable creep of the quartz spring, whose length under tension increases. For one STD measurement, several occupations of the surface and depth sites are done to constrain

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Topography surrounding BEAU site. Contour lines are every 1.5 m. Elevations are in metres. Blue triangle, surface measurement site; red circle, depth measurement site. The pothole and chamber are represented in grey. (b) Cross-section along line AB: $z_s$, surface site elevation; $z_d$, depth site elevation; $h$, elevation difference between the surface and depth sites, measuring some 60 m.</figcaption></figure>
<table><caption>**Table 1.** Results of the least square inversion for the different time periods.</caption><thead><tr><th>Date</th><th>Surface occup.</th><th>Depth occup.</th><th>Calibration correction factor</th><th>Δ<sub>STDg</sub> (mGal)</th><th>σ<sub>STD</sub> (mGal)</th><th>RMS residuals (mGal)</th></tr></thead><tbody><tr><td>t<sub>0</sub>:19/09/06</td><td>4</td><td>3</td><td>1.00042</td><td>-5.8924</td><td>0.0015</td><td>0.0031</td></tr><tr><td>t<sub>1</sub>:07/11/06</td><td>5</td><td>4</td><td>1.00051</td><td>-5.8664</td><td>0.0015</td><td>0.0037</td></tr><tr><td>t<sub>2</sub>: 26/01/07</td><td>4</td><td>3</td><td>1.00065</td><td>-5.8732</td><td>0.0012</td><td>0.0021</td></tr><tr><td>t<sub>3</sub>: 07/09/07</td><td>5</td><td>4</td><td>1.00000</td><td>-5.8865</td><td>0.0024</td><td>0.0074</td></tr><tr><td>t<sub>4</sub>: 09/11/07</td><td>4</td><td>3</td><td>1.00030</td><td>-5.8835</td><td>0.0012</td><td>0.0019</td></tr><tr><td>t<sub>5</sub>: 04/02/08</td><td>4</td><td>3</td><td>1.00040</td><td>-5.8608</td><td>0.0012</td><td>0.0018</td></tr></tbody></table>

this drift. Typically, the surface site was occupied between 4 and 5 times and the depth site 3 to 4 times for one STD measurement (see Table 1). For each occupation of a site, five sets of 90 s at 6 Hz sampling measurements were performed.

Special attention was given so that the CG-5 gravity sensor was at a fixed height and location for all STD measurements. This was achieved by fixing the height of the instrument's tripod with a brass ring and carving the rock at the surface and depth sites to fix the tripod location.

## 3.2 Surface to depth gravity difference evaluation

### 3.2.1 Data corrections

The gravity measurements are corrected for Earth tides using ET-GTAB software (Wenzel 1996) with the Tamura tidal potential development (Tamura 1987), ocean loading effects with FES2004 ocean tide model (Lyard et al. 2006) and atmospheric pressure loading using an admittance value of -0.3 µGal hPa⁻¹. Atmospheric pressure is measured at SALV station (see Fig. 1) at a 15 mn rate with a 0.1 hPa precision. Polar motion effects are not corrected because they are constant for the surface and depth site over the time span of one STD measurement (~4 hr), and we are looking at the gravity difference between these two sites.

### 3.2.2 Instrument calibration

As demonstrated by Budetta & Carbonne (1997), Scintrex relative gravimeters need to be regularly calibrated when used to detect small gravity variations over extended periods of time. The calibration correction factor with respect to the instrument calibration constant is constrained by performing calibration line measurements, where large known differences in gravity are compared with those measured by the relative gravimeter (Debeglia & Dupont 2002). Fig. 4 shows the calibration factor change for CG-5 #167, calculated from various calibration lines. The Scintrex calibration line is known as the Orangeville calibration line and is located in Ontario, Canada. It spans 119 mGals between five stations over a 70 km distance (Scintrex limited 2007). The accuracy on the calibration factor is 10⁻⁴. The Aigoual calibration line is composed of three stations between Montpellier and Mont Aigoual, south of France, spanning more than 300 mGal. The accuracy on the calibration is also 10⁻⁴. The Larzac calibration line is between the three AG stations CANA, BLAQ and SALV, spanning some 20 mGals. The calibration factor is calculated for gravity surveys composed of 40 stations, including the three AG stations BLAQ, CANA and SALV. Its accuracy is also 10⁻⁴, due to the strong ties between the stations (Jacob, unpublished results). Following Budetta & Carbonne (1997), a second-degree polynomial was fitted to the data. The evolution of the calibration factor is of one part in a thousand over a 2 yr period.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of the calibration factor for CG-5 no. 167; grey lines represent the BEAU measurement dates, purple dashed lines represent the 95 per cent confidence interval of the fit.</figcaption></figure>

The polynomial fit is used to obtain the calibration correction factor during the BEAU site measurement periods. The interpolated calibration factor for each measurement period is shown in Table 1. Taking into account the calibration factor evolution yields µGal order differences on the STD gravity difference ΔSTDg. Furthermore, the error of the calibration factor change does not significantly affect the corrected ΔSTDg. Because no calibration data is available for CG5 #323, we consider a calibration correction factor of 1.

### 3.2.3 Least-square adjustment

The gravity difference between the surface and depth ΔSTDg for one time period is obtained in the following manner. The instrumental drift is assumed linear due to the short time span of the measurements (~4 hr). The surface and depth measurements can be expressed as the sum of a linear drift and a constant, forming the following observation equations:

$$C_f(m_s)_i + v_i = Dt_i + Ks, \quad (1)$$

$$C_f(m_d)_j + v_j = Dt_j + Kd, \quad (2)$$

where C<sub>f</sub> is the calibration correction factor, (m<sub>s</sub>)<sub>i</sub> the surface gravity reading at time i in mGal, (m<sub>d</sub>)<sub>j</sub> the depth reading at time j in mGal, v<sub>i</sub> and v<sub>j</sub> the residuals, D the drift in mGal s⁻¹ and K<sub>s</sub> and K<sub>d</sub> constants in mGal. These constants equate to the surface and depth gravity value g<sub>s</sub> and g<sub>d</sub> plus the same constant. K<sub>s</sub> - K<sub>d</sub> therefore yields the STD gravity difference ΔSTDg.
The matrix representation of the observation equations for n gravity readings is

$$L + V = AX, \qquad (3)$$

where L is a $n \times 1$ vector of relative gravity readings with a weight matrix P given by the inverse of the variance of the measurements, V is the $n \times 1$ matrix of residuals, A is the design matrix and X is a $3 \times 1$ vector of unknowns, that is, the drift D and the constants Ks and Kd. The variance for one 90-s-gravity reading is the square of its standard error, which reflects the amount of microseismicity during the reading. Additional errors come from inaccurate gravity corrections and instrument setup.

The STD gravity difference $\Delta_{STDg}$ equates to the difference between Ks and Kd. Because Ks and Kd are not independent but linked by the drift variable D, the standard deviation $\sigma_{STD}$ of the STD gravity difference is

$$\sigma_{STD} = \sqrt{\sigma_{Ks}^2 + \sigma_{Kd}^2 - 2\sigma_{Ks-Kd}^2} \qquad (4)$$

where $\sigma_{Ks}^2$ and $\sigma_{Kd}^2$ are the variances of Ks and Kd and $\sigma_{Ks-Kd}^2$ is the covariance between Ks and Kd, given by the a posteriori covariance matrix from the least-square adjustment. To account for the calibration correction factor error of $10^{-4}$, 0.001 mGal is quadratically added to $\sigma_{STD}$. Table 1 summarizes the results.

$\sigma_{STD}$ for measurement at time $t_3$ is high. This is due to the fact that the repeatability of CG-5 #323 is not as good as that of CG-5 #167. Fig. 5 shows the gravimeter readings and the least-square adjusted linear trends as well as the residuals from these trends for measurement period $t_4$. The small dispersion of the residuals (less than 0.005 mGals) illustrates the good quality of the gravity readings and hence the robustness of the STD gravity difference determination.

The $\Delta_{STDg}$ values show significant variations ranging from -5.8608 to -5.8924 mGal. These variations are now interpreted.

## 3.3 Interpretation

### 3.3.1 Vertical gravity differences from a homogeneous layered model

The following hypotheses are made to interpret the $\Delta_{STDg}$ variations. We assume the sedimentary formations between the two measurement sites to be horizontal and of uniform density $\rho_{app}$, and that water storage and its variations occur homogeneously within these beds. The density structure of the rock mass below the depth site is assumed to be uniform of density $\rho$. Let $g_s$ and $g_d$ be the gravity value respectively at the surface and at depth at heights $z_s$ and $z_d$, respectively (see Fig. 3). We define $h$ as the height difference between the two sites.

The surface and depth gravity $g_s$ and $g_d$ corrected for all known tidal, polar motion and atmospheric effects can be expressed in the following manner:

$$g_s = 2\pi G\rho_{app}h + 2\pi G\rho z_d + T_s + z_s \text{grad}(g_0) + g_0(\phi_s) + \Delta g_B(z_s) + g_{LW}^s, \qquad (5)$$

$$g_d = -2\pi G\rho_{app}h + 2\pi G\rho z_d + T_d + z_d \text{grad}(g_0) + g_0(\phi_d) + \Delta g_B(z_d) + g_{LW}^d, \qquad (6)$$

where G is the universal gravitational constant, $\rho$ the density of the beds, $T_s$ and $T_d$ the terrain effects for the surface and depth sites, grad($g_0$) the vertical normal gravity gradient known as the free-air gradient, $g_0(\phi_s)$ and $g_0(\phi_d)$ are the normal gravity for the surface and depth sites at latitudes $\phi_s$ and $\phi_d$, $\Delta g_B$ the Bouguer anomaly and $g_{LW}^d$ and $g_{LW}^s$ are the long wavelength effect of global hydrology.

The Bouguer gravity anomaly is caused by spatial variations in the gravity field mainly induced by lateral variations in the density and thickness of the Earth's crust. The long-wavelength effect of global hydrology is dominated by surface deformation induced by

<figure><img src="image_6.png" /><figcaption>Figure 5. (a) Gravimeter readings corrected for the classic corrections (see text) at time $t_4$; dotted lines mark the linear trends fitted to the data, (b) Residuals of the fit for this measurement. Error bars represent the standard errors of the 90 s gravity readings.</figcaption></figure>
hydrological continental loading. Over the pothole depth, the vertical deformation induced by hydrological loading is constant. We therefore consider the long wavelength gravity effects due to global hydrology to be constant over the pothole depth. The STD gravity difference $\Delta_{STDg}$ can therefore be expressed as follows:

$$ \Delta_{STDg} = 4\pi G \rho_{app} h + \Delta_{STD} T + h \text{grad}(g_0) + \Delta_{STDg_0}(\phi) + \Delta_{STD} \Delta g_B, \quad (7) $$

where $\Delta_{STD} T$ is the difference in terrain effect between surface and depth sites, $\Delta_{STDg_0}(\phi)$ is the change in gravity due to the latitude difference and $\Delta_{STD} \Delta g_B$ is the difference in Bouguer anomaly between the two sites. These terms are explicitly described and determined in Sections 5.1.2–5.1.4.

Let us now consider the time-evolution of the terms of eq. 7. The free-air gradient and normal gravity are constant with time. The height difference $h$ can also be considered time independent. The time evolution of $\Delta_{STD} T$ is negligible for two reasons. First, topography does not evolve over the studied time period. Second, taking into account the density evolution due to water storage variation (see Table 4) yields a negligible time evolution of $\Delta_{STD} T$ (<0.001 mGal). Finally, the time evolution of $\Delta_{STD} \Delta g_B$ can also considered negligible, as homogenous water storage variations within the plate are invoked.

Therefore, the evolution of $\Delta_{STDg}$ over time period $\delta t$ is

$$ \Delta_{\delta t} \Delta_{STDg} = 4\pi G \Delta_{\delta t} \rho_{app} h, \quad (8) $$

where $\Delta_{\delta t} \rho_{app}$ is the apparent density change over time $\delta t$ within height $h$.

Hence $\Delta_{\delta t} \Delta_{STDg}$ equates to twice the Bouguer attraction of a plate of density $\Delta_{\delta t} \rho_{app}$ and height $h$. The time evolution of the density between the two measurement sites can therefore be monitored.

Under the hypothesis that water storage and its variation is uniformly spread out, the evolution of plate density with time can only be attributed to variations in the amount of stored water within that plate. Processes that have an effect on the density of carbonates, such as dissolution and erosion, can be discarded because they occur over much longer time periods.

Let $\Delta_{\delta t} \Delta_{STDl}$ be the equivalent water level change over time $\delta t$ within height $h$. $\Delta_{\delta t} \Delta_{STDl}$ induces the density change $\Delta_{\delta t} \rho_{STD}$. The time-evolution of $\Delta_{STDg}$ can therefore also be expressed in the following manner:

$$ \Delta_{\delta t} \Delta_{STDg} = 4\pi G \rho_w \Delta_{\delta t} \Delta_{STDl}, \quad (9) $$

where $\rho_w$ is the water density and $\Delta_{\delta t} \Delta_{STDl}$ the equivalent water level change over time $\delta t$ within height $h$. The time lapse STD gravity difference monitors the equivalent water level change within the known depth $h$. The detection of water storage change is facilitated because the time evolution of STD gravity differences senses twice the Bouguer plate effect of these water storage changes. Furthermore, 90 per cent of the gravity effect felt by $\Delta_{\delta t} \Delta_{STDg}$ is derived from within five times the vertical separation between measurement sites (McCulloh 1965). This means that 90 per cent of the gravity effect comes from a cylinder with ~300 m radius and height $h$.

## 4 COMPARISON TO SURFACE ABSOLUTE GRAVITY MEASUREMENTS

AG is monitored on the karst system at three sites with monthly measurements: CANA, BLAQ and SALV (see Fig. 1) at, respectively, 693, 773 and 749 m elevation. A detailed description of these measurements is done in (Jacob et al. 2008). To obtain a gravity signal related exclusively to local hydrological changes, the regional or long wavelength contribution of hydrology has to be removed. Regional hydrology-induced gravity changes are corrected for using the European Center for Medium-range Weather Forecasts (ECMWF; available at www.ecmwf.int) soil moisture model and snow model (Viterbo & Beljaars 1995). This is done by convolving surface mass distribution associated with soil moisture with the Green's functions associated with Newtonian and deformation effects. The amplitude of the corrected signal is 2 µGals at most; as shown in Jacob et al. (2008). A mass balance model taking into account rainfall, actual evapotranspiration and spring discharge yields the average water stock in the karst system for daily time steps. Precipitation is measured at BLAQ station at quarter-hourly time steps, daily actual evapotranspiration is calculated using Penmann-Monteith's formula (Penman 1948) for daily potential evapotranspiration, scaled by Turc's yearly total actual evapotranspiration (Turc 1961; also see Jacob et al. 2008). This water stock is then converted into a gravimetric water stock (GWS) using a Bouguer plate approximation, which can then be compared with the AG variations. The GWS plotted on Fig. 5 reasonably accounts for the observed AG variations. However, SALV gravity is systematically higher than that of the other AG sites and the GWS, this site is therefore interpreted as a preferential water storage zone (Jacob et al. 2008). The AG variations at the three karst sites are shown in Fig. 6.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Absolute gravity variations corrected for the regional contribution of hydrology calculated from ECMWF model at the three karst sites. The gravimetric water stock (GWS) from the mass balance modelling is represented. Dates of the surface to depth experiments are marked from $t_0$ to $t_5$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Surface to depth gravity differences compared to AG values and GWS for the same time periods. AG and GWS are shifted along the y-axis for legibility. Dashed lines, linear fits to the data.</figcaption></figure>

<table><caption>**Table 2.** Fit parameters between $\Delta_{STDg}$ and AG values and GWS. $R^2$ is the coefficient of determination.</caption><thead><tr><th></th><th>Slope s</th><th>σ slope</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>CANA</td><td>0.273</td><td>0.013</td><td>0.95</td></tr><tr><td>BLAQ</td><td>0.481</td><td>0.037</td><td>0.89</td></tr><tr><td>SALV</td><td>0.392</td><td>0.087</td><td>0.49</td></tr><tr><td>GWS</td><td>0.537</td><td>0.031</td><td>0.93</td></tr></tbody></table>

Within the homogeneous layered model (see Section 3.3.1), a surface AG change over time $\delta t$ corrected for global hydrology is

$$ \Delta_{\delta t} g_{AG} = 2\pi G \rho_w \Delta_{\delta t} l_{AG}, \quad (10) $$

where $\Delta_{\delta t}l_{AG}$ is the equivalent water height change occurring beneath the gravimeter over an undetermined depth interval. Surface time-lapse AG therefore yields information on the amplitude of water storage variation; it can, however, not determine at which depth these variations occur.

STD gravity differences are compared with the AG values and GWS for the same time periods (see Fig. 7). AG values for the STD time periods $t_0$ to $t_5$ are interpolated using a piecewise cubic Hermite interpolating polynomial shown in Fig. 6. Error bars on the estimated AG values are evaluated to 3 $\mu$Gals. Correlation between $\Delta_{STDg}$ and both AG measurements and GWS is now examined. Linear trends are least-square adjusted to the data sets (see Fig. 7). The fit parameters are shown in Table 2.

The best correlation between AG sites and STD measurements occur at sites CANA, followed by BLAQ. Moderate correlation is observed between SALV and $\Delta_{STDg}$. GWS and $\Delta_{STDg}$ show good correlation. The slope $s$ of the linear trends between AG, GWS and $\Delta_{STDg}$ and is of prime interest. Indeed, this slope is the ratio between eqs 10 and 9:

$$ s = 0.5 \Delta_{\delta t} l_{AG} / \Delta_{\delta t} \Delta_{STD} l. \quad (11) $$

Therefore $s$ equates to half the ratio of the equivalent water level change at an AG site over an undetermined depth to the equivalent water level change at BEAU within height $h$. The slope $s$ is smaller than 0.5 for CANA AG, therefore $\Delta_{\delta t}l_{CANA} < \Delta_{\delta t}\Delta_{STD}l$. Water

storage variations are therefore lower at CANA than at BEAU site. The slope $s$ is smaller than 0.5 for SALV AG, however, poor data set correlation precludes to infer a simple water storage relationship between BEAU and SALV sites. For BLAQ AG and GWS, $s \sim 0.5$, therefore $\Delta_{\delta t}\Delta_{STD}l \sim \Delta_{\delta t}l_{BLAQ}$ and $\Delta_{\delta t}\Delta_{STD}l \sim \Delta_{\delta t}l_{GWS}$. Water storage variations at BEAU site are similar to those at BLAQ site and also to those of the GWS. The implications of these findings will be discussed in Section 6.2.

# 5 MEAN RESERVOIR PROPERTIES FROM STD GRAVITY AND CORE SAMPLE MEASUREMENTS

The aim of this section is to determine mean porosities and saturation changes for the investigated depth from STD gravity measurements and core sample measurements. All physical properties derived from STD gravity will be referred to as 'apparent' properties. Apparent densities are first determined, then core samples are analysed to obtain the grain density. The apparent porosity is obtained from the comparison of apparent and grain density.

## 5.1 Apparent density determination

The gravity difference between two depths can be used to determine the apparent density structure between these depths, assuming that the homogeneous layered approximation discussed in Section 3.3.1 is valid. This technique has been used in mine shafts (Hussain et al. 1981) and in boreholes (LaFehr 1983; Kasameyer & Hearst 1988). The apparent density derived from $\Delta_{STDg}$ is obtained from eq. 7, using the fact that the terrain effect is a linear function of density:

$$ T(\rho) = \rho T(\rho = 1), \quad (12) $$

$$ \rho_{\text{app}}(t) = \frac{(\Delta_{\text{STD}}g(t) - h \operatorname{grad}(g_0) - \Delta g_0(\phi) - \Delta_{\text{STD}}g_B)}{(4\pi G h + \Delta_{\text{STD}} T(\rho = 1))}. \quad (13) $$

This apparent density is that of the slab between the two measurement sites, therefore corresponding to a mean density representative of a large volume of rock. The apparent density determination needs the precise evaluation of following terms (see eq. 13): the height $h$ between the two sites, the difference in topographic effects $\Delta_{STD}T$ between surface and depth sites, the latitude dependent term $\Delta g_0(\phi)$ and the Bouguer term, which will be expressed in terms of gravity gradient.

### 5.1.1 Height determination

The difference in height between the surface and depth site was acquired using a combination of three techniques: geometrical levelling; measuring tapes and water tubes. We measure the height difference $h$ between the surface and depth site to be 63.685 m, with an estimated measurement error of 0.2 m.

### 5.1.2 Terrain effect calculation

The terrain effects $T_s$ and $T_d$ from eqs 5 and 6 need to be calculated to relate to a Bouguer plate of thickness $h$. The stations surrounding topography has a critical effect on the terrain effects (Nowel 1999), and therefore a digital elevation model (DEM) based on a real time kinematic GPS survey is used to calculate terrain effects. This DEM covers an area 100 m around the stations with 1 m grid spacing. The topography further away is taken from the Shuttle Radar Topography Mission (SRTM) data (Farr & Kobrick 2000; available
<table><caption>**Table 3.** Results for the terrain effects and chamber and pothole effects on both sites.</caption><thead><tr><th></th><th>Terrain effects (mGal)</th><th>Chamber and pothole (mGal)</th><th>S. D. chamber</th><th>Total (mGal)</th><th>Difference (mGal)</th><th>Error (mGal)</th></tr></thead><tbody><tr><td>Surface</td><td>-0.231</td><td>-0.027</td><td>0.029</td><td>-0.258</td><td></td><td></td></tr><tr><td>Depth</td><td>-0.648</td><td>0.439</td><td>0.057</td><td>-0.209</td><td>-0.050</td><td>0.058</td></tr></tbody></table>

at http://www.jpl.nasa.gov/srtm/) with a 3" grid. The pothole and a chamber adjacent to the pothole (see Fig. 3) contribute to the local terrain effects. The dimensions of the chamber were measured using a compass and a measuring tape with a 2 per cent precision, and its average height was evaluated to 6 m. The topography, fault and chamber were triangulated, and the gravitational effects were calculated using Okabe's analytical solution (Okabe 1979) for triangular prisms. The evaluation of the errors on the terrain effects due to chamber and pothole dimension and shape was done by adding a normally distributed random noise to these structures' boundaries. One thousand computations were performed to determine an average error. 2 and 1 m standard deviation noises were added to the chamber and pothole boundaries, respectively. Results are shown in Table 3.

The terrain effect was calculated to within 22 km of the sites; terrain effects further away were not calculated because they are constant for the surface and depth sites. The following table summarizes the different contributions of the topographic effects calculated with a 2600 kg m⁻³ density. Due to the linearity of terrain effects with density (eq. 12), these effects can be obtained for any density.

Terrain effects are larger for the depth site because the topography above the surface site induces a larger attraction on the depth site than on the surface site. Indeed, being on a plateau, the surface site is at the same level with the surrounding topography, hence there is little deviation from the Bouguer plate hypothesis. By contrast, the depth site does sense the deviation from the upper boundary of the Bouguer plate much better, leading to a larger correction. The chamber and pothole have a very large effect on the depth site because it is located within the bottom of the pothole. The standard deviation associated to the pothole and chamber dimension uncertainties are consequently large for the depth site (see Table 2). Therefore, for a density of 2600 kg m⁻³, the difference between the surface and depth sites for terrain correction equates to -0.050 ± 0.058 mGal.

### 5.1.3 Latitude dependent term

The surface and depth sites are 20.32 m apart in the north-south direction (see Fig. 3a), which equates to a normal gravity difference of 0.0166 mGal, calculated with the Gravity Formula for Geodetic Reference System 1980.

### 5.1.4 Bouguer anomaly term

Following Kasameyer and Hearst (1988), we evaluate the Bouguer anomaly term in terms of anomalous gravity gradient. We define the regional masses to be those distant enough that their effect on the gradient is constant over the depth of the pothole. This regional gradient is none other than the vertical gradient derived from Bouguer anomalies. The Bouguer anomalies and associated gravity gradients were obtained from the Bureau de Recherches Géologiques et Minières (BRGM; Martelet et al. 2002). The studied area is located in a broad negative Bouguer anomaly, and therefore the associated gradient is very weak, it equates to -0.00036 mGal m⁻¹. Hence, the gradient at BEAU site equates to the free-air gradient

<table><caption>**Table 4.** Calculated apparent densities as a function of time.</caption><thead><tr><th>Date</th><th>Δ<sub>STDg</sub><br>(mGal)</th><th>ρ<sub>app</sub><br>(kg m<sup>-3</sup>)</th><th>ρ<sub>app no<br>pothole</sub><br>(kg m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>t<sub>0</sub>:19/09/06</td><td>-5.8924</td><td>2589</td><td>2498</td></tr><tr><td>t<sub>1</sub>:07/11/06</td><td>-5.8664</td><td>2593</td><td>2503</td></tr><tr><td>t<sub>2</sub>: 26/01/07</td><td>-5.8732</td><td>2592</td><td>2501</td></tr><tr><td>t<sub>3</sub>: 07/09/07</td><td>-5.8865</td><td>2590</td><td>2499</td></tr><tr><td>t<sub>4</sub>: 09/11/07</td><td>-5.8835</td><td>2590</td><td>2500</td></tr><tr><td>t<sub>5</sub>: 04/02/08</td><td>-5.8608</td><td>2594</td><td>2504</td></tr><tr><td>Mean</td><td>-5.8771</td><td>2591</td><td>2501</td></tr></tbody></table>

-0.30860 mGal m⁻¹ and the Bouguer anomalies derived gradient, yielding a vertical gradient of -0.30896 mGal m⁻¹.

### 5.1.5 Results

The apparent density for each ΔSTDg measurement is calculated taking into account all corrections (eq. 13); it is also calculated discarding the pothole and chamber effects. Results are summarized in Table 4.

The apparent densities range from 2589 to 2594 kg m⁻³, with a mean apparent density of 2591 kg m⁻³ when taking into account all corrections. When the pothole and chamber effects are discarded, the mean density is 2501 kg m⁻³. It is normal that the apparent density excluding pothole and chamber corrections is smaller than that including all corrections. Indeed, the air-filled volumes of the pothole and chamber contribute to lowering the overall density structure between surface and depth sites. By correcting for pothole and chamber effect, the apparent density is representative of the rock mass excluding these known air-filled cavities.

With the uncertainties described in Sections 5.1–5.4, the error budget on the apparent density calculation is evaluated to 22 kg m⁻³, mainly due to the error on the terrain effects. Indeed, for the apparent densities excluding pothole and chamber, the error falls to 8 kg m⁻³.

As discussed in Sections 3.3.1 and 3.3.2, the time-lapse gravity evolution senses the water storage variations within the investigated depth. These water storage variations have a direct impact on the evolution of the apparent density (Table 3).

## 5.2 Core sample porosity and density measurements

Fourteen Bathonian dolomite core samples from the surface and 12 from depth were analysed and their density and porosity measured. The dry mass $M_d$, the saturated mass $M_s$ and the immersed mass $M_{im}$ were weighted with a 10⁻³g precision. From these three weight measurements, the porosity $\omega$ and the grain density $\rho_g$ can be calculated.

$$ \omega = (M_s - M_d)/(M_s - M_{im}), \qquad (14) $$

$$ \rho_g = \rho_f M_d/(M_d - M_{im}), \qquad (15) $$

where $\rho_f$ is the fluid density, evaluated to be 980 kg m⁻³.

Knowing the porosity and the grain density, one can determine the bulk and dry densities $\rho_{bulk}$ and $\rho_{dry}$ of the core samples, which
<table><caption>Table 5. Mean porosities, grain densities, bulk densities and dry densities for the surface, depth and total samples.</caption><thead><tr><th></th><th>Samples</th><th>Mean poro.</th><th>Std poro.</th><th>Mean ρ<sub>g</sub></th><th>Std ρ<sub>g</sub></th><th>Mean ρ<sub>bulk</sub></th><th>Std ρ<sub>bulk</sub></th><th>Mean ρ<sub>dry</sub></th><th>Std ρ<sub>dry</sub></th></tr></thead><tbody><tr><td>Surface</td><td>14</td><td>9.2</td><td>2.7</td><td>2729.2</td><td>29.2</td><td>2569.7</td><td>54.4</td><td>2477.4</td><td>79.2</td></tr><tr><td>Depth</td><td>12</td><td>5.5</td><td>2.3</td><td>2707.9</td><td>14.9</td><td>2613.4</td><td>49.8</td><td>2558.0</td><td>72.6</td></tr><tr><td>Total</td><td>26</td><td>7.5</td><td>3.1</td><td>2718.5</td><td>25.7</td><td>2589.9</td><td>55.9</td><td>2514.6</td><td>25.7</td></tr></tbody></table>

are, respectively, the densities of the saturated and dry samples. Table 5 shows the results.

The surface samples are characterized by a mean grain density and porosity of, respectively, 2729.2 ± 29.2 kg m⁻³ and 9.2 ± 2.7 per cent, and the depth samples have a mean grain density and porosity of 2707.9 ± 14.9 kg m⁻³ and 5.5 ± 2.3 per cent, respectively. Each weight measurement is known to ±2.10⁻³g, the error on the porosity therefore ranges between 0.01 and 0.02 per cent and the error on the densities is inferior to 1 kg m⁻³. The surface samples are characterized by a higher porosity compared with the depth samples. The higher porosity can be attributed to a more effective weathering at the surface. The mean grain density being larger than the mean apparent density, this demonstrates that an apparent porosity ωₐpp does exist and can be evaluated.

## 5.3 Apparent porosity and saturation determination

### 5.3.1 Overview of the measured and calculated densities

The measured core densities and calculated apparent densities are plotted in Fig. 8.

The apparent gravimetric density is included into the bulk surface and depth samples density range; this demonstrates that the calculated apparent density is coherent. The grain density for both depth and surface samples is logically some 200 kg m⁻³ larger than the apparent density. The apparent density is closest to bulk density

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of core samples grain densities ρ<sub>g</sub>, bulk densities ρ<sub>bulk</sub>, dry densities ρ<sub>dry</sub> and calculated apparent densities ρ<sub>app</sub>. Density ranges are represented.</figcaption></figure>

of the depth samples. The mean grain density is now compared with the apparent density to estimate an apparent porosity.

### 5.3.2 Apparent porosity and saturation change

The apparent density for every time period is compared with the mean grain density measured from rock samples, yielding apparent porosities and apparent saturation changes. Indeed, the apparent density ρ<sub>app</sub> can be expressed as a function of the grain density ρ<sub>g</sub>, the apparent porosity ω<sub>app</sub> and the void density ρ<sub>v</sub>:

$$ \rho_{\text{app}}(t) = \rho_{\text{g}}(1 - \omega_{\text{app}}) + \rho_{\text{v}}(t)\omega_{\text{app}}. \quad (16) $$

The void density is the product of saturation S(t) and water density ρ<sub>w</sub>:

$$ \rho_{\text{v}}(t) = S(t)\rho_{\text{w}}. \quad (17) $$

Combining eqs 16 and 17, the apparent porosity ω<sub>app</sub> can therefore be expressed as follows:

$$ \omega_{\text{app}}(t) = [\rho_{\text{app}}(t) - \rho_{\text{g}}]/[S(t)\rho_{\text{w}} - \rho_{\text{g}}]. \quad (18) $$

The saturation S(t) is unknown; however, its range, between 0 and 1, can be useful to determine a range for the porosity. For an acceptable porosity range ω<sup>r</sup><sub>app</sub>, an apparent saturation change range Δ<sub>δt</sub>S<sup>r</sup><sub>app</sub> over time period δt can thus be retrieved:

$$ \Delta_{\delta t} S_{\text{app}}^{r} = (\Delta_{\delta t} \rho_{\text{app}} - \rho_{\text{g}}(1 - \omega_{\text{app}}^{r})) / \rho_{\text{w}} \omega_{\text{app}}^{r} \quad (19) $$

where Δ<sub>δt</sub>ρ<sub>app</sub> is the apparent density change over time period δt.

The apparent porosity and the degree of pore saturation are unknowns that need to be constrained from eq. 18. The pore saturation can range between 0 and 1, whether the rock mass is dry or saturated. Table 6 shows the ranges in porosities calculated from eq. 18 for the different time periods for dry and saturated pores.

Because the porosity does not vary during the studied time period, taking the maximum value from the dry porosities and the minimum value from the saturated porosities yields the acceptable porosity range ω<sup>r</sup><sub>app</sub>. We find ω<sup>r</sup><sub>app</sub> to be in the range from 4.8 to 7.3 per cent for the investigated volume. The errors on these porosity estimates are 1.3–2 per cent for, respectively, the dry and wet porosities. The apparent porosity range is compared to the depth and surface samples porosities obtained from eq. 14 in Fig. 9.

The apparent porosity is an integrative value for a 60 m thick volume of rocks, where voids of all scales such as pores, fissures and conduits may have an impact on it. The measured porosities from the

<table><caption>**Table 6.** Ranges for the calculated porosities for the different time periods, considering dry or saturated pores.</caption><thead><tr><th>Time period</th><th>Porosity (dry pores)</th><th>Porosity (saturated pores)</th></tr></thead><tbody><tr><td>t<sub>0</sub></td><td>4.8</td><td>7.6</td></tr><tr><td>t<sub>1</sub></td><td>4.6</td><td>7.3</td></tr><tr><td>t<sub>2</sub></td><td>4.7</td><td>7.4</td></tr><tr><td>t<sub>3</sub></td><td>4.8</td><td>7.5</td></tr><tr><td>t<sub>4</sub></td><td>4.8</td><td>7.5</td></tr><tr><td>t<sub>5</sub></td><td>4.6</td><td>7.3</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of depth and surface samples porosity ranges to the apparent porosity range. Error bars for the minimum and maximum values of the apparent porosity are represented.</figcaption></figure>

core samples are only representative of the sampled rock volumes where pores make up most of the porosity. The sample porosities and the apparent porosity range are not very different (see Fig. 9), this could mean that there are no large unknown open chambers or potholes within 300 m of the measurement sites. Each $\Delta_{STDg}$ measurement yields a different porosity to saturation relationship (eq. 18). In Fig. 10(a), saturation is represented with respect to apparent porosity for the various time periods. For any given apparent porosity, the saturation change between different time periods can be retrieved over the acceptable porosity range (eq. 19). This is illustrated in Fig. 10(b) for consecutive time periods.

Over the acceptable porosity range, the saturation change is between 6 and 9 per cent for recharge periods ($t_0-t_1$ and $t_4-t_5$), and can be as low as -3.5 to -5.5 per cent for discharge period $t_2-t_3$ (see Fig. 10).

The comparison of the apparent density to the measured grain density yields a likely range of porosities for the investigated volume and a range of saturation change over different time periods.

However, the aforementioned physical properties derived from STD gravity are mean values for the investigated depth. We now discuss how the presence of an epikarst overlying an infiltration zone influences the porosity and saturation change within the investigated depth.

# 6 DISCUSSION

## 6.1 Porosity and depth of the epikarst

Physical properties derived from the STD gravity method are mean values considered homogeneous over the investigated depth. However, epikarst structures overlying an infiltration zone are bound to exist within the investigated depth (Mangin 1975; Williams 1985). The phreatic zone is thought to lie more than 160 m below the surface, as attested by vertical cave systems reaching it at this depth, in the vicinity of BEAU site. The apparent physical properties derived from STD gravity measurements are therefore representative of both the epikarst and the uppermost part of the infiltration zone. The former homogeneous approximation could therefore not be representative of the aquifer structure; it is worthwhile to consider a layered case of an epikarst overlying an infiltration zone (see Fig. 11).

The epikarst, describing a horizon at the top of the vadose zone of a karst aquifer, is characterized by enhanced storage capacity due to its specific soil and rock properties. Indeed, its porosity is thought to be at least one order of magnitude higher than that of the underlying infiltration zone (Klimchouk 2004; Williams 2008). Estimates of overall epikarst porosity are few and range between 5 and 10 per cent (Williams 1985) and as much as 10–30 per cent (Williams 2008). However, these values are only speculative and may largely vary following rock type, climatic and other environmental factors. Due to its heterogeneous nature, no direct porosity measurement representative of the epikarst can generally be made. The lower limit of the epikarst also depends on the lithology and geomorphological history (Klimchouk 2004) and can range between 3 to 30 m or more (Williams 2008).

Using apparent properties derived from the STD gravity experiment, epikarst and infiltration zone heights, porosities and saturation changes are now discussed.

The apparent porosity is first used to yield information on epikarst porosity $\omega_{ep}$ and thickness $h_{ep}$.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) Saturation with respect to apparent porosity for the various time periods. (b) Saturation change over consecutive time periods with respect to apparent porosity.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Representation of the investigated depth in a homogeneous and in a layered case. In the homogeneous case, apparent porosity $\omega_{app}$ describes the media (see Section 4.3.2). The epikarst and infiltration zone make up the layered case, each having their specific thicknesses and porosities.</figcaption></figure>

The total void volume within the investigated depth is the sum of epikarst and infiltration zone void volumes, therefore yielding the following equations:

$$h\omega_{\text{app}} = h_{\text{ep}}\omega_{\text{ep}} + h_{\text{inf}}\omega_{\text{inf}}, \quad (20)$$

$$h = h_{\text{ep}} + h_{\text{inf}}, \quad (21)$$

where $h_{inf}$ and $\omega_{inf}$ are, respectively, the thickness and porosity of the infiltration zone within height $h$ (see Fig. 11). Combining eqs 20 and 21 yields the epikarst height $h_{ep}$ as a function of the normalized porosities $\Omega_{inf}$ between $\omega_{inf}$ and $\omega_{app}$ and $\Omega_{ep}$ between $\omega_{ep}$ and $\omega_{app}$:

$$h_{\text{ep}} = h(1 - \Omega_{\text{inf}})/(\Omega_{\text{ep}} - \Omega_{\text{inf}}). \quad (22)$$

We impose that $\Omega_{inf}$ be smaller than unity and that $\Omega_{ep}$ be larger than unity, that is, that the infiltration zone and epikarst porosity be, respectively, smaller and larger than the apparent porosity. This condition accounts for the a priori hydrological knowledge of the epikarst having a higher porosity than that of the infiltration zone (Klimchouk 2004; Williams 2008). Epikarst thickness is plotted in Fig. 12 with respect to $\Omega_{inf}$ and $\Omega_{ep}$.

Let us consider that the apparent porosity for the investigated depth is 6 per cent, the mean of the acceptable porosity range. If we consider an epikarst height $h_{ep}$ of 15 m and an infiltration zone porosity $\omega_{inf}$ of 1.5 per cent, that is, $\Omega_{inf} = 0.25$, then $\Omega_{ep} \sim 3.5$ (see Fig. 12), therefore the average epikarst porosity is $\sim 21$ per cent.

Let us now consider that the surface and depth rock samples yield representative properties for, respectively, the epikarst and the infiltration zone. Taking the mean values of the measured porosities (see Table 4) yields $\Omega_{inf} = 0.92$ and $\Omega_{ep} = 1.53$, which results in an epikarst thickness equal to 8.4 m (see Fig. 12). This thickness is not incoherent with thicknesses from the literature. However, the surface and depth sample porosities cannot be considered representative of the epikarst and infiltration zone mean porosities. The epikarst is a very weathered zone, where voids, dissolution widened fissures and high porosity detritic material exist (Bruxelles 2001b; Williams 2008). The porosity from the surface samples can therefore be considered as a lower bound for epikarst porosity. The depth cores were sampled within a pothole at the atmosphere-rock interface, where weathering occurs. They may be more porous than

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Epikarst height $h_{ep}$ as a function of the porosity ratios $\Omega_{ep} = \omega_{ep}/\omega_{app}$ and $\Omega_{inf} = \omega_{inf}/\omega_{app}$. Heights are in metres, and contour lines are every 2.5 m. Black lines represent $\Omega_{ep} = 1.53$ and $\Omega_{inf} = 0.92$, and grey shaded area mark the possible epikarst heights, see text for explanations.</figcaption></figure>

the average rock mass in the infiltration zone. Because of this, the depth samples porosity may be considered as an upper bound for infiltration zone porosity. The aforementioned conditions imply that acceptable epikarst depths lie within the domain defined by $\Omega_{ep} > 1.53$ and $\Omega_{inf} < 0.92$. This is represented in Fig. 12 as grey shaded area. Epikarst depths cannot be larger than 42.5 m ($\Omega_{ep} < 1.53$) and lower than 1 m ($\Omega_{inf} > 0.92$).

STD gravity measurements yield apparent densities, which when compared to measured sample grain densities, allow for the determination of an apparent porosity range. Coherent epikarst properties such as porosity and depth are then estimated with success.

## 6.2 Implications on water storage dynamics

We show that water storage variation within the first 60 m of the karst at BEAU site is very similar to that of BLAQ site over an undetermined depth. We shall now argue that the water storage variation at BLAQ site is likely to occur over the same depth interval as that of BEAU site.

Two main arguments are considered:

(1) BEAU and BLAQ sites lie in the same central zone of the karst system, separated by less than 2 km (see Fig. 1).

(2) Surface lithology is unchanged between the two sites, consisting of Bathonian dolomites having an alteration layer at the surface.

These arguments suggest that karst structure is unchanged between the two sites, and that similar storage properties exist at both sites.

If we assume the existence of a high porosity epikarst overly-ing a low porosity infiltration zone, it is tempting to conjecture that epikarst is the main storage entity. Indeed, water storage varia-tions and karst features being nearly identical between the two sites, it is plausible that gravity variation at BLAQ site may be driven by epikarst water storage variations. Unfortunately, potholes in the vicinity of the other AG sites do not exist, and therefore the results
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Conceptual model of Durzon karst functioning inferred from gravimetric measurements. A.D.S represents annex to drain systems (see text for explanations).</figcaption></figure>

found for BLAQ site cannot be confirmed for CANA and SALV sites with certainty. However, epikarst water storage could also be invoked to explain AG variations at CANA and SALV sites. Indeed, AG variations show the same trends at each site, a gravity rise after precipitation and a gravity decline during drier periods, due to evapotranspiration and deep infiltration (Jacob *et al.* 2008). Amplitudes of the aforementioned mechanisms are, however, different for each site (see Fig. 7). This may, in part, be due to the different storage and drainage properties of the epikarst at each site, yet also to different evapotranspirated amounts of water at each site.

Accounting for the global water budget on the aquifer (including rainfall, evapotranspiration and spring discharge), the GWS is representative of the time evolution of the average water storage for the karst system as a whole (Jacob *et al.* 2008). We demonstrate that water storage variation at BEAU site is very similar to that of the GWS (see Section 4). This may indicate that water storage variation on the Durzon karst occurs mainly within the STD experiment depth. In other words, our results suggest that the main dynamic storage entity accounting for seasonal scale water storage variations lies within the first 60 m of the karst. It is plausible that the epikarst plays the main storage function, as discussed above. This implies that the phreatic zone does not exhibit significant water storage variations on a seasonal timescale, and that it, along with the infiltration zone, has a transmissive function rather than a storage function. Water seems therefore mostly stored in the epikarst, and once it reaches the base of the epikarst by leakage, it goes through the infiltration zone to the phreatic zone and ultimately to the spring in a short time. Vertical shafts and fissures in the infiltration zone and drains in the phreatic zone make a fast transfer from the base of the epikarst to the spring possible. Epikarst leakage occurs continuously as long as epikarst storage is not depleted, and it contributes to the slow component of the Durzon discharge (base flow). This is materialized as 'continuous infiltration' in Fig. 13.

After important precipitation events, fast infiltration occurs, which leads to flood events. These events are interpreted by invoking the overflow of the epikarst (Bakalowicz 2005). As a threshold storage value is reached, the epikarst overflows into vertical shafts of the infiltration zone. This is materialized as 'fast infiltration' in Figs 2 and 13.

Some authors invoke karst cavities with high head loss connections to the phreatic conduits as important storage entities (Mangin 1975; Bakalowicz 2005). These are called annex-to-drain systems (see Fig. 13). Such entities may well exist within the Durzon karst. However, since epikarst storage variation accounts for average storage variation in the karst as a whole, storage variation within these low permeability volumes is small. The storage function of the epikarst is illustrated in other studies, as the one in the Milandre test site, Switzerland, where a study of isotopic data leads to the demonstration of dominantly epikarstic storage compared with phreatic storage (Perrin *et al.* 2003).

## 7 CONCLUSION

Vertical time-lapse relative gravity monitoring is a complimentary technique to surface AG monitoring. Under the assumption that water storage is spatially homogeneous, the time-lapse gravity difference within an investigated depth is representative of water storage variation within this depth. We compare the AG and STD gravity differences and find that they are very similar. Specifically, we find that the STD gravity differences and BLAQ AG variations very closely correspond. The STD gravity differences are used to determine the apparent density of the investigated zone. With the appropriate corrections, apparent densities ranging from 2589 to 2594 kg m⁻³ are calculated for the different measurements. Comparing these densities to measured grain densities from rock samples gives an insight on the average porosity of the first 60 m of the karst. The acceptable
porosity range found is from 4.8 to 7.3 per cent. With this porosity range, the pore saturation change for consecutive measurements is calculated. Over the studied time period, strong recharge is materialized by saturation increases of 6–9 per cent, whereas discharge periods show negative saturation changes as low as –5.5 per cent.

We argue that epikarst is an efficient storage entity due to its high porosity. As a consequence, we propose that STD gravity and AG mostly sense epikarst water storage. We claim that surface AG measurements on a karst system are an effective tool to quantify epikarst water storage variation. Furthermore, the GWS, representative of the karst system behaviour as a whole, compares favourably to both the AG time-series and the STD gravity differences. We come to the conclusion that on this karst system, the main water storage variation entity is the epikarst, and that the saturated zone only plays a transmissive function, with little water storage variation.

Time-lapse absolute and STD gravity monitoring of a karst system therefore provide original data, which enhance the understanding of karst aquifer functioning.

However, it is important to remind that any physical property derived from vertical gravity differences is an average value between the two gravity measurements sites. Therefore, depth variations of density, porosity and saturation cannot be determined with certainty. There is no way to discriminate if water storage occurs homogeneously within the investigated depth, or if it is concentrated within the first 20 or even 10 m. To address this matter, several measurement sites at different depths in the pothole are needed, every 10 m or so. This way, the depth interval presenting the highest gravity change should indicate where the water storage variations occur. This could lead to the determination of storage variations with depth. This strategy is similar to that of borehole gravimeters. Indeed, within a borehole, gravity differences can be determined at any depths, and the gravity-derived physical properties can be calculated over much shorter depth intervals, if not continuously, as a function of depth.

## ACKNOWLEDGMENTS

We would like to thank Sebastien Deroussi and Michel Diament for their advice with the Scintrex relative gravimeters provided by Institut National des Sciences de l'Univers (INSU). This project is part of the program 'Ecosphère Continentale' (ECCO) funded by the 'Agence Nationale de la Recherche' (ANR). We would like to thank Gilles Barrau for his precious help in getting the gravimeters at the bottom of the pothole and back to the surface. We also would like to express our gratitude towards Didier Loggia and Jeff Ritz for their help on the field and in the laboratory.

## REFERENCES

* Abe, M. et al., 2006. Hydrological effects on the superconducting gravimeter observation in Bandung, *J. Geodyn.*, **41**(1–3), 288–295.

* Amelung, F., Galloway, D.L., Bell, J.W., Zebker, H.A. & Laczniak, R.J., 1999. Sensing the ups and downs of Las Vegas-InSAR reveals structural control of land subsidence and aquifer-system deformation, *Geology*, **27**, 483–486.

* Bakalowicz, M., 2005. Karst groundwater: a challenge for new resources, *Hydrogeol. J.*, **13**(1), 148–160.

* Ballu, V., Diament, M., Briole, P. & Ruegg, J.-C., 2003. 1985–1999 gravity field variations across the Asal Rift; insights on vertical movements and mass transfer. *Earth planet. Sci. Lett.*, **208**(1–2), 41–49.

* Bonvalot, S., Remy, D., Deplus, C., Diament, M. & Gabalda, G., 2008. Insights on the March 1998 eruption at Piton de la Fournaise volcano

(La Reunion) from microgravity monitoring. *J. geophys. Res.*—*Solid Earth*, **113**, B5.

* Bower, D.R. & Courtier, N., 1998. Precipitation effects on gravity measurements at the Canadian absolute gravity site. *Phys. Earth planet. Inter.*, **106**, 353–369.

* Bruxelles, L., 2001a. Dépôts et altérites des plateaux du Larzac central: causes de l'Hospitalet et de Campestre (Aveyron, Gard, Hérault). Evolution morphogénique, conséquences géologiques et implications pour l'aménagement, *Ph.D thesis*. Université de Provence, Aix-en-Provence, 266 pp.

* Bruxelles, L., 2001b. Reconstitution morphologique du Causse du Larzac (Larzac central, Aveyron, France) Rôle des formations superficielles dans la morphogenèse karstique, *Karstologia*, **38**, 25–40.

* Budetta, G. & Carbone, D., 1997. Potential application of the Scintrex CG-3M gravimeter for monitoring volcanic activity; results of field trials on Mt Etna, Sicily, *J. Volcan. Geotherm. Res.*, **76**, 199–214.

* Burbey, T.J., 2003. Use of time-subsidence data during pumping to characterize specific storage and hydraulic conductivity of semi-confining units, *J. Hydrol.*, **281**, 3–22.

* Dal Moro, G. & Zadro, M., 1998. Subsurface deformations induced by rainfall and atmospheric pressure: tilt/strain measurements in the NE-Italy seismic area, *Earth planet. Sci. Lett.*, **164**, 193–203.

* Debglia, N. & Dupont, F., 2002. Some critical factors for engineering and environmental microgravity investigations, *J. applied Geophysics*, **50**(4), 435–454.

* Farr, T.G. & Kobrick, M., 2000. Shuttle Radar Topography Mission produces a wealth of data, *Am. geophys. Un., EOS*, **81**, 583–585.

* Goodkind, J.M., 1999. The superconducting gravimeter, *Rev. Sci. Instrum.*, **70**(11), 4131–4152.

* Harnisch, G. & Harnisch, M., 2006. Hydrological influences in long gravitational data series, *J. Geodyn.*, **41**, 276–287.

* Hinderer, J. & Crossley, D., 2000. Time variations in gravity and inferences of the Earth's structure and dynamics, *Surv. Geophys.*, **21**, 1–45.

* Hussain, A., Walach, G. & Weber, F., 1981. Underground gravity survey in Alpine regions, *Geophys. Prospect.*, **29**(3), 407–425.

* Imanishi, Y., Kokubo, K. & Tatehata, H., 2006. Effect of underground water on gravity observation at Matsushiro, Japan, *J. Geodyn.*, **41**, 221–226.

* Jacob, T. et al., 2008. Absolute gravity monitoring of water storage variation in a karst aquifer on the Larzac plateau (Southern France). *J. Hydrol.*, **359**(1–2), 105–117, doi:10.1016/j.jhydrol.2008.06.020.

* Jousset, P., Dwipa, S., Beauducel, F., Duquesnoy, T. & Diament, M., 2000. Temporal gravity at Merapi during the 1993–1995 crisis; an insight into the dynamical behaviour of volcanoes, *J. Volcanol. Geotherm. Res.*, **100**(1–4), 289–320.

* Kasameyer, P. & Hearst, J., 1988. Borehole gravity measurements in the Salton Sea scientific drilling project, *J. Geophys. Res.*, **93**(B11), 13 037–13 045.

* King, N.E. et al., 2007. Space geodetic observation of expansion of the San Gabriel Valley, California, aquifer system, during heavy rainfall in winter 2004–2005, *J. geophys. Res.*, **112**, B03409.

* Klimchouk, A., 2004. Towards defining, delimiting and classifying epikarst: its origin, processes and variants of geomorphic evolution, *Epikarst*, **9**(1), 23–25, (in *Proceedings of the Symposium held 2003 October 1–4*, Sheperdstown, West Virginia, USA, Karst Water Institute special publication).

* Kroner, C., 2001. Hydrological effects on gravity data of the geodynamic observatory Moxa, *J. Geod. Soc. Japan*, **47**(1), 353–358.

* Kümpel, H.-J., Lehmann, K., Fabian, K. & Mentes, G., 2001. Point stability at shallow depths: experience from tilt measurements in the Lower Rhine Embayment, Germany, and implications for high-resolution GPS and gravity recordings, *Geophys. J. Int.*, **146**, 699–713.

* LaFehr, T.R., 1983. Rock density from borehole gravity surveys, *Geo-physics*, **48**(3), 341–356.

* Lyard, F., Lefevre, F., Letellier, T. & Francis, O., 2006. Modelling the global ocean tides: modern insights from FES2004, *Ocean Dyn.*, **56**(5–6), 394–415.

* Mangin, A., 1975. Contribution à l'étude hydrodynamique des aquifères karstiques, *Ph.D thesis*. Université de Dijon, 124 pp.
Martelet, G., Debeglia, N. & Truffert, C., 2002. Homogénéisation et validation des corrections de terrain gravimétriques jusqu'à la distance de 167 km sur l'ensemble de la France, *C. R. Geosci.*, **334**, 449-454.

McCulloh, T.H., 1965. A confirmation by gravity measurements of an underground density profile on core densities, *Geophysics*, **30**, 1108-1132.

Merlet, S., Kopaev, A., Diament, M., Geneves, G., Landragin, A. & Pereira Dos Santos, F., 2008. Micro-gravity investigations for the LNE watt balance project, *Metrologia*, **45**, 265-274, doi:10.1088/0026-1394/45/3/002.

Naujoks, M. et al., 2007. Detection of small hydrological variations in gravity by repeated observations with relative gravimeters, *J. Geod.*, **82**, 543-553, doi:510.1007/s00190-007-0202-9.

Naujoks, M., Weise, A., Kroner, C. & Jahr, T., 2007. Detection of small hydrological variations in gravity by repeated observations with relative gravimeters, *J. Geod.*, doi:10.1007/s00190-007-0202-9.

Niebauer, T.M., Sasagawa, G.S., Faller, J.E., Hilt, R. & Klopping, F., 1995. A new generation of absolute gravimeters, *Metrologia*, **32**, 159-180.

Nowel, D.A.G., 1999. Gravity terrain corrections; an overview, *J. appl. Geophys.*, **42**(2), 117-137.

Okabe, 1979. Analytical expression for gravity anomalies due to homoge-
neous polyhedral bodies and translations into magnetic anomalies, *Geo-
physics*, **44**(4), 730–741.

Penman, H.L., 1948. Natural evaporation from open water, bare soil and grass, *Proc. R. Soc. Lond.*, **193**, 120-145.

Perrin, J., Jeannin, P.-Y. & Zwahlen, F., 2003. Epikarst storage in a karst aquifer: a conceptual model based on isotopic data, Milandre test site, Switzerland, *J. Hydrol.*, **279**, 106–124.

Plagnes, V., 1997. Structure et fonctionnement des aquifères karstiques. Cartérisation par la chimie des eaux. *Ph. D thesis*. University of Montpellier II, Document du BRGM n°294, 372 pp.

Ricard, J. & Bakalowicz, M., 1996. Connaissance, aménagement et pro-

tection des ressources en eau du Larzac septentrional, Aveyron (France),
Report R38953, BRGM, Orléans.

Scintrex, 2006. *CG-5 Scintrex autograv system Operation Manual*, Scintrex Limited, Concord, Ontario.

Scintrex limited, 2007. *CG5 Scintrex autograv system Operation Manual*, Concord, Ontario.

Seo, K.-W., Wilson, C.R., Famiglietti, J.S., Chen, J.L. & Rodell, M., 2006. Terrestrial water mass load changes from Gravity Recovery and Climate Experiment (GRACE). *Water Resour. Res.*, **42**, W05417, doi:10.1029/2005WR004255.

Takemoto, S. et al., 2002. Effect of groundwater changes on SG observations in Kyoto and Bandung. *Bulletin d'Informations des Marées Terrestres*: 10 839–10 848.

Tamura, Y., 1987. A harmonic development of the tide generating potential. *Bull. d'Inf. Marées Terr., 99*.

Turc, L., 1961. Evaluation des besoins en eau d'irrigation, évapotranspiration potentielle. *Ann. Agron.*, **12**(1), 13-49.

Van Camp, M. et al., 2006. Hydrogeological investigations at the Mem-
bach station, Belgium, and application to correct long periodic gravity
variations. *J. geophys. Res.*, **111**, B10403.

Viterbo, P. & Beljaars, A.C.M., 1995. An improved land surface parame-
terization scheme in the ECMWF model and its validation. *J. Clim.*, **8**,
2716–2748.

Wenzel, H.-G., 1996. The Nanogal software: earth tide data processing package ETERNA 3.30, *Bulletin d'Informations des Marées Terrestres*, **124**, 9425-9439.

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

Williams, P.W., 1985. Subcutaneous hydrology and the development of do-
line and cockpit karst, *Zeitschrift für Geomorphologie*, **29**, 463–482.

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