O Bulliard-Sauret
email: odin.bulliard-sauret@yncrea.fr

J Berindei

S Ferrouillat

L Vignal

A Memponteil

C Poncet

J M Leveque

N Gondrexon

Heat transfer intensification by low or high frequency ultrasound: thermal and hydrodynamic phenomenological analysis

Keywords: Heat transfer enhancement, ultrasound, forced convection, turbulence Highlights

The aim of this work is to quantitatively demonstrate the intensification of heat transfer in forced convection by mean of ultrasonic irradiation at low (25 kHz) or high (2 MHz) frequency. High frequency ultrasound induces convective acoustic streaming while low frequency ultrasonic waves produce mainly cavitation effects. These hydrodynamic phenomena are at the origin of strongly different observations in terms of flow pattern modification and thereby in Nusselt number values. A link is tentatively established between hydrodynamic behaviors at both ultrasonic frequencies and corresponding thermal results.

Hydrodynamic approach was performed with a 2D-2C PIV device while thermal one was carried out under uniform heat flux conditions. It seems that thermal enhancement effect of acoustic streaming (2 MHz) decreases as flow rate increases. This behavior is consistent with the decrease of turbulent kinetic energy produced by acoustic streaming in the same

Introduction

This study compares heat transfer enhancement induced by low and high frequency.ultrasound. A link is proposed between the observed intensification to the hydrodynamic behavior produced by ultrasound in flowing fluid. Low and high frequency ultrasound are at the origin of hydrodynamic phenomena of different nature. In addition, it is investigated to what extent heat transfer intensification by ultrasound can be linked to the nature of this induced effect. It is therefore expected to give recommendations on the choice of ultrasonic frequency in the design of future vibrating heat exchangers depending on parameters such as flow regime, Reynolds number and erosion issues.

State of the art

Ultrasound is a mechanical wave with frequencies above 18 -20 kHz. Acoustic cavitation and acoustic streaming are two hydrodynamic effects induced by the ultrasonic irradiation of a liquid [START_REF] Mason | Applied sonochemistry: Uses of power ultrasound in chemistry and processing[END_REF]. Their intrinsic characteristics are strongly dependent on the frequency of the irradiative waves. A schematic description of these both effects is shown in figure 1.

Fig. 1: Illustration of ultrasonic cavitation and acoustic streaming effects as presented in [START_REF] Legay | Enhancement of heat transfer by ultrasound: Review and recent advances[END_REF] Acoustic cavitation can be defined as a growing, oscillating and collapsing phenomenon of gas and vapor cavities induced by ultrasound propagation within a liquid. Potential uses of acoustic cavitation for industrial process intensification are wide [START_REF] Mason | Applied sonochemistry: Uses of power ultrasound in chemistry and processing[END_REF] (heat exchangers, chemical reactors, antifouling/defrosting system…). Nevertheless, strong technological development is still required in terms of reactors design and optimization [START_REF] Peshkovsky | Acoustic cavitation theory and equipment design principles for industrial applications of high-intensity ultrasound[END_REF]. Acoustic cavitation is able to disturb strongly fluid environment from a hydrodynamic point of view [START_REF] Nomura | Flow pattern in a channel during application of ultrasounic vibration[END_REF] thanks to three different mechanisms:

-At first, micro-streaming, induced by the oscillation of the cavitation bubble wall, produces fluid recirculation around it on a distance scale known as acoustic boundary layer. The thickness of this layer relies on both the kinematic viscosity of the medium and the pulsation of the acoustic wave [START_REF] Landau | Physique théorique -Mécanique des fluides[END_REF]. As an example, in water at 20 °C, acoustic boundary layer ranges from 8 microns with 2 MHz ultrasound to 90 microns with 20 kHz ultrasound [START_REF] Leighton | The acoustic bubble[END_REF].

-The second hydrodynamic phenomenon is the acoustic bubbles displacement. Indeed, these latter can be displaced across the fluid by radiation strength because of the absorption and scattering of incident acoustic waves. This is true for both progressive and stationary waves fields. In the second case, the time averaged radiation pressure is equal to zero only at nodes or anti-nodes, driving then bubbles to stable locations [START_REF] Leighton | The acoustic bubble[END_REF]. For strong acoustic intensity, void ratio of the liquid increases because of the growing number of bubbles producing an impedance wall which attenuates the acoustic intensity and then the acoustic driven bubbles phenomenon [START_REF] Leighton | The acoustic bubble[END_REF][START_REF] Mandroyan | Modification of the ultrasound induced activity by the presence of an electrode in a sonoreactor working at two low frequencies (20 an 40 kHz). Part I: Active zone visualisation by laser tomography[END_REF].

-The third phenomenon is linked to the collapse of the acoustic bubbles which produces strong shock waves within the fluid leading to liquid micro-jets when occurring near a solid wall [START_REF] Leighton | The acoustic bubble[END_REF]. In this latter case, collapse leads generally to surface erosion issues, advantageously used in cleaning processes [START_REF] Mazue | Influence of a Perpendicular Liquid Flow on a Cleaning Process using 20 kHz ultrasoundcharacterisation of the agitation at vicinity of the surface opposite to the transducer[END_REF][START_REF] Bru | System and method of ultrasonic cleaning of a soiled surface, comprising an ultrasonic wave generation device[END_REF].

Nevertheless, the erosion phenomenon is an issue to be solved in order to avoid material damage.

For all these reasons, it has been widely shown that acoustic cavitation is able to trigger or to increase turbulence rate within a liquid flow under laminar or transitional hydrodynamic conditions [START_REF] Nomura | Effects of Turbulence by Ultrasonic Vibration on Fluid Flow in a Rectangular Channel[END_REF][START_REF] Bulliard-Sauret | PhD thesis: Étude expérimentale de l'intensification des transferts thermiques par les ultrasons en convection forcée[END_REF][START_REF] Bulliard-Sauret | Experimental study of heat transfer enhancement using ultrasound on a flat plate in forced convection[END_REF].

Acoustic streaming is a generic term defining a liquid or a gas flow driven by acoustic waves.

There are two kinds of acoustic streaming depending on conditions at which it is enabled.

Thanks to a null net mass flow characteristic, acoustic streaming is able to produce recirculating flows with strong global mixing properties.

-The first one is named Eckart's streaming [START_REF] Eckart | Vortices and streams caused by sound waves[END_REF] and appears because of the spatial attenuation of acoustic wave intensity due to liquid absorption, diffraction of the wave leading to energy spreading, or because of the acoustic beam enlargement. In anyway, all these phenomena lead to energy attenuation in the direction of the acoustic wave propagation inducing pressure gradient and then fluid movement.

Eckart's streaming usually occurs under propagative wave conditions.

-The second type of acoustic streaming occurs under stationary wave conditions [START_REF] Wiklund | Acoustofluidics 14: Applications of acoustic streaming in micro fluidic devices[END_REF]. In this case, the streaming is governed by the dissipation phenomenon occurring within the acoustic boundary layer, i.e. near the oscillating wall of the transducer or near a solid wall exposed to an acoustic wave. Such as for oscillating acoustic bubbles, this leads to recirculating flows within the acoustic boundary layer giving rise to others recirculation modes into the bulk by propagating the movement from one cell to the next one. Near the walls, mixing effects combined with strong shearing strength, are able to lyse cells walls [START_REF] Wiliams | Ultrasonic dispersal of aggregates of Sarcina lutea[END_REF].

Both types of streaming are able to induce turbulent flows [START_REF] Atkas | Heat Transfer Enhancement by acoustic streaming in an enclosure[END_REF][START_REF] Moudjed | Oscillating acoustic streaming jet[END_REF]. Mixing properties of ultrasound were evidenced within sonochemical reactors by Residence Time Distribution (RTD) measurements [START_REF] Gondrexon | Experimental study of the hydrodynamic behaviour of a high frequency ultrasonic reactor[END_REF][START_REF] Monnier | Influence of ultrasound on mixing on the molecular scale for water and viscous liquids[END_REF]. According to the literature, acoustic streaming observed at observed with mechanically full-stirred reactors [START_REF] Gondrexon | Experimental study of the hydrodynamic behaviour of a high frequency ultrasonic reactor[END_REF].

The main effect limiting the convective heat transfer is the presence of a laminar boundary layer (also called viscous sub-layer) near the wall. This is usually defined as a flow region where the thermal diffusion is predominant, contrary to thermal convective effects that prevail in the bulk. A decrease of its thickness allows to reduce the thermal resistance of this interface and thereby to increase heat transfer rate. Stirring, mixing and turbulence properties of both acoustic cavitation and acoustic streaming render consequently ultrasound as an attractive and active method to enhance heat transfer. By this way, heat transfer intensification could also help to achieve two main objectives:

-thermal transfer increase between two fluids separated by a solid interface with fixed parameters (case of heat exchangers) -temperature difference decrease keeping the same heat flux when sensitive fluids are processed (case of liquid food pasteurization)

Studies on heat transfer enhancement by ultrasound began in the middle of the 20 th century. It can be noticed that most of these studies are focused on phase change, like boiling [START_REF] Boziuk | Enhanced Boiling Heat Transfer on Micromachined Surfaces using acoustic actuation[END_REF][START_REF] Zheng | Ultrasonic heat transfer enhancement on different structural tubes in LiBr solution[END_REF], or on natural convection [START_REF] Uhlenwinkel | Investigation of heat transfer from circular cylinders in high power 10 kHz and 20 kHz acoustic resonant fields[END_REF]. Some studies dealt also with heat transfer in forced convection or in the field of heat exchangers but in a less extend. Most of the recent studies have for main purpose fouling reduction [START_REF] Chen | Influence of ultrasound to convectional heat transfer with fouling of cooling water[END_REF]. However, some works since 2000 are devoted to heat transfer enhancement in forced convection by using ultrasound:

-Most of them use low frequency ultrasound (20 kHz to 100 kHz [START_REF] Legay | Enhancement of heat transfer by ultrasound: Review and recent advances[END_REF][START_REF] Gondrexon | Intensification of heat transfer process: Improvement of shell-and-tube heat exchanger performances by means of ultrasound[END_REF][START_REF] Tam | Experimental study of the ultrasonic effect on heat transfer inside a horizontal mini-tube in the laminar region[END_REF]).

Intensification of heat transfer up to 270 % has been reported using a 35 kHz vibrating shell and tube heat exchanger [START_REF] Tisseau | Conception et étude préliminaire d'un échangeur de chaleur tubes et calandre assisté par ultrasons, 18ème Congrès Français de Mécanique conference[END_REF].

-High frequency effect of ultrasound on heat transfer has also been highlighted through recent works showing how acoustically induced streaming can enhance convective heat transfer coefficients [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF].

-Studies on cooling process have also been successfully performed [START_REF] Kurbanov | Use of ultrasound for intensification of heat transfer process in heat exchangers[END_REF][START_REF] Duan | Experimental study on the influence of ultrasonic vibration on heat transfer and pressure drop in heat exchanger tubes[END_REF][START_REF] Yao | Experimental Study on Heat Transfer Enhancement of Waterwater Shell-and-Tube Heat Exchanger Assisted by Power Ultrasonic[END_REF] leading to a 27% heat transfer enhancement.

One might here consider that a positive effect on heat transfer has always been reported in all studies involving ultrasound. A comprehensive overview can be found in two publications [START_REF] Gondrexon | Intensification of heat transfer process: Improvement of shell-and-tube heat exchanger performances by means of ultrasound[END_REF][START_REF] Yao | Research and applications of ultrasound in HVAC field: A review[END_REF]. However, it appears that most of the studies deal with heat transfer enhancement by using ultrasound in a very general way. Only a few of them have been carried out in a phenomenological way especially with the impact of ultrasound on the boundary layer or on the turbulence rate within the liquid flow [START_REF] Nomura | Flow pattern in a channel during application of ultrasounic vibration[END_REF][START_REF] Mazue | Influence of a Perpendicular Liquid Flow on a Cleaning Process using 20 kHz ultrasoundcharacterisation of the agitation at vicinity of the surface opposite to the transducer[END_REF][START_REF] Nomura | Effects of Turbulence by Ultrasonic Vibration on Fluid Flow in a Rectangular Channel[END_REF][START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF][START_REF] Barthès | Characterization of the activity of ultrasound emitted in a perpendicular liquid flow using Particle Image Velocimetry (PIV) and electrochemical mass tranfer measurements[END_REF].

According to the literature, it is reasonable to identify four main parameters that clearly influence heat transfer enhancement by ultrasound. Whereas two of them do concern geometry configuration and fluid properties, the two others relate to ultrasonic wave characteristics:

-Geometry of the system: the distance between the exchange surface and the ultrasonic transducer is an important parameter to be taken into account as it acts directly on acoustic intensity especially in the presence of cavitation [START_REF] Nomura | Streaming induced by ultrasonic vibration in a water vessel[END_REF][START_REF] Lee | Smart cooling technology utilizing acoustic streaming[END_REF].

Moreover, the angle of radiation has a significant impact on thermal results for natural convection in the presence of acoustic cavitation [START_REF] Yukawa | Effect of ultrasonic vibration on free convection heat transfer from an inclinate plate in water[END_REF].

-Fluid used may have different physicochemical properties affecting ultrasonic effects previously described. For example, bubble behavior is driven by saturation pressure and surface tension [START_REF] Leighton | The acoustic bubble[END_REF]. Viscous fluids are less sensitive to cavitation [START_REF] Leighton | The acoustic bubble[END_REF].

On the other hand, increasing the fluid viscosity (bulk and dynamic) leads to an increase of the acoustic absorption coefficient. As an example, this coefficient is of prime importance into acoustic forces induced by a plane acoustic wave responsible for acoustic streaming phenomenon [START_REF] Moudjed | Near-field acoustic streaming jet[END_REF].

-Ultrasound frequency appears to be a very important parameter that determines the nature of induced-hydrodynamic phenomena and thereby the heat transfer intensification trend [START_REF] Bulliard-Sauret | PhD thesis: Étude expérimentale de l'intensification des transferts thermiques par les ultrasons en convection forcée[END_REF]. Indeed, hydrodynamics phenomena produced by ultrasound (cavitation or streaming) are different with respect to frequency leading to different thermal behaviors.

-Ultrasonic power: thermal transfer enhancement caused by ultrasound seems to be interrelated with delivered ultrasonic powers [START_REF] Legay | Enhancement of heat transfer by ultrasound: Review and recent advances[END_REF]. As an example, Nusselt number is directly linked to the supplied ultrasonic power in the case of heat transfer intensification due to the presence of acoustic streaming at high frequency [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF].

As aforesaid, very few articles deal with enhancement of heat transfer in forced convection using ultrasound. Main experimental set-ups use shell and tube heat exchanger geometry and low ultrasound frequency (20 kHz to 100 kHz). Furthermore, most of them were limited to cavitation effect on heat transfer intensification and reported the influence of liquid flow rate on enhancement values [START_REF] Gondrexon | Intensification of heat transfer process: Improvement of shell-and-tube heat exchanger performances by means of ultrasound[END_REF][START_REF] Yao | Experimental Study on Heat Transfer Enhancement of Waterwater Shell-and-Tube Heat Exchanger Assisted by Power Ultrasonic[END_REF][START_REF] Legay | Improvement of heat transfer by means of ultrasound: Application to a double tube heat exchanger[END_REF].

Modeling of thermo-hydraulic effects generated by ultrasound is of prime importance for the integration of this technology to optimize many industrial processes using heat exchangers. Some mathematical [START_REF] Moudjed | Scaling and dimensional analysis of acoustic streaming jet[END_REF] and numerical [START_REF] Louisnard | A viable method to predict acoustic streaming in presence of cavitation[END_REF][START_REF] Talebi | Numerical investigation of natural convection heat transfer in a cylindrical enclosure due to ultrasonic vibrations[END_REF][START_REF] Kiani | Experimental analysis and modeling of ultrasound assisted freezing of potato spheres[END_REF] models dealing with hydrodynamic effects of ultrasound have been released. These models have been designed for low [START_REF] Louisnard | A viable method to predict acoustic streaming in presence of cavitation[END_REF][START_REF] Talebi | Numerical investigation of natural convection heat transfer in a cylindrical enclosure due to ultrasonic vibrations[END_REF][START_REF] Kiani | Experimental analysis and modeling of ultrasound assisted freezing of potato spheres[END_REF] and high [START_REF] Moudjed | Scaling and dimensional analysis of acoustic streaming jet[END_REF] frequency range and specific conditions of use of ultrasound. It is therefore necessary to investigate, identify and qualify parameters involved on thermo-hydraulic effects induced by ultrasound in a larger way in order to feed these models. Over all, it is still phenomenologically unclear how ultrasonically induced phenomena (cavitation and streaming) can modify a flow.

Objectives

The objectives of this work are thus to carry out both hydrodynamic and thermal approaches in the presence of intense acoustic cavitation (25 kHz) or strong acoustic streaming (2 MHz) in order to propose a relevant phenomenological explanation of the frequency effect.

Motivation for using high frequency ultrasound comes from strong erosion issues that could occur in heat exchangers by using low frequency ultrasound even if antifouling action at this frequency has been proved [START_REF] Legay | Experimental investigations of fouling reduction in an ultrasonically-assisted heat exchanger[END_REF]. In the aim of designing future optimized vibrating heat exchangers, it is highly desired to compare thermal performances in relation with phenomena implied at low and high frequency ultrasound. Benefits and drawbacks for both frequencies with respect to specific requirements and operating conditions of these future heat exchangers will then be analyzed.

Experimental apparatus and data reduction 2.1 General presentation of the test section

An experimental setup was designed and built (Fig. 2) in order to investigate heat transfer intensification in forced convection by ultrasound. An overall description of the bench test used for this study has been detailed in a previous article in which preliminary results of heat transfer enhancement using 2 MHz ultrasound were exposed [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF]. The canal is made of a block composed of an insulating material that encompass five heating plate blocks (respectively 5 and 6) and by the ultrasonic transducer [START_REF] Mazue | Influence of a Perpendicular Liquid Flow on a Cleaning Process using 20 kHz ultrasoundcharacterisation of the agitation at vicinity of the surface opposite to the transducer[END_REF]. The heating plate ( 6) is made of five stacked heating elements, thermally insulated from each other and individually power supplied [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF]. Each element is equipped with temperature sensors in order to determine local convective heat transfer coefficients. The thickness of the canal (xaxis) can be changed by using the device 8. In this present study, the canal size is 30 * 100 mm (x-z plan) and the total heat plate (6) length is 100 mm (y axis). Finally, two perpendicular visualization windows (9 and 10) allow measuring velocity and turbulent fields in the x-y plan thanks to a 2D-2C Particle Image Velocimetry (PIV) method.

Local heat transfer coefficients in front of each heating element can be determined following the procedure described in [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF]. It is thereby possible to calculate Nusselt number as follows:

= ℎ (1)
with ℎ [W/m 2 /K] the local convective heat transfer coefficient in front of the i th heating element,

[m] the hydraulic diameter and [W/m/K] the thermal conductivity of the fluid.

Heat Transfer Enhancement Factor ( ) can then be defined as:

= - (2)
where and [-] are the Nusselt numbers with and without ultrasound respectively.

Characterization of ultrasound transducers

In this work, two different ultrasonic frequencies were used. in order to compare specific effects of acoustic streaming and of acoustic cavitation on convective heat transfer in forced convection. To ensure comparison, the ultrasonic power delivered to the liquid flowing throughout the canal was estimated by means of a classical calorimetric method and kept identical for both transducers. The 25 kHz transducer was specifically designed by Sinaptec with an active surface size identical to that of the 2 MHz Sonosys transducer. Both transducers were purposely characterized in order to evaluate the presence or not of acoustic cavitation and of acoustic streaming using 25 kHz or 2 MHz. It was confirmed that 2 MHz ultrasound irradiation mainly results in strong acoustic streaming. Acoustic cavitation was not revealed by chemical (iodide oxidation) and aluminum foil tests for this frequency (2 MHz).

On the other hand, none significant chemical activity was observed using iodide oxidation test at 25 kHz. Nevertheless, strong mechanical effects due to intense cavitation were clearly illustrated by the aluminum foil erosion test [START_REF] Yukawa | Effect of ultrasonic vibration on free convection heat transfer from an inclinate plate in water[END_REF] as it can be observed in Fig. 3.

According to these results, it was assumed that the low-frequency transducer (25 kHz) only produces acoustic cavitation with strong mechanical effects whereas the high frequency one

Particle Image Velocimetry (PIV) set-up

Hydrodynamic modifications in the flow induced by ultrasound were investigated by PIV technique using the set-up described in [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF]. instantaneous velocity fields can then be calculated thanks to the following equation:

= 3 4 (3)
with [m/s] the root mean square velocity as calculated by the following equation:

where V [m/s] is the value of the modulus of the instantaneous velocity of the i th image on a total of 2000 images and V [m/s] is the modulus of the average velocity. It is assumed that the velocity component V , which is not measured because of the x-y plan orientation of the laser sheet, contains a similar contribution as the one contained within the values of V and V ! [START_REF] Lavision | Product Manual DaVis 8.3 Software[END_REF].

This assumption implies isotropy of the turbulence even if it is questionable near the heating plate. Nevertheless, the correction of 3 / 4 instead of ½, classically used for turbulence calculated from planar velocity fields, does not change the observed phenomenology. It only affects the absolute values. Since this study is mainly phenomenological and aims at

= # ∑ %V -V & ' () %N -1& (4)
comparing the effects at two different ultrasonic frequencies, the bias will therefore be identical for all measurements without being detrimental to the analysis.

Operating parameters and experimental protocol

Operating hydrodynamic parameters used in this study are displayed in Tab. Canal thickness (x-axis) was kept constant at 30 mm in this study whereas canal width (zaxis) and canal height (y-axis) were both set at 100 mm. Thermal measurements reported in this paper were done in front of the fifth heating block, i.e. the upper block of the heating plate, which is 20 mm height (see Fig. 1). This choice was made because of the fully developed hydrodynamic and thermal conditions in front of the fifth block as previously

shown [START_REF] Bulliard-Sauret | Heat transfert enhancement using 2 MHz ultrasound[END_REF]. The hydraulic diameter is equal to 0.046 m for this series of measurements. To ensure a fair comparison, the ultrasonic power was kept constant to 110 W for both 25 kHz and 2 MHz with the assumption of uniform distribution of 22 W per heating element.

Moreover, the total heating power being equal to 450 W, each heating element is assumed to deliver uniformly 90 W. Nusselt increases as Reynolds increases. This can be fully explained by the decrease of the laminar boundary layer thickness near the heating wall leading thereby to the decrease of the corresponding thermal resistance [START_REF] Kakaç | Handbook of convective heat transfer[END_REF]. However, it is important to note that the laminar/turbulent transition classically observed for Reynolds numbers between 2000 and 3000 in internal forced convection is not experimentally observed here. This can be explained by the presence of a 90° turn at the fluid inlet. This unusual configuration seems to generate hydrodynamic disturbances, causing a gradual enhancement of heat transfer for Reynolds numbers lower than 2000. This behavior is contrary to theory, which predicts constant Nusselt in laminar regime. As stated before, the presence of an ultrasonic field enhances heat transfer but more interestingly, these results highlight its dependency towards the ultrasonic frequency. Indeed, at high frequency (2 MHz), an increase of the flow velocity leads to a slow decrease of the observed heat transfer intensification. The reverse occurs at 25 kHz where an increase of flow velocity leads to an increase of the Nusselt number reaching a plateau for the highest tried flow velocities and Reynolds numbers. Surprisingly, at the lowest investigated flow velocity in this study, the intensification of heat transfer is stronger at 2 MHz than at 25 kHz; whereas beyond flow velocities of 0.033 m/s (i.e. Re = 1500) heat transfer enhancement is higher using low frequency ultrasound. In order to understand and to propose relevant explanations for the difference highlighted by thermal results despite the same supplied ultrasonic and thermal power, a hydrodynamic study using PIV was then performed.

Results and discussions

Hydrodynamic study results under 25 kHz and 2 MHz sonication

Preliminary observations without ultrasound

The average velocity field observed without ultrasound is given in Fig. 5 As mentioned in the state of the art, the relationship between acoustic streaming and induced turbulences within a liquid or a flow has been already established [START_REF] Moudjed | Oscillating acoustic streaming jet[END_REF]. The hydrodynamic behavior analysis of the average velocity fields presented in Fig. 6 can thus explain those presented in Fig. 7 wherein turbulent kinetic energy fields are reported. Turbulent kinetic energy (TKE) field with a mean flow velocity of 0.020 m.s -1 (Re = 900) reveals that positions of TKE highest values (12.5 mm < y < 20 mm) do correspond to the ones where the average velocity field is the most disturbed. Attenuation of velocity field modification by 2 MHz ultrasound, with respect to the mean flow velocity increase (and corresponding Reynolds number) illustrated in Fig. 6, can be directly correlated to the decrease of the TKE values in the same experimental conditions as reported in Fig. 7. Similarly, the decrease of Turbulent Kinetic Energy produced by 2 MHz ultrasound with respect to the mean flow velocity increase can be linked to the decrease of Nusselt number as detailed in Fig. 4 under the same conditions.

The present study is not able to determine real causes of this extinguishing phenomenon of the streaming strength when liquid flowrate is increased. Nevertheless, two hypotheses can be drawn at this point. At first, the flow observed within the canal for 2 MHz sonication may be produced by the combination of the pumping pressure due to the mean flow on the y-axis direction and radiation pressure due to the secondary flow on the x-axis direction (acoustic streaming). Consequently, the higher the pumping pressure in comparison with the radiation pressure (which is constant because of ultrasound power kept constant in this study), the lesser the influence of acoustic streaming is. The second hypothesis could be that the pressure field within the canal is disturbed when the flow becomes more turbulent because of pressure fluctuations induced by turbulences leading possibly to a masking effect, thus hindering acoustic wave propagation and then acoustic streaming production. A future Reynolds similarity study could be helpful to answer to this question.

Experimental results in the presence of 25 kHz ultrasound

Average velocity fields for 25 kHz sonication are displayed in disturb the flow velocity within a disturbance area, which is of bigger size than the cluster bubble itself (black doted arrow on image (c) of Fig. 10). Two hypotheses can be drawn to explain this phenomenon. In the one hand, the violent bubble wall move induces fluid entrainment away from it. This may be supported by the fact that interrogation windows all around the cluster area display also higher velocity values than 0.020 m/s, which is the mean flow velocity. The second hypothesis is a bias of the velocity calculation. Indeed, the velocity is calculated on a 16 * 16 pixel interrogation window whose each pixel measures 16.5 * 16.5 μm 2 . Consequently, PIV measurements are not able to clearly distinguish phenomena inducing effects within an area smaller than 264 *264 μm 2 . In this study, the used interrogation window size gives the most accurate measurement technically possible. As already aforementioned, some previous studies have established a clear relationship not only between acoustic cavitation and turbulences increase [START_REF] Nomura | Flow pattern in a channel during application of ultrasounic vibration[END_REF][START_REF] Bulliard-Sauret | PhD thesis: Étude expérimentale de l'intensification des transferts thermiques par les ultrasons en convection forcée[END_REF] but also between bubble collapse and recirculation flows around bubbles walls [START_REF] Leighton | The acoustic bubble[END_REF]. Consequently, it seems reasonable to link the increase of turbulent kinetic energy in the presence of 25 kHz ultrasound (Fig. 11)

to the collapses and oscillations of bubbles and bubble clusters (Fig. 10). Whatever the set mean flow velocity in the presence of 25 kHz ultrasound, turbulent kinetic energy increases as compared to silent regime. Nevertheless, highest values are measured close to the heating wall and for the highest flow velocity where the turbulent kinetic energy can reach locally up to 0.15 m 2 /s 2 (Fig. 11). Turbulent kinetic energy mean values within the bulk are overall similar for all experimental flow conditions and are equal to 0.014 m 2 /s 2 (see Fig. 11). This means that acoustic cavitation can induce strongest turbulence rates near the heating wall. This observation is consistent with the fact that acoustic cavitation bubble collapse is known to be more violent when occurring near a solid wall [START_REF] Leighton | The acoustic bubble[END_REF]. Bubbles behavior is driven by the ultrasonic frequency and power, liquid properties such as viscosity, surface tension and vapor pressure and by the presence or not of solid wall. Any change in these parameters leads to the modification of bubbles size and collapsing diameter, oscillating frequency and strength of the collapse. Consequently, turbulent kinetic energy within the bulk remains constant whatever mean flow velocity. This observation is consistent with the unchanged operating parameters influencing bubbles behaviors. is dependent on several aforesaid parameters but not on the flow regime in contrary to the boundary layer thickness (. in Fig. 12) which decreases as mean flow velocity and Reynolds number increase. The velocity spots size (-) distribution is ranging from 250 micrometers to 500 micrometers, i.e. one or two PIV interrogation windows (squares in Fig. 12). This size distribution was observed for all mean flow velocities tested in this study. The flow velocity within the canal does not influence the velocity spots size nor their spatial distribution. Because the boundary layer thickness decreases as flow velocity and associated

Reynolds increase, the ratio -. / increases for the same conditions. Consequently, the disturbance effect within the boundary layer due to bubbles collapse increases as mean flow velocity and Reynolds number increase. Indeed, violent collapse of bubbles literally disrupts the boundary layer, laminar flow being turned into a very turbulent one over a thickness defined by -. This value is driven by the collapse strength of the bubble and by viscous effects whose attenuate disturbance propagation. The smaller the thickness of the boundary layer is (i.e. the higher flow velocity and Reynolds number are), the easier bubbles are able to disrupt it over all its thickness It seems relevant to say that this phenomenon leads to the increase of the heat transfer enhancement with respect to flow rate increase.

Thermo-hydraulic study between 25 kHz and 2 MHz sonication

As already shown in this paper, thermal and hydrodynamic behaviors at 25 kHz and 2 MHz sonication are different. Moreover, in Fig. 13, average turbulent kinetic energy values are plotted for both frequencies with mean flow velocities equal to 0.020 and 0.122 m/s (Re = 900 and 5500 respectively) as a function of canal position. These calculated values give the mean value for each interrogation window column (y-axis) for a fixed row (x-axis) in order to obtain an average turbulent kinetic energy profile with respect to canal thickness (x-axis).

This leads to quantitative values required to compare hydrodynamic effects for each ultrasound frequency.

As illustrated in Fig. 13, for 2 MHz ultrasound the lower the mean flow velocity is, the higher the average turbulent kinetic energy is. This observation is in agreement with the Nusselt number behavior, which experiences a decrease when flow velocity and thereby Reynolds number increase (Fig. 4). This naturally leads to a decrease of the Heat Transfer Enhancement Factor (HTEF) as plotted in Fig. 14.

With 25 kHz sonication, average turbulent kinetic energy values of Fig. 13 are also in good agreement with TKE fields presented in Fig. 11. The highest values are obtained close to the heating wall whatever the mean flow velocity. The higher the mean flow velocity is, the higher the turbulent kinetic energy near the wall is. Nevertheless, whatever the considered flow velocity, the turbulent kinetic energy profile remains unchanged. Indeed, when going far from the heating wall, i.e. at a distance higher than 4 mm (x-axis), a nearly constant value of the average TKE is measured within the bulk. Consequently, 25 kHz ultrasound modifies bulk turbulences whatever the mean flow velocity value, allowing a constant mixing effect within the bulk of the liquid. This could be of great interest for processes that require constant mixing effects of the fluid within the bulk on large flow velocity range.

However, the thermal enhancement increase, with respect to the mean flow velocity increase observed for 25 kHz sonication (Fig. 14), is also consistent with the -. / ratio increase with same operating conditions. Indeed, the increase of the mean flow velocity and associated

Reynolds number induces a decrease of the boundary layer thickness. This phenomenon increases the disruption effect due to cavitation bubbles on the laminar boundary layer. By consequence, heat transfer from the heating wall to the bulk increases leading to a Nusselt number increase. Therefore, it can be expected that when the ratio -. / tends towards one, Nusselt number would not increase anymore with respect to the mean flow velocity. This could explain the asymptotical behavior for experimental curves observed for both Nusselt number (see Fig. 4) and HTEF profiles (see Fig. 14). As a conclusion, it is of great importance to highlight that using 25 kHz ultrasound will be preferable to enable optimum intensification of the transfer process under transitory-turbulent flows for processes involving parietal transfer such as heat exchangers or membrane processes.

One can remark that the laminar boundary layer thickness is also decreasing with respect to Reynolds number increase for 2 MHz sonication. So why is this behavior different from the one observed with 25 kHz sonication at high flow velocity? It is probable that increasing the pressure fluctuations or the flow velocity can affect the wave propagation and then the streaming production at 2 MHz. Thus, the acoustic streaming is weaker for high Reynolds numbers. This induces less intense convective effects and turbulence kinetic energy values for both bulk and laminar boundary layer. When average turbulent kinetic energy curves are compared with respect to frequency for fixed hydrodynamic conditions, it appears that 2 MHz ultrasound produces a higher level of turbulence at low mean flow velocities near the heating wall (Fig. 13). The contrary occurs at high flow velocities. This observation is consistent with thermal results both expressed in terms of Nusselt numbers (Fig. 4) or in terms of Heat Transfer Enhancement Factor (Fig. 14). It is therefore reasonable to make the hypothesis that for a fixed ultrasonic power, acoustic streaming is more efficient to enhance parietal transfer under laminar flow regimes while acoustic cavitation is more relevant under transitory-turbulent flows.

Conclusions

The aim of this work, which belongs to a heat exchangers optimization project, was to study mechanisms involved and responsible of the intensification of heat transfer with two different ultrasound frequencies. One of the objectives was to establish a link between observed acoustically induced heat transfer enhancements and related induced hydrodynamic phenomena depending on the experimental conditions. This work focused on a comparison between the differences highlighted using low (25 kHz) or high (2 MHz) frequency from both thermal and hydrodynamic points of view.

Thermal enhancement by ultrasound, i.e. Nusselt number increase compared to silent condition, is observed for both frequencies. On the one hand, 2 MHz ultrasound frequency produces strong convective effect usually named acoustic streaming that induces turbulence enhancement and global mixing effect. On the second hand, acoustic cavitation produced by 25 kHz ultrasound leads to an increase of turbulences both in the bulk and within the boundary layer. Highest values on turbulent kinetic energy are obtained near the heating wall (i.e. within the boundary layer) which is consistent with the fact that acoustic bubble collapse is more violent near a solid wall. The higher the mean flow velocity is, the higher the thermal enhancement is. This stronger action of 25 kHz ultrasound at high flow velocity could be explained by the enhancement of the effect of cavitation bubbles collapse on the boundary layer. Indeed, as mean flow velocity increases, laminar boundary layer thickness decreases.

This leads to an easier and a more efficient disruption of the boundary layer by bubble collapse and a reduction of the thermal resistance. In order to quantify this influence, a ratio defined as the size of the fluid disturbance area -[m] on the laminar boundary layer thickness .

[m] can be estimated. When this ratio -. / is equal to 1, it means that the length of disturbance area produced by bubbles collapse is equal to boundary layer thickness. This allows to directly transfer physical quantities (heat, concentration…) from the wall to the bulk flow. It is then expected that when the ratio -.

/ is equal or superior to 1, Nusselt number may stop to increase with respect to mean flow velocity for 25 kHz ultrasound regime.

Further studies will prospect the boundary layer disruption explanation proposed here with experiments performed for higher flow velocities and Reynolds numbers. At least, results reported in this paper provide a better understanding of thermal intensification by using low (25 kHz) or high (2 MHz) frequency ultrasound according to main hydrodynamic phenomena induced by acoustic cavitation and acoustic streaming. It is hoped that this work will be helpful to the way to the design of future vibrating heat exchangers.

Fig. 2 :

 2 Fig. 2: Side view of the test section ALTO used for this study.

( 2 Fig. 3 :

 23 Fig. 3: Erosion pit clustering (black part) of the aluminum foil covering the entire 5 th heating

3. 1

 1 Thermal results for 25 kHz and 2 MHz sonicationThermal results obtained in silent conditions and for both ultrasound configurations (25 kHz and 2 MHz) are plotted on Fig.4. Nusselt numbers are given with respect to the fluid flow velocity and corresponding Reynolds number (constant canal thickness). Experimental results clearly demonstrate the thermal enhancement induced by ultrasound when Nusselt numbers are compared to those of silent regime, whatever the incident irradiative frequency. On silent mode, the evolution of the Nusselt number is consistent with the fact that the heat transfer coefficient increases when the liquid flow rate increases and thereby corresponding flow velocity. These results confirm classical results in forced convection for silent regime, i.e.

Fig. 4 :

 4 Fig. 4: Nusselt number vs mean flow velocity and Reynolds number for silent conditions,

Fig. 5 : 3 . 2 . 2

 5322 Fig. 5: Average velocity fields observed under silent conditions, heating block n°5, canal

Fig. 6 :

 6 Fig. 6: Average velocity field observed with 2 MHz irradiation, heating block n°5, canal

Fig. 7 :

 7 Fig. 7: Turbulent kinetic energy fields estimated in front of the 5 th heating block for 6 mean

Fig. 8 .Fig. 8 :Fig. 9 Fig. 9 :

 8899 Fig. 8: Average velocity fields observed with 25 kHz ultrasound sonication, mean flow

Fig. 10 :

 10 Fig. 10: PIV doublet recording a bubble cluster collapse and resulting instantaneous velocity

Fig. 11 also

 11 Fig. 11 also shows the increase of the turbulent kinetic energy values close to the heating wall

Fig. 11 :

 11 Fig. 11: Details of Turbulent Kinetic Energy (TKE) field near the 5 th heating block from 0 to

Fig. 12 :

 12 Fig. 12: Illustration of the influence area of cavitation bubbles on the laminar boundary layer. The theoretical viscous boundary layer δ for silent regime is calculated from [45] (f = 25 kHz, e = 30 mm, Q v = 1075 L/h, Mean flow velocity = 0.099 m/s, Re = 4500, δ = 1.76 mm, P us =

Fig. 13 :

 13 Fig. 13: Average turbulent kinetic energy profile vs canal position on the first 10 millimeters from the 5 th heating block to the bulk (x-axis) under 25 kHz and 2 MHz

  Indeed, the highest Nusselt Number and HTEF values at the lowest liquid flow rate are obtained with 2 MHz ultrasound while the contrary is observed at highest flow velocity under 25 kHz sonication.

Fig. 14 :

 14 Fig. 14: Heat Transfer Enhancement Factor vs Reynolds number and mean flow velocity

Fig. 15 is

 15 Fig. 15 is finally suggested as a qualitative illustration of the relationship between

Fig. 15 :

 15 Fig. 15 : Illustration of the phenomenon involved in heat transfer enhancement by ultrasound

  Additional work will aim at determining causes of the turbulence attenuation with respect to flow velocity and Reynolds number increase in presence of 2 MHz ultrasound observed in this study. The objective would be to determine if this phenomenon is reliable to the increase of the Reynolds number, i.e. to the natural turbulence in the flow in silent regime, or to the increase of the mean flow velocity. A Reynolds similarity study, including canal thickness modification for fixed Reynolds values, and then different flow velocities should help to answer to this question. Finally, a study dealing with global energy balance of thermal enhancement by ultrasound should be done. A comparison between energy efficiency of such a device and passive methods could be of great interest. Nevertheless, this kind of study needs an optimization of the ultrasound technology in accordance with thermal enhancement goals.

  The laser is a Nd-Yag (Litron provided by equals to 16.5 * 16.5 µm 2 . Interrogation windows are 16 * 16 pixels leading to a square of 0.264 * 0.264 mm 2 . Consequently, PIV measurement provides measurement fields of the entire canal thickness (x-axis) on the entire i th heating block height (y-axis).

	These PIV measurements allow us to link thermal phenomenon to hydrodynamic ones that
	take place in front of this i th heating block.
	PIV device gives instantaneous velocity fields from which it is possible to calculate average
	velocity fields from 2000 instantaneous velocity fields. RMS velocity [m/s] fields and
	turbulent kinetic energy (	) [m 2 /s 2 ] on the two dimensions x-y plan for the 2000
	Dantec Dynamics) with two cavities generating 200 mJ of maximum energy for 6-9 ns pulse
	time. The camera is an ImagerProX2 M (provided by LaVision) with a 1600 * 1200 pixels
	spatial resolution. Each of the pixels of the camera measures 7.4 * 7.4 µm 2 . Average fields are
	obtained from 2000 instantaneous fields.

The measurement field size is equal to the size of the canal in front of the i th heating block (30 mm thick by 20 mm high on x-y plan). The camera lens is a Nikon 150 mm, leading to PIV field pixels


Benjamin Fores

Cédric Champollion

Nolwenn Lesparre

Sylvain Pasquet

Aurélie Martin

Frédéric Nguyen

Variability of the water stock dynamics in karst: insights from surface-to-tunnel geophysics

Keywords: Karst, Heterogeneity, Geophysical methods, Unsaturated zone

A hydrogeophysical field experiment was conducted on a karst hydrosystem in the south of France to investigate groundwater transfer and storage variability at a scale of a few hundred meters. A 200-m-long N/S tunnel going through limestone provided the unique opportunity to set up measurements with original configurations inside the unsaturated zone. Three geophysical methods were used: gravimetry, electrical, and seismic. Two-dimensional electrical resistivity and seismic velocity images were retrieved by surrounding the medium with electrodes and geophones, both at the surface and inside the tunnel to improve the sensitivity in depth. This gave information about the weathering state but also about the limestone content and associated porosity characteristics, as the methods are sensitive to distinct properties with different resolution patterns. A time-lapse gravity surfaceto-tunnel profile supplied information on the seasonal water mass changes and its variations along the tunnel. Besides, tracers were injected on each side of the profile from the surface and the restitution was sampled in the tunnel drip flows. A contrasting hydrological behavior was evidenced on each side of the tunnel from temporal gravity measurements and tracing tests. The analysis of the whole dataset allowed for better interpretation of the imaged structures, with different hydrological functioning. This study demonstrates the variability of the karst behavior at the scale of a few hundred meters and the benefits of a multimethod approach coupling hydrological and geophysical measurements. This kind of experiment provides fundamental understanding of systems that cannot be directly observed.

Introduction

A quarter of the world relies on carbonate aquifers for its freshwater needs [START_REF] Ford | Karst hydrogeology and geomorphology[END_REF]. Those aquifers, subject to complex karstification processes, are highly heterogeneous from the pore-scale to the basin scale. They are often seen as compartmentalized, with a buffering epikarst close to the surface [START_REF] Klimchouk | Towards defining, delimiting and classifying epikarst: its origin, processes and variants of geomorphic evolution[END_REF]). The latter can be considered as a more weathered, more porous and more homogeneous zone compared to an underlying massive medium crossed by fractures and/or karstic channels. Given its higher porosity, the epikarst is more permeable and is considered as an active zone favoring infiltration. Karstic aquifers are also laterally heterogeneous with sinkholes acting as infiltration points, spatially variable weathering, localized fractures and extremely transmissive network of drains. Each of these features are inherited from the karst evolution history, which depends on the landscape, climate, vegetation, geological settings, carbonates nature, preexisting fractures (e.g. ghost rocks karstification: [START_REF] Dubois | The process of ghost-rock karstification and its role in the formation of cave systems[END_REF]epikarst weathering: Klimchouk 2004). Such heterogeneity is challenging to assess and to account for in models for both scientific or groundwater resources management studies [START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: a brief review[END_REF], making the comparison and combination of data at different scales very difficult [START_REF] Dausse | Multi-scale assessment of hydrodynamic properties in a karst aquifer (Lez, France)[END_REF]. In karstic environments, boreholes provide very localized information which usually cannot be interpolated at the scale of the catchment, nor a few meters away (e.g. [START_REF] Noushabadi | Influence of the observation scale on permeability estimation at local and regional scales through well tests in a fractured and karstic aquifer (Lez aquifer, southern France)[END_REF]. At the opposite scale, springs studies (geochemical, hydrogram analyses; e.g. [START_REF] Bicalho | A conceptual model for groundwater circulation using isotopes and geochemical tracers coupled with hydrodynamics: a case study of the Lez karst system, France[END_REF][START_REF] Mazzilli | KarstMod: a modelling platform for rainfalldischarge analysis and modelling dedicated to karst systems[END_REF]) integrate information about structures and processes taking place across the entire catchment. Between those two end members, the water storage capacity is likely to strongly fluctuate in space over the watershed [START_REF] Vouillamoz | Localization of saturated karst aquifer with magnetic resonance sounding and resistivity imagery[END_REF][START_REF] Jacob | Time-lapse microgravity surveys reveal water storage heterogeneity of a karst aquifer[END_REF][START_REF] Mazzilli | Contribution of magnetic resonance soundings for characterizing water storage in the unsaturated zone of karst aquifers[END_REF].

Spatialized models at the catchment scale [START_REF] Hartmann | Karst water resources in a changing world: review of hydrological modeling approaches[END_REF] require distributed data to inform the medium permeability and the dynamics of water circulation. In this context, geophysical imaging methods are pertinent tools to identify hydrofacies distribution [START_REF] Al-Fares | Analysis of the karst aquifer by means of a ground penetrating radar (GPR): example of the Lamalou area (Hérault, France)[END_REF][START_REF] Valois | Assessment of doline geometry using geophysics on the Quercy plateau karst (South France)[END_REF][START_REF] Chalikakis | Contribution of geophysical methods to karst-system exploration: an overview[END_REF]Carrière et al. 2013;[START_REF] Pasquet | Geophysical imaging of shallow degassing in a Yellowstone hydrothermal system[END_REF][START_REF] Novitsky | Mapping inherited fractures in the critical zone using seismic anisotropy from circular surveys[END_REF]. Geophysical properties often reflect many properties of the medium such as lithology, the degree of weathering/fracturing, salinity or water content. To grasp a better understanding of complex systems such as karsts, multiple geophysical methods are often colocalized and jointly interpreted (e.g. [START_REF] Gallardo | Characterization of heterogeneous nearsurface materials by joint 2D inversion of dc resistivity and seismic data[END_REF][START_REF] Robert | Assessing the contribution of electrical resistivity tomography (ERT) and self-potential (SP) methods for a water well drilling program in fractured/karstified limestones[END_REF][START_REF] Valois | Assessment of doline geometry using geophysics on the Quercy plateau karst (South France)[END_REF]. Whereas different geophysical methods are often complementary, one still requires additional information to constraint joint interpretation. This information can take the form of prior knowledge in a joint inversion process (e.g. [START_REF] Gallardo | Characterization of heterogeneous nearsurface materials by joint 2D inversion of dc resistivity and seismic data[END_REF], invasive sampling or more recently time-lapse surveys (e.g. [START_REF] Robert | A salt tracer test monitored with surface ERT to detect preferential flow and transport paths in fractured/karstified limestones[END_REF][START_REF] Valois | Application of combined time-lapse seismic refraction and electrical resistivity tomography to the analysis of infiltration and dissolution processes in the epikarst of the Causse du Larzac (France)[END_REF][START_REF] Watlet | Imaging groundwater infiltration dynamics in the karst vadose zone with long-term ERT monitoring[END_REF]. Joint inversion typically requires one to constrain the inverted parameters together using structural prior or petrophysical relationships. Such constraints are in practice difficult to achieve successfully in complex systems such as karst given the variability of the physical properties. In some sites, human-made or natural structures allow one to make multidirectional measurements, with sensors surrounding the medium to better characterize it. They often require specific methodological approaches to account for the uncommon survey geometries, e.g. surface-to-borehole [START_REF] Tsourlos | Measurement and inversion schemes for single borehole-to-surface electrical resistivity tomography surveys[END_REF], surface-to-tunnel [START_REF] Gritto | Surface-to-tunnel seismic tomography studies at Yucca Mountain, Nevada[END_REF][START_REF] Simyrdanis | Surface-to-tunnel electrical resistance tomography measurements[END_REF][START_REF] Lesparre | Electrical resistivity imaging in transmission between surface and underground tunnel for fault characterization[END_REF]), surface-to-cliff (Dussauge-Peisser et al. 2003;[START_REF] Jeannin | Multiconfiguration GPR measurements for geometric fracture characterization in limestone cliffs (Alps)[END_REF]. In cases where existing tunnels are available, they provide a singular surface-to-depth geometry for gravity measurements leading to a direct and unique insight on the water mass variations in the medium in between the surface and depth acquisition points [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF][START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF]. Such configurations allow a considerable improvement of the spatial resolution of the geophysical images in depth [START_REF] Maufroy | Travel time inversion from ground level to gallery: protocol for the characterization of P-wave seismic signature in a fractured-porous Urgonian platform at hectometric scale[END_REF][START_REF] Lesparre | Electrical resistivity imaging in transmission between surface and underground tunnel for fault characterization[END_REF], leading to more reliable interpretation of geophysical images. In karstic systems, time-lapse micro-gravity experiments complete characterization studies by shedding light on the heterogeneity of an aquifer in terms of water storage dynamics [START_REF] Jacob | Time-lapse microgravity surveys reveal water storage heterogeneity of a karst aquifer[END_REF][START_REF] Hinderer | Hybrid gravimetry as a tool to monitor surface and underground mass changes[END_REF]). The latter is informed without ambiguity through time-lapse gravity within a radius of tens to hundreds of meters around the instrument (e.g. [START_REF] Gehman | Estimating specific yield and storage change in an unconfined aquifer using temporal gravity surveys[END_REF][START_REF] Creutzfeldt | Measuring the effect of local water storage changes on in situ gravity observations: case study of the geodetic observatory Wettzell, Germany[END_REF][START_REF] Pfeffer | Evaluating surface and subsurface water storage variations at small time and space scales from relative gravity measurements in semiarid Niger[END_REF][START_REF] Piccolroaz | On the use of spatially distributed, time-lapse microgravity surveys to inform hydrological modeling[END_REF].

In this study, an innovative hydrogeophysical experiment is performed to explore the water storage dynamics heterogeneity at the intermediate scale of a few hundreds of meters in a shallow karst on the Larzac plateau (south of France). This unique experimental dataset provides a rare opportunity to better understand how epikarst heterogeneity and fractures affect water storage dynamics. A tracing experiment is combined with time-lapse gravimetry, electrical resistivity tomography (ERT) and seismic refraction tomography (SRT) to understand how the structural information validates the hypothesis inferred from the observed temporal water changes and vice-versa. The presence of a 200-m-long tunnel allows implementing surface-to-tunnel (S2T) measurements. Observations given by the different data types are analyzed together in light of the geological and hydrological field observations across the watershed and along the tunnel walls.

In the following, the geological and hydro-meteorological context of the experiment are presented. Then, the methodology of tracing tests and geophysical experiments are described: time-lapse S2T gravimetry for storage changes and static S2T SRT and ERT images of the underground structures. Finally, the results are described with a discussion about how the observed water storage dynamic spatial variations might be explained by the presence of structures with different hydrogeological properties. A conceptual model details the role of the identified compartments on the different water transfer processes involved. Finally, the benefits and resolution improvements of such S2T configurations are discussed.

Site description

The experiment takes place along and above the tunnel of Saint-Ferréol, located on the south part of the Campestre karstic plateau (a part of the Larzac karstic plateau) in the south of France. The Campestre plateau's elevation varies between 600 and 900 m asl and has a small area of 36 km 2 covered by fields and pine forests.

Hydrogeological context

The Campestre plateau is mainly made of limestones, dolostones and marls deposited during marine regressions around -150 Myr. The climate is Mediterranean with dry summers, wet winters and springs, while extreme rainfall events happen regularly in the fall. Annual precipitation strongly varies between 1,050 and 1,650 mm a year (Météo-France, ALZON station, years 2011 to 2015). The plateau is highly karstified, which explains the absence of perennial surface stream. An active network of karst conduits (e.g. [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF]) is expected to conduct the water falling on the Campestre plateau towards the "Foux de la Vis" karstic spring, 8 km south-east [START_REF] Bruxelles | Dépôts et altérites des plateaux du Larzac central: Causses de l'Hospitalet et de Campestre (Aveyron, Gard, Hérault)-evolution morphogénétique, conséquences géologiques et implications pour l'aménagement [Deposits and alterites of plateau of Larzac central: causses from l'Hospitalet and de Campestre (Aveyron, Gard, Hérault)-morphogenetic evolution, geological consequences and implications for development[END_REF]). The soil is thin and disappears towards the south where limestone outcrops are common. The volume investigated with geophysical methods (from the tunnel to the surface) is made entirely of Thitonian limestone whose thickness may be greater than 150 m [START_REF] Bruxelles | Dépôts et altérites des plateaux du Larzac central: Causses de l'Hospitalet et de Campestre (Aveyron, Gard, Hérault)-evolution morphogénétique, conséquences géologiques et implications pour l'aménagement [Deposits and alterites of plateau of Larzac central: causses from l'Hospitalet and de Campestre (Aveyron, Gard, Hérault)-morphogenetic evolution, geological consequences and implications for development[END_REF]. Stratification in the tunnel shows a small and constant 10°dip to the east. Fractures observed on the tunnel walls are mainly ~N140-180 with 50-90°dip to the west, which is consistent with the geological map [START_REF] Gèze | Notice explicative de la feuille de Nant à 1/50 000[END_REF].

Tunnel configuration

The tunnel is 220 m long and is oriented N-S, with a section of about 2 m × 2 m (Fig. 1). It is entirely manmade and was dug for commercial purposes at the end of the nineteenth century to connect the bottom of a sinkhole (called a "doline" in this region), easily reachable, and a large cavity at the bottom of a 60 m depth abyss (Fig. 1), to store and mature Roquefort cheese. The entrance of the tunnel is located 30 m beneath the plateau (which is 700 m asl) on the sinkhole side (south). The tunnel presents a descending slope varying between 10 and 15% towards the abyss (north) and leading to the cavity 60 m beneath the plateau. Surface topography shows a slight dip to the SW, with pines to the east and fields to the west.

The tunnel section (see Fig. 2) near the sinkhole has a very regular shape dug in massive limestone. Approaching the abyss, the limestone becomes more and more fractured with marks of runoff on the walls and on the roof. A pile of rubble is followed by a reinforced masonry vault between 160 and 190 m (Fig. 2), highlighting an unstable section. Observations of drops falling in the tunnel provide insights on the water transfers occurring above the tunnel. The following qualitative observations were made during gravity measurements and water sampling for tracer tests. Over the entire distance of the tunnel, the relationship between heavy rainfalls and drip flows is pronounced; however, in the first 50 m from the tunnel entrance located on the sinkhole side, drip flows dry up after a few weeks without important rainfall events. On the other hand, in the last 50 m of the tunnel (abyss side), a different behavior is observed with perennial drip flows subsisting after several months without precipitation at the surface, even though flows are then weaker (Fig. 2). One can also note the existence of a small spring on the wall near the abyss at about 200 m from the tunnel entrance. This spring has a sustained output and shows important concretions. Finally, a few clay pockets (mainly kaolinite) are observed on the tunnel roof and walls, on the second half of the tunnel only and especially close to the abyss [START_REF] Bruxelles | Dépôts et altérites des plateaux du Larzac central: Causses de l'Hospitalet et de Campestre (Aveyron, Gard, Hérault)-evolution morphogénétique, conséquences géologiques et implications pour l'aménagement [Deposits and alterites of plateau of Larzac central: causses from l'Hospitalet and de Campestre (Aveyron, Gard, Hérault)-morphogenetic evolution, geological consequences and implications for development[END_REF].

Methods

Electrical resistivity tomography (ERT)

Electrical resistivity tomography (ERT) was carried in November 2016 with 48 electrodes in the tunnel and 48 electrodes at the surface, with inter-distances of 4 and 5 m, respectively (Fig. 2). Electrical resistivity measurements were acquired using an IRIS Syscal-Pro connected to a switch allowing current injections and measurements in between the tunnel and the surface. Such a configuration helps in constraining the medium properties at depth, notably close to the tunnel [START_REF] Simyrdanis | Surface-to-tunnel electrical resistance tomography measurements[END_REF][START_REF] Lesparre | Electrical resistivity imaging in transmission between surface and underground tunnel for fault characterization[END_REF]. Classical dipole-dipole and gradient protocols were used separately from the surface and from the tunnel. A third protocol was added by injecting current between the surface and the tunnel electrodes following a dipole-dipole AM-BN configuration with 1 injecting and 1 measuring electrode, both in the tunnel and at the surface [START_REF] Bing | Cross-hole resistivity tomography using different electrode configurations[END_REF]. Such a configuration was favored since it provided a higher sensitivity than regular AB-MN dipole-dipole (both injecting electrodes in the tunnel and both measuring electrodes at surface or the contrary) or pole-dipole configurations (one electrode at the infinite; [START_REF] Demirel | Two-dimensional joint inversions of cross-hole resistivity data and resolution analysis of combined arrays[END_REF]. In all, 12,682 data were acquired using the different protocols and 5,156 additional data were collected in a reciprocal mode to assess the data quality. Reciprocal measurements corresponded to a swap of the electrodes used for current injection and voltage measurements during the "normal" acquisition.

All the data sets (dipole-dipole, gradients and transmission) are gathered for the ERT inversion. Only data presenting a repeatability error lower than 1% and a normal to reciprocal difference lower than 5% were selected. Additional filters are applied to keep only data with apparent resistivity in the range [0-40000] Ω.m, leading to a final set of 9,526 data that fulfilled the selection. The complete workflow for data selection, error analysis and weighting is provided in the electronic supplementary material (ESM). The medium is represented using a 3D finite element model built using Netgen and Gmsh softwares [START_REF] Schöberl | NETGEN an advancing front 2D/3D-mesh generator based on abstract rules[END_REF][START_REF] Geuzaine | A three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The mesh of the forward model is unstructured and refined close to the electrodes for a better accuracy of the resistance estimation, while elements are coarser further away from the electrodes allowing a reduced computation time. The mesh includes the superficial topography as well as the abyss and the tunnel. The geometry of the investigated region is defined with a two-dimensional (2D) regular coarse mesh following a vertical plane between the surface and the tunnel. That second mesh is used to define the geometry of the resistivity values sought during the inversion. Thus, the sought resistivity values are defined in 2D and extruded in the threedimensional (3D) mesh for the forward model computations. This allows the drastic reduction of the number of sought values in the inversion, while the estimate of resistivity remains accurate since the forward problem includes the 3D geometry of the site. The measured resistances are inverted using the EIDORS software to obtain an image of the electrical resistivity distribution above the tunnel [START_REF] Polydorides | A MATLAb toolkit for threedimensional electrical impedance tomography: a contribution to the electrical impedance and diffuse optical reconstruction software project[END_REF]. Resistance values for a given distribution of electrical resistivity can then be estimated from the computation of the forward model with resistivity values allocated to each element. A cumulative sensitivity was used as a transparent filter to quantify uncertainty [START_REF] Kemna | Tomographic inversion of complex resistivity: theory and application[END_REF][START_REF] Caterina | A comparison study of different image appraisal tools for electrical resistivity tomography[END_REF].

Seismic refraction tomography (SRT)

Seismic data were acquired in October 2018 using four Geode seismic recorders (Geometrics) and 96 14-Hz vertical-component geophones. A total of 48 geophones were installed at the surface and 48 geophones in the tunnel, both with 4 m spacing. A total of 48 shots were performed about every 8 m, both in the tunnel and at the surface, with a 5-kg sledgehammer swung onto a metal plate. For each shot, seismic data were recorded simultaneously on all 96 geophones with a sampling rate of 0.125 ms and a recording time of 800 ms, so as to obtain 4,608 traces. This acquisition setup allowed for the estimation of travel times in four different configurations:

(1) surface to surface, (2) surface to tunnel, (3) tunnel to tunnel, and (4) tunnel to surface. Here again, the set-up of sensors and the possibility to shoot at depth allows increasing the image quality deeper, down to the tunnel [START_REF] Maufroy | Travel time inversion from ground level to gallery: protocol for the characterization of P-wave seismic signature in a fractured-porous Urgonian platform at hectometric scale[END_REF].

A set of 3,367 travel times were picked manually on the recorded shot gathers, where signal-to-noise ratio is high enough to confidently identify first arrivals. These traveltime observations were then inversed for subsurface pressure-wave velocity structure (Vp) using pyGIMLi refraction tomography inversion module [START_REF] Rücker | pyGIMLi: an open-source library for modelling and inversion in geophysics[END_REF]. In pyGIMLi, the inversion domain is parameterized in 2D as a sheared mesh of constant velocity triangles in which rays are traced using a shortest path algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF][START_REF] Moser | Shortest path calculation of seismic rays[END_REF]. Velocity updates are found in each of the 7,473 cells of the mesh using a generalized Gauss-Newton inversion framework. The program starts with an initial model consisting of a velocity field that increases linearly with depth, and then finds an appropriately smooth update to the model that reduces the difference between predicted and observed travel times. To quantify uncertainty and model the sensitivity, inversions were performed with 135 combinations of starting models and regularization parameters [START_REF] Pasquet | Geophysical imaging of shallow degassing in a Yellowstone hydrothermal system[END_REF]. Top velocity of the starting model ranged between 250 and 750 m/s; bottom velocity varied between 1,000 and 5,000 m/s. All the models presenting a satisfactory fit to the data are then eventually used to create an average final model. The final model has a travel-time mean absolute residuals of 5.3% for a root mean square error of 1.79 ms.

Time-lapse micro-gravimetry

Nine couples of gravity stations were set along the tunnel, starting 50 m after the tunnel entrance with a spacing of 20 m (Fig. 2). For each station in the tunnel, a paired station was set up vertically at the surface on the nearest limestone solid outcrop. Precise positions and heights were obtained using a differential GPS (Trimble R8) outdoors, coupled to laser measurements (Trimble 5600) in the tunnel. Most of the 18 gravity stations (9 at the surface and 9 in the tunnel) were measured at least 3 times a year from December 2013 to September 2015 to sample the hydrological cycle: one measurement after the wet season (December/January), one at the end of the dry season (August/September) and one in between. The repetition of the measurements at the exact same positions and heights were ensured by setting bolts in the limestone outcrops or in the tunnel floor; and by blocking one foot of the CG5 tripod with a brass ring [START_REF] Jacob | Time-lapse microgravity surveys reveal water storage heterogeneity of a karst aquifer[END_REF]. The same CG5 gravimeter (Scintrex Limited, 2006) was used for all measurements. This instrument has a standard field repeatability of 50 nm s -2 corresponding to a water slab variation of 120 mm. The calibration factor of this CG5 instrument has always been kept up to date [START_REF] Jacob | Time-lapse microgravity surveys reveal water storage heterogeneity of a karst aquifer[END_REF][START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF].

All the recommendations to minimize errors associated with gravimeter transport were applied [START_REF] Lederer | Accuracy of the relative gravity measurement[END_REF], i.e. gently handling the gravimeter and settling it 2 h before starting the measurements after its transport by car to the site. CG5 measurements were programmed every 95 s (90-s reading, 5-s pause), following the recommendations of [START_REF] Merlet | Micro-gravity investigations for the LNE watt balance project[END_REF]. The gravity value at one station is the mean of about several measurements of 90 s, with variations depending on the repeatability (up to 15 measurements). All the surface measurements could not be performed each time because of a logging event and field mishaps. Indeed, an N/S strip of pines has been cleared during spring 2014 after the first gravity survey, right at the vertical of the tunnel (Fig. 1), leading to a few data gaps for stations pairs 6 and 7.

The instrumental drift was determined and corrected by repeating the first station of the survey three times a day, following the measurement strategy of [START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF]. The removal of all unwanted gravity temporal changes, including tides, atmospheric pressure or any other nonhydrological temporal variation over 1 day of survey is performed by subtracting the continuous and precise gravity monitoring of a superconducting gravimeter (SG), located 14 km away to the West [START_REF] Fores | Assessing the precision of the iGrav superconducting gravimeter for hydrological models and karstic hydrological process identification[END_REF]. The correction of the topography and karstic voids effects were done using a precise digital elevation model of the site and a scan of the cave.

As stated in the introduction, making gravity differences between two stations vertically aligned provides the density of the slice between the two stations, assuming a onedimensional (1D) model [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF]). After the corrections of earth tides and atmospheric pressure changes (and assuming that there are no surface deformation), variations of S2T gravity over time are expected to be related only to water mass changes above the tunnel. Therefore, there is a straightforward and linear relationship of between gravity S2T changes over time (also called double-difference) and equivalent water height changes in between both measuring points [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF][START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF]. The complete set of equations are available in the supplementary materials, as the novelty of this study lies in the possibility to investigate the lateral variability of S2T measurements thanks to the tunnel.

Tracer experiment

In order to study solute tracers' transport and their spatial variations in the unsaturated zone, a multitracer experiment was carried out on 16 March 2018, after a period of rain. Two tracers were injected simultaneously: 750 g of Fluorescein was injected on the sinkhole side in a rather clayey surface; and 250 g of Sulforhodamine B was injected on the abyss side (Fig. 2). These two tracers are easily detectable (around the μg L -1 ), conservative, little adsorbed by the soil but sensitive to light [START_REF] Käss | Hydrological tracing practice on underground contaminations[END_REF]). Both tracers were diluted in 2.5 L of water and 15 L of water were added to push the tracers in depth below the soil layer, down to the epikarst, and limit their adsorption by the soil. The excitation and emission wavelengths of the Fluorescein are 490 and 515 nm, respectively, and 565 and 586 nm for Sulforhodamine B.

Groundwater that percolated in the tunnel roof was sampled manually by holding vials to gather the drops every 10-20 m (Fig. 2). Samplings were made at specific dropping points or on a roof section, depending on the flow rate. The day of the injection, samplings were performed at least every hour to monitor the background noise level of tracer for concentration analyses. From March 19-22, a daily sampling was carried out. Then, the temporal sampling rate slowed down, with measurements on April 7, April 13 and May 5. All samples were kept in the dark until May 23rd, when they were analyzed using a PerkinElmer LS45 spectrofluorometer.

Results

Electrical resistivity tomography

The electrical resistivity image is displayed only for regions with a high sensitivity (Fig. 3). The electrical resistivity varies on a wide range spanning 6 orders of magnitude, the extreme values being located along the tunnel. At the surface, a region of high resistivity is visible below the sinkhole, where electrodes are located in a scree slope, just above the tunnel entrance. Then, outside of the sinkhole, the shallow subsurface is characterized by low resistivity values (<500 Ω.m) down to a few meters depth, between 50 and 160 m from the tunnel entrance. This conductive region extends down to 20 m below the surface at the vertical around the 70 m landmark. On the north side, the resistivity is higher near the surface, locally above 7,500 Ω.m. However the depth of that highly resistive medium cannot be accurately estimated due to the poor constraint in that region as illustrated by the low sensitivity values.

At depth, the south part of the tunnel shows very high resistivity values, globally above 10 4 Ω.m and reaching locally 10 6 Ω.m. Those high resistive patches are inserted in less resistive zones. Open fractures filled with conductive materials dividing rock in blocks might produce such patterns. However, the distance between the electrodes and the smooth regularization used in the inversion prevent the distinction of such thin and sharp structures. Moreover, the high resistivity values induce a low sensitivity in these regions since the electrical current circulation is very small. Along the north side of the tunnel, the resistivity values are lower than 1,000 Ω.m. A restricted region with a very low resistivity below 50 Ω.m coincides with clay observed on the tunnel walls. Some higher resistivity patches might be related to the region less mechanically stable at the level of the masonry (160-195 m, Fig. 2). On the abyss side, 10 m above the tunnel, the resistivity presents values below 100 Ω.m; nevertheless, the low sensitivity in that region prevents from assessing the dimension of that conductive zone. Thus, it remains difficult to delimit the boundary between the resistive zone below the surface from the conductive one above the tunnel on the north side of the tunnel.

Seismic refraction tomography

The pressure-wave velocity model obtained with SRT is represented on Fig. 4. The S2T configuration allows one to identify and pick travel times from a large number of sources and receivers pairs, thus considerably improving the resolution and coverage of the estimated Vp model. In this model, velocities vary between 500 and 5,000 m/s and highlight three distinct zones within the investigated medium. A first zone, characterized by low velocities (<1,500 m/s), covers the whole profile in the first 2-3 m below the surface. Underneath, intermediate velocities (>1,500 and <3,000 m/s) extend down to 5 m at both ends of the profile, and down to about 25 m in its thickest part at about 110 m away from the sinkhole. Higher velocities (>3,000 m/s) occur deeper and extend down to the

Density profile

The apparent density along the tunnel for each pair of station is computed from the mean of all gravity surveys (Fig. 5) with the error bars showing the minimum and maximum values. The apparent density appears to be constant at 2,530 kg m -3 from the sinkhole to 160 m (masonry in the tunnel) with a sharp and significant increase of +25 kg m -3 to reach 2,555 kg m -3 after 180 m. The porosity range is obtained directly from the apparent density profile assuming a homogeneous grain density of 2,700 kg m -3 (theoretical calcite value). The minimum and maximum porosity is evaluated by assuming the medium saturation to be 0 and 100% respectively. The porosity ranges from 5.5 to 10.0% with an average of 7.8 ± 1.8%, which is realistic and coherent with previous estimations on other karstic basins on the Larzac plateau [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF][START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF].

Water storage changes

The S2T gravity difference (Δ z g) is represented over time for each couple of gravity stations alongside meteorological data (Fig. 6). Gravity variations are coherent with meteorological data as S2T gravity increases after rainfall events and decreases after dry periods because of evapotranspiration but also because of water infiltration beneath the tunnel. The most important rain events take place between August 2014 and January 2015. The Alzon Météo-France station measures more than 800 mm of precipitation between these two surveys (Fig. 6). Globally, Δ z g increases after rainy periods (August 2014 to January 2015; September 2015 to January 2016) and decreases after the dry season (December 2014 to August 2014; January 2015 to May 2015).

A strong spatial variability is also evidenced with much more important changes for stations 1 to 5 (50-130 m from the entrance) than for stations 6 to 9 (150-210 m from the entrance). To highlight the spatial heterogeneity along the tunnel, Fig. 7 shows the gravity double differences ( Δ t z g À Á ) along the tunnel, which correspond to the differences of Δ z g between two successive surveys [START_REF] Jacob | Time-lapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF][START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-todepth gravity measurements[END_REF]. Water storage changes are important on the southernmost point, near the sinkhole, and smoothly decreases until 150 m from the sinkhole entrance. From 150 m to the abyss, there are no significant water storage changes according to the gravity data. Such observations confirm the existence of strong lateral variations in water storage behavior between the surface and the tunnel.

Tracer restitution

The maximum restitution happens about 2 months (67 days) after the injections (last samplings on May 22, 2018) for both tracers, and for both tracers below their injection points (Fig. 8a): 74 m from the sinkhole for the Fluorescein and 162-170 m for B-Sulforhodamine. No significant response is observed in the drips sampled elsewhere in the tunnel. A fluorescein response were also detected a few days after the injection in the sampling points close to the vertical of the injection point of this tracer (30-74 m from the sinkhole entrance, Fig. 8b), while there is no B-Sulforhodamine response at all before May 22nd.

Thereby slow and vertical transfers of about 1 m/day can be estimated. Rapid and vertical transfers, close to the sinkhole only, are suggested by a small fluorescence peak on March 19, 3 days after the injection (Fig. 8b). One cannot completely exclude lateral N/S flows from the tracings, but this is coherent with the fractures observed in the tunnel which dip vertically to 50°W but never N/S. One also cannot exclude that a part of the water flow circumvents the tunnel or is directed westward by fractures.

Discussion

Joint interpretation of the structural images

The density profile, the SRT and ERT images show common features of the sounded karstic area and attest its Fig. 6 a Surface-to-tunnel gravity monitoring for each vertical couple of stations (couple 1 is the southernmost and couple 9 the northernmost). The mean of each profile has been subtracted. b Potential evapotranspiration and precipitation from the Alzon Météo-France station (7 km to the north-east) Fig. 7 Gravity double differences (changes in the surface-to-tunnel gravity difference between each consecutive survey) and equivalent water height changes for each couple of stations vertically aligned heterogeneous nature. Near the topographic surface, resistivity correlates with seismic velocity, especially close to the sinkhole and the abyss. From the sinkhole to ~150 m northward, the low resistivity strip beneath the surface (250 and 500 Ω.m, Fig. 3) coincides seismic velocities lower than 1,000 m/s in the very first meters and 2,000 m/s in the first 10-20 m (Fig. 4). This shallow layer is about 5 m thick on the south edge and reaches 20 m deep 70 m away from the sinkhole. Close to the abyss, higher surface resistivity and surface velocities are present. Such a positive correlation between resistivity and seismic velocity can be related to different weathering stage of the carbonates. Massive limestone presents high P waves velocities and high resistivity. The alteration increases the porosity that reduces P wave velocity, and may be filled with more conductive material (water, organic matter, clay), lowering the resistivity. The thickness variation of this weathered layer is also coherent with the density increase northward according to the S2T gravimetry (Fig. 5). Since the tunnel is dipping progressively to the abyss, the gravity density sharp rise to the north might correspond to a larger proportion of healthy limestone. Following the classical vertical structure of the karst [START_REF] Klimchouk | Towards defining, delimiting and classifying epikarst: its origin, processes and variants of geomorphic evolution[END_REF]), this less dense, porous, weathered subsurface layer could be identified as the so-called epikarst.

Below the epikarst, the limestone is generally expected to be healthy with some fractures or much localized altered areas (the "infiltration zone"). Close to the tunnel roof, seismic velocities are high and homogeneous (3,000-5,000 m/s) with weathering decreasing with depth (Fig. 4). However, ERT shows lateral variations, with a resistive south and a more conductive north (Fig. 3) and several small pockets of extreme resistivity values (very high or very low). Note that positive and negative targets influence differently both geophysical methods, which in addition have different resolution patterns. Highly conductive objects like small fractures filled by clay could have a large impact on the global resistivity even far from the electrodes (notion of short-circuit) but have little influence on the P wave velocity if the volume of clay is small compared to the sampled volume. In a more general way, anomalies larger than the methods resolution are correlated, but smaller anomalies can be uncorrelated. Note that the resolution depends on the wavelengths for the SRT, and electrodes spacing and protocol for the ERT. Then, high values of resistivity might reflect the presence of open fractures in the south, while spots of low resistivity might show the occurrence of fractures clogged by clay in the north. The observation of clay pockets on the tunnel walls around the abyss is another argument supporting this interpretation. Around 160-180 m close to the abyss a small anomaly of low velocity (~2,000 m/s) and very low resistivity (less than 100 Ω.m) is seen just above the tunnel. At this exact location, the tunnel is reinforced with masonry because of damages and instabilities which explains those geophysical observations. The studied vertical profile exhibits two zones divided at about 130 m with different features that might induce the observed temporal dynamics from the tracer test, the drip flows and the gravity. All the complementary information gathered on that particular site provides clues to sketch the distinct hydrological processes susceptible to occur so closely in karstic environments.

Hydrological processes on the southern side S2T gravity measurements (Figs. 6 and7) show important water mass changes on the southern side. The water storage increases during the "wet" seasons (up to 450 nm s -2 , i.e. 500 mm of water). A porous reservoir, that may be an intermittent perched aquifer, must exist to support this groundwater transient storage. Such a reservoir might correspond to the weathered shallow region revealed by the geophysical images. This subsurface layer, identified as the epikarst is thicker in the south (between 50-90 m, Fig. 3), where gravity changes are the highest (Fig. 7).

On this south part, the drip flows observed in the tunnel depend on the climatic context, drying up in summer and immediately restarting after significant rainfalls. The tracing tests have demonstrated the presence of vertical flows with two residence times: a very fast one (a few days) and a slow transfer (2 months). Then one can infer that the south epikarst reservoir acts as a buffer distributing the water by slow water percolation through the underlying less permeable infiltration zone (5-20 m deep to the tunnel and below). The presence of open fractures is deduced from the ERT and SRT images, which is coherent with the drip reactivity to rainfalls and the tracing test. In terms of processes, seasonal storage is located in a shallow weathered porous limestone with a transfer through vertical fractures. Fast transfer is enhanced when the reservoir reaches water saturation Hydrological processes on the northern side Perennial flows are observed on the northern side, with a variable rate depending on the season. The tracing shows only vertical and slow transfer (more than one month). Seasonal water storage from gravity observations is almost insignificant (±100 mm of water) over 3 years. This is in agreement with the geophysical images showing almost no weathered epikarst and an infiltration zone composed of healthy limestone with fractures partially clogged with clay. The saturation in the infiltration zone must remain constant and near the saturation to allow perennial flow and the absence of transient storage. Evidence of high saturation in the infiltration zone has already been described in the Larzac plateau from magnetic resonance soundings experiment [START_REF] Mazzilli | Contribution of magnetic resonance soundings for characterizing water storage in the unsaturated zone of karst aquifers[END_REF].

In terms of processes, both storage and transfers differ from the south side. In the north a vertical transfer is assumed, thanks to a piston mechanism [START_REF] Aquilina | Water storage and transfer in the epikarst of karstic systems during high flow periods[END_REF]) that explains the absence of gravity changes (all the incoming water pushing an equivalent quantity), the reactivity of drip flows to rainfall and the slow tracer restitution. Clays in fractures reduce the available space for water and could explain the triggering of this piston effect at low water contents. Clays may also prevent a good drainage, so the reactivity to rain is lower than on the sinkhole side and flows are perennial. A summarizing scheme of the identified hydrological processes and structures conciliates all observations gathered (Fig. 9).

Conclusions

The site configuration with a surface-to-tunnel geometry provides a unique way to observe and understand karst heterogeneity at the scale of few hundreds of meters. Geophysical images reveals that the epikarst, an equivalent porous subsurface layer acting like a reservoir, has a variable thickness. Time-lapse S2T gravity suggests that this reservoir can store up to 80% of the rainfall, distributing the water through open fractures in the underlying healthy limestone. The surface reservoir thins up and almost disappears where the water storage is divided by a factor of 5. In depth, limestone is massive as shown by the high seismic velocities, while the varying resistivity is interpreted by distinct clay content. While open fractures seem to dominate below the epikarst reservoir, a substantial clay content may clog the vertical fractures where the reservoir thins. This is also consistent with tracer tests showing only vertical flows, with a fast and slow transport underneath the thick epikarst, and only a slow transport underneath the massive limestone where clay is imaged.

S-SE Sinkhole

N-NW Abyss

In terms of methodology, the three geophysical methods provide complementary information for a joint interpretation of the underground structures. The different properties imaged (electrical resistivity, p wave velocities and density) are often consistent. This is the case for the weathered epikarst with a high porosity filled with water, and consequently low density, p wave velocity and resistivity. However, geophysical properties might reveal apparent contradictions where the resolution power or the sensitivity of the methods strongly differ. That is why a multimethod approach is valuable but needs to account for the specific resolution loss and uncertainty of each method. This study also demonstrates the importance of combining geophysical imaging with time-lapse measurements and tracer tests relating the observed structures, to solute transport and water storage processes to provide tangible data on actual conceptual models that would remain hypothesis otherwise. This type of experiment, which is seldomly available, should be encouraged and repeated where possible in order to validate conceptual models and to provide fundamental understanding of systems inherently difficult if not impossible to observe directly.

Fig. 1

 1 Fig. 1 Aerial image of the site located on the Campestre plateau showing the tunnel, abyss and sinkhole positions. The colorscale represents the elevation of the tunnel and abyss. Topographic isolines have a height interval of 5 m

  Fig. 2 S-N vertical cross-section of the tunnel connecting the sinkhole and the abyss, showing the configuration of all the geophysical experiments (gravimetric stations, electrodes, geophones) and hydrological experiments (tracings and drip flows observations)

Fig. 3 a

 3 Fig. 3 a Electrical resistivity tomography and b sensitivity analysis. The black crosses indicate the electrodes and the red circles record the positions of the 9 gravity pairs of stations. The cross section is blanked where sensitivity is considered too low to be reliable

Fig. 4 a

 4 Fig. 4 a Seismic refraction tomography and b travel-time residuals. Black and blue crosses indicate geophones and sources, respectively. Contour lines at 1,500 and 30,00 m/s are represented with thick and thin black lines, respectively

Fig. 8 a

 8 Fig. 8 a Normalized concentrations of Fluorescein and B-Sulforhodamine along the tunnel on 22 May 2018, which is the day of maximum restitution. b Fluorescein concentration over time for the sampling point vertically aligned with the injection (74 m from the sinkhole entrance) along with daily rainfalls

Fig. 9

 9 Fig. 9 Interpretation of the site based on all the geophysical and hydrological experiments. The karstic compartments (epikarst/ infiltration zone) and the clay content are deduced from the relations between the ERT and SRT images. The reservoir role and flows path are deduced from the time-lapse gravity and tracing experiments

Acknowledgements The Scintrex CG5 relative gravimeter was loaned by the Gravity Mobile facility (GMOB) of RESIF-INSU. We thank N. Denchik, S. Baudin, E. Doerflinger, J. Chéry, O. Khairoun, S. Furst, A. Maréchal, G. Mainsant, A-K. Cooke and A. Fort for their help during the field experiments. We are grateful to Eric Cotteux (ITAP IRSTEA) for the Spectrometer and tracings equipment, and to Thomas Kremer for Magnetic Resonance Soundings attempts. The seismic equipment was provided by the CRITEX ANR-11-EQPX-0011 project and the METIS laboratory at Sorbonne Université. We would also like to acknowledge the farmer who let us access the cave. All the data used in this study (geophysics, tracers and meteorological data) are available on the OSU OREME observatory website (http://data.oreme.org/gek/home).

Funding Gravity surveys were part of the OSU OREME observatory funded by the Institut de Recherche pour le Développement (IRD) and the Institut National des Sciences de l'Univers (INSU) and of the National observation system (SNO) H+. N. Lesparre benefited from a BEcome a WAlloon REsearcher fellowship fund co-financed by the Department of Research Programs of the Federation Wallonia -Brussels and the COFUND program of the European Union that both founded the project SUITE4D.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1007/s10040-021-02365-5.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.