Sandeep Kumar

Séverine Rosat
email: severine.rosat@unistra.fr

Jacques Hinderer

Maxime Mouyen

Jean-Paul Boy

M Israil

J-P Boy

Delineation of aquifer boundary by two vertical superconducting gravimeters

Keywords: hydro-gravimetry, superconducting gravimeters, hydrogeology, gravity timeseries, karst aquifer

Introduction

Keeping in mind global warming and growing water demand, the importance of characterizing and monitoring the hydrological cycle has rapidly increased. The scientific community is continuously developing more robust monitoring techniques in order to assess the depletion of water resources. Karst systems represent an essential resource of fresh water but are heterogeneous matrices with fractures and open voids revealing specific hydrodynamic behaviours. Classical hydrogeological methods (e.g., piezometers) provide direct measurements of aquifer properties, but they are invasive and characterize a local region of a larger heterogeneous karstic aquifer system. Other methods, like flow rate and chemical measurements at the spring, provide an overall hydrologic information at the basin scale. Noninvasive geophysical methods appear relevant for investigating the spatial variability of karst systems covering a few square kilometers [START_REF] Chalikakis | Contribution of geophysical methods to karst-system exploration: an overview[END_REF]. Among the geophysical methods, gravimetry provides a direct quantification of the water-mass fluctuations at the catchment scale [START_REF] Kroner | Hydrological signals in gravity -foe or friend? Dynamic Planet[END_REF][START_REF] Hasan | Evaluating catchment-scale hydrological modelling by means of terrestrial gravity observations[END_REF][START_REF] Jacob | Time-lapse microgravity surveys reveal water storage heterogeneity of a karst aquifer[END_REF][START_REF] Hector | Gravity effect of water storage changes in a weathered hard-rock aquifer in West Africa: results from joint absolute gravity, hydrological monitoring and geophysical prospection[END_REF][START_REF] Hemmings | Investigating hydrological contributions to volcano monitoring signals: A time-lapse gravity example[END_REF][START_REF] Imanishi | Effect of underground water on gravity observation Matsushiro, Japan[END_REF][START_REF] Jacob | Timelapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF]Pool & Eychaner, 1995;[START_REF] Camp | Hydrogeological investigations at the Membach station, Belgium, and application to correct long periodic gravity variations[END_REF][START_REF] Wilson | Field Test of the Superconducting Gravimeter as a Hydrologic Sensor[END_REF]. Besides, hydrological models require calibration that can be achieved by assimilating gravity measurements [START_REF] Naujoks | Evaluating local hydrological modelling by temporal gravity observations and a gravimetric three-dimensional model[END_REF][START_REF] Creutzfeldt | Total water storage dynamics in response to climate variability and extremes: Inference from long-term terrestrial gravity measurement[END_REF][START_REF] Chaffaut | Continuous Monitoring with a Superconducting Gravimeter As a Proxy for Water Storage Changes in a Mountain Catchment[END_REF]. Moreover, gravity gradiometry or multiple gravimeters across the local area can further improve the assessment of temporal water mass changes [START_REF] Naujoks | Detection of small hydrological variations in gravity by repeated observations with relative gravimeters[END_REF][START_REF] Kennedy | Direct measurement of subsurface mass change using the variable baseline gravity gradient method[END_REF], by mitigating the common non-uniqueness issues in singlegravimeter surveys.

The low noise underground laboratory, the so-called LSBB ("Laboratoire Souterrain à Bas Bruit") of Rustrel (France) is located in a karstic massif at the southern edge of the Fontaine de Vaucluse impluvium. The LSBB constitutes a typical experimental site for underground karstic water-fractured hydrosystems [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF]. The LSBB site is well-equipped with various instrumentations, as for our concern, with two Superconducting Gravimeters (SGs): an observatory SG, the iOSG24, recording data since September 2015 and an iGrav instrument [START_REF] Warburton | Initial results with the new GWR iGrav TM superconducting gravity meter[END_REF], the iGrav31 continuously recording since February 2019. SGs are relative instruments that record temporal gravity variations at the Earth's surface and have a wide range of geodynamics applications [START_REF] Hinderer | Superconducting Gravimetry in Treatise on Geophysics[END_REF]. SGs have precision and longterm stability appropriate for monitoring water-storage changes (e.g., [START_REF] Kroner | Sensitivity of superconducting gravimeters in central Europe on variations in regional river and drainage basins[END_REF][START_REF] Hector | Hydro-gravimetry in West-Africa: first results from the Djougou (Benin) superconducting gravimeter[END_REF][START_REF] Fores | Assessing the precision of the iGrav superconducting gravimeter for hydrological models and karstic hydrological process identification[END_REF][START_REF] Güntner | Landscape-scale water balance monitoring with an iGrav superconducting gravimeter in a field enclosure[END_REF][START_REF] Chaffaut | Continuous Monitoring with a Superconducting Gravimeter As a Proxy for Water Storage Changes in a Mountain Catchment[END_REF]. The set-up of both SGs at the LSBB site is such that iGrav31 is 520 m nearly vertically above the iOSG24, with horizontal offset of about 90 m. This specific configuration of SGs is unique in the world and constrains the local water masses influencing SG measurements [START_REF] Mouyen | Monitoring groundwater redistribution in a karst aquifer using a superconducting gravimeter[END_REF]. This is the first time that hydrological processes are investigated when the baseline configuration of two SGs is vertical. Previously, several significant findings were made based on signals from horizontal baselines [START_REF] Naujoks | Detection of small hydrological variations in gravity by repeated observations with relative gravimeters[END_REF][START_REF] Kennedy | Direct measurement of subsurface mass change using the variable baseline gravity gradient method[END_REF]. Here, the vertical configuration of two gravimeters has the advantage to increase the sensitivity to water mass changes in the area in between, by a factor two when considering a simple Bouguer plate approximation. Indeed, the gravity change associated with an equivalent liquid water thickness SGs is that a single gravimeter is insensitive to water mass changes at the same height of the instrument while the vertical configuration of two SGs solves this issue (see e.g., Fig. S2 in Carrière et al., 2022).

h is ∆g = 2 ℎ
In this paper, we first describe the hydrogeological setting of the LSBB site and the unique configuration of two vertically arranged SGs. We then explain the SG data processing used to retrieve gravity time-residuals compared with two global hydrological models. A hydro-gravimetric modelling based on integration of rectangular prisms is used to compute hydrological admittances (i.e., gravity response to a 10-cm water layer thickness). Then, we use the obtained admittances to calculate the gravity effects induced by a homogeneous water layer spread over the topography with water head changes inferred from the global hydrological model. Finally, the sensitivity of the differential gravity signal to the extent and depth of the main water storage change is discussed.

Experimental Site: Geological and Hydrological Setting

Fontaine de Vaucluse Hydrosystem

The Fontaine de Vaucluse (FdV) is a karst hydrosystem located in the southeast of France.

Groundwater discharges from the system at the FdV spring with an average flow-rate of 23 m 3 /s [START_REF] Cognard-Plancq | Historical monthly rainfall-runoff database on Fontaine de Vaucluse karst system: review and lessons[END_REF]. The impluvium spreads over an area of 1,115 km 2 , and the average thickness of its unsaturated zone (UZ) is around 800 m [START_REF] Puig | Le système karstique de la Fontaine de Vaucluse [The karst system of the Fontaine de Vaucluse[END_REF]. This karst system consists of massive and continuous lower Cretaceous limestone and its base extends to a depth of 1500 m [START_REF] Masse | Contribution à l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF].

The reef limestone part of this carbonate platform has thickness of around 450 m in the study area. The Urgonian facies present in the reef limestone has high total matrix porosity up to 25-30 percent among limestone and is present in half of the FdV [START_REF] Masse | Contribution à l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF][START_REF] Masse | Drowning discontinuities and stratigraphic correlation in platform carbonates: the Late Barremian-Early Aptian record of southeast France[END_REF]. The stable low flow rate, the thick UZ, and some small permanent outlets in the limestone plateau indicate that the low flow discharge occurs mainly in the UZ [START_REF] Garry | Etude des processus d'écoulement de la zone non saturée pour la modélisation des aquifères karstiques[END_REF].

The base of this karst system lies on Valanginian and -Upper Hauterivian impermeable marls.

The presence of dolines and dry valleys cutting through the plateau indicates strong karstification of the FdV hydrosystem due to limestones dissolution [START_REF] Blavoux | The karst system of the Fontaine de Vaucluse (Southeastern France)[END_REF].

Low Noise Underground Laboratory

The LSBB is a 3.8 km long, nearly horizontal underground tunnel in the north of Rustrel village, France (Fig. 1). It is a former component of the French nuclear missile system, which has been turned into a cross disciplinary underground research laboratory. The LSBB is a unique low noise underground set-up because of its initial conception, aiming to overcome nuclear blast, and its location far from anthropic noise [START_REF] Waysand | The Laboratoire Souterrain Bas Bruit (lsbb) In Rustrel-pays D'apt (france): A Unique Opportunity For Low-noise Underground Science[END_REF]. This tunnel passes through the massif and randomly intersects karstic features. It also goes through some channels in the UZ. The maximum thickness of rock cover over the tunnel reaches 519 m. The saturated zone (SZ) in this catchment lies approximately 400 m beneath the LSBB tunnel. The fracturing and karstification are very diverse throughout the LSBB gallery. [START_REF] Garry | Etude des processus d'écoulement de la zone non saturée pour la modélisation des aquifères karstiques[END_REF] carried out several hydro-chemical and hydrodynamic investigations in this site and developed different hydrodynamic models of these flow points. [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF] demonstrated the effects of the water in the rock using an integrated hydro-geophysical approach.

SG Observations and Global Hydrological Models

Gravity data processing

The iOSG24 has been recording data since September 2015 [START_REF] Rosat | A two-year analysis of the iOSG24 superconducting gravimeter at the low noise underground laboratory (LSBB URL) of Rustrel, France: environmental noise estimate[END_REF], whereas the iGrav31 was set up more recently and started recording data since May 2019. The recorded SG data correspond to the feedback voltage that balances the position of the levitating sphere at a 1 Hz sampling frequency. The pre-processing of SG data includes the change of voltage to acceleration through a calibration factor using FG5 parallel measurements, signal processing, and drift correction. Many studies related to SG calibration using parallel absolute gravity measurements can be cited (e.g., [START_REF] Francis | Calibration of the C021 superconducting gravimeter in Membach (Belgium) using 47 days of absolute gravity measurements[END_REF][START_REF] Imanishi | Calibration of the superconducting gravimeter T011 by parallel observation with the absolute gravimeterFG5#210 -a Bayesian approach[END_REF]. SG calibrations are known to be very stable over time (e.g., [START_REF] Goodkind | The Superconducting Gravimeters. principles of operation, current performance, and future prospects[END_REF][START_REF] Calvo | Time stability of spring and superconducting gravimeters through the analysis of very long gravity record[END_REF]. Two absolute gravity campaigns were carried out during the considered time epoch, in October 2019 and in September 2020 for both instruments. The obtained calibrations factors are -451 ± 3 nm/s²/V for iOSG24 [START_REF] Rosat | A two-year analysis of the iOSG24 superconducting gravimeter at the low noise underground laboratory (LSBB URL) of Rustrel, France: environmental noise estimate[END_REF] and -851 ± 7 nm/s²/V for the iGrav31, that is an accuracy better than 1 %. Please note that the uncertainty of the scale factor for iGrav31 is larger because of its installation at the surface, where environmental noise is larger than inside the tunnel.

The instrumental drift was removed using a linear trend since we considered data recorded several months after the installation, when the drift has become linear. Discussions about SG drift were published in Van Camp and Francis (1997) and more recently in [START_REF] Hinderer | Intercomparing Superconducting Gravimeter Records in a Dense Meter-Scale Network at the J9 Gravimetric Observatory of Strasbourg, France[END_REF]. Contributions from signals that are not due to local hydrology are finally removed.

These signals are listed below:

• Solid and oceanic tides: they were both removed using a local tidal model resulting from a least-squares adjustment to SG data so that it contains both solid and oceanic tides. This tidal analysis was performed with ETERNA3.4 [START_REF] Wenzel | The nanogal software: earth tide data processing package[END_REF] software based on [START_REF] Hartmann | The HW95 tidal potential catalog[END_REF] potential catalogue (as in [START_REF] Rosat | A two-year analysis of the iOSG24 superconducting gravimeter at the low noise underground laboratory (LSBB URL) of Rustrel, France: environmental noise estimate[END_REF]) and for which semi-annual and annual tides were set to the nominal 1.16 gravimetric factors (else part of the seasonal hydrological signal would be artificially reduced).

• Atmospheric pressure masses: these effects were calculated using atmospheric pressure data from the second Modern-Era Retrospective analysis for Research and Applications (MERRA2) assuming an inverted barometer response of the oceans. This atmospheric loading (including non-tidal oceanic loading) is calculated using a Green's function formalism as in [START_REF] Boy | Reduction of surface gravity data from global atmospheric pressure loading[END_REF]; the atmospheric loading consists of two contributions: a local contribution corresponding to an integration within 0.10° around the gravimeter and converted to a gravity signal using a barometric value of -2.21047 nm/s 2 /hPa , a non-local contribution resulting from the Green function loading integrating deformation beyond the 0.10° radius.

• Polar motion: the gravity effect for pole tides is calculated using Earth's rotation parameters taken from IERS (ftp://hpiers.obspm.fr/iers/eop/eopc04/) using elastic gravimetric factor as in [START_REF] Hinderer | Superconducting Gravimetry in Treatise on Geophysics[END_REF].

• Non-local hydrological loading: the "non-local" contribution of hydrology is computed using the Green's function [START_REF] Farrell | Deformation of the Earth by surface loads[END_REF] as in [START_REF] Boy | Study of the seasonal gravity signal in superconducting gravimeter data[END_REF]. The non-local part corresponds to a region outside 20-30 km radius from the station, and its amplitudes are about 15 percent of that of the local signal. Values are calculated using either MERRA2 (Reichle et al., 2017a) or ERA5 [START_REF] Hersbach | The ERA5 global reanalysis[END_REF] hydrological models (soil moisture, and snow components).

The hence obtained gravity residuals are plotted on Fig. 2.In the following, we assume that air mass changes were perfectly reduced from the gravity time-series, so that the remaining gravity signals can be interpreted directly in terms of water storage changes.

Comparison of gravity data with global hydrological models

In this section, we compare gravity residuals of the iOSG24 (∆g(iOSG24)) and iGrav31

(∆g(iGrav31)) (Fig. 3a) between May 2019 and June 2020 with corresponding local contributions from MERRA2 and ERA5 hydrological models. The local hydrological effect was calculated using a Bouguer plate approximation with a nominal admittance value of ±4.2677 nm.s -2 / cm (the sign depends on the location of the station with respect to the ground:

it is minus when the station is underground). We recall that gravity residuals should only contain local hydrological effects.

The MERRA2 reanalysis is a data product of the National Aeronautics and Space Administration (NASA) Global Modelling and Assimilation Office [START_REF] Gelaro | The modern-era retrospective analysis for research and applications, version-2 (MERRA-2)[END_REF]. The MERRA2 contains hourly-sampled data from 1980 to the present with a horizontal resolution of approximately 50 km. In MERRA2, the land surface water budget was estimated using observation-based precipitation (Reichle et al. 2017b). This analysis is calculated on a latitudelongitude grid at the exact spatial resolution as the atmospheric model.

ERA5, a fifth-generation reanalysis data product, is the most recent release by the European Centre for Medium-Range Weather Forecasts (ECMWF). Several significant improvements were incorporated compared to its previous ERA-Interim [START_REF] Hersbach | The ERA5 global reanalysis[END_REF]. This is a 1hourly sampled data and was computed using a more advanced 4D VAR assimilation algorithm. The horizontal resolution of ERA5 is approximately 30 km. The ERA5 contains data from 1979 to the present.

The temperature and precipitation biases of ERA5 are found to be more accurate relative to MERRA2 for hydrological modelling. Local hydrological contribution to gravity loading at the LSBB site is slightly smaller for MERRA2 than for ERA5 model (Fig. 3).

In Fig. 3(a), we plot gravity residuals recorded by both SGs with corresponding local effects from both hydrological models ERA5 and MERRA2. We see that the ERA5 and MERRA2 models predict well the main features of gravity changes at the surface (iGrav31) while for iOSG24, after taking into account the sign inversion, the predicted variability is overestimated.

Particularly between November 2019 and January 2020, the predicted gravity fluctuations reach ~75 nm/s² while the observed ones are less than 20 nm/s² underground and larger than 100 nm/s² at the surface.

In Fig. 3 SGs. The overall agreement between the hydrological models and the differential gravity timeseries is quantified by a RMS of 44.3 nm/s² while for each time series taken separately the RMS is of 72.4 nm/s² for iOSG24 and 28.5 nm/s² for iGrav31. This suggests that the hydrological models better represent the groundwater storage variation from the surface down to the UZ when the differential ∆g(iGrav31)-∆g(iOSG24) gravity signal is considered. Indeed, under Bouguer plate approximation, the contribution of the water below both gravimeters cancels out when we compute the differential gravity. Since the Bouguer plate is insensitive to depth, the water-mass contribution is the same on both gravimeters.

However, misfits between local hydrological models and differential gravity time-series are still visible, highlighting the effects of lateral flows and the role of topography, which mitigate the Bouguer plate approximation and the adequateness of the global hydrological models in such a mountainous karstic system.

Subsequently, we zoomed in the residual gravity time-series from Dec 1, 2019 to Jan 10, 2020 (Fig. 4) and plotted the observations-based precipitation as provided within MERRA2 products (Reichle et al. 2017b). We marked five major timings (numbered 1 to 5 in Fig. 4) where the residual gravity trends of both SGs change abruptly. An important criterion for our analysis is whether the gravity variations measured at each gravimeter at a given time, have the same sign or are opposite. If they are opposite, it means that water storage variation occurred at a depth between the two gravimeters. If they have the same sign, then the water storage variation occurred below both gravimeters. For timing 1 corresponding to a large rain event, ∆g(iGrav31)

increases significantly whereas ∆g(iOSG24) decreases slightly, because of the SGs configuration. However, for timing 2, which does not correspond to any precipitation event, the gravity trends are reversed indicating that all the surface water was seeped into the zone below both the SGs. This transition from rainfall to infiltration has taken approximately 6 days, which is very quick in terms of hydrological events. This shows that the UZ is highly porous and permeable. Timing 3 illustrates another rain event smaller than event 1, with an amount of gravity increase seen by iGrav31 smaller than the amount of gravity decrease seen by iOSG24.

Similarly, timing 4 marks another precipitation occurrence whereas timing 5 is for infiltration.

We also see that the slope of ∆g(iGrav31) during infiltration periods is much larger than observed water storage changes (∆g(iOSG24)), indicating faster water mass transfers as seen by gravimetry.

Forward Gravity Modelling

Methodology

Here we present a rectangular prism method developed to compute time-lapse gravity signals from the hydrological models. In this forward problem, the hydrological model cells are assumed to be incremental prismatic mass storages, and corresponding gravity effects are computed by integrating over all cells within a given a priori radius. However, we consider a homogenous water layer change spread over the topography, so our modelling is not 3D since we do not consider the possibility to have different water head changes in different cells. We used the RGE ALTI digital elevation model (DEM) provided by the National Institute of Geographical and Forestry Information in France (IGN-F) on a regular grid with 1 m resolution but decimated to 5 m in order to save computation time. This method comprises three formulas:

the prism formula, the MacMillan formula, and the point-mass formula. These three equations were derived from Newton's law of gravity with different approximations depending on the distance to the mass change (see [START_REF] Leiriao | Calculation of the temporal gravity variation from spatially variable water storage change in soils and aquifers[END_REF]. They are used to calculate gravity changes due to an unconfined aquifer with varying hydraulic heads. From previous studies, we know that a single prism causes gravity changes only if its vertical length varies. The DEM (Fig. 5) is discretized into rectangular prisms of sizes 5 m x 5 m x 5 m to incorporate the site topography. We consider the lateral offset between both SGs in the modelling by computing the forward gravity effects at the exact locations of each SG.

The gravity variation is then calculated using all three formulas depending upon the radial distance to the location of the SG [START_REF] Forsberg | A study of Terrain Reductions Density Anomalies and Geophysical Inversion Methods in Gravity Field Modelling[END_REF]. The entire computation process involves four steps comprising applying the prism formulas and summing up the obtained results. The Forsberg formula is used when the SG is close to the mass change:

∆ = log + + log + -tan !" ,
where the vertical lines indicate that the integration of the mass element is performed in the x,y and z directions between the two points of coordinates that define the prism. x, y and z are the components of the instrument-prism vector (see [START_REF] Leiriao | Calculation of the temporal gravity variation from spatially variable water storage change in soils and aquifers[END_REF] for further details).

For intermediate radial distances, the spatial detail in the discretization decreases. This leads to the application of the MacMillan formula [START_REF] Macmillan | The Theory of Potential[END_REF]:

∆ = - ∆ ∆ ∆ #-$ - 5 24 & ' + ( ' + ) ' * + 1 12 )
, -

For farther radial distances, it is assumed that the whole mass is in the centre of the cell.

Consequently, the 3D point-mass formula is used in the remaining model grid:

. = -∆ . . ℎ -/ ' + ' + ℎ -/ ' $/'

Finally, the total gravity change is a cumulative sum of the contributions from applications of all three formulas (Forsberg, MacMillan, and point-mass).

The formulation of the rectangular prism method is sketched in Fig. 6 where d=r i.e., the distance between the prism and the instrument, and ∆h are hydraulic heads representing the vertical changes in water layer thickness. Fig. 6 is basically a schematic design indicating manifestations of the Newton's forward gravity equation.

The total gravity response for a hydro-system is the cumulative effect of all the prisms.

Therefore, it is crucial to decide which of the three formulas to use when a prism is located randomly away from the SGs location. The normalized switching factor has been defined using the size of the prism and its relative distance to include all three formulas. The switching factor f 2 is defined below:

1 ' = ' . '
Where ' = ' + ' + ' is the distance between the prism and the instrument, . ' = . ' + . ' + . ' is the square of length of the diagonal line of a prism, and other symbols are described in [START_REF] Leiriao | Calculation of the temporal gravity variation from spatially variable water storage change in soils and aquifers[END_REF]. Here in this work, the prism formula is used for f 2 below 25, the MacMillan formula is used for f 2 between 25 and 36, and the point-mass formula for f 2 above 36 as done in [START_REF] Leiriao | Calculation of the temporal gravity variation from spatially variable water storage change in soils and aquifers[END_REF] and [START_REF] Chaffaut | Continuous Monitoring with a Superconducting Gravimeter As a Proxy for Water Storage Changes in a Mountain Catchment[END_REF].

Hydrological Admittances

The hydrological admittance usually refers to an increment in the observed gravitational acceleration due to an infinitely extended horizontal layer of water (infinite Bouguer plate),

with the SG at any distance above the layer. It is approximately equal to 4.2677 nms ' cm for the density of water. For practical reasons, if the distance between the sensor and the layer is less than or equal to one-tenth of the horizontal size, the layer is considered infinite. The admittance value enables quantifying the contribution from local hydrology and is obtained using a simple Bouguer approximation in its infinite limit. Using the prism method described in section 4.1, we computed admittances for 10-cm thick layers of water that follow the topographic model. These admittances are 2D transfer functions that depend on the radius of integration and depth of the water layer. In the infinite limit for the radius of integration, they correspond to the classical admittance. In the following, we denote this 2D transfer function as "admittance".

The obtained hydrological admittances (Fig. 7) show that the sensitivity of iOSG24 to waterstorage changes tends to a Bouguer plate with the asymptotic value of -3.4 nm s ' cm for radii larger than 4 km, while for iGrav31, the admittance tends to an asymptotic 3.86 nm s ' cm for radii larger than 6 km. Please note that for depths of 450 m and larger, the admittance computed at the iOSG24 location becomes positive since, upon incorporating the topography height at this depth, the water layer is located below the SG.

Relative contribution of Unsaturated and Saturated Zones

In the previous section, we performed the integration for all depths independently to the fact the water layer was located in the UZ or SZ. However, SG observations integrate gravity signals from changes in water distribution in both SZ and UZ, in particular with different sensitivities for iGrav31 and iOSG24 due to their respective distances to the SZ. In this section, we want to check the importance of the SZ w.r.t. the UZ to gravity changes. For that, we divided the forward gravity modelling in two steps: one for the integration in the UZ and one for the integration in the SZ. The temporal gravity variation in the UZ is computed by dividing the soil into several layers between 0 and 900 m depth, and then contributions from each layer are summed up together. A similar integration is performed for the SZ by summing up the contributions of water layers located between 900 and 1500 m deep.

The time-variability of the hydraulic head is simulated by using the local gravimetric effect computed for ERA5 model divided by the nominal admittance of 4.2677 nm s² cm -1 to convert it into water-head changes. The obtained water layer height changes are then multiplied by the respective hydrological 2D admittances as computed in the section 4.1 for various integration radii. Please remember that we call admittance not the limit for an infinite Bouguer plate (large radii) but the Newtonian gravity effect integrated along various radii and depths. The admittance is then a 2D function of radius and depth and can be directly multiplied with the time-varying water-level changes to simulate the time-dependent gravity effect. An integration radius of 150 m predicts gravity changes that best fit ∆g(iGrav31) (Fig. 8). For iOSG24, finding a best-fit solution is more complicated since the influence of the SZ becomes larger and the UZ is located partly below and partly above the instrument. A specific inverse problem could be developed but is beyond the scope of this paper, which is to focus on the use of the gravity difference ∆g(iGrav31)-∆g(iOSG24), indeed, to mitigate the uncertainty linked to the use of a unique gravimeter.

The water storage within the SZ in karst conduits plays a significant role in regulating flow in hydro-systems [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF]. We see that the water mass fluctuations in the SZ are, as expected, much less than the effects of the UZ on the iGrav31, 96% of the gravity contribution coming from the UZ (Fig. 8). As seen in most hydro-systems, main contributions to the gravity signals are from the time-varying soil moisture content in the UZ, between both SGs. This result matches with Mangin's schema [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], which claimed the presence of an epikarst zone in this site. This zone lies at the top of the UZ, especially above deep infiltration zone. Therefore, the simulated gravity signals in the UZ are corresponding to the epikarst and the infiltration zone. The hydrological models are also essentially informing on the epikarst in this area, given the agreement obtained with the differential gravity signal ∆g(iGrav31)-∆g(iOSG24) (section 3.2).

Sensitivity of the differential gravity

Our aim here is to delineate major water storage in the unsaturated zone using differential gravity data. The underlined idea is to find the zone where the water-mass changes best explain observed gravity fluctuations. Therefore, we computed the differential gravity effect ∆g(iGrav31)-∆g(iOSG24) for various depths and radii of integration (Fig. 9(a)) that we compare with the observed ∆g(iGrav31)-∆g(iOSG24). The RMS of the difference between predicted and observed ∆g(iGrav31)-∆g(iOSG24) shows the depths and integration radii for which the RMS is minimum (dark blue zone in Fig. 9(b)). Some examples for specific radii of integration and depths of the water layer are illustrated in Fig. 10. We can get a few major findings upon thoroughly analysing Fig. 9 and Fig. 10. Firstly, the largest RMS values above iOSG24 suggests that residual data of both SGs have significant differences. This is mostly because the water present in this zone (i.e., from surface to 500 m depth) lies below the iGrav31 and above the iOSG24. The differential gravity signal in this upper part of the UZ quickly extends to a Bouguer plate approximation when the radius of integration exceeds 400 m (Fig.

9(a)

) that is closer than for an individual gravimeter as inferred from the computed admittances (Fig. 7). At depths between 0 and 400 m (above iOSG24), the admittance ranges from 2 to 4 nms 2 cm -1 indicating that the surface iGrav31 fluctuations dominate in this zone. The largest RMS difference of observed differential gravity with simulated one is located in the zone located at depths between 300 and 500 m and radii between 150 and 500 m (Fig. 9(b)) where the admittance is closer to zero (Fig. 9(a)) indicating a lack of sensitivity of the differential gravity to water-storage changes in this zone. We indeed see that an integration radius of 200 m at a depth of 300 m is not sufficient to explain most of the observed ∆g(iGrav31)-∆g(iOSG24) (Fig. 10(a)). A radius between 1000 and 1500 m would better explain the amplitude of the gravity fluctuations. When considering a water-layer below iOSG24, at a depth of 600 m (Fig. 10(b)), the integration radius needs to be larger than 1500 m to fit the amplitude of the observed differential gravity fluctuations. When the integration radius is too small (less than 600 m), whatever the depth of the water layer, the amplitude of the simulated ∆g(iGrav31)-∆g(iOSG24) is too small to explain the observed one (Fig. 10(a) and Fig. 10(c)).

When the depth of the water layer is close to the depth of iOSG24 (between 400 and 500 m), the predicted differential gravity signal is close to zero (Fig. 10(d)) when the integration radius is large enough (here 1000 m). This is because we are closer to iOSG24 and there is a cancellation of the gravity effect due to the water layer, which is partly above and partly below iOSG24.

Discussion

Comparison of gravity observations with local contributions from global hydrological models

We observe that both local hydrological contributions of ERA5 and MERRA2 show a few trends of overestimating and underestimating gravity effects (Fig. 3(a)) compared to the residual gravity time-series recorded by the iGrav31. This could be due to the following reasons:

i)

The excessive lateral fluxes of groundwater may be inhibiting the vertical infiltration of rainwater into the ground.

ii)

The SGs are located at the south border of the impluvium, and hence, the water mass distributions outside it may partly cancel the gravity effect.

iii)

The time lags between observed rainfalls at the site, residual gravity and the hydrological loading model may be the cause for this misfit.

iv)

The water storage variations given by the hydrological models are not strictly localised at the surface (where the iGrav31 is) but integrate a greater thickness of ground However, when we compare the ERA5 and MERRA2 signals with differential gravity timeseries ∆g(iGrav31)-∆g(iOSG24), we could see that those hydrological models are better mimicking the observed differential gravity fluctuations (Fig. 3b). Due to the above reason, for further investigations, we mostly used differential signals (∆g(iGrav31)-∆g(iOSG24)).

Comparison of surface gravity changes with gravity changes inside the tunnel after rainfalls By comparing the two gravity time series at the surface and inside the tunnel, we can estimate the transfer time of the water from the surface to the SZ after major rain events. Indeed, when we have a sign concordance between the two gravimetric series, after a rain event, that implies that the water is below both the instruments. For instance, during the period between Dec 3, 2019, and Dec 22, 2019 (Fig. 4), we can infer that the peak of flow inside the tunnel occurs as approximately 6 days after the major rain as already observed in past studies (e.g. [START_REF] Garry | Etude des processus d'écoulement de la zone non saturée pour la modélisation des aquifères karstiques[END_REF].

Please note that since the amplitude of the gravity effect is different at the surface and at the depth of iOSG24, the differential gravity (Fig. 3(b)) cannot show such a transfer simply as a negative anomaly. A few days of transfer indicates the presence of a high porosity infiltration zone corresponding to the epikarst overlying the infiltration zone. This zone in the LSBB tunnel was already identified as a fast circulation flow point with a dripping within the LSBB tunnel occurring a few days after a rain event [START_REF] Garry | Contribution of artificial galleries to knowledge of karstic system behaviour in addition to natural cavern data[END_REF] but other flow points were identified representative either of a karstic discharge or of a fractured circulation. The spatial sensitivity of both SGs raises the problem of the representatively of such gravity data at the system scale.

Besides, this time is not an infiltration or transfer time since it corresponds to a mass transfer without taking into account the identity of these water masses. In other words, geochemistry measurements are required to identify the origin of this water mass reaching iOSG24 and resulting in a dripping inside the LSBB tunnel (e.g. [START_REF] Blondel | Punctual and continuous estimation of transit time from dissolved organic matter fluorescence properties in karst aquifers, application to groundwaters of 'Fontaine de Vaucluse' experimental basin (SE France)[END_REF].

Sensitivity of the gravity signals to water storage changes

In section 4.2, we have computed the 2D hydrological admittance using a prism method for a uniform water layer of 10 cm as a function of the radius of integration and depth of the water layer (Fig. 9). This modelling shows that the area of largest sensitivity for the differential gravity lies from 450 m to 700 m deep and in the epikarst zone between 0 and 400 m deep approximately (Fig. 9(a)). This result is very crucial in the estimation of the total water storage in the LSBB site. When we further compare the observed and simulated differential gravity time-series for different depths and radii of integration, we can make the following remarks i)

The simulated signals best fit the observed gravity data for a water-layer depth below 500 m that is below iOSG24 (Fig. 9(b) in dark blue and Fig. 10).

ii)

The RMS of the difference between predicted and observed ∆g(iGrav31)-∆g(iOSG24) is largest in the UZ between 400 and 500 m depth and for radii smaller than 500 m (Fig. 10(b)) where the admittance is close to zero (Fig. 10(a)). In this zone, induced gravity effects on iGrav31 and on iOSG24 cancel out.

iii)

We note that as we go deeper below iOSG24, the predicted differential gravity signal fits better the observed one when the integration radius increases (Fig. 9(b)).

Conclusions

In this study, we put complete emphasis on investigating the hydrological processes occurring at the LSBB site. These processes were inferred from the observed differential gravity timeseries from a vertical dipole of two superconducting gravimeters. The observed mass transfer time of 6 days gives an estimate of the vertical hydraulic conductivity of the infiltration zone.

The residual and differential gravity time series of the SG vertical dipole clearly show the groundwater redistribution at seasonal time-scales at the LSBB site. The simulated gravity responses from a global hydrological model indicate that most groundwater is stored in the unsaturated zone between both SGs. This paper also presented the relative contribution from the saturated and unsaturated zones separately. Firstly, most of the contributions in gravity changes are due to water in the unsaturated zone. Secondly, since the saturated zone lies farther away from the SGs, the gravity contributions are very low. We also computed true and asymptotic admittance values using a rectangular prism method matching the topography. We started with simple comparisons of recorded gravity time-series with the gravity effects of global hydrological models and then investigated the following observations further. Finally, we mapped the entire aquifer zone by minimizing the misfit between the observed and simulated differential gravity signals from both saturated and unsaturated zones. We could conclude that a part of the aquifer lies above the iOSG24, whereas the remaining part lies below iOSG24 in the unsaturated zone. Even though we successfully delineated the groundwater boundary in this catchment, further work is needed to improve the simplified hydro-gravimetric model proposed in this study. Mainly, the mask effect due to the shelter at the iGrav31 location and the mask effect due to the tunnel around iOSG24 should be considered in the gravity modelling by removing prism cells around the instruments (e.g., [START_REF] Chaffaut | Continuous Monitoring with a Superconducting Gravimeter As a Proxy for Water Storage Changes in a Mountain Catchment[END_REF]. Mask effects could have important consequences on the parameters of the unsaturated zone, which is shallower.

The main hypothesis made in this study is that all other contributions, particularly air mass changes, were completely and accurately reduced from gravity observations. However, it is well-known that the dynamics of atmospheric mass changes is also complex with turbulent effects, particularly during heavy rain events, that are not considered in the loading models used in this study (see, for example, [START_REF] Neumeyer | Gravity reduction with threedimensional atmospheric pressure data for precise ground gravity measurements[END_REF][START_REF] Gitlein | Modeling of Atmospheric Gravity Effects for High-Precision Observations[END_REF]. Some small atmospheric contribution in the gravity residuals may remain.

We have considered in this paper the problem of determining water layer depth and extent in the unsaturated zone from the vertical gradient of time-varying gravity changes using preimposed water head changes with a uniform layer fitting the topography. Using the hydrogeological knowledge that has been accumulated among the years at the LSBB karstic site by the numerous previous studies, we should in the future extend the prism modelling into 3D by considering heterogeneous water content field [START_REF] Chaffaut | Hybrid Gravimetry to Map Water Storage Dynamics in a Mountain Catchment[END_REF]. Finally, we have shown how the vertical gradient of gravity could infer properties on water storage content within a mountainous karstic environment and the potential improvements that could be done from this initial work.

  where is the volumetric mass of water and G is universal gravitational constant. The gravimeter at the surface measures ∆g(iGrav31) = 2 ℎ while the one below the water layer measures ∆g(iOSG24) = -2 ℎ since the water masses are located above the sensor. The differential gravity measured by these two SGs ∆g(iGrav31)-∆g(iOSG24) is then 4 ℎ, twice the Bouguer plate effect. Another advantage of the vertical configuration of two

  (b), we compare differential gravity time-series (∆g(iGrav31)-∆g(iOSG24))) with twice the corresponding values from the ERA5 and the MERRA2. Note that we doubled ERA5 and MERRA2 local hydrological contribution because, as we explained in the introduction, the vertical configuration of two gravimeters increases the gravitational effect of the water located between the gravimeters by a factor two, when considering a simple Bouguer plate approximation and assuming a homogeneous mass distribution.When the difference ∆g(iGrav31)-∆g(iOSG24) is negative (Fig.3(b)), we can say that most of the water storage change is located below the iOSG24. While when ∆g(iGrav31)-∆g(iOSG24) is positive, most of the water mass change is temporarily stored in the UZ between the two

Figure 1 A

 1 Figure 1 A simplistic model of the Fontaine de Vaucluse hydrosystem with the LSBB site (modified from Carrière et al., 2016) in the south-east of France. The location of the two Superconducting Gravimeters (iGrav31 and iOSG24) is indicated: iOSG24 is located within the LSBB tunnel while iGrav31 is at the surface. The separation between the unsaturated and saturated zones is at around 800 m depth.

Figure 2

 2 Figure 2 Time-varying gravity recorded by iGrav31. (a) Original data showing mainly tidal variations. Gravity residuals after subtraction of (b) tides and polar motion effects. (c) Tides, polar motion and air mass effects (d) Tides, polar motion, atmospheric and non-local hydrological loading effects.

Figure 3

 3 Figure 3 (a) Residual gravity time-series recorded by iGrav31 (black) and by iOSG24 (gray) compared to the gravity effect induced by the local contribution calculated from ERA5 (blue) and MERRA2 (cyan) hydrological models; (b) Differential gravity signal (red) compared to twice the gravity effect computed for iGrav31 corresponding to the local hydrological contributions from ERA5 (blue) and MERRA2 (cyan) global models.

Figure 4

 4 Figure 4 Zoom-in view of the residual gravity time-series recorded by iGrav31 (black) and by iOSG24 (gray) between December 2019 to January 2020 compared with the local gravity effect induced by hydrology computed from ERA5 (blue) and MERRA2 (cyan) global models. Observed total precipitation (mm) and water storage (kg/m²) used in MERRA2 modeling are also plotted in cyan and green.

Figure 5

 5 Figure 5 Digital elevation model of the LSBB site centered on iOSG24 (black dot). The black circle indicates the radius of integration of a supposed water layer around iOSG24.

Figure 6

 6 Figure 6 Schematic view of the parameters used in the forward gravity modelling of the gravity effect due to water storage change in an aquifer.

Figure 7

 7 Figure 7 Gravity effects of the 10-cm uniform water layer w.r.t radii, for various depths (one color per depth). Red boxes mark the asymptotic gravity value for layers above/below the SGs.

Figure 8

 8 Figure 8 Residual gravity time-series recorded by iGrav31 and the simulated gravity effects due to saturated (SZ), unsaturated (UZ), or sum of the two zones (UZ+SZ).

Figure 10

 10 Figure 10 Observed (black) and predicted (colored) ∆g(iGrav31)-∆g(iOSG24) gravity signal for various depths of the water layer and integration radii. (a) For a depth of 300 m and integration radii of 50, 100 and 200 m; (b) For a depth of 600 m and integration radii of 1000, 2000 and 3000 m; (c) For an integration radius of 600 m and a water of layer at depths 100, 400 and 800 m; (d) For an integration radius of 1000 m and a water layer at depths 100, 400 and 800 m.

Acknowledgments

Surface loading models based on MERRA2 and ERA5 are available through the EOST loading service (http://loading.ustrasbg.fr). S.R. thanks N. Mazzilli (University of Avignon, France) for some discussions on this work. We are grateful to F. Littel, D. Boyer, J.-B. Decitre and S.

Gaffet for installing and maintaining the SGs at the LSBB site. The digital elevation model provided by the "Institut national de l'information géographique et forestière" (IGN-F) was downloaded from https://geoservices.ign.fr/rgealti

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Author Contributions