Sylvain Pasquet
email: sylvain.pasquet@upmc.fr

Ludovic Bodet

Amine Dhemaied

Amer Mouhri

Quentin Vitale

Fayçal Rejiba

Nicolas Flipo

Roger Guérin

Detecting different water table levels in a shallow aquifer with combined P-, surface and SH-wave surveys: insights from V P /V S or Poisson's ratios

Keywords: hydrogeology, seismic methods, surface waves, P-wave, shear-wave, water table, V P /V S ratio, Poisson's ratio

When applied to hydrogeology, seismic methods are generally confined to the characterisation of aquifers geometry. The joint study of pressure-(P) and shear-(S) wave velocities (V P and V S ) can provide supplementary information and improve the understanding of aquifer systems. This approach is proposed here with the estimation of V P /V S ratios in a stratified aquifer system characterised by tabular layers, well-delineated thanks to electrical resistivity tomography, log and piezometer data. We carried out seismic surveys under two hydrological conditions (high and low flow regimes) to retrieve V S from both surface-wave dispersion inversion and SH-wave refraction interpretation, while V P were obtained from P-wave refraction interpretation. P-wave first arrivals provided 1D V P structures in very good agreement with the stratification and the water table level.

Both V S models are similar and remain consistent with the stratification. The theoretical dispersion curves computed from both V S models present a good fit with the maxima of dispersion images, even in areas where dispersion curves could not be picked. Furthermore, V P /V S and Poisson's ratios computed with V S models obtained from both methods show a strong contrast for both flow regimes at depths consistent with the water table level, with distinct values corresponding to partially and fully saturated sediments.

A C C E P T E D

M A N U S C R I P T

Introduction

Characterisation and monitoring of groundwater resources and associated flow and transport processes mainly rely on the implementation of wells (piezometers). The interpretation of hydrogeological observations is however limited by the variety of scales at which these processes occur and by their variability in time. In such a context, using geophysical (mostly electromagnetic and electrical) methods often improves the very low spatial resolution of borehole data and limit their destructive nature [START_REF] Guérin | Borehole and surface-based hydrogeophysics[END_REF][START_REF] Hubbard | Hydrogeophysics[END_REF]. These methods regularly help to characterise the geometry of the basement [START_REF] Mouhri | Designing a multi-scale sampling system of stream-aquifer interfaces in a sedimentary basin[END_REF], identify and assess the physical and environmental parameters affecting the associated flow and transport processes [START_REF] Mcclymont | Investigating groundwater flow paths within proglacial moraine using multiple geophysical methods[END_REF], and possibly follow the evolution of these parameters over time [START_REF] Michot | Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography[END_REF][START_REF] Gaines | Detecting perched water bodies using surface-seismic time-lapse traveltime tomography[END_REF]. They also tend to be proposed to support the implantation of dense hydrological monitoring networks [START_REF] Mouhri | Designing a multi-scale sampling system of stream-aquifer interfaces in a sedimentary basin[END_REF].

Among the geophysical tools applied to hydrogeology, seismic methods are commonly used at different scales, but remain mainly confined to the characterisation of the aquifer geometry. With dense acquisition setups and sophisticated workflows and processing techniques, seismic reflection produce detailed images of the basement with the resolution depending on the wavelength (Haeni, 1986a;[START_REF] Juhlin | High-resolution reflection seismics applied to detection of groundwater resources in glacial deposits, Sweden[END_REF][START_REF] Bradford | Depth characterization of shallow aquifers with seismic reflection, Part I-The failure of NMO velocity analysis and quantitative error prediction[END_REF][START_REF] Bradford | Depth characterization of shallow aquifers with seismic reflection, Part II-Prestack depth migration and field examples[END_REF][START_REF] Haines | Hydrogeologic structure underlying a recharge pond delineated with shear-wave seismic reflection and cone penetrometer data[END_REF][START_REF] Kaiser | Ultrahigh-resolution seismic reflection imaging of the Alpine Fault, New Zealand[END_REF]. These images are routinely used to describe the stratigraphy in the presence of strong impedance contrasts, but do not allow for distinguishing variations of a specific property [START_REF] Pride | Relationships between seismic and hydrological properties[END_REF][START_REF] Hubbard | Hydrogeophysics[END_REF]. From these images, hydrogeologists are able to retrieve the geometry of aquifer systems, and allocate a lithology to the different layers with the help of borehole data [START_REF] Paillet | Integrating surface geophysics, well logs and hydraulic test data in the characterization of heterogeneous aquifers[END_REF][START_REF] Guérin | Borehole and surface-based hydrogeophysics[END_REF].

Surface refraction seismic provides records from which it is possible to extract the propagation velocities of seismic body waves. This method has the advantage of being relatively inexpensive and quick to implement, and is easily carried out with a 1D to 3D coverage [START_REF] Galibert | Seismic study of the low-permeability volume in southern France karst systems[END_REF].

It is frequently chosen to determine the depth of the water table when the piezometric surface is considered as an interface inside the medium (i.e. free aquifer) [START_REF] Wallace | Some limitations of seismic refraction methods in geohydrological surveys of deep alluvial basins[END_REF]Haeni, 1986b[START_REF] Haeni | Application of seismic-refraction techniques to hydrologic studies[END_REF][START_REF] Paillet | Integrating surface geophysics, well logs and hydraulic test data in the characterization of heterogeneous aquifers[END_REF][START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, Part I: Water table, fluid flow, and saturation[END_REF]. But the seismic response in the presence of such interfaces, and more generally in the context of aquifer characterisation, remains complex [START_REF] Ghasemzadeh | Effect of subsurface hydrological properties on velocity and attenuation of compressional and shear wave in fluid-saturated viscoelastic porous media[END_REF]. The interpretation of the estimated velocities is often difficult because

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

their variability mainly depend on the "dry" properties of the constituting porous media. In these conditions, borehole seismic (up-hole, down-hole, cross-hole, etc.) are regularly used to constraint velocity models in depth, though they remain destructive and laterally limited [START_REF] Haeni | Application of seismic-refraction techniques to hydrologic studies[END_REF][START_REF] Sheriff | Exploration seismology[END_REF][START_REF] Liberty | Surface and borehole seismic characterization of the boise hydrogeophysical research site[END_REF][START_REF] Steeples | Shallow seismic methods[END_REF][START_REF] Dal Moro | Unambiguous determination of the Vs profile via joint analysis of multi-component active and passive seismic data[END_REF].

Geophysicists seek to overcome these limitations, especially through the joint study of compression (P-) and shear (S-) wave velocities (V P and V S , respectively), whose evolution is by definition highly decoupled in the presence of fluids (Biot, 1956a,b). The effect of saturation and pore fluids on body wave velocities in consolidated media has been subject to many theoretical studies [START_REF] Berryman | Origin of Gassmann's equations[END_REF][START_REF] Lee | Modified Biot-Gassmann theory for calculating elastic velocities for unconsolidated and consolidated sediments[END_REF][START_REF] Dvorkin | Yet another Vs equation[END_REF] and experimental developments [START_REF] Wyllie | Elastic wave velocities in heterogeneous and porous media[END_REF][START_REF] King | Wave velocities in rocks as a function of change in overburden pressure and pore fluid saturants[END_REF][START_REF] Nur | The effect of saturation on velocity in low porosity rocks[END_REF][START_REF] Domenico | Effect of water saturation on seismic reflectivity of sand reservoirs encased in shale[END_REF][START_REF] Gregory | Fluid saturation effects on dynamic elastic properties of sedimentary rocks[END_REF][START_REF] Domenico | Elastic properties of unconsolidated porous sand reservoirs[END_REF][START_REF] Murphy | Effects of partial water saturation on attenuation in Massilon sandstone and Vycor porous glass[END_REF][START_REF] Dvorkin | Acoustic signatures of patchy saturation[END_REF][START_REF] Foti | Porosity of fluid-saturated porous media from measured seismic wave velocities[END_REF][START_REF] Prasad | Acoustic measurements in unconsolidated sands at low effective pressure and overpressure detection[END_REF][START_REF] Adam | Gassmann's fluid substitution and shear modulus variability in carbonates at laboratory seismic and ultrasonic frequencies[END_REF][START_REF] Uyanık | The porosity of saturated shallow sediments from seismic compressional and shear wave velocities[END_REF], especially in the fields of geomechanics and hydrocarbon exploration. From a theoretical point of view, this approach proves suitable for the characterisation of aquifer systems, especially by estimating V P /V S or Poisson's ratios [START_REF] Stümpel | The use of seismic shear waves and compressional waves for lithological problems of shallow sediments[END_REF][START_REF] Castagna | Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks[END_REF][START_REF] Bates | Studies in P-wave and S-wave seismics[END_REF][START_REF] Bachrach | Seismic velocities and Poisson's ratio of shallow unconsolidated sands[END_REF]. Recent studies show that the evaluation of these ratios, or derived parameters more sensitive to changes in saturation of the medium, can be systematically carried out with seismic refraction tomography using both P and SH (shear-horizontal) waves [START_REF] Turesson | A comparison of methods for the analysis of compressional, shear, and surface wave seismic data, and determination of the shear modulus[END_REF][START_REF] Grelle | Seismic refraction methodology for groundwater level determination: "Water seismic index[END_REF][START_REF] Mota | 2D sections of porosity and water saturation from integrated resistivity and seismic surveys[END_REF].

The estimation of the V P /V S ratio with refraction tomography requires to carry out two separate acquisitions for V P and V S . While P-wave seismic methods are generally considered well-established, measurements of V S remain delicate because of well-known shear-wave generation and picking issues in SH-wave refraction seismic methods [START_REF] Sheriff | Exploration seismology[END_REF][START_REF] Jongmans | The importance of surface waves in vibration study and the use of Rayleigh waves for estimating the dynamic characteristics of soils[END_REF][START_REF] Xia | A pitfall in shallow shear-wave refraction surveying[END_REF][START_REF] Haines | A hammer-impact, aluminium, shear-wave seismic source[END_REF]. Indirect estimation of V S is commonly achieved in a relative straightforward manner by using surface-wave prospecting methods, as an alternative to SH-wave refraction tomography (e.g. [START_REF] Gabriels | In situ measurements of shear-wave velocity in sediments with higher-mode Rayleigh waves[END_REF][START_REF] Jongmans | The importance of surface waves in vibration study and the use of Rayleigh waves for estimating the dynamic characteristics of soils[END_REF][START_REF] Park | Multichannel analysis of surface waves[END_REF][START_REF] Socco | Surface-wave method for near-surface characterization: a tutorial[END_REF][START_REF] Socco | Surface-wave analysis for building near-surface velocity models -Established approaches and new perspectives[END_REF]. Such approach has recently been proposed for geotechnical [START_REF] Heitor | Characterising compacted fills at Penrith Lakes development site using shear wave velocity and matric suction[END_REF] and hydrological applications in sandy aquifers [START_REF] Cameron | A new approach to predict hydrogeological parameters using shear waves from multichannel analysis of surface waves method[END_REF][START_REF] Konstantaki | Determining hydrological and soil mechanical parameters from multichannel surface-wave analysis across the Alpine Fault at Inchbonnie, New Zealand[END_REF][START_REF] Fabien-Ouellet | Using all seismic arrivals in shallow seismic investigations[END_REF]. [START_REF] Konstantaki | Determining hydrological and soil mechanical parameters from multichannel surface-wave analysis across the Alpine Fault at Inchbonnie, New Zealand[END_REF] highlighted major variations of V P /V S and Poisson's ratios that was correlated with the water table level. Retrieving V P and V S from a single acquisition setup thus appears attractive in terms of time

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

and equipment costs, even if SH-wave methods provide high quality results in reflection seismic [START_REF] Hunter | Surface and downhole shear wave seismic methods for thick soil site investigations[END_REF][START_REF] Guy | High-resolution SH-wave seismic reflection investigations near a coal mine-related roadway collapse feature[END_REF][START_REF] Haines | Shear-wave seismic reflection studies of unconsolidated sediments in the near surface[END_REF][START_REF] Ghose | Signature of fault zone deformation in near-surface soil visible in shear wave seismic reflections[END_REF]. Moreover, [START_REF] Pasquet | Surface-wave dispersion inversion versus SH-wave refraction tomography in saturated and poorly dispersive quick clays[END_REF] recently evaluated the applicability of the combined use of SH-wave refraction tomography and surface-wave dispersion inversion for the characterisation of V S .

In order to address such issues in more complex aquifer systems (e.g. unconsolidated, heterogeneous or low permeability media), we performed high spatial resolution P-, surface-and SH-wave seismic surveys in the Orgeval experimental basin (70 km east from Paris, France) under two distinct hydrological conditions. This basin is a part of a research observatory managed by the OR-ACLE network (http://bdoracle.irstea.fr/) and has been studied for the last 50 years, with particular focuses on water and pollutant transfer processes occurring at different scales throughout the basin [START_REF] Flipo | Caractérisation des fonctionnements hydrologique et hydrogéologiques du bassin de l'Orgeval[END_REF]. The basin drains a stratified aquifer system characterised by tabular layers, well-delineated all over the basin by [START_REF] Mouhri | Designing a multi-scale sampling system of stream-aquifer interfaces in a sedimentary basin[END_REF] thanks to extensive geological and geophysical surveys including Electrical Resistivity Tomography (ERT), Electrical Soundings (ES), Time Domain ElectroMagnetic (TDEM) soundings and borehole core sampling. The hydrogeological behaviour of the Orgeval watershed is influenced by the aquifer system, which is composed of two main geological units: the Oligocene sand and limestone (Brie formation on Fig. 1b) and the Middle Eocene limestone (Champigny formation on Fig. 1b) [START_REF] Mouhri | Designing a multi-scale sampling system of stream-aquifer interfaces in a sedimentary basin[END_REF]. These two aquifer units are separated by a clayey aquitard composed of green clay and marl (Fig. 1b). Most of the basin is covered with table-land loess of about 2-5 m in thickness, essentially composed of sand and loam lenses of low permeability. These unconsolidated deposits seem to be connected to the Oligocene sand and limestone, forming a single aquifer unit. This upper aquifer is monitored by a dense network of piezometers (Fig. 1a) [START_REF] Mouhri | Designing a multi-scale sampling system of stream-aquifer interfaces in a sedimentary basin[END_REF] which have allowed for establishing maps of the piezometric level for high and low water regimes in 2009 and 2011 [START_REF] Kurtulus | Hydraulic head interpolation in an aquifer unit using ANFIS and Ordinary Kriging[END_REF][START_REF] Kurtulus | Hydraulic head interpolation using anfis-model selection and sensitivity analysis[END_REF]. It thus offers an ideal framework for the study of the V P /V S ratio through the combined analysis of P-wave refraction, SH-wave refraction and surface-wave dispersion data. Measurements were carried out under two distinct hydrological conditions in order to evaluate the ability of this approach to detect variations of the water table level, and assess its practical limitations.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

2. Location of the experimentation and acquisition strategy

Choice of the site

The experiment location has been selected in a plateau area, where the upper layers of the aquifer system are known to be the most tabular. The site is located in the southeast part of the Orgeval basin, at 70 km east from Paris, near the locality of Les Granges (black square Fig. 1a).

A piezometer (PZ3 on Fig. 1a) with its water window in the Brie aquifer is situated in the middle of a trail crossing the survey area in the southeast-northwest direction. Thanks to the ORACLE facilities, the piezometric head level in the upper aquifer is continuously recorded in PZ3 on an hourly basis (Fig. 2a). Two acquisition campaigns were carried out in the site under two distinct hydrological conditions. The first campaign took place between March 12 th and March 14 th 2013 during a high flow regime (i.e. high water level or HW on Fig. 2a), with a piezometric head level measured at 1.15 m. The second campaign was conducted between August 26 th and August 28 th 2013 during a low flow regime (i.e. low water level or LW on Fig. 2a), with a recorded piezometric head level of 2.72 m. During both HW and LW campaigns, the piezometric head level was measured from ground level at the base of PZ3.

Electrical Resistivity Tomography was performed during both HW and LW campaigns to accurately describe the stratigraphy in the upper aquifer unit and confirm the tabularity required for our experiment. We used a multi-channel resistivimeter with a 96-electrode Wenner-Schlumberger array (Fig. 2b). ERT profiles were implanted on the side of the trail and centred on PZ3 (Fig. 1a), 1 m away from the piezometer and 0.25 m below, respectively. Electrodes were spaced with 0.5 m to obtain 41.5-m long profiles. The inversion was performed using the RES2DINV commercial software [START_REF] Loke | Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method[END_REF]. The origin of the depth axis in Fig. 2 and in figures hereafter was chosen at ground level in the centre of the line (i.e. the water table level is 0.25 m higher than recorded in PZ3). The ORACLE experimental facilities provided soil and air temperatures during both campaigns thanks to probes installed near the survey area. At HW, air temperature was below 0 • C and soil temperature was increasing from 6.3 • C at 0.5 m in depth to 6.5 • C at 1 m in depth. In comparison, air temperature was around 22 • C at LW, with a soil temperature varying from 18.5 • C at 0.5 m in depth to 18 • C at 1 m in depth. With such fluctuations between both campaigns, the variation of ground resistivity due to temperature cannot be neglected. To account for those effects, [START_REF] Campbell | Change of electrical conductivity with temperature and the relation of osmotic pressure to electrical conductivity and ion concentration for soil extracts[END_REF] proposed an approximation stating that an increase of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 1 • C in
temperature causes a decrease of 2 % in resistivity. We used this approximation to correct resistivity values obtained at HW from the temperature differences observed between HW and LW periods, after extrapolating both temperature profiles in depth with an exponential trend [START_REF] Oke | Boundary layer climates[END_REF]. The comparison of the corrected HW ERT profile with the LW ERT profile shows no significant variation of the resistivity values and clearly depicts the stratigraphy with three distinct tabular layers (Fig. 2c) that are consistent with those observed at the basin scale (Fig. 1b). The most superficial layer has a thickness of 0.2 to 0.25 m and an electrical resistivity (ρ) of about 30 Ω.m. This thin layer, corresponding to the agricultural soil, was not observed at the basin scale.

It presents higher resistivity values at LW that can be explained by lower water content at the surface. The second layer, associated with the table-land loess, is characterised by lower electrical resistivity values (around 12 Ω.m), with a thickness of about 3.5 m. The semi-infinite layer has higher electrical resistivity values (around 35 Ω.m), and can be related to the Brie limestone layer. ERT and log results offer a fine description of the site stratigraphy. These results, combined with piezometric head level records, provide valuable a priori information for the interpretation of seismic data.

Seismic acquisition 2.2.1. Acquisition setup

An identical seismic acquisition setup was deployed during both HW and LW campaigns. It consisted in a simultaneous P-and surface-wave acquisition followed by a SH-wave acquisition along the same line. The seismic line was centred on PZ3 (Fig. 1a) along the ERT profile, with the origin of the x-axis being identical to the one used for ERT (Fig. 2b). While a small receiver spacing is required to detect thin layers with seismic refraction, a long spread is needed for surface-wave analysis in order to increase spectral resolution and investigation depth. To meet both requirements, we used a dense multifold acquisition setup with 72 geophones and a 0.5 m receiver spacing to obtain a 35.5-m long profile (Fig. 3). We carried out a topographic leveling using a tacheometer to measure the relative position and elevation of each geophone. The maximum difference of elevation along the profile is around 0.5 m which represent a slope of less than 1.5 %. A 72-channel seismic recorder was used with 72 14-Hz vertical component geophones for the P-wave profile, and 72 14-Hz horizontal component geophones for the S-wave profile. First shot location was one half receiver spacing away from first trace, and move up between shots was one receiver interval. 73 shots were

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

recorded along each profile for a total number of 5256 active traces.

The P-wave source consisted in an aluminium plate hit vertically by a 7-kg sledgehammer. The plate was hit 6 times at each position to increase signal-to-noise ratio. The SH-waves were generated with a manual source consisting of a heavy metal frame hit laterally by a 7-kg sledgehammer. The SH-wave source was hit 8 times at each position. For both P-and SH-wave acquisitions, the sampling rate was 1 ms and the recording time was 2 s (anticipating low propagation velocities).

A delay of -0.02 s was kept before the beginning of each record to prevent early triggering issues (i.e. time shift between the recording starting time and the actual beginning of the seismic signal).

Recorded seismograms

The collected data presented on Fig. 4 are of good quality with low noise level, and did not require specific processing other than basic trace normalisation. P-wave seismograms recorded during both HW (Fig. 4a) and LW (Fig. 4c) campaigns present similar characteristics. P-wave first arrivals are clearly visible before 0.04 s (P on Fig. 4a and4c), with three different apparent velocities visible at first glance: 200 m/s for the first two traces, then 800 m/s for the next 7

to 10 traces, and around 2000 m/s for the farthest traces. They are followed by the air wave, characterised by higher frequencies and a velocity of 340 m/s (A on Fig. 4a and4c). At last come P-SV waves (or Rayleigh waves), corresponding to a high-amplitude and low-frequency wave-train with an apparent velocity of about 150 m/s (R on Fig. 4a and4c). SH-wave shots records obtained during both HW (Fig. 4b) and LW (Fig. 4d) campaigns also show similar features. They contain lower frequency signal, with coherent events consistent with SH-wave first arrivals (SH on Fig. 4b and4d). These first arrivals have three distinct apparent velocities (around 70 m/s for the first two traces, 175 m/s for the next 30 traces, and 450 m/s for the farthest traces). SH-wave first arrivals are directly followed by Love waves (L on Fig. 4b and4d), which present an apparent velocity of about 175 m/s. Early P-wave arrivals are visible on horizontal geophones records, especially on Fig. 4b between 15 and 20 m and before 0.1 s. Even under such excellent experimental conditions, it is always challenging to guarantee the horizontality of geophones. These early events are one of the main features that make first arrival picking delicate when carrying out SH-wave surveys.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

3. Processing and results

Body waves

For both HW and LW, P-and SH-wave traveltimes were easily identified and picked in the raw data from near to long offsets. First arrivals of 5 shots (1 direct shot, 1 reverse shot and 3 evenly spaced split-spread shots) were interpreted as simple 2D models with tabular dipping layers [START_REF] Wyrobek | Application of delay and intercept times in the interpretation of multilayer refraction time distance curves[END_REF][START_REF] Dobrin | Introduction to geophysical prospecting[END_REF]. Traveltimes corresponding to the interpreted models were computed and represented along with observed traveltimes. In the absence of a proper estimation of the traveltimes relative errors and in order to propose an estimate of the accuracy of the interpreted models, we introduced a perturbation of ±5 % on interpreted models (+5 % on velocities and -5 % on thicknesses for the lower model, and -5 % on velocities and +5 % on thicknesses for the upper model), and calculated the corresponding theoretical traveltimes. For the sake of readability, only direct and reverse shots traveltimes were represented Fig. 5 along with ±5 % perturbations. 1D models corresponding to the centre of the profile (i.e. the position of PZ3) were extracted and represented with the corresponding ±5 % perturbation (Fig. 5).

P-wave first arrivals picked for the HW campaign (Fig. 5a) were interpreted as a 3-layer model, with interfaces between layers slightly dipping southeast (less than 1 %). These three layers have P-wave velocities from surface to depth of 250 ± 12.5 m/s, 750 ± 37.5 m/s and 2000 ± 100 m/s, respectively. The two upper layers have thicknesses at the centre of the profile of 0.85±0.043 m and 3±0.15 m, respectively (Fig. 5c). P-wave first arrivals observed for the LW campaign (Fig. 5d) were interpreted with 4 layers presenting slightly dipping interfaces towards southeast (less than 0.5 %).

The corresponding velocities are 170 ± 8.5 m/s, 300 ± 15 m/s, 825 ± 41.25 m/s and 2000 ± 100 m/s from top to bottom. The thicknesses of the three upper layers at the centre of the model are 0.15 ± 0.008 m, 1.2 ± 0.06 m and 2.65 ± 0.133 m, respectively (Fig. 5f). The first layer observed during the LW campaign is missing in the interpretation of first arrivals of the HW campaign.

Indeed, early triggering issues prevented us from picking first arrivals corresponding to this thin layer.

SH-wave first arrivals picked for both HW (Fig. 5b) and LW (Fig. 5e) campaigns were interpreted as 3-layer models, with interfaces slightly dipping southeast (less than 0.25 %). For HW,

ACCEPTED MANUSCRIPT

3.65 ± 0.183 m thick, respectively (Fig. 5c). As for LW, the V S model at the centre of the profile is composed of a low velocity (65 ± 3.25 m/s) and thin (0.3 ± 0.015 m) layer in surface, a 3.5 ± 0.175 m thick layer with a velocity of 170 ± 8.5 m/s, and a semi-infinite layer with a velocity of 425 ± 21.25 m/s (Fig. 5f).

Despite known limitations of the refraction interpretation technique (e.g. in presence of low velocity layers, velocity gradients, etc.), the interpreted velocity models are highly satisfying and provide a description of the stratigraphy in very good agreement with ERT and log results. When

V S show 3 layers corresponding to this stratigraphy, V P present a fourth layer that is consistent with the observed water table level, especially for HW (Fig. 5c). These velocity models are quite stable in depth, as demonstrated by the ±5 % error bars displayed on Fig. 5. Furthermore, the calculated residuals between observed and calculated traveltimes remain mostly below 5 %, with only a few over 10 %, and Root Mean Square (RMS) errors calculated for direct and reverse shots are around 2-2.5 % (Fig. 6). These low values point out the good consistency of the estimated velocity models and reinforce the confidence in our interpretations.

P-SV waves

Extraction of dispersion

Surface-wave dispersion images were obtained from P-wave shot gathers for both HW and LW campaigns (Fig. 7). After correction for geometrical spreading, the wavefield was basically transformed to the frequency-phase velocity (fc) domain in which maxima should correspond to Rayleigh-wave propagation modes [START_REF] Russel | Multi-channel processing of dispersed surface waves[END_REF][START_REF] Mokhtar | Seismic velocity and Q model for the shallow structure of the Arabian Shield from short-period Rayleigh waves[END_REF]. Anticipating slight shallow lateral variations, we used the entire spread to analyse surface waves. A 70-trace extraction window (34.5-m wide) was actually used in order to be roughly centred on PZ3 (x = 24.25 m). For both flow regimes, we obtained dispersion images from direct (Fig. 7a, HW and 7d, LW) and reverse (Fig. 7b, HW and 7e, LW) shots on each side of the window. The comparison of both single dispersion images presented only slight differences, confirming the validity of the 1D approximation [START_REF] Jongmans | Geophysical investigation of a large landslide in glaciolacustrine clays in the Trièves area (French Alps)[END_REF]. These images were thus stacked in order to increase the signal-to-noise ratio (Fig. 7c, HW and 7f, LW). The stacking was achieved by summing the frequency-phase velocity spectra of windowed data (e.g. [START_REF] O'neill | Full-waveform P-SV reflectivity inversion of surface waves for shallow engineering applications[END_REF], which clearly enhanced the maxima.

The dispersion data present a strong "effective character", which aspects are for instance discussed by Forbriger (2003a,b) and O' Neill and Matsuoka (2005). In shallow seismic data, large

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

velocity contrasts and/or velocity gradients often generate wavefields with dominant higher modes.

Guided waves may also appear with large amplitudes at high frequencies and phase velocities. In that case, the identification of different propagation modes and the picking of dispersion curves is challenging and requires a thorough analysis of the observed dispersion images, or alternative inversion approaches (e.g. [START_REF] Maraschini | A new misfit function for multimodal inversion of surface waves[END_REF][START_REF] Boiero | Surface-and guided-wave inversion for near-surface modeling in land and shallow marine seismic data[END_REF]. To facilitate mode identification, we relied on preliminary picking and inversions along with trial and error forward modelling based on a priori geological knowledge and results from refraction analysis. Such approach actually highlighted a "mode-jump" occurring around 35 Hz on each dispersion image, confirming the presence of overlapping modes. Some maxima yet remained hard to identify as propagation modes in the extracted dispersion images, either because they could be seen as secondary lobes of the wavefield transform, or because they were too close to other maxima. To prevent from including "misidentified modes" in dispersion data, maxima were not picked in those areas.

On each dispersion image, coherent maxima were finally extracted with an estimated standard error in phase velocity defined according to the workflow described in O'Neill (2003). Corresponding error bars are not presented on Fig. 7 to keep images readable. Four propagation modes were observed and identified as fundamental (0), first (1), second (2) and third (3) higher modes (Fig. 7).

The apparent phase velocity of the fundamental mode increases with decreasing frequency (from 175 to 350 m/s). As recommended by [START_REF] Bodet | Limites théoriques et expérimentales de l'interprétation de la dispersion des ondes de Rayleigh : apport de la modélisation numérique et physique[END_REF] and [START_REF] Bodet | Near-offset effects on Rayleigh-wave dispersion measurements: Physical modeling[END_REF], we limited dispersion curves down to frequencies (f lim ) where the spectral amplitude of the seismogram became too low (15 Hz on Fig. 7), thus defining the maximum observed wavelength λ max (∼ 22.5 m on Fig. 7).

Inversion

Assuming a 1D tabular medium below each extraction window, we performed a 1D inversion of dispersion data obtained during both HW and LW campaigns. We used the Neighbourhood Algorithm (NA) developed by [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space[END_REF] and implemented for near-surface applications by [START_REF] Wathelet | Surface-wave inversion using a direct search algorithm and its application to ambient vibration measurements[END_REF] and [START_REF] Wathelet | An improved neighborhood algorithm: Parameter conditions and dynamic scaling[END_REF]. Theoretical dispersion curves were computed from the elastic parameters using the Thomson-Haskell matrix propagator technique [START_REF] Thomson | Transmission of elastic waves through a stratified solid medium[END_REF][START_REF] Haskell | The dispersion of surface waves on multilayered media[END_REF]. NA performs a stochastic search of a pre-defined parameter space (namely V P , V S ,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

density and thickness of each layer), using the following misfit function (MF ):

M F = N f i=1 (V cal i -V obs i ) 2 N f σ i 2 , ( 1
)
with V cal i and V obs i , the calculated and observed phase velocities at each frequency f i ; N f , the number of frequency samples and σ i , the phase velocity measurement error at each frequency f i .

Based on site a priori geological knowledge and results from refraction analysis, we used a parametrisation with a stack of three layers (soil, partially saturated loess and fully saturated loess) with an uniform velocity distribution overlaying the half-space (Brie limestone layer). An appropriate choice of these parameters is considered as a fundamental issue for the successful application of inversion [START_REF] Socco | Surface-wave method for near-surface characterization: a tutorial[END_REF][START_REF] Renalier | Influence of parameterization on inversion of surface wave dispersion curves and definition of an inversion strategy for sites with a strong contrast[END_REF]. The thickness of the soil layer was allowed for ranging from 0.05 to 1 m, while the thicknesses of the partially and fully saturated loess could vary between 0.5 and 3.5 m. The half-space depth (HSD), of great importance since it depends on the poorly known depth of investigation of the method, was fixed to about 40 % of the maximum observed wavelength (8 m) as recommended by O'Neill ( 2003)

and [START_REF] Bodet | Surface-wave inversion limitations from laser-Doppler physical modeling[END_REF][START_REF] Bodet | Near-offset effects on Rayleigh-wave dispersion measurements: Physical modeling[END_REF]. The valid parameter range for sampling velocity models was 1 to 750 m/s for V S (based on dispersion observations and refraction analysis). Anticipating a decrease of V S in the saturated zone, we did not constraint velocities to increase with depth in the two layers corresponding to the partially and fully saturated loess, as it is usually done in surface-wave methods [START_REF] Wathelet | An improved neighborhood algorithm: Parameter conditions and dynamic scaling[END_REF]. P-wave velocity being of weak contraint on surface-wave dispersion, only S-wave velocity profile can be interpreted. V P however remain part of the actual parameter space and were generated in the range 10 to 2500 m/s. Density was set as uniform (1800 kg/m 3 ).

A total of 75300 models were generated with NA (Fig. 8a, HW and 8c, LW). Models matching the observed data within the error bars were selected, as suggested by [START_REF] Endrun | S velocity structure and radial anisotropy in the Aegean region from surface wave dispersion[END_REF]. The accepted models were used to build a final average velocity model associated with the centre of the extraction window (dashed line Fig. 8b, HW and 8d, LW). Thickness and velocity accuracy was estimated with the envelope containing the accepted models.

For both HW (Fig. 8b) and LW (Fig. 8d) campaigns, the inversion led to very similar 4-layer V S models. While velocities in the second and third layers were not constrained to increase with depth, neither final V S model presents decreasing velocities. These two models are characterised by the same very thin low velocity layer in surface (around 0.052 ± 0.025 m in thickness with a

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

S-wave velocity of 8 ± 3 m/s). The second layer is slightly thicker for LW (0.67 ± 0.14 m) than for HW (0.56 ± 0.11 m), and has higher V S values for LW (86 ± 15 m/s) than for HW (79 ± 10 m/s).

The third layer has identical thickness for both flow regimes (3.47 ± 0.25 m), but V S is slighlty higher for LW (179 ± 10 m/s) than for HW (169 ± 5 m/s). The half-space is also characterised by very similar velocities for both flow regimes, with 459 m/s for HW (between 430 and 570 m/s), and 464 m/s for LW (between 380 and 740 m/s). Dispersion curves being less well defined at low frequencies, a larger variability (i.e. larger error bars) of half-space velocities is observed, especially for LW.

This first layer is actually very thin and "slow" but was identified on the field and corresponds to a "mode-jump" in the fundamental mode at about 35 Hz. The high frequency part of this mode could not be picked on dispersion images (Fig. 7) due to stronger higher modes above that frequency, and was thus not included as a priori information in our parameterisation. Using only the fundamental mode in the inversion would obviously have given different results, with theoretical dispersion curves not necessarily presenting this "mode-jump". The incorporation of higher modes in the inversion process allowed us to constrain the fundamental mode behavior at high frequency, even though we could not identify it above 35 Hz. Indeed, all models included within the error bars (Fig. 8) present the same "mode-jump" at frequencies higher than 35 Hz, leading to velocity models with a thin low velocity layer at the surface.

Cross-validation of V S models

Models obtained from surface-wave dispersion inversion (in red, Fig. 9a for HW and 9c for LW) are remarkably similar to the models obtained from SH-wave refraction interpretation (in green, Fig. 9a for HW and 9c for LW), and are thus very consistent with the stratigraphy observed on ERT and log results (Fig. 2). V S obtained through surface-wave dispersion inversion are however characterised for both flow regimes by a very thin and low velocity layer in surface that is not observed with SH-wave refraction interpretation. The error bars of V S models retrieved from refraction analysis were estimated by introducing a perturbation of ±5 % on the central model parameters (in green, Fig. 9a for HW and 9c for LW). As for error bars of V S models retrieved from surface-wave dispersion inversion (in red, Fig. 9a for HW and 9c for LW), they correspond to the envelope of accepted models for each hydrological regime (i.e. fitting the error bars in Fig. 8).

As a final quality control of inversion results, forward modelling was performed using the 1D

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

V S average models obtained from both surface-wave dispersion inversion and SH-wave refraction interpretation. While models obtained from both methods are remarkably similar, the theoretical dispersion curves computed from surface-wave dispersion inversion results (in red, Fig. 9b for HW and 9d for LW) provide the best fit with the coherent maxima observed on measured dispersion images. The theoretical modes are consistent with the picked dispersion curves, and are well-separated from each other while they looked like a unique and strong mode at first glance.

Interestingly, theoretical dispersion curves calculated from refraction models (in green, Fig. 9a for HW and 9c for LW) are clearly following this effective dispersion which remains representative of the stratigraphy since models from both methods are in good agreement. There is however no evidence of water table level detection, though several authors noticed a significant V S velocity decrease in the saturated zone [START_REF] O'neill | Dominant higher surface-wave modes and possible inversion pitfalls[END_REF][START_REF] Heitor | Characterising compacted fills at Penrith Lakes development site using shear wave velocity and matric suction[END_REF].

Discussion and conclusions

When studying aquifer systems, hydrogeologists mainly rely on piezometric and log data to estimate the spatial variations of water table level and lithology. However, these data provide only local information and require the implantation of boreholes which remain expensive and destructive.

Geophysical methods are increasingly proposed to interpolate this piezometric and lithological information between boreholes and build high resolution hydrological models. If electrical and electromagnetic methods have shown their efficiency for the fine characterisation of the lithology, they remained nonetheless unable to detect the water table level in clayey geological formations such as loess. In order to assess the ability of seismic methods to retrieve water table level variations, we carried out seismic measurements in a site characterised by a tabular aquifer system, welldelineated thanks to ERT, log and piezometer data. Measurements were completed under two distinct hydrological conditions (HW and LW). A simultaneous P-and surface-wave survey was achieved with a single acquisition setup, followed by a SH-wave acquisition along the same line.

A simple refraction interpretation of P-and SH-wave first arrivals provided quasi-1D V P and V S models in conformity with the stratigraphy depicted by ERT and logs during both campaigns. V S models obtained through surface-wave dispersion inversion are matching those obtained with SHwave refraction interpretation, except for a thin low velocity layer in surface, which has only been identified in surface-wave dispersion inversion results. The recomputation of theoretical dispersion

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
curves provided results that are very consistent with the measured dispersion images and proved to be a reliable tool for validating the 1D V S models obtained from SH-wave refraction interpretation and surface-wave dispersion inversion.

While V S remains constant in partially and fully saturated loess, V P exhibits a strong increase at a depth consistent with the observed water table level, especially for HW. This correlation is yet not so obvious for LW. Furthermore, V P values observed in the saturated loess remain lower (around 800 m/s) than the expected values in fully saturated sediments (usually around 1500-1600 m/s). It is however quite hard to find in the literature a range of typical V P values that should be expected in various partially and fully saturated sediments. Most of the existing studies present V P values in saturated sands, where the relationship between V P and water saturation remains quite simple and is thoroughly described by many authors (e.g. [START_REF] Bachrach | Seismic velocities and Poisson's ratio of shallow unconsolidated sands[END_REF][START_REF] Foti | Porosity of fluid-saturated porous media from measured seismic wave velocities[END_REF][START_REF] Prasad | Acoustic measurements in unconsolidated sands at low effective pressure and overpressure detection[END_REF]Zimmer et al., 2007a,b). With more complex mixtures (e.g. containing a significant proportion of clays), the behavior of V P with the saturation becomes more complicated [START_REF] Fratta | Combined TDR and P-wave velocity measurements for the determination of in situ soil density-experimental study[END_REF]. V P values around 800 m/s have already been observed in saturated loess by [START_REF] Danneels | Geophysical investigation and dynamic modelling of unstable slopes: case-study of Kainama (Kyrgyzstan)[END_REF] when studying unstable slopes in Kyrgyzstan. In such low permeability materials, full saturation can be hard to reach (due to an irreducible fraction of air in the pores), thus limiting the maximum V P velocity [START_REF] Lu | Effects of soil water potential and moisture content on sound speed[END_REF][START_REF] Lorenzo | Benchmark hydrogeophysical data from a physical seismic model[END_REF]. The study of V P alone thus remains insufficient to lead back to hydrological information. In order to cope with this limitation, V P /V S (Fig. 10a for HW, 10c for LW) and Poisson's ratios (Fig. 10b for HW, 10d for LW) were computed with V S models retrieved from SH-wave refraction interpretation (in green) and surface-wave dispersion inversion (in red). In any case, V P /V S and Poisson's ratios were computed with V P retrieved from P-wave refraction interpretation.

For HW, V P /V S ratio (Fig. 10a) is around 4 in the soil layer, and Poisson's ratio (Fig. 10b) ranges between 0.45 and 0.48. These values are typical of saturated soils [START_REF] Uyanık | The porosity of saturated shallow sediments from seismic compressional and shear wave velocities[END_REF], and may be explained by the presence of a melting snow cover on the site during the acquisition. Directly down the soil, the loess layer is characterised down to 0.75-0.85-m deep by V P /V S ratio values of 1.5 and Poisson's ratio values of 0.1. These values are unusually low, even for non-saturated sediments, and might be explained by the presence of a frozen layer [START_REF] Wang | Application of ultrasonic technology for physical-mechanical properties of frozen soils[END_REF]. At this depth, consistent with the water table level (0.9 m), V P /V S and Poisson's ratios values increase to 4.5 and 0.47-0.48, respectively. This kind of contrast in a single lithological unit is typical of a transition between

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

partially saturated (low V P /V S and Poisson's ratios) and fully saturated sediments (high V P /V S and Poisson's ratios). V P /V S and Poisson's ratios remain constant in the deepest part of loess and in the Brie limestone layer, reinforcing the assumption of a continuously saturated aquifer. A similar contrast is visible for LW on V P /V S (Fig. 10c) and Poisson's (Fig. 10d) ratios. The depth of this contrast (between 1.25 and 1.40 m) is not in very good agreement with the water table level

(2.47 m), but yet do not correspond to any stratigraphic limit. The low V P /V S and Poisson's ratios values (around 1.7 and 0.24, respectively) in the upper part of the loess support the assumption of a partially saturated area, while the high values of these ratios (around 4.5 and 0.48, respectively) computed in the deepest part of the loess and in the Brie limestone layer are consitent with a fully saturated porous medium.

These results are supported by water content measurements performed on auger sounding samples collected during the LW campaign (soil samples could not be collected during the HW campaign due to unfavorable weather conditions). As can be observed on Fig. 10e, the water content decreases between the soil and the upper part of the loess, and reaches a minimum around 0.8-0.9 m. Between 1.2 and 1.5 m, a small peak of moisture is observed, probably corresponding to a rainfall event that occurred 24 hours before the sounding (pluviometry data are available at http://bdoracle.irstea.fr/). This peak is followed by a progressive increase of water content that reaches a maximum at a depth corresponding to the water table level. Auger refusal was encountered at 2.70 m, thus limiting the number of measurements in the saturated zone. The differences observed for LW between the water table level and the depth of the contrast of V P /V S and Poisson's ratios can be explained by several mechanisms. In near-surface sediments, capillary forces create a saturated zone above the water table [START_REF] Lu | Unsaturated soil mechanics[END_REF][START_REF] Lorenzo | Benchmark hydrogeophysical data from a physical seismic model[END_REF] that can reach up to 60 cm in such silty sediments [START_REF] Lu | Unsaturated soil mechanics[END_REF]. Refraction probably occurs above the water table on this capillarity fringe. The rainfall event observed on Fig. 10 might have a similar effect, since the depth of the peak of moisture corresponds to the depth at which the V P /V S contrast occurs. The decrease of water content between the rainfall peak and the water table probably creates a low velocity zone that alters the first arrivals interpretation (irrespective of the acquisition configuration). The relevance of this tabular interpretation might be called into question if the studied medium is characterised by continuously varying properties with velocities increasing progressively from the partially saturated area to the fully saturated

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
area [START_REF] Cho | Unsaturated particulate materials -particle-level studies[END_REF]. Despite an advanced and thorough analysis of surface-wave dispersion, no decrease of V S is detected in the fully saturated zone. This is probably due to very weak variations of water content between the partially and fully saturated areas (Fig. 10e), which do not produce a significant decrease of V S in such material [START_REF] Dhemaied | Effet de l'état hydrique sur la raideur mécanique[END_REF]. Such issues have to be adressed thanks to laboratory experiments by combining analogue modelling and ultrasonic techniques [START_REF] Bergamo | Physical modelling of a surface-wave survey over a laterally varying granular medium with property contrasts and velocity gradients[END_REF][START_REF] Bodet | Small-scale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] on water saturated porous media [START_REF] Pasquet | Apport des méthodes sismiques à l'hydrogéophysique : importance du rapport Vp/Vs et contribution des ondes de surface[END_REF]. Despite these theoretical issues, our approach provided encouraging results that call for more experimental validation. Furthermore, the use of single acquisition setup to retrieve both V P and V S from refraction interpretation and surface-wave analysis appears promising in terms of acquisition time and costs. Associated with existing piezometric data, seismic measurements could be carried out at a wider scale throughout the entire basin to build high resolution maps of the piezometric level. Its application in more complex (e.g. 2D) cases should also provide valuable information for the study of stream-aquifer interactions. A C C E P T E D M A N U S C R I P T

Acknowledgments

ACCEPTED MANUSCRIPT

x Seismic methods were proposed to assess piezometric level variations x We worked on a well-constrained experimental site x A single acquisition setup were used to retrieve V P and V S x We retrieved V S from surface-wave analysis and V P from P-wave refraction x V P /V S ratios show strong contrasts at depths consistent with piezometric levels

  these three layers are characterised from top to bottom by SH-wave velocities of 50 ± 2.5 m/s, 165 ± 8.25 m/s and 400 ± 20 m/s, respectively. The two upper layers are 0.35 ± 0.018 m and A C C E P T E D M A N U S C R I P T

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 7 :Figure 8 :

 123478 Figure 1: (a) Situation of the Orgeval experimental basin, and location of the experiment. (b) Geological log interpreted at PZ3.

Figure 9 :

 9 Figure 9: Comparison of 1D VS models obtained from SH-wave refraction interpretation (in green) and surfacewave dispersion inversion (in red) for HW (a) and LW (c), with corresponding error bars. The error bars of VS models retrieved from refraction analysis were estimated by introducing a perturbation of ±5 % on the central model parameters. As for error bars of VS models retrieved from surface-wave dispersion inversion, they correspond to the envelope of accepted models for each hydrological regime (i.e. fitting the error bars in Fig. 8). Dispersion curves calculated from both surface-wave dispersion inversion (in red) and refraction interpretation (in green) models are superimposed on the stacked dispersion image obtained for HW (b) and LW (d).

A C C E P T E D

M A N U S C R I P T

ACCEPTED MANUSCRIPT