Marion Gaudeau

Mengxue Zhang

Michaël Tatoulian

Camille Lescot

Stéphanie Ognier

Michael Tatoulian

Fast carbonylation reaction from CO 2 using plasma gas/liquid microreactors for radiolabeling applications

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Positron emission tomography (PET) is a powerful non-invasive imaging technique that uses radiolabeled compounds as molecular probes to image biological processes in vivo, for research and clinical diagnostic 1 . These imaging probes are labeled with short-lived radioisotopes (e.g., 18 F, 11 C, 13 N, 15 O). Its rising use in recent years has increased the demand for the development of new PET tracers and that of new synthetic methods for these tracers. Carbon-11 is a valuable positron-emitting radionuclide, because carbon is the main constituent of all biological molecules 2 . However, due to the very short half-life time of 11 C (t1/2 = 20.4 min), the synthesis of 11 C-labeled radiotracers is not without challenge: chemists must find rapid synthesis routes so that the labeled molecule can be produced and isolated within two halflives, i.e. 40 min for 11 C.

As carbon-11 mainly comes from 11 CO2 generated in biomedical cyclotrons [START_REF] Rotstein | [END_REF] , it is crucial to have a method allowing the rapid radiolabeling of molecules directly from 11 CO2. However, considering the high stability of CO2, it is necessary to envisage the transformation in situ of 11 CO2 into a more reactive species. Among the existing carbon-11 radiolabeling methods, carbonylation reactions with 11 CO are of great interest since it leads to a wide range of functional groups (amides, ketones, ureas, esters, lactones, carboxylic acids, …) having interesting pharmaceutical properties [4][5][6][7][8][9][10][11][12] . These functional groups are thus present in many biologically interesting compounds, and more specifically the amide function obtained by aminocarbonylation [13][14][15][16] .

Our objective is to develop an original process of introducing a carbon atom through aminocarbonylation using CO2 as carbon source. The development will be performed with carbon-12 with the vision of future use in radiolabeling with carbon-11, so taking into account the time constraint and the security issues related to the radioactive nature of carbon-11.

Considering these parameters, carbonylation reactions usually require the use of intensified processes to ensure a high reaction rate. As for the gas-liquid reaction, high pressure autoclave reactors generating high partial pressures of CO have proven to be effective for 11 C radiolabeling 17 , but these devices are expensive and require special safety precautions. Another way to intensify the process is to use microreactors [18][19][20][21][22][23][24][25][26][27][28][29][30][31] . Indeed, the large surface area to volume ratio of microfluidic systems improves the gas-liquid mass transfer between the flow of gas containing CO and the liquid phase reagent, leading to efficient carbonylation reactions with a short residence time 32 . Improved transfer in microreactors also induces milder reaction conditions compared to conventional methods, and their small volume is perfectly suited to the manipulation of submicromolar amounts of radioactive reagents.

Concerning the reduction of CO2 to CO, carbon dioxide is kinetically and thermodynamically stable [START_REF]Chemical Bonds and Bonds Energy -2nd Edition[END_REF] , and therefore requires significant energy or powerful catalysts for its chemical transformation. Several methods can be used for the splitting of CO2 into CO (chemical, electrochemical, photochemical) 10-15,34- 45 but plasma has been demonstrated as an effective and clean way [START_REF] Fridman | Plasma Chemistry[END_REF][START_REF] Aerts | [END_REF] . Indeed, under electric discharges, the highly energetic species (excited neutrals, radicals, metastables, photons…) and electrons generated (with energy of [1-10] eV [START_REF] Boenig | Fundamentals of plasma chemistry and technology[END_REF] ) in the plasma can break the C-O bond (2.9 eV/molecule at 300 K) and transform CO2 molecules. Compared to classical routes, the conversion of CO2 by plasma does not generate waste and is an "atom-economical" approach without the use of catalysts and reducing agents. 16,[START_REF] Laitar | [END_REF][50][51] Despite the complex plasma chemistry of the reaction, the splitting of CO2 can be summarized as follows:

𝐶𝑂 2 → 𝐶𝑂 + 1 2 𝑂 2 (1)
∆𝐻 𝑅 = 2.9 𝑒𝑉 𝑚𝑜𝑙𝑒𝑐𝑢𝑙𝑒 ⁄ 𝑎𝑡 300𝐾 = 280 𝑘𝐽/𝑚𝑜𝑙 Among the different plasma systems, non-thermal (or cold) atmospheric plasma appears as a convenient and innovative method, since it is easy-to-handle and is operated under mild conditions (room temperature, atmospheric pressures) 52 . One of the most exploited non-thermal plasma configurations for CO2 splitting is the dielectric barrier discharge (DBD), for which one or two dielectric layers are inserted between the electrodes. It prevents the formation of electric arcs detrimental to the device integrity and allows homogenous discharges in the discharge gap. Simple to use and to scale up, the DBD configuration is a promising technique for the CO2 conversion to CO at industrial scales 53 .

Thus, the combination of microfluidics and plasma seems to be an ideal solution for carrying out the rapid carbonylation reactions required for 11 C radiolabeling for PET.

Here we propose a new method to synthesize carbonylated molecules starting from CO2 directly, associating microfluidics and plasma, within a very short residence time (several minutes) and under mild conditions. The proposed system is composed of two microreactors in a continuous 2-step process for the carbonylation reaction starting from CO2: a first plasma microreactor for the CO2 splitting and a second gas-liquid microreactor for the carbonylation reaction. To evaluate the efficiency of this innovative process combining plasma and microfluidics, aminocarbonylation of p-iodoanisole and nhexylamine in the presence of a palladium catalyst to form Nhexyl-4-methoxybenzamide (1) has been chosen as proof-ofconcept because of its high selectivity and efficiency to trap CO. Furthermore the reaction mixture is easy to analyse by 1 H NMR 16 . This process would subsequently be applied to radiolabeling with carbon-11. To the best of our knowledge, such a fast and efficient process, starting from CO2 leading to amide in mild conditions, has never been studied. Scheme 1 presents the general reaction scheme.

Scheme 1: General reaction scheme

Experimental

The study of the decomposition of CO2 into CO by plasma and the study of aminocarbonylation reaction have firstly been carried out separately in two separate microreactors, then the two microreactors were coupled in line for the 2 steps synthesis.

Microreactors

Two glass microreactors with identical geometry have been used in this study: one with metallic electrodes for the plasma reactor and the other without electrodes for the carbonylation reaction. The patented reactor geometry has been previously developed in our group [START_REF] Wengler | Biphasic Plasma Microreactor and Method of Using the Same[END_REF] and demonstrated efficient for gasliquid reactions [START_REF] Wengler | [END_REF][56][57] . The reactor has been designed to provide a stabilized gas-liquid parallel flow in 1-m long channel. The gas flows in the central channel, while the liquid phase flows in the two side channels along the gas phase (Figure 1a and Figure 3). The internal volume of the gas phase is 10 µL and that of the liquid phase is 25 µL.

Aminocarbonylation Material

The liquid phase was prepared by mixing in the following order: 0.0125 mmol of Pd(dba)2 [Bis(dibenzylideneacetone)palladium(0)] (5 mol%), 0.025 mmol of Xantphos [4,5-bis(diphenylphosphino)-9,9dimethylxanthene] (10 mol%), 0.25 mmol of p-iodoanisole (1 eq), 3 mL of Dioxane, 0.5 mmol of n-hexylamine (2 eq), 0.5 mmol of trimethylamine (2 eq).

Fluidic set up

A biflow reactor was used for the aminocarbonylation reaction (Figure 1b). Details about the connectors assembling are given in Electronic Supplementary Information (ESI). Pure CO from commercial CO cylinder has been used for this part and was managed by means of mass flow controller (Bronkhorst EL-FLOW PRESTIGE FG-201CV) in the range of 0.014 -4 sccm. The liquid flows were controlled by a syringe pump (kdScientific Legato 180) in the range of 7.5 -20 µL/min. The liquid phase was introduced through the liquid inlet and the gas phase through the gas inlet. In both cases, the respective flowrate of each phase determines its residence time in the chip. One of the reactor outlets was sealed while the other outlet of the microreactor was connected to a vial for liquid collection.

TEMPERATURE CONTROL

An infrared lamp (Bioseb Lab Instruments) was used to heat the microreactor, and the temperature was controlled by two thermocouples: one was fixed below the microreactor (Figure 2b), and the other above. The temperature inside the microreactor was determined supposing that the temperature changes linearly between the top and the bottom of the reactor. During the experiments, the typical temperature difference between the bottom and the top of the reactor was about 10 °C. The operating temperature was controlled in the range of 20 -80 °C by adapting the distance from the IR lamp.

Flow characterisation

The chip was placed under a macroscope (Leica Z16 APO, Germany) onto which was mounted a CCD camera (PixeLINK PL-B781U) for flow visualisation, with the chip being backilluminated by a diffuse LED lamp. Figure 3 shows the gas-liquid flow inside microchannel through the CCD camera. A stable gasliquid flow was observed in all the tested conditions.

Liquid analysis

The liquid was collected at the outlet in a vial during 45 min. Conversion of p-iodoanisole into the amide 1 was determined by crude 1 H NMR analysis (CDCl3) without internal standard, by comparing to the remained p-iodoanisole quantity * (see also ESI).

Each experiment was conducted at least twice, and the presented results correspond to the average of the measured substrate conversions. An error of 5% is appreciated on NMR results. The limit of detection of NMR spectroscopy is 5-10%, that is why the result >90% was given when no more trace of starting material was detected.

Cleaning

Prior to the reaction, the reactor was initially rinsed by a stable dioxane-CO flow during 30 min with the reactor temperature well controlled. Then, the dioxane syringe was replaced by a syringe containing the reaction liquid phase, and the reaction was carried out, collecting the solution from the outlet of the microreactor in a new vial (synthesized product vial). After the reaction, the reactor was thoroughly rinsed again with a stable dioxane-CO flow for 15 minutes. for Ar, to one of the inlets of the reactor, while the other was sealed. One of the reactor outlets was sealed while the other outlet was connected to a gas sampling bulb. The outlet of the bulb was connected to a bubble flowmeter for a regular flowrate control during the experiment, ensuring the absence of leakage.

Electrical set up

A sine wave signal (AC) was sent from a function generator (ELC, GF467F, 5MHz) to a voltage amplifier (TREK 10/40A high voltage amplifier ×1000) which multiplies the input voltage by 1000. This high voltage was applied to the microreactor and a 3.11 nF capacitor in series and was monitored with an oscilloscope (PicoScope 5000 Series) through the 1000:1 outlet of the amplifier. The capacitor was used to obtain discharge electric charge necessary to calculate the discharge power, and the capacitance voltage was measured via a low voltage proble (Teledyne LeCroy PP024 500 MHz 10:1). The frequency used was between 1000 and 3000 Hz, and the high voltage was between 9 and 18 kVpp. Figure 5 shows the diagram of the experimental setup.

Power assessment

The measurement of the discharge power supplied to the microreactor was determined using the Lissajous method and is described in ESI. The corresponding standard deviation of the calculated power lies below 3%.

Gas analysis

After 20 minutes' reaction (necessary time for the system to achieve the steady state), the gas leaving the microreactor was collected from the bulb with a syringe and was injected into an Agilent Technologies 490 Micro GC equipped with two columns (column PoraPLOT U, 10 meters long with heated injector, and column Molesieve 5A with backflush capability and heated injector) and micro-machined Thermal Conductivity Detector (TCD). Calibrations were made for CO/Ar and CO2/Ar mixtures distinctly with mixtures of known compositions. Each experiment was conducted three times.

Temperature measurement

A thermocouple was positioned under the reactor (blue cable on Figure 6) to measure the temperature at the surface of the bottom PDMS layer. The temperature inside microchannels was supposed to be equal to this temperature, thus neglecting the thermal resistances of the PDMS coating and the glass slide relatively to the convection boundary layer. The temperature was recorded once stabilized.

Coupling

The two microreactors previously described were set in line to carry out the aminocarbonylation from CO2. In the first microreactor, the decomposition of CO2 into CO by plasma was achieved as explained above, and the outlet gas mixture was introduced directly into the second microreactor along with the liquid phase for aminocarbonylation reaction.

The second microreactor was heated to the appropriate temperature with IR lamp and thermocouples, and the gasliquid flow was observed by the CCD camera, as explained before. The cleaning protocol was similar to that used for aminocarbonylation reaction prior and after the reaction. The liquid was collected at the outlet of the second reactor during 45 min. The collected liquid was analysed by 1 H NMR spectroscopy.

Results & discussion

Aminocarbonylation from commercial CO Aminocarbonylation reaction was first studied using commercially available CO in order to find the optimized conditions leading to a maximum conversion of substrate into amide 1.

Optimization of the temperature

In a first time, the reaction temperature was optimized, starting from 80°C which was the reaction temperature in batch process according to C. Lescot et al 16 . All the reactions were performed with a liquid flowrate of 10 µL/min, corresponding to a liquid residence time of 60 seconds, and a CO flowrate of 1.5 sccm. Under these conditions, the CO was in excess, the molar ratio between CO and the substrate being 80. Results are shown in Table 1. A full conversion was obtained for temperatures above 50°C (Table 1, entries 1 -4). Reducing the reaction temperature to 40°C and 30 °C resulted in conversions of 81% and 76% respectively (Table 1, entries 5 and 6). Finally, 56% conversion was obtained at room temperature (Table 1, entry 7).

The optimized temperature was then found to be 50°C with a full conversion after 60 s residence time (Table 1, entry 4). To compare with the batch process (details in ESI), only 41% conversion was obtained for a reaction time of 45 min at 50°C. The better result obtained in microfluidic conditions can be explained by the intensified gas/liquid transfer thanks to the very large surface area to volume ratio. In all that follows, the temperature was set to 50°C.

Optimization of the liquid residence time

The liquid residence time is a key parameter in this study, as the aminocarbonylation reaction takes place in the liquid phase. The liquid residence time is calculated by dividing the volume of the liquid phase in the microreactor by the liquid flowrate. The residence time of the liquid was varied by changing the liquid flowrate. The CO flowrate was adjusted to keep the stable gas/liquid flow inside microchannels.

Starting from previously optimized parameters (Table 2, entry 1), the liquid residence time was reduced to 40 s by increasing the liquid flowrate to 15 µL/min, and a full conversion was still obtained (Table 2, entry 2). Reducing further the liquid residence time to 30 s resulted in a 62% conversion (Table 2, entry 3). This low conversion may be explained by the residence time which is too short to allow the reaction to reach completion. Hence, the optimized liquid residence time was 40 s corresponding to a liquid flowrate of 15 µL/min, and has been kept thereafter.

The CO flowrate was slightly reduced to 2.5 sccm, leading to full conversion. With these conditions, a 95% isolated yield of amide 1 was obtained.

Study of the CO volume concentration in the gas phase

Then, the influence of CO volume concentration in the gas phase was studied in order to be both in conditions close to the actual coupling conditions where the conversion of CO2 to CO by plasma is low, but also close to radiolabeling conditions where 11 CO2 is diluted in low concentration in a carrier gas. Ar/CO mixtures of various CO volume concentrations in the gas phase (100%, 5.6%, 1.1%, and 0.6%) were injected while the total gas flowrate was maintained at 2.5 sccm and the liquid flowrate was 15 µL/min. These CO volume concentrations corresponded to CO to substrate molar flow ratio FCO/FSUBSTRATE of 89, 5, 1 and 0.5 respectively. Calculation details are given in ESI. The results are shown in Table 3. The reaction with previous optimized parameters corresponds to a CO volume concentration of 100% and a CO to substrate molar flow ratio FCO/FSUBSTRATE of 89 (Table 3, entry 1). The results indicate that the conversion decreases by decreasing the CO volume concentration, because of the decreasing quantity of CO entering the liquid phase. Despite the CO excess, the conversion for FCO/FSUBSTRATE = 5 was not complete, and can be explained by slower kinetics, either because the mass transfer is slower or because the intrinsic kinetics of the reaction are slower. On the other hand, the result at 1.1% CO volume concentration (Table 3, entry 3) was very promising for the future, because despite the low conversion of CO2 into CO by plasma, a good conversion of the substrate into the amide 1 can be expected for coupling.

Influence of the liquid residence time in equimolar conditions

Keeping a CO volume concentration of 1.1% (FCO/FSUBSTRATE = 1), the influence of the liquid residence time on the substrate conversion was studied by varying the liquid flowrate (3 rd column, Table 4). The total gas flowrate was also adjusted in order to ensure a stable gas-liquid flow (4 th column, Table 4), while the CO volume concentration was fixed at 1.1% in the CO/Ar mixture. CO was supplemented with argon (6 th column, Table 4) to reach the total gas flowrate. Calculation details about flowrates, and the estimation of pressure drop (5 th column, Table 4) and CO partial pressure (6 th column, Table 4) are given in ESI.

Contrary to what we could expect, conversion decreases with the liquid residence time. This unexpected result can be explained by the decrease, even moderate, of the CO partial pressure. Indeed, to increase the liquid residence time from 40 to 80 seconds, the gas and liquid flowrates were decreased, leading to a decrease of the inlet pressure, and therefore the CO partial pressure at the reactor inlet. For example, the partial pressure of CO at the reactor inlet decreases by roughly 25% when the liquid residence time was multiplied by a factor 2. Hence, this decrease of CO partial pressure could explain the decrease of substrate conversion, as less CO was transferred to the liquid phase. Those results suggest that the CO diffusion at the gas/liquid interface can be a hindering step for reactions under low pressure conditions, contrary to the work of Jensen et al 60 Comparing Table 2 and Table 4, it can be concluded that the influence of the liquid residence time on the substrate conversion depends on the percentage of CO in the gas mixture: if the CO is in excess in regard to the substrate, the conversion increases with the residence time (Table 12), on the other hand, if the CO is not abundant supply, the increase of the residence time can negatively affect the conversion because of the influence of the gas flowrate on the CO inlet partial pressure (Table 4).

Aminocarbonylation from CO resulting from the splitting of CO2 by DBD plasma

Determination of plasma parameters

In the previous part, it has been shown that, using a gas-liquid microreactor, it is possible to obtain good conversions with CO volume concentrations in the gas phase as low as 1%. The objective of this second part is to find operating conditions allowing the production of such CO volume concentrations from CO2 splitting using a plasma DBD microreactor. In this part, the liquid and gas flowrates are fixed at 15 µL/min and 2.5 sccm respectively.

The decomposition of CO2 into CO by DBD plasma was studied with CO2/Ar mixtures of 10/90 and 5/95 volume ratios. Indeed, Ramakers and al. showed that addition of Ar led to higher CO2 conversions, and the effect was more pronounced for Ar fractions above 70% 58 .

For each gas mixture, various plasma parameters (voltage, frequency) were used to vary Specific Energy Input (SEI), defined as the ratio of the discharge power to the gas flow rate:

𝑆𝐸𝐼( 𝐽 𝑚𝐿 ) = 𝑃𝑜𝑤𝑒𝑟 (𝑊) 𝑡𝑜𝑡𝑎𝑙 𝑓𝑙𝑜𝑤𝑟𝑎𝑡𝑒 ( 𝑚𝐿 𝑚𝑖𝑛 ) × 60( 𝑠 𝑚𝑖𝑛 )
The specific energy input is generally considered as the characteristic parameter of the plasma discharge.

The influence of the SEI on the reactor temperature, the CO volume concentrations at the outlet and the CO2 conversion were studied for each gas composition.

Temperature

The reactor temperature was monitored with the thermocouple. It can be observed that the reactor temperature increased linearly with SEI for a given gas composition, while the temperature never exceeded 60 °C. The temperature is moderate because of fast dissipation of the low amount of heat generated by the plasma thanks to the high surface to volume ratio of the microreactor.

CO2 conversion

The CO2 conversion is defined as the change in CO2 concentration over the initial CO2 concentration.

𝑋 𝐶𝑂2 𝑎𝑏𝑠 (%) = %𝐶𝑂 2 𝐼𝑁𝐿𝐸𝑇 -%𝐶𝑂 2 𝑂𝑈𝑇𝐿𝐸𝑇 %𝐶𝑂 2 𝐼𝑁𝐿𝐸𝑇 × 100
Where %CO2 INLET and %CO2 OUTLET represent the volume concentration of CO2 (in %) before applying plasma discharge and the volume concentration of CO2 (in %) at the outlet of the reactor after 20 min run, respectively. %CO OUTLET represents the volume concentration (in %) of CO produced at the outlet of the reactor.

In Figure 8, CO2 conversion is plotted according to SEI. For both gas mixtures, increasing SEI led to an increase in the conversion of CO2. Increasing input power could effectively enhance the electrical field, electron temperature, electron density and gas temperature in the discharge, which contributes to the improvement of CO2 conversion. However, the CO2 conversion with SEI was higher for the 5% CO2 mixture in comparison with the 10% CO2 mixture. This is what is generally reported in the literature. 58 There are 2 possible explanations: firstly, the more argon there is, the more the dielectric strength of the gas mixture decreases, which makes it possible to benefit a greater amount of energy available for the dissociation of CO2; secondly, the charge and energy transferred from excited Ar species to CO2 molecules increase with a higher Ar content and there was a higher chance for the CO2 molecule to collide with radicals or excited Ar.

CO volume concentration

In Figure 9, the CO volume concentration (in %) produced by the decomposition of CO2 is plotted as a function of the SEI. We can observe that the CO volume concentration increased when SEI increased for the gas mixture with 10% CO2 in Ar, but remained almost constant for the gas mixture with 5% CO2 in Ar. Although the increase of the Ar ratio has a positive impact on the conversion of CO2 as we observed previously, it can also have a limiting effect as we can see here. Indeed, we can suppose that the charge transfer process between the Ar + or Ar 2+ ions and CO2 (Ar + + CO2 → Ar + CO2 + ), followed by the dissociative electron-ion recombination of the CO2 + ions (CO2 + + e -→ CO + O), is limited if the amount of CO2 decreases, and so Ar + or Ar 2+ ions are less likely to meet CO2 molecules to produce CO. A maximum percentage of CO of 1.6% ±0.4%, corresponding to a CO yield of 16%, was obtained for 10% CO2 in Ar with a frequency of 3 kHz and a voltage of 12 kVpp. In these conditions, the conversion of CO2 was equal to 17.0% ±5.7%, that is approximately equal to the CO yield and assume that the production of other carbon containing products is negligible.

The temperature was around 42°C, and the discharge power was around 300 mW. The CO volume concentration was determined by means of a calibration previously done on the µGC using very low CO flowrates to have a wide range of CO volume concentrations, which can imply uncertainties on the flowrates. This could explain the high uncertainties of the results here.

In the previous study of aminocarbonylation from commercial CO, we obtained 53% conversion for a CO volume concentration of 1.1%. So, the maximum CO volume concentration of 1.6% ±0.4% obtained from the decomposition of CO2 by plasma should be enough to carry out successfully the aminocarbonylation reaction in line.

Splitting of CO2 by DBD plasma coupled with aminocarbonylation in line

The two microreactors were assembled (Figure 10) to carry out the aminocarbonylation reaction with CO synthesized in situ from decomposition of CO2 by plasma. The first microreactor was dedicated to the decomposition of CO2 into CO by plasma, and the experimental conditions were a total gas flowrate of 2.5 sccm composed by CO2/Ar mixture of 10/90 volume ratio, 3 kHz and 12 kVpp. The out-flowing gas moved into the second microreactor, in which the aminocarbonylation takes place, with a liquid flowrate of 15 µL/min at 50°C. This reaction was performed twice, and 70% conversion of the substrate into amide 1 was obtained each time. This result is closed to conversions obtained using commercial CO at similar low CO concentrations (Table 3, entry 2 & 3), the liquid residence time being 40 s for both cases. For the first time, aminocarbonylation was performed from CO2 using a 2-step microfluidic synthesis, within a total residence time shorter than 1 min.

Conclusion and outlooks

In this work, amide 1 was synthesized from CO2 with high yield, under mild conditions, thanks to an innovative process combining microfluidics and plasma. This innovative process could also be used to carry out other carbonylation reactions, and could be transferred in carbon-11 chemistry for carbon-11 radiolabeling for PET. The combination of plasma and microfluidics allows to reach a very short total process time and the radiolabeling process is no longer limited by the very short half-life of carbon 11, while keeping milder conditions compared to conventional methods. The 11 C-labeledcarbonylative molecules could be synthesized, purified and analysed in less than 40 min without making any compromise between radiochemical yield (RCY) and radioactive decay. Such a process of carbonylation starting from CO2, leading to this high substrate conversion with such mild conditions, is unprecedented in the literature and exceeds the limitations encountered so far. In addition to the radiolabeling application, the synthesis in-situ of CO is a real advantage for safety constraints since the CO is not being directly handled.

Figure 1 :

 1 Figure 1: a) Top-view of the channel. b) Picture of the glass microreactor

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: a) Picture of the entire installation; b) Microreactor being back-illuminated by a diffuse LED lamp, with thermocouple below (thermocouple above does not appear on this picture)

Figure 5 :Figure 6 :

 56 Figure 5: Experimental set-up for CO2 splitting reactions

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Graph of the temperature as a function of the SEI

Figure 10 :Table 5 :

 105 Figure 10: Experimental set-up for coupling synthesis

Table 1 :

 1 where reactions are conducted with high CO pressures (2.8 -13.8 bars) that suggest that the oxidative addition was the rate-limiting step below 120°C. Optimization of temperature, for a 10 µL/min liquid flowrate and 1.5 sccm CO flowrate † Standard deviation calculated from the different measured substrate conversions obtained for the same experimental conditions

	Entry	Temperature (°C) [error +/-1°C]	Conversion of p-iodoanisole into amide 1 (%)	Standard deviation (%) †
	1	80		>90		0
	2	70		>90		0
	3	60		>90		0
	4	50		>90		0
	5	40		81		± 14
	6	30		76		± 8
	7	20		56		± 8
	Entry	Liquid flowrate (µL/min)	Liquid residence time (s)	CO flowrate (sccm)	Conversion of p-iodoanisole into amide 1 (%)	Standard deviation (%)
	1	10	60	1.5	>90		0
	2	15	40	3.0	>90		0
	3	20	30	4.0	62		± 4

Table 2 :

 2 Optimization of the liquid residence time, at 50°C

	Entry	CO volume concentration (%)		Conversion (%)	Standard deviation (%)
	1	100	89	90	± 7
	2	5.6	5	61	± 7
	3	1.1	1	53	± 10
	4	0.6	0.5	14	± 3

Table 3 :

 3

	Entry	Liquid residence time (s)	Liquid flowrate (µL/min)	Total gas flowrate (sccm)	Inlet experimental partial pressure of CO (mbar)	Conversion (%)	Standard deviation (%)
	1	40	15	2.5	21	53	± 10
	2	60	10	1.8	17	33	± 3
	3	80	7.5	1.3	16	30	± 11

Study of CO volume concentration with a total gas flowrate of 2.5 sccm gas residence time < 1s), and 15 µL/min as liquid flowrate (liquid residence time = 40 s), at 50°C

Table 4 :

 4 Influence of the liquid residence time in equimolar conditions at 50°C, CO volume concentration = 1.1%

Acknowledgements

This work has received the support of "Institut Pierre-Gilles de Gennes" (laboratoire d'excellence, "Investissements d'avenir" program ANR-10-IDEX-0001-02 PSL and ANR-10-LABX-31; équipement d'excellence, "Investissements d'avenir" program ANR-10-EQPX-34).

Conflicts of interest

There are no conflicts to declare.

Notes and references

where A is the area of the peak of methoxy group of amide 1 (at 3.81 ppm) and B is the area of the peak of methoxy group of amide 1 (at 3.81 ppm) + the area of the peak of methoxy group of p-iodoanisole (at 3.75 ppm)


Simone Di Prima
email: sdiprima@uniss.it

Thierry Winiarski

Rafael Angulo-Jaramillo

Ryan D Stewart

Mirko Castellini

Majdi R Abou Najm

Domenico Ventrella

Mario Pirastru

Filippo Giadrossich

Giorgio Capello

Marcella Biddoccu

Laurent Lassabatere

Detecting infiltrated water and preferential flow pathways through time-lapse groundpenetrating radar surveys

Keywords: nonuniform flow, GPR, water infiltration, wetting zone, infiltrometer, non-Newtonian fluid

 We combined time-lapse GPR surveys with automatized single-ring infiltration tests.  Preferential flow (PF) paths were detected from the 360-degree view of the wetted zones.  Results revealed a link between lithological heterogeneity and the occurrence of PF paths.  Results revealed a link between shrub root system and the occurrence of PF paths.  The approach offers a simplified way to detect preferential flow pathways in the field.

Introduction

Infiltration basins are increasingly used in urban areas for flood protection, groundwater recharge, and stormwater disposal [START_REF] Fletcher | SUDS, LID, BMPs, WSUD and more -The evolution and application of terminology surrounding urban drainage[END_REF]. To ensure sufficient infiltration capacity, these structures are usually established over highly permeable and strongly heterogeneous soils.

Water treatment can also occur when mobile contaminants become filtered and sorbed by the soil.

However, the filtration capacity of infiltration basins can be strongly affected by the occurrence of preferential flow (PF) pathways [START_REF] Ben Slimene | The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flowa numerical study[END_REF], making it important to detect the presence of PF, its effect on infiltration, and related risks for groundwater contamination.

Developing non-invasive and easily replicable procedures to study infiltration processes and the occurrence of PF pathways is crucial for assessing hydraulic functioning of infiltration basins and other soil-based green infrastructure.

Currently, PF pathways are identified either by using invasive techniques (e.g., Abou [START_REF] Abou Najm | New method for the characterization of three-dimensional preferential flow paths in the field[END_REF][START_REF] Stewart | Nondestructive Quantification of Macropore Volume using Shear-Thinning Fluid[END_REF] or by direct observation of disturbed soil profiles. Dye-staining, using compounds such as uranine, brilliant blue FCF, or methylene blue, can reveal PF pathways (e.g., [START_REF] Alaoui | A Preferential Flow Effects on Infiltration and Runoff in Grassland and Forest Soils[END_REF][START_REF] Kodešová | Using Dye Tracer for Visualization of Preferential Flow at Macro-and Microscales[END_REF][START_REF] Kung | Improved Solute-Sampling Protocol in a Sandy Vadose Zone Using Ground-Penetrating Radar[END_REF][START_REF] Luo | Influence of Root Distribution on Preferential Flow in Deciduous and Coniferous Forest Soils[END_REF][START_REF] Weiler | An experimental tracer study of the role of macropores in infiltration in grassland soils[END_REF]. The dyes can be sprinkled onto the soil surface [START_REF] Gerke | Preferential flow mechanisms identified from staining experiments in forested hillslopes: Preferential Flow Mechanisms Identified from Staining Experiments[END_REF] or injected during infiltration experiments [START_REF] Rockware | Preferential Flow Patterns in Paddy Fields Using a Dye Tracer[END_REF], for instance with the help of infiltrometers [START_REF] Cey | Field study of macropore flow processes using tension infiltration of a dye tracer in partially saturated soils[END_REF]. However, an adequate observation of deep dye patterns requires sectioning the soil along different vertical or horizontal planes, which necessitates extensive effort and destroys the soil in the process. Moreover, dye-staining often only reveals preferential pathways that are directly connected to the source [START_REF] Beven | Macropores and water flow in soils revisited: REVIEW[END_REF]. Other investigations have used different tracers, such as 18 O-enriched water or bromide, to get an insight on PF (e.g., [START_REF] Ahuja | Macropore Transport of a Surface-Applied Bromide Tracer: Model Evaluation and Refinement[END_REF][START_REF] Angulo-Jaramillo | Field measurement of soil surface hydraulic properties by disc and ring infiltrometers: A review and recent developments[END_REF][START_REF] Jardine | Transport of Inorganic and Natural Organic Tracers Through an Isolated Pedon in a Forest Watershed[END_REF][START_REF] Köhne | Spatial and temporal dynamics of preferential bromide movement towards a tile drain[END_REF]. These approaches require the collection of multiple samples in time and space to get information on tracer distribution within the porous medium, thus providing only discrete information about the tracer concentrations.

Ground-penetrating radar (GPR) surveys constitute a valid alternative to these traditional methods for characterizing PF paths. For the purposes of GPR, the soil may be thought of as an ensemble of multiple homogeneous regions with contrasting electromagnetic properties (e.g., dielectric permittivity, electrical conductivity), that are separated by interfaces. A GPR system sends into the soil a pulse of high-frequency electromagnetic waves (radiowaves) which oscillate near a particular frequency. When radiowaves come into contact with an interface, portions of the incoming radiowaves are reflected, transmitted and refracted. A GPR system thus allows the user to identify these interfaces and thereby gain information about subsurface structures.

Water flow within the unsaturated zone causes a variation of the dielectric contrast between interfaces that alters the reflection coefficient [START_REF] Truss | Imaging rainfall drainage within the Miami oolitic limestone using high-resolution time-lapse ground-penetrating radar: IMAGING DRAINAGE USING TIME-LAPSE GPR[END_REF]. These variations are manifested as amplitude changes. Based on this principle, time-lapse three-dimensional (3D) groundpenetrating radar surveys can be used to obtain differenced radargrams from pre-and post-wetting surveys [START_REF] Birken | Use of four-dimensional ground penetrating radar and advanced visualization methods to determine subsurface fluid migration[END_REF], which has been used to investigate links between non-uniform flow and heterogeneous soil structures. For instance, [START_REF] Kung | Improved Solute-Sampling Protocol in a Sandy Vadose Zone Using Ground-Penetrating Radar[END_REF] first used the GPR method to detect deep soil layers with textural discontinuities that triggered PF paths on a sandy soil. This technology allowed [START_REF] Harari | Ground-penetrating radar (GPR) for imaging stratigraphic features and groundwater in sand dunes[END_REF] to observe non-uniform downward water movement along more active pathways in sand dunes. [START_REF] Truss | Imaging rainfall drainage within the Miami oolitic limestone using high-resolution time-lapse ground-penetrating radar: IMAGING DRAINAGE USING TIME-LAPSE GPR[END_REF] collected two-dimensional (2D) timelapse GPR surveys before, during, and after rainfall events. The 2D time-lapse surveys were subsequently used to create a 3D GPR dataset from which these authors derived the geometry of PF paths and wetting fronts. [START_REF] Guo | Subsurface lateral preferential flow network revealed by time-lapse groundpenetrating radar in a hillslope[END_REF] combined trench infiltration tests with time-lapse GPR surveys along parallel transects to study the subsurface hydrologic dynamic of a complex lateral PF network in a hillslope. [START_REF] Allroggen | 4D ground-penetrating radar during a plot scale dye tracer experiment[END_REF] studied the potential of time-lapse GPR surveys in conjunction with rainfall simulations and brilliant blue dye injection to improve the interpretation of soil cross-sections and wetted patterns photographed after excavating the site. [START_REF] Guo | Pairing geophysical techniques improves understanding of the near-surface Critical Zone: Visualization of preferential routing of stemflow along coarse roots[END_REF] combined time-lapse GPR surveys with another non-invasive geophysical technique, electrical resistivity tomography (ERT), to trace stemflow infiltration and identify root-derived preferential flow beneath an American beech tree.

Similar approaches were adopted by other investigations for non-invasive monitoring of temporal and spatial distribution of water flow and PF detection (e.g., [START_REF] Birken | Use of four-dimensional ground penetrating radar and advanced visualization methods to determine subsurface fluid migration[END_REF][START_REF] Guo | Preferential flow through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology[END_REF][START_REF] Jackisch | Form and function in hillslope hydrology: in situ imaging and characterization of flow-relevant structures[END_REF][START_REF] Klenk | Monitoring infiltration processes with high-resolution surface-based Ground-Penetrating Radar[END_REF][START_REF] Trinks | Monitoring water flow in the unsaturated zone using georadar[END_REF]. Most of these investigations were conducted at the hillslope scale along transect survey lines, given that PF incidence and connectivity need to be characterized at appropriate scales. However, non-uniform and preferential flow phenomena may be induced by various types of soil heterogeneity (e.g., sealed soil, multi-porosity and permeability systems, porous media of contrasting materials), and can emerge at the point scale [START_REF] Angulo-Jaramillo | Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales[END_REF]. In such cases, infiltration measurements combined with time-lapse GPR surveys may allow for 3D visualization of active preferential flowpaths. However, obtaining precise 3D imaging of soil wetting at small spatial scale still remains challenging, as it requires high accuracy in GPR position during the repeated surveys [START_REF] Allroggen | 4D ground-penetrating radar during a plot scale dye tracer experiment[END_REF].

In this investigation, we combined time-lapse GPR surveys with automatized single-ring infiltration experiments to gain insight on the occurrence of PF at two experimental sites located in Lyon, France. The sites were chosen to represent two different conditions that lead to the establishment of preferential flows: lithological heterogeneity and deeply developed root systems.

The present study aimed to provide experimental evidence on PF processes, with the specific objective of evaluating the ability of time-lapse GPR surveys to identify, at a small-plot scale (i.e., 100 × 100 cm), the incidence and extent of preferential flow. To achieve this objective, we used a non-invasive and easily replicable procedure to visualize the 3D wetting behaviors under ponding experiments.

Material and methods

Experimental areas and sampling scheme

The first field site was a stormwater infiltration basin, named Django Reinhardt basin, located in Chassieu in the eastern suburbs of Lyon, France. A detailed description of the experimental area can be found in [START_REF] Goutaland | Hydrostratigraphic Characterization of Glaciofluvial Deposits Underlying an Infiltration Basin Using Ground Penetrating Radar[END_REF] and [START_REF] Winiarski | The impact of stormwater on a soil profile in an infiltration basin[END_REF]. The infiltration basin was constructed above a heterogeneous glaciofluvial deposit by mixing together the upper 50-80 cm of the soil profile. The glaciofluvial deposit was composed of four main lithofacies: i) the uppermost mixture of the soil matrix and gravel, ii) a mixture of the soil matrix and gravel with a bimodal particle size distribution that formed most of the deposit below 50-80 cm, iii) large lenses of sand, and iv) smaller lenses of matrix-free gravel [START_REF] Ben Slimene | The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flowa numerical study[END_REF][START_REF] Goutaland | Sedimentary and hydraulic characterization of a heterogeneous glaciofluvial deposit: Application to the modeling of unsaturated flow[END_REF].

Previous investigations identified how PF pathways affect groundwater quality at the studied basin, especially with regards to the transfer of pollutants (e.g., heavy metals, organics) carried by infiltrating water [START_REF] Coutinho | Vadose Zone Heterogeneity Effect on Unsaturated Water Flow Modeling at Meso-Scale[END_REF][START_REF] Lassabatere | Effect of the settlement of sediments on water infiltration in two urban infiltration basins[END_REF][START_REF] Lassabatere | Concomitant Zn-Cd and Pb retention in a carbonated fluvio-glacial deposit under both static and dynamic conditions[END_REF][START_REF] Lassabatere | Retention of Three Heavy Metals (Zn, Pb, and Cd) in a Calcareous Soil Controlled by the Modification of Flow with Geotextiles[END_REF]. Previous modeling investigations have reported that lithological heterogeneities in the basin may induce PF, thus increasing water flux rates by more than a factor of ten compared to surface infiltration rates [START_REF] Ben Slimene | The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flowa numerical study[END_REF].

For the second site, the trials were carried out in an open-air flume located in an experimental garden within the La DOUA scientific campus in the municipality of Villeurbanne (France). The sampled area was colonized by ray grass and a few shrub plants with more developed root systems.

The upper horizon was a mixture of matrix and gravel, 20 cm thick, under which lay a dense horizon acting as restrictive layer. The two horizons appear quite homogeneous without any lithological heterogeneity apart from layering. At the DOUA site, we determined the soil bulk density from 24 undisturbed soil cores (~ 100 cm 3 ) collected at different depths from 0 to 50 cm (Table 1). In contrast to the first site, which was dedicated to the study of the effect of lithological heterogeneity on flow pattern, this second site was used to detect preferential flow due to plant roots in a layered profile.

Two locations (referred to as A and B) were sampled at the Django Reinhardt infiltration basin and one (referred to as DOUA) at the experimental site of "Parc de La DOUA". At each site, a GPR grid (100 cm × 100 cm), consisting of six horizontal (1-6) and six vertical (A-F) parallel survey lines with 20 cm intervals between them, was established using white/red signaling tape (Figure 1a-b). The vertical lines were oriented along the north-south axis. Within each grid, an infiltration test was conducted. Two GPR surveys were carried out just before and then 20 min after the infiltration test. A total of 72 (3 sites × 2 GPR surveys × 12 survey lines) time-depth cross-sections (radargrams) were collected by moving the antenna along the survey lines following the sampling scheme reported in Figure 1d. Radar depth penetration and spatial resolution both vary with antenna frequency [START_REF] Guo | Preferential flow through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology[END_REF], so we elected to use a GSSI (Geophysical Survey System Inc., Salem, NH) SIR 3000 system with a 900-MHz antenna as the best compromise between high resolution and adequate penetration for the spatial scale of our wetting experiments.

Automated ponding experiments

In this investigation, the use of dyes at the Django Reinhardt basin was restricted for operational reasons. We therefore infiltrated only water at this location. In contrast, we carried out our infiltration experiment at the DOUA site using a solution that was prepared by adding 1 g L -1 of dry Xanthan gum powder and 1 g L -1 of brilliant blue dye (E133) to beakers of deionized water and allowing them to mix with magnetic stirrers. The xanthan gum created a shear-thinning viscous solution that was expected to fill preferential pathways with limited infiltration into the soil matrix, and thus reveal complex geometries or macropore networks in highly heterogeneous soils [START_REF] Abou Najm | Non-Newtonian Fluids in Action: Revisiting Hydraulic Conductivity and Pore Size Distribution of Porous Media[END_REF][START_REF] Atallah | Characterization of synthetic porous media using non-Newtonian fluids: experimental evidence[END_REF][START_REF] Stewart | Nondestructive Quantification of Macropore Volume using Shear-Thinning Fluid[END_REF]. The Xanthan solution was used to boost preferential flow in the root system of the shrub plants.

At each site we used the automated single-ring infiltrometer proposed by [START_REF] Concialdi | An open-source instrumentation package for intensive soil hydraulic characterization[END_REF] to infiltrate a volume of water corresponding to a total cumulative infiltration of 28 cm, based on an infiltration surface of 94 cm 2 (Figure 1a-b). We used a stainless steel ring with a 15-cm inner diameter inserted shallowly into the soil according to the Beerkan procedure for single-ring infiltration experiments [START_REF] Lassabatere | Beerkan estimation of soil transfer parameters through infiltration experiments-BEST[END_REF]. To setup the experiment, we followed the procedure described in Di Prima (2015) and Di [START_REF] Di Prima | Testing a new automated single ring infiltrometer for Beerkan infiltration experiments[END_REF]. We firstly positioned a plastic film on the soil surface inside the ring and applied a small water head of 0.2-0.3 cm, depending on the surface roughness. The infiltrometer was positioned inside the ring and regulated in height so that the base was in contact with the ponded water. The Mariotte bottle was filled with the fluid, and then activated through lifting the piston, at which time data acquisition began. Finally, the infiltration experiment started when the plastic film was removed. At the DOUA site, the infiltration surface was established around the root system of hawthorn shrub after removing its shoot system (Figure 1c). At the end of the infiltration experiment, we excavated the soil to expose the dyed patterns.

The automated procedure proposed by [START_REF] Concialdi | An open-source instrumentation package for intensive soil hydraulic characterization[END_REF] to treat the transducer output was subsequently applied to determine the cumulative infiltration curve. A video tutorial showing the field procedure and the data processing can be viewed online [START_REF] Di Prima | An open source instrumentation package for intensive soil hydraulic characterization[END_REF]. The use of the automatic infiltrometers in this protocol simplified the field setup, while also reducing the amount of effort needed to characterize water infiltration. During the infiltration tests the GPR antenna was maintained in a static position close to the ring with the SIR 3000 sets for point data collection (Figure 1a-b). We recorded a single scan (i.e., one complete reflected wave from transmission to reception) every one minute in order to monitor the time evolution of infiltrated water. The first ten scans, representative of the pre-wetting condition, were recorded before starting the infiltration process. The acquisition was stopped twenty minutes after the end of the test, and a picture of the infiltration bulb at the infiltration surface was taken. Then we carried out the GPR survey on the established grid. These radargrams were combined with the infiltrometer data. This coupling of information allowed us to monitor movement of water from the infiltration surface in real-time.

GPR data processing and solid modeling

A GPR display typically consists of a radargram, namely a collection of time-series returns, giving a time-depth cross-sections [START_REF] Al-Nuaimy | Automatic target detection in GPR data[END_REF]. We processed each collected radargram using the GSSI Radan 7 software. The process consisted of distance normalization, a static time shift to align direct ground wave arrival to 0 ns, and exponential scaling to compensate GPR energy attenuation with propagation depth. In accordance with [START_REF] Truss | Imaging rainfall drainage within the Miami oolitic limestone using high-resolution time-lapse ground-penetrating radar: IMAGING DRAINAGE USING TIME-LAPSE GPR[END_REF], the same gain curve (established using the pre-wetting data) was applied to all data acquired within a timelapse series. The distance normalization allowed us to obtain an equal distance between the 20 cm marks taken along the survey line intersections.

To determine the appropriate velocity to convert two-way travel times into actual depths we relied on previous GPR-based sedimentological studies performed on the site [START_REF] Goutaland | Sedimentary and hydraulic characterization of a heterogeneous glaciofluvial deposit: Application to the modeling of unsaturated flow[END_REF][START_REF] Goutaland | Hydrostratigraphic Characterization of Glaciofluvial Deposits Underlying an Infiltration Basin Using Ground Penetrating Radar[END_REF]. Those authors performed a Common-Mid-Point (CMP) acquisition to determine velocity profiles at two nearby profiles. The velocity profiles were generally uniform, with values on the order of 0.1 m ns -1 (Figure 2a). In addition, the authors confirmed their estimated velocities by excavating two trenches and measuring the real distances between soil surface and several interfaces detected by the GPR, which resulted in mean values of 0.100 and 0.103 m ns -1 for the two trenches [START_REF] Goutaland | Hydrostratigraphic Characterization of Glaciofluvial Deposits Underlying an Infiltration Basin Using Ground Penetrating Radar[END_REF]. Based on these results, we assumed a uniform wave velocity in our analysis.

For each grid area, the amplitude values, G, from the 12 radargrams were combined to build a unique XYZG dataset, with the X, Y, and Z values representing location coordinates (easting, northing, and elevation). The elevation, Z, was determined from the wave time, given by GPR acquisition and considering the value of wave velocity mentioned above. We therefore obtained three pre-wetting and three post-wetting XYZG datasets, one for each site (Figure 3).

Next, we created three other datasets based on absolute differences between pre-and postwetting G values. Indeed, water flowing within the unsaturated zone caused a variation of the dielectric contrast between layers that altered the reflection coefficient [START_REF] Truss | Imaging rainfall drainage within the Miami oolitic limestone using high-resolution time-lapse ground-penetrating radar: IMAGING DRAINAGE USING TIME-LAPSE GPR[END_REF]. These variations were manifested as amplitude changes. Therefore, the differenced datasets were constructed to highlight the amplitude fluctuations between repeated GPR radargrams collected over the same survey lines before and after the infiltration tests (e.g., [START_REF] Guo | Preferential flow through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology[END_REF][START_REF] Holden | Hydrological connectivity of soil pipes determined by ground-penetrating radar tracer detection[END_REF]. We assumed that the absolute difference between the pre-and post-wetting amplitude values were negatively correlated with the soil water pressure head [START_REF] Birken | Use of four-dimensional ground penetrating radar and advanced visualization methods to determine subsurface fluid migration[END_REF]. We therefore expected strongest reflection differences in correspondence to the saturated soil bulb that extended outwards from the infiltration surface, and decreasing reflection differences in the partially saturated (transmission) zone that surrounded the wetting bulb. In this zone, the water pressure head decreases as the wetting front moves away due to soil capillarity from the saturated bulb [START_REF] Angulo-Jaramillo | Infiltration Measurements for Soil Hydraulic Characterization[END_REF]. For each site, we considered all differential amplitude values greater than 1.5 standard deviations of the mean absolute amplitude for each differenced XYZG dataset to represent substantial changes in reflection due to infiltrated water. A similar criterion considering amplitudes greater than two standard deviations from the mean was used by [START_REF] Guo | Preferential flow through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology[END_REF] to identify preferential flow through saprock. Note that we used a less strict criterion because we were interested in identifying the entire wetted soil volume.

The interpretation of differenced radargrams may be complicated by the occurrence of spurious amplitude differences due to geometrical mismatches between measurements. In addition, infiltrated water causes an increase in bulk dielectric constant, which in turn increases the travel time and lowers the reflection position, resulting in net time shifts [START_REF] Truss | Imaging rainfall drainage within the Miami oolitic limestone using high-resolution time-lapse ground-penetrating radar: IMAGING DRAINAGE USING TIME-LAPSE GPR[END_REF]. To alleviate these problems, we smoothed the radargrams by applying a moving-average filter with a window of 7 vertical time samples and 5 horizontal traces. This procedure allowed us to reduce the geometric mismatch between time-lapse survey lines, yielding lower differences between corresponding preand post-wetting amplitude values. A similar approach was also adopted by [START_REF] Guo | Subsurface lateral preferential flow network revealed by time-lapse groundpenetrating radar in a hillslope[END_REF], who corrected the mismatch in the Z direction by stacking successive samples on each trace using their mean value.

In this investigation, we checked the efficiency of the moving-average filter by comparing the absolute differences between the pre-and post-wetting amplitudes of the processed and unprocessed data. Figure 2b depicts an example involving two survey lines at Site A. In this example, only Survey line 2 crosses the infiltration bulb. Consequently, for this survey line the amplitude values substantially increased after the infiltration of water. That is the reason why the magnitudes scored approximately 0.0318 for Survey Line 2 (impacted by the infiltration bulb) versus 0.0185 for Survey Line 1 (less impacted by the infiltration bulb). However, in both cases, the moving-average filter reduced the difference in amplitude, with a shift from 0.0328 to 0.0318 for Survey Line 2, and 0.0188 to 0.0185 for Survey Line 1 (Figure 2b). In both cases, the mismatch correction was effective, yielding smaller differences between pre-and post-wetting. Based on this finding, we were confident that this process reduced the difference in amplitude that resulted from the time shift induced by the change of velocity with wetting, while keeping the difference in amplitude that resulted from the change in water content. We used the RockWorks 17 software (RockWare, Inc., 2015) and the inverse-distance anisotropic modeling method to perform a 3D interpolation of the differenced datasets. This type of directional search can improve the interpolation of voxel (volumetric picture element) values that lie between data point clusters, and it has been applied in different geophysical investigations (e.g., [START_REF] Attwa | An integrative approach for preliminary environmental engineering investigations amidst reclaiming desert-land: a case study at East Nile Delta, Egypt[END_REF][START_REF] Lange-Athinodorou | The sacred canals of the Temple of Bastet at Bubastis (Egypt): New findings from geomorphological investigations and Electrical Resistivity Tomography (ERT)[END_REF][START_REF] Longo | Prospecting for clay minerals within volcanic successions: Application of electrical resistivity tomography to characterise bentonite deposits in northern Sardinia (Italy)[END_REF]. Here, the method assigned a voxel node value on the basis of the weighted average of neighboring data points looking in each 90-degree sector around the node. The amplitude value of each data point was weighted according the inverse of its distance from the voxel node, taken to a power of 2 (inverse-distance squared).

Results

Determination of the 3D representation of wetted zones

The GPR determination of the wetted zone is illustrated for the Django Reinhardt site. The absolute differences between the pre-and post-wetting amplitude values ranged from 0 to 0.394 for Site A and from 0 to 0.293 for Site B. According to our stated criterion (i.e., amplitude values greater than 1.5 standard deviations), the estimated thresholds used to filter the differenced XYZG datasets were equal to 0.04 (Site A) and 0.025 (Site B), which corresponded to cumulative frequencies of 0.92 and 0.88 (Figure 4). Thus, only 8 to 12% of the amplitude values measured by the radargrams showed substantial changes between pre-and post-wetting.

After filtering, 3D interpolation was used to visualize and quantify the wetted soil volumes associated with each infiltration experiment. The 3D diagrams clearly demarcated the dimensions and shapes of the wetted zones created during both experiments (Figures 5a and5f). Then, we compared these demarcated wetted zones with the observations of the infiltration bulbs at the soil surface at the end of the experiment. Figure 6 illustrates the comparison between the observed infiltration bulb at the surface (Figure 6a, delineated on the picture) and the expansion of the wetted zone in the 3D diagram at the soil surface (Figure 6b), based on results from Site A. Similar results were obtained for Site B. The correspondence between the two shapes constituted a first encouraging signal of the validity of the 3D interpolation procedure to determine the wetted pattern at the soil surface. However, dye-staining results obtained at the DOUA site better supported the use of GPR to identify 3D wetting patterns, as discussed in the next section.

Both experiments had hemispherical-shaped wetted zones in the upper soil horizon that extended to depths of approximately 20 cm (Figure 5: Site A in subpanels c-d and Site B in subpanels h-j).

In this hemispherical zone, infiltration was mainly controlled by the matrix, resulting in a uniform distribution of water. Below this upper layer, the wetting patterns became vertically elongated and irregularly shaped. As discussed in the next section, these patterns were assumed to be associated with PF paths [START_REF] Regalado | Simplified Method to Estimate the Green-Ampt Wetting Front Suction and Soil Sorptivity with the Philip-Dunne Falling-Head Permeameter[END_REF]. The estimated volumes of the wetted zones were 48.1 L (Site A) and 72.8 L (Site B), with the second volume characterized by lower reflection changes (Figure 7). Given that identical volumes of water infiltrated during the tests at both locations, we argue that the expansion of the wetting front due to capillary forces was more relevant at Site B, yielding a larger volume of wetted soil with a smaller degree of wetting [START_REF] Moncef | Approach for predicting the wetting front depth beneath a surface point source: theory and numerical aspect[END_REF].

Detection of preferential flow pathways

At the DOUA site, dyed zones were observed after excavating the soil (Figure 8). The use of time-lapse GPR surveys in conjunction with staining-dye tracer allowed us to verify the 3D interpolation procedure. The GPR data resolved water perching above a shallow restrictive layer (demarcated zone Z1 in Figure 8a-c), which was confirmed by the concentration of dye at that interface (Figure 8c). Below this layer, the 3D diagram showed elongated and irregularly shaped wetting patterns (demarcated zones Z2 and Z3 in Figure 8a-b), which corresponded with roots present at 30 and 40 cm depths (demarcated zones Z2 and Z3 in Figure 8c). The GPR data and dye staining patterns therefore both indicated that water infiltrated through this restrictive layer mainly following the plant root system, revealing that roots act as important flow pathways at the studied site.

At the Django Reinhardt basin, the use of dyes was restricted for operational reasons; thus, the main PF paths were firstly identified from the 360-degree view of the wetted zones. Then, horizontal and vertical slices were extracted from the 3D diagrams to identify zones with strong reflection differences corresponding to patches of wetted soil located below the main wetting front (Figure 5: Site A in subpanels c-e and Site B in subpanels h-j; also shown in a video included as supplementary material). The slices were then overlapped with the pre-wetting radargrams to identify the source of spatial heterogeneity that triggered PF. In both sites, subsurface coarse gravels with low capillarity, and thus low unsaturated hydraulic conductivity [START_REF] Goutaland | Sedimentary and hydraulic characterization of a heterogeneous glaciofluvial deposit: Application to the modeling of unsaturated flow[END_REF], acted as capillary barriers that concentrated flow into relatively narrow pathways (Figure 9), in a process commonly referred to as funneled flow [START_REF] Kung | Improved Solute-Sampling Protocol in a Sandy Vadose Zone Using Ground-Penetrating Radar[END_REF]. In Site A, the 3D diagram showed that flow channeling occurred at a depth of 22 cm. In the underlying horizons the soil had finer soil textures and likely greater capillary forces, which encouraged lateral spreading of the wetting front. In Site B, the 3D diagram showed that flow channeling occurred at a depth of 43 cm. The wetting front exceeded the maximum soil depth sampled by GPR surveys in Site B (70 cm). As a result, the wetted soil volume was likely underestimated at this location.

Real time monitoring of infiltration at surface (GPR punctual acquisition)

We also analyzed the ability of GPR punctual acquisition to provide real-time monitoring of infiltration processes from the soil surface. Our first finding here was that no clear information is provided when the wetting bulb is too small to be detected by the antenna. That was the case of the location B of the Django Reinhardt site. The wetted zone here expanded mainly in the southeast direction, far from the position of the antenna (which was located northwest side of the ring).

Consequently, the point GPR measurements collected during the infiltration test at that site did not provide clear evidence of water flow. However, we measured substantial reflection changes below the antenna position during the infiltration test at Site A, where initial infiltration rates were relatively low before starting to increase ~2 minutes after the start of the experiment (Figure 10). The resultant cumulative infiltration curve exhibited a convex shape typically associated with soil water repellency [START_REF] Alagna | Alternative analysis of transient infiltration experiment to estimate soil water repellency[END_REF]Di Prima et al., 2019;Lassabatere et al., 2019a). This infiltration impedance caused a delayed response of reflection changes, which were initially detected 3 minutes after infiltration started (Figure 10a). In addition, the antenna was positioned at the northwest edge of the infiltrometer, so detection was possible only after sufficient lateral expansion of the wetted zone. At the end of the infiltration test (i.e., at time = 18 min), a strong increase in signal amplitude appeared at approximate depths of 50 and 70 cm, indicating the wetting zone reached the lower boundary of 70 cm (Figure 10a). After that, drainage and water redistribution rates diminished and the water profile approached hydrostatic conditions, as signaled by constant amplitudes for the final 15 minutes of the measurement. Because the wetted zone did not change in shape and dimension at the end of the infiltration run, we assume that water drainage and redistribution was slow enough to allow proper acquisition of GPR over the complete grid survey. In other words, the last GPR acquisition (last Survey line) described the same soil moisture condition than the first survey line.

Consequently, all survey lines described the same wetted zone and can be concatenated to form a 3D representation of the entire wetted zone.

GPR punctual acquisition and time-lapse GPR (3D block) were crossed-checked. At the end of the measurement, we identified the depth of the stronger signal amplitudes in the final scans of the radargram (Figure 10a). We subsequently extracted horizontal slices from 3D diagram at the corresponding depths (Figure 10b). There was a perfect agreement between stronger signal amplitudes and the occurrence of deep wetter patches positioned vertically below the antenna.

These observations increased our confidence in the 3D interpolation procedure that was previously performed to determine the wetted zone volume.

Discussion

The use of time-lapse GPR surveys in conjunction with staining-dye tracer at the DOUA site allowed us to verify the 3D interpolation procedure. Here the comparison confirmed that the GPR was able to accurately resolve perched water above a shallow restrictive layer along with deeper preferential flow associated with a shrub plant root system (Figure 8). That pattern is in line with the expected water infiltration in layered systems and with the potential role of leaks by the plant root systems in layered soils [START_REF] Gerke | Preferential flow mechanisms identified from staining experiments in forested hillslopes: Preferential Flow Mechanisms Identified from Staining Experiments[END_REF]. However, the wetted zones detected by the GPR in the vicinity of the root system were larger than the zones visibly stained with dye. Further investigations are thus needed to better understand the dimensions of flow pathways detectable with GPR, including whether the technology is able to resolve very small fingers or pores.

For the Django Reinhardt site, we observed a uniform distribution of water with hemisphericalshaped wetted zones in the upper layer where the soil matrix controlled water flow, and vertically elongated and irregularly shaped wetting patterns corresponding to the coarse gravels below.

Therefore, the 3D observation of the wetted zones using GPR allowed the detection of preferential funneled flow due to lithological heterogeneities. This result can be viewed as an experimental confirmation of previous modeling investigations carried out at the same infiltration basin by [START_REF] Coutinho | Vadose Zone Heterogeneity Effect on Unsaturated Water Flow Modeling at Meso-Scale[END_REF] and Ben [START_REF] Ben Slimene | The role of heterogeneous lithology in a glaciofluvial deposit on unsaturated preferential flowa numerical study[END_REF]. The latter study speculated that lithological heterogeneities in the basin (e.g., the subsurface gravel lenses) may induce funneled flow, so the present study provides experimental evidence showing that this process does commonly occur (as we observed PF at both tested locations). Moreover, the 3D diagrams provided additional insight into spatial differences in PF processes. Here, the depth of funneled flow varied between locations in the same site, as did the subsequent shape, volume, and depth of wetting in the finer-textured underlying soil. The two infiltration tests introduced nearly identical volumes of water into the soil, so the GPR scans were crucial to identify actual differences in infiltration and wetting that occurred.

Our findings clearly show that GPR can be used for the detection for preferential flow and to get a picture of water pathways for both cases: lithological heterogeneity and root systems. These data have notable practical applications. For example, infiltration patterns measured at the surface of a stratified medium are often influenced by hydraulic behaviors of any restrictive layers, yet infiltration data alone typically do not contain sufficient information to characterize different layers. This limitation was identified by [START_REF] Lassabatere | Effect of the settlement of sediments on water infiltration in two urban infiltration basins[END_REF] for the case of the sediment settlement in the Django Reinhardt basin, as well as by many authors in other locations (e.g., [START_REF] Alagna | Using Beerkan experiments to estimate hydraulic conductivity of a crusted loamy soil in a Mediterranean vineyard[END_REF][START_REF] Di Prima | Laboratory testing of Beerkan infiltration experiments for assessing the role of soil sealing on water infiltration[END_REF][START_REF] Yilmaz | Influence of Carbonation on the Microstructure and Hydraulic Properties of a Basic Oxygen Furnace Slag[END_REF][START_REF] Yilmaz | Hydrodynamic Characterization of Basic Oxygen Furnace Slag through an Adapted BEST Method[END_REF]. Similarly, when water flows simultaneously through the soil matrix and a macropore network, both flow domains contribute to the hydraulic properties of the soil (Lassabatere et al., 2019b[START_REF] Lassabatere | New Analytical Model for Cumulative Infiltration into Dual-Permeability Soils[END_REF][START_REF] Stewart | A Generalized Analytical Solution for Preferential Infiltration and Wetting[END_REF]. Our analysis showed that non-uniform wetting patterns are detectable using GPR, so this information could be used to improve the parameterization of hydraulic properties for both homogeneous and heterogeneous soils.

One potential limitation of the applied protocol is the uncertainty associated with spurious amplitude differences. This problem depends in part on small differences of the antenna position between successive surveys. Also, the lowering of the reflection position after the increase of the water content may induce a time shift and consequently a mismatch in the vertical direction. In this investigation, we alleviated these problems, at least in part, by smoothing the radargrams both in the horizontal and vertical directions. The smoothing process proved satisfactory for removing these potential artifacts while retaining sufficient resolution to demarcate the extent and degree of wetting during localized infiltration processes. Therefore, we encourage other practitioners to employ similar techniques when using GPR-based measurements to analyze preferential wetting behaviors.

Summary and conclusions

In this study, we combined time-lapse GPR surveys with automated single-ring infiltration experiments. The proposed procedure used non-invasive monitoring of the spatial distribution of infiltrated water to create 3D representations of infiltrated water. A dye-staining experiment conducted at one of the experimental sites provided verification that the applied procedure was able to identify preferential flow caused by shrub roots. In all experimental sites, the 3D diagrams created using the proposed procedure clearly demarcated the wetting zone dimensions and shapes, and facilitated detection of preferential flow pathways and funneled flow in the subsurface. In the future such information may serve to verify hydraulic characterization procedures using water infiltration experiments, and may be particularly useful in highly heterogeneous soils, including sealed soils, multi-porosity/multi-permeability systems, and stratified profiles. Indeed, this study can be viewed as a step towards time-lapse GPR surveys becoming more widely applied in the field as a reliable and non-invasive way for investigating linkages between non-uniform flow and heterogeneous soil structures.

List of Figures

Figure 1 .

 1 Figure 1. Experimental setup at the (a) Django Reinhardt infiltration basin and (b) Parc de La

  Figure 1. Experimental setup at the (a) Django Reinhardt infiltration basin and (b) Parc de La

Figure 2 .

 2 Figure 2. (a) Common-Mid-Point (CMP) determination of the velocity profile over depth (acquired

Figure 3 .

 3 Figure 3. Flowchart illustrating the procedure to obtain the 3D diagram of the wetting zone from

Figure 4 .

 4 Figure 4. Cumulative frequency distributions of the absolute difference between pre-and post-

Figure 5 .

 5 Figure 5. (a), (f): Three-dimensional representation of the wetting zones at Sites A and B. (b), (e),

Figure 6 .

 6 Figure 6. Comparison between (a) the infiltration bulb observed at the end of the experiment

Figure 7 .

 7 Figure 7. Volumes of the wetted zones for Sites A and B (left Y-axes) and associated cumulative

Figure 8 .

 8 Figure 8. (a): Three-dimensional representation of the wetting zones at the DOUA site. (b):

Figure 9 .

 9 Figure 9. Overlapping of selected vertical slices extracted from 3D diagrams with the pre-wetting

Figure 10 .

 10 Figure 10. (a) Overlapping of the cumulative infiltration curve with the radargram measured at the

  Figure 1 Click here to download high resolution image

Table 1 .

 1 Bulk density values (g cm -3 ) measured at the DOUA site.

	644
	Depth (cm)	Sample size	Min	Mean	Max
	0-10	9	0.928	1.074	1.173
	10-20	3	1.046	1.090	1.168
	20-30	6	1.349	1.479	1.648
	30-40	3	1.587	1.658	1.738
	40-50	3	1.677	1.716	1.763

Acknowledgements:

The authors thank Sara Puijalon and the Lyon city Field Observatory for Urban Water Management (OTHU), for technical and scientific support, and Jean-Philippe Bedell for advices on plant root systems.

Funding: This work was supported through i) the INFILTRON Project (ANR-17-CE04-0010, Package for assessing infiltration & filtration functions of urban soils in stormwater management; https://infiltron.org/), ii) the European Regional Development Fund (ERDF) and the Italian Ministry of Education, University and Research (MIUR) through the "Programma Operativo Nazionale (PON) Ricerca e Innovazione 2014-2020 (Linea 1 -Mobilità dei ricercatori, AIM1853149, CUP: J54I18000120001), and iii) the Short Term Mobility (STM) Program 2019 of CNR (Research Programme: "Innovative techniques for hydraulic characterization of soil").

Conflicts of Interest:

The authors declare no conflict of interest.


A superhydrophobic capillary tube immersed in water and brought in contact with the bath surface will be invaded by air, owing to its aerophilicity. We discuss this phenomenon where the ingredients of classical capillary rise are inverted, which leads to noticeable dynamical features. [START_REF] De Laplace | Supplément à la théorie de l'action capillaire[END_REF] The main regime of air invasion is linear in time, due to the viscous resistance of water. [START_REF] Bell | The flow of liquids through capillary spaces[END_REF] Menisci in tubes with millimetre-size radii strongly oscillate before reaching their equilibrium depth, a consequence of inertia. On the whole, capillary descent provides a broad variety of dynamics where capillary effects, viscous friction and liquid inertia all play a role.

Capillary rise occurs when a narrow tube is brought in contact with a wetting liquid, as first reported by Leonardo and quantified by Hauksbee and Jurin at the beginning of the 18 th century. Despite the efforts of Newton, the phenomenon was only understood one century later [START_REF] De Laplace | Supplément à la théorie de l'action capillaire[END_REF]. The meniscus is driven by a negative Laplace pressure ∆P = -2γcosθ/r, denoting γ as the surface tension of the liquid, ߠ < 90° as the contact angle and r as the tube radius (assumed to be smaller than the millimeter-size capillary length κ -1 ). The rise stops at the height H where ∆P balances the hydrostatic pressure loss -ρgH (with ρ the liquid density, and g the acceleration of gravity), which classically yields:

ܪ = 2ߛ cos ߠ ݎ݃ߩ (1)
Considering water and assuming ߠ = 0° (as observed on clean glass), H is expected (and found) to be typically 10 cm for r = 100 µm -a spectacular distance responsible for the early interest on this question: confinement amplifies the natural meniscus height κ -1 = (γ/ρg) 1/2 by a factor κ -1 /r. The rise dynamics was described by [START_REF] Bell | The flow of liquids through capillary spaces[END_REF], followed by [START_REF] Lucas | Ueber das Zeitgesetz des kapillaren Aufstiegs von Flüssigkeiten[END_REF] and [START_REF] Washburn | The dynamics of capillary flow[END_REF] [START_REF] Bell | The flow of liquids through capillary spaces[END_REF][START_REF] Lucas | Ueber das Zeitgesetz des kapillaren Aufstiegs von Flüssigkeiten[END_REF][START_REF] Washburn | The dynamics of capillary flow[END_REF]. As long as the meniscus position z is small compared to H, the meniscus speed ݖሶ generally results from a balance between the capillary force ܨ ఊ = 2ߨߛݎcosߠ and the Poiseuille friction ܨ ఎ = ሶݖݖߟߨ8 [START_REF] Guyon | Physical Hydrodynamics[END_REF],

denoting η as the liquid viscosity and defining ݖሶ as the fluid flux per cross sectional area.

This leads to the Washburn law:

ݖ ଶ ሺݐሻ = ݎߛ cos ߠ 2ߟ ݐ (2)
The meniscus slows down as it progresses ݖ( ሶ ~ ,)ݖ/1 a consequence of viscous friction.

Later, as the final height H is approached, gravity also opposes the rise and z exponentially relaxes to H with a characteristic time ηH/ρgr 2 . Interestingly, the typical time ߟH 2 /ߛrcosߠ ~ ηH/ρgr 2 needed to reach the height H following ሺ2ሻ obeys the same scaling as the relaxation, showing that the two phases of the rise (parabolic regime, exponential relaxation) have comparable timescales. The scaling in t 1/2 resists many configurations. Hele-Shaw cells, grooves or rough solids are generally invaded with such dynamics [START_REF] Deng | Characterization of capillary rise dynamics in parallel micro V-grooves[END_REF][START_REF] Shen | A model for capillary rise in nano-channels with inherent surface roughness[END_REF], which is also observed when complex porous media (fabrics, sand, sponges, papers) contact a wetting liquid [START_REF] Perwuelz | Experimental study of capillary flow in yarns[END_REF][START_REF] Lago | Capillary rise in porous media[END_REF][START_REF] Marmur | Kinetics of penetration into uniform porous media: Testing the equivalentcapillary concept[END_REF][START_REF] Nia | Theoretical Analysis of Capillary Rise in Porous Media[END_REF]. This universality was exploited to determine either contact angles or characteristic radii in porous media [START_REF] Marmur | Characterization of porous media by the kinetics of liquid penetration: the vertical capillaries model[END_REF][START_REF] Bachmann | Extended methodology for determining wetting properties of porous media[END_REF].

3

Conversely, a few uncommon situations that escape Washburn law were recently evidenced. (i) Staple (2001), [START_REF] Reyssat | Imbibition in geometries with axial variations[END_REF] and [START_REF] Figliuzzi | Rise in optimized capillary channels[END_REF] analytically and experimentally studied the progression of liquids in capillary tubes with non-uniform cross-sections [START_REF] Staples | Wicking flow in irregular capillaries[END_REF][START_REF] Reyssat | Imbibition in geometries with axial variations[END_REF][START_REF] Figliuzzi | Rise in optimized capillary channels[END_REF]. They reported scaling laws of the type ݖሺݐሻ ~ ݐ ఉ , where β was found to depend on the tube shape and to vary between 0 and ½. (ii) V-shaped channels (for which there is no more characteristic distance such as its radius for a tube) also generate special dynamics [START_REF] Romero | Flow in an open channel capillary[END_REF][START_REF] Weislogel | Compound capillary rise[END_REF][START_REF] Higuera | Capillary rise of a liquid between two vertical plates making a small angle[END_REF], such as z(t) ~ t 1/3 [START_REF] Ponomarenko | A universal law for capillary rise in corners[END_REF]. (iii) As recently shown by [START_REF] Kim | Capillary rise of non-aqueous liquids in cellulose sponges[END_REF], porous media with largely-distributed pores also have unusual dynamics of impregnation (as t 1/4 ), due to the coexistence of characteristic length scales [START_REF] Kim | Capillary rise of non-aqueous liquids in cellulose sponges[END_REF]. (iv) A final case concerns inertial effects at short times, especially for wide tubes: the invading liquid is then resisted by inertia and the balance of the capillary force ܨ ఊ with the momentum variation dሺߨݎ ଶ ሶݖݖߩ ሻ/dݐ leads to a surprisingly simple law, ݖሺݐሻ = ,ݐܿ where c = ሺ2ߛcosߠ/ߩݎሻ ଵ/ଶ is called the inertio-capillary speed [START_REF] Quéré | Inertial capillarity[END_REF][START_REF] Fries | The transition from inertial to viscous flow in capillary rise[END_REF].

For non-wetting liquids (θ > 90°), Equation 1predicts negative values for the equilibrium position of the meniscus, a situation that can be generated by bringing an "aerophilic" tube immersed in water in contact with air. Then we expect the meniscus separating air and water to be pulled down to the depth H -as observed by Jurin for quicksilver in a glass tube [START_REF] Jurin | An account of some new experiments relating to the action of glass tubes upon water and quicksilver[END_REF]. Aerophilic tubes are obtained by coating the inner walls of tube with a hydrophobic texture, which leads to the well-known situation of superhydrophobicity, where solids are wet by air rather than by water [START_REF] De Maleprade | Spreading of bubbles after contacting the Lower side of an aerophilic slide immersed in water[END_REF][START_REF] Wang | Air bubble bursting effect of lotus leaf[END_REF]. In our experiments, we use glass capillary tubes with radius r and total length L respectively ranging between 0.51 ± 0.01 mm and 1.50 ± 0.01 mm, and between 5.5 ± 0.1 cm and 13.5 ± 0.1 cm. The tubes are treated by filling them with isopropanol containing silanized silica nanobeads with diameter 30 nm (Glaco, Soft99). After removing the solvent with compressed air, the nanocoating is consolidated at 180°C for 30 min. We repeat the treatment three times to make the aerophilic layer more robust. The resulting tubes are transparent, and their walls have a typical roughness of 100 nm, much smaller than the tube radius, as shown in the Supplementary Information (figure SI-1). The advancing and receding contact angles of water on a flat surface with the same coating are found to be 165 ± 1° and 160 ± 2°, which provides for air on immersed solids the complementary values θr = 15 ± 1° and θa = 20 ± 2°. The coated tubes are then vertically immersed at room temperature (between 20°C and 22°C) inside a pool of deionised water (surface tension γ = 72.0 ± 0.5 mN/m, density ρ = 998 ± 1 kg/m 3 , viscosity η = 0.98 ± 0.03 mPa.s). By connecting the outer extremity of the capillary tube to a small funnel, we pump the air trapped inside with a syringe, and we visually control that no air remains.

Then tubes are slowly brought up to the bath surface, using a micrometric screw (figure 1a). Contact with air defines the origin of time, from which we follow the position of the air/water interface with a high-speed video camera (Phantom v7). As seen in figures 1b and 1c, the meniscus moves downward and follows two types of dynamics: either the descent is roughly linear in time before relaxing to equilibrium (thin tubes, figure 1b), or the meniscus oscillates before reaching equilibrium (wide tubes, figure 1c). We also observe that the final position is deeper when confinement is stronger. We first discuss the centimetre-size equilibrium depth reached by the meniscus. To that end, we extract from the movies the evolution z(t) of the position of the menisci for four tubes having the same length yet different radii. We plot the collection of graphs in figure 2 where dotted lines indicate the final depth expected from equation (1). The agreement is good, and it can be improved by noticing that the final depth slightly depends on the dynamics. (i) At small r (r = 0.51 mm), equilibrium is reached monotonously, so that we must compare data with the value ܪሺθ ୟ ሻ calculated in equation 1 with the advancing angle θ ୟ (black dotted line). The observed depth is smaller by 4% than expected, which might be due to coating defects. Superhydrophobic treatment at small r is challenging because it must remain of homogeneous quality despite both the confinement and the large distance z. Coating heterogeneities make tubes less aerophilic, and thus limit the depth reached by air. Conversely, our experiment provides a simple test for checking the quality of aerophilic treatments in confined geometries. (ii) At large r, capillary descent is accompanied by oscillations, so that the contact angle switches between its advancing and receding values. Thus, the equilibrium depth is ill-defined between ܪሺθ ୟ ሻ and ܪሺθ ୰ ሻ, whose theoretical values are drawn in black and red, respectively. We observe that oscillating menisci systematically reach a position between ܪሺθ ୟ ሻ and ܪሺθ ୰ ሻ, as expected from this discussion. The uncertainty on H remains modest, owing to the small angle hysteresis ∆ߠ = ߠa -ߠr in our tubes (about 5°): the relative difference between ܪሺθ ୟ ሻ and ܪሺθ ୰ ሻ is found in equation 1

(at small ߠ) to be ∆ߠ.ߠ ̅ (noting ߠ ̅ the mean contact angle), about 3% in our experiments.

We now discuss the dynamics of the descent, and distinguish the two kinds of behaviour (monotonic, oscillatory). As seen in figure 2 (light blue data), the meniscus in a thin tube is first accelerated (on a short time, smaller than 50 ms) until it reaches a constant velocity of the order of 10 cm/s. Later, the meniscus relaxes to its equilibrium position H ≈ -26 mm in a timescale of a fraction of a second, comparable to that of the descent.

A regime of constant speed markedly differs from what is observed for a liquid rising in a tube (equation 2) [START_REF] Bell | The flow of liquids through capillary spaces[END_REF][START_REF] Lucas | Ueber das Zeitgesetz des kapillaren Aufstiegs von Flüssigkeiten[END_REF][START_REF] Washburn | The dynamics of capillary flow[END_REF]. While rising columns of liquid gradually increase in height, capillary descent consists of replacing a column of water by air of negligible viscosity.

Viscous effects should mainly concern the water column in which the time necessary to establish a Poiseuille flow is ߬ ≈ ݎߩ ଶ /8ߟ, that is, ߬ ≈ 30 ms for r = 0.51 mm. After this short inertial phase, but still far from equilibrium (z < H), the force ܨ ఊ = ߛݎߨ2 cos θ driving the meniscus is balanced by the Poiseuille friction in water ܨ ఎ = ሶݖܮߟߨ8 , assuming z << L -a condition fulfilled for our tubes, all longer than H. Hence we deduce that the meniscus progresses at a constant velocity given by:

ܸ = ߛ cos θ 4ߟ ݎ ܮ (3)
The speed of descent is determined on the one hand by the materials properties (surface tension ߛ, angle ߠ, viscosity ߟ), and on the other hand by the tube dimensions (radius ,ݎ length L). For water, r ≈ 1 mm and L ≈ 10 cm, we expect V to be on the order of 10 cm/s, as indeed observed in figures 1 and 2. The characteristic time of descent

߬ ௗ = H/V is ܮߟ8 ݎ݃ߩ ଶ ⁄ .
For the thinnest tube in figure 2, this time is 0.24 s, in excellent agreement with the observations. Since V characterizes the beginning of the descent, it is independent of g, which we checked by inclining the immersed tubes (figure SI-2).

We test this original visco-capillary regime by varying the aspect ratio r/L of the tubes.

To that end, we measure the descent velocity V at the inflection point in the curves z(t), that is, between the acceleration and relaxation phases. We display in figure 3a The relaxation following the linear descent is quite similar to that observed in capillary rise [START_REF] Zhmud | Dynamics of capillary rise[END_REF], and we provide some details about it in the SI (figure SI-3). Its duration is obtained after introducing the hydrostatic thrust ݎߨ- ଶ ݖ݃ߩ in the model, which yields an exponential relaxation of characteristic time ߬ = 8ߟ(L-H)/ߩgr 2 . This time is expected to be ~0.2 s for the thinnest tubes, in agreement with the observations in figure 2. For long tubes (H << L), ߬ is similar to the time of descent H/V = 8ߟL/ߩgr 2 deduced from equations 1 and 3, a common feature with Washburn law, as discussed after equation 2.

For wide tubes (r > 1 mm), inertial effects cannot be neglected, as seen in figure 2 where the meniscus significantly overcomes the equilibrium depth and then oscillates, a behaviour also observed in capillary rise [START_REF] Quéré | Inertial capillarity[END_REF][START_REF] Zhmud | Dynamics of capillary rise[END_REF]. Close to equilibrium (z = H + ߝ, ߝ << H), viscous friction, gravity force and inertia of water respectively scale as ߟ(L-H)ߝሶ , ߩgr 2 ߝ and ߩr 2 (L-H)ߝሷ . Prefactors in these expressions can be easily calculated and the balance of forces writes: ߩgr 2 ߝ + 8ߟ(L-H)ߝሶ + ߩr 2 (L-H)ߝሷ = 0. This second order differential equation results in oscillatory regimes for tube radius r larger than (16ߟ 2 (L-H)/ߩ 2 g) 1/4a quantity of ~1 mm for water and decimeter-long tubes, as observed in figure 2. This condition can also be expressed through a critical tube length:

ܮ = ݃ ቆ ݎߩ ଶ 4ߟ ቇ ଶ + ܪሺݎሻ (5)
For L > ܮ , the descent is expected to be monotonous while it should be oscillatory for L < ܮ . Equation 5 predicts ܮ ~ 10 cm for millimetric tubes, a distance that rapidly decreases with the viscosity. But ܮ is even more sensitive to the tube radius found to be the main parameter controlling the descent dynamics. While the critical length is several meters for r = 1.5 mm, it is only a few centimetres for r = 0.5 mm.

In order to establish a dynamical phase diagram of capillary descent, we perform experiments with a collection of fifteen vertical tubes with different lengths and radii

(5 < L (cm) < 13, 0.5 < r (mm) < 1.5). We report in figure 4 whether the meniscus oscillates (green data) or not (red data), and compare our observations with equation 5.

The agreement is satisfactory, even if the transition to the oscillatory regime is slightly overestimated, which might reflect the existence of other sources of damping such as the singular pressure loss at the tube exit (if the liquid descends) or entrance (if the liquid rises), which we discuss in the SI (figure SI-4).

The regime of monotonous descent is only observed for long and thin tubes, but the latter condition is not sufficient. As seen in figure 4, oscillations may also occur in thin tubes, provided that L is smaller than typically 5 cm, corresponding to the minimum of the function Lc(r) in equation 5. Hence the minimal length required to observe monotonous regimes is Lmin ≈ ߢ -4/5 D 1/5 where D = ߛ 2 /gߟ 2 is a kilometric distance for water, making Lmin a few centimetres, as reported in figure 4.

In summary, capillary descent provides a broad variety of dynamics, where the different liquid parameters (surface tension, viscosity, density) all play a role, owing to the combination of capillarity, viscous effects and inertia at stake in these experiments. It would be interesting to extend our results to the case of immersed hydrophobic porous media to see how our findings resist the reduction of scale and the complexity encountered in such materials. More generally, it is worth realizing that the inversion of the two fluids present in capillary rise generates similar equilibrium conditions, but quite different dynamics, due to the contrast in viscosity between water and air.

  de Maleprade, Christophe Clanet & David Quéré Physique & Mécanique des Milieux Hétérogènes, UMR 7636 du CNRS, ESPCI, 75005 Paris, France. and LadHyX, UMR 7646 du CNRS, École polytechnique, 91128 Palaiseau, France.

Figure 1 .

 1 Figure 1. (a)Sketch of an experiment of capillary descent in an aerophilic tube. Blue and white represent water and air, respectively. The origin of time is taken when the tube contacts air. (b) Meniscus position at t = 50 ms, 90 ms and 250 ms for a tube with r = 0.51 mm and L = 7.5 cm. Water appears in black, and air in white. The meniscus dynamics can be reconstructed from such images, and its depth z is linear in time (after a short regime of acceleration) before relaxing to equilibrium. (c) Same experiment with r = 1.50 mm and L = 7.4 cm. The three pictures are taken at t = 40 ms, 280 ms and 1600 ms, and the dynamics is now found to be non monotonous, with strong oscillations before equilibrium is reached.

Figure 2 .

 2 Figure 2.Meniscus dynamics in a vertical superhydrophobic tube immersed in water, after its contact to the surface of the bath. Data correspond to four tubes with the same length (L = 7.5 ± 0.1 cm) but different radii between r = 0.51 mm (light blue) and r = 1.50 mm (black). The dotted lines show the equilibrium depths given by equation 1, taking either the advancing angle of air in water θa = 20° (black lines) or the receding angle θr = 15° (red lines).

  these velocities as a function of the tube length L, at fixed r. As predicted by equation 3, V decreases in a hyperbolic way (thin lines) with L. Data with various r and L are collected in figure 3b, where we plot the descent velocity as a function of the aspect ratio .ܮ/ݎ All data collapse on a common straight line of slope 16.0 ± 0.5 m/s, in good agreement with the value expected from equation 3, γcosߠ/4η ≈ 17.3 m/s.

Figure 3 .

 3 Figure 3.Visco-capillary regime of constant descent velocity V. In (a), V is plotted as a function of the total length L of the immersed tube, for two radii (r = 0.51 mm, red data; r = 0.70 mm, blue data). Lines show the hyperbolic behaviour predicted by equation 3. In (b), V is plotted as a function of the tube aspect ratio r/L. The linear fit has a slope of 16 m/s, close to γcosθ/4η ≈ 17 m/s, the value expected from equation 3.

Figure 4 .

 4 Figure 4. Dynamical phase diagram of capillary descent in the (L, r) space. The green and red regions respectively indicate whether oscillations are expected or not from equation 5. The hatched region stands for tubes smaller than the equilibrium depth, so that the meniscus pins at the tube end. Dots are the data (green dots: oscillating behaviour, red dots: no oscillation).

Published on 21 May 2018. Downloaded by Kaohsiung Medical University on 23/05/2018 13:53:23. View Article Online DOI: 10.1039/C8SM00453F