Xianjue Chen
email: xianjue.chen@unsw.edu.au

Marc Dubois
email: marc.dubois@uca.fr

Silvana Radescu

Aditya Rawal

Chuan Zhao

Liquid-phase exfoliation of F-diamane-like nanosheets

Keywords: diamane, graphite fluoride, fluorination, exfoliation, 2D materials

Introduction

Atomically thin diamond ("diamane") has recently emerged as a new two-dimensional carbon allotrope. Theoretical studies first predicted the conversion of AB-stacked few-layer graphene into diamane films through sufficient formation of carbon-fluorine (C-F) or carbon-hydrogen (C-H) bonds on the two free surfaces [START_REF] Chernozatonskii | Diamond-Like C2H Nanolayer, Diamane: Simulation of the Structure and Properties[END_REF][START_REF] Ribas | Patterning Nanoroads and Quantum Dots on Fluorinated Graphene[END_REF][START_REF] Kvashnin | Phase Diagram of Quasi-Two-Dimensional Carbon, From Graphene to Diamond[END_REF][START_REF] Antipina | Converting Chemically Functionalized Few-Layer Graphene to Diamond Films: A Computational Study[END_REF]. Odkhuu et al. showed that when a one-side surface of a bilayer graphene grown on a transition metal surface was fluorinated, the energetics for the conversion to sp 3 -bonded films were significantly lower than that for the free-standing bilayer [START_REF] Odkhuu | Conversion of multilayer graphene into continuous ultrathin sp 3 -bonded carbon films on metal surfaces[END_REF]. Fluorination and hydrogenation of graphene have been experimentally investigated [START_REF] Elias | Control of Graphene's Properties by Reversible Hydrogenation: Evidence for Graphane[END_REF][START_REF] Nair | Fluorographene: A Two-Dimensional Counterpart of Teflon[END_REF].

Recent advances in chemical vapor deposition (CVD) growth of graphene towards precise quality control (number of layers, crystallinity, coverage and uniformity) have indicated pathways for converting few layers of graphene into "diamane" films [START_REF] Huang | Large-Area Single-Crystal AB-Bilayer and ABA-Trilayer Graphene Grown on a Cu/Ni(111) Foil[END_REF]. Bakharev et al. [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF] first demonstrated that fluorine chemisorption on CVD-grown, AB-stacked, bilayer graphene led to a fluorinated diamond monolayer ("F-diamane"). Rajasekaran et al. [START_REF] Rajasekaran | Interlayer Carbon Bond Formation Induced by Hydrogen Adsorption in Few-Layer Supported Graphene[END_REF] reported experimental evidence that hydrogen adsorption could induce partial phase transition of few-layer graphene to an sp 3 -bonded structure on a Pt(111) surface. Diamane-like structures have also been achieved by compressing trilayer graphene to ultrahigh-pressure (>20 GPa) at room temperature [START_REF] Ke | Synthesis of Atomically Thin Hexagonal Diamond with Compression[END_REF], or exposing bilayer graphene to hydrogen radicals generated by a hot filament process [START_REF] Piazza | Raman Evidence for the Successful Synthesis of Diamane[END_REF]. The stability, mechanical, electronic, optical and thermal properties of diamane structures stabilized with H, F or Cl atoms have been predicted using DFT simulations [START_REF] Mortazavi | First-Principles Investigation of Mechanical, Electronic and Optical Properties of H-, F-and Cl-Diamane[END_REF][START_REF] Raeisi | High Thermal Conductivity in Semiconducting Janus and Non-Janus Diamanes[END_REF].

Fluorination of graphite was first reported in 1934, where graphite was exposed to fluorine gas at moderate temperatures (400-600 °C) [START_REF] Ruff | Die Reaktionsprodukte der verschiedenen Kohlenstoffformen mit Fluor II (Kohlenstoff-monofluorid)[END_REF]. The product is a solid fluorocarbon compound that has been loosely classified as a graphite intercalation compound due to its layered structures derived from graphite [START_REF] Nakajima | Fluorine-Carbon and Fluorine-Carbon Materials[END_REF]. By controlling the fluorination process, two forms of graphite fluoride can be obtained: "Stage-1" (CF)n, fluorine atoms are bonded above and under graphene planes, where every carbon atom is covalently bonded with a fluorine atom, producing a chair-type structure; "Stage-2" (C2F)n, fluorine atoms are inserted into every second layer of graphite where the carbon atoms establish an sp 3 orbital with a double-decked structure. The structural model proposed in previous studies [START_REF] Kita | Chemical Composition and Crystal Structure of Graphite Fluoride[END_REF][START_REF] Watanabe | Two Types of Graphite Fluorides, (CF)n and (C2F)n, and Discharge Characteristics and Mechanisms of Electrodes of (CF)n and (C2F)n in Lithium Batteries[END_REF][START_REF] Touhara | On the Structure of Graphite Fluoride[END_REF] is illustrated in Fig. 1a,b. Note that the (C2F)n is essentially composed of stacked layers of "F-diamane". The low surface energy caused by the short, strong C-F bonds results in a low friction coefficient between the layers [START_REF] Tuteja | Designing Superoleophobic Surfaces[END_REF], and thus the "F-diamane" layers could be exfoliated under shear.

Liquid-phase exfoliation has been widely used in the preparation of single-layer and fewlayer platelets from bulk flakes of layered materials including graphite, hexagonal boron nitride, black phosphorus, transition metal dichalcogenides, and more [START_REF] Hernandez | High-Yield Production of Graphene by Liquid-Phase Exfoliation of Graphite[END_REF][START_REF] Chen | Vortex Fluidic Exfoliation of Graphite and Boron Nitride[END_REF][START_REF] Chen | p-Phosphonic Acid Calix[8]arene Assisted Exfoliation and Stabilization of 2D Materials in Water[END_REF][START_REF] Brent | Production of Few-Layer Phosphorene by Liquid Exfoliation of Black Phosphorus[END_REF]. In the case of (CF)n, many different "top-down" exfoliation methods have been developed [START_REF] Zbořil | Graphene Fluoride: A Stable Stoichiometric Graphene Derivative and its Chemical Conversion to Graphene[END_REF][START_REF] Yang | One-Step Preparation of Fluorographene: A Highly Efficient, Low-Cost, and Large-Scale Approach of Exfoliating Fluorographite[END_REF][START_REF] Zhang | Melamine Assisted Solid Exfoliation Approach for the Synthesis of Few-Layered Fluorinated Graphene Nanosheets[END_REF][START_REF] Liu | Ammonia Borane Assisted Solid Exfoliation of Graphite Fluoride for Facile Preparation of Fluorinated Graphene Nanosheets[END_REF][START_REF] Dubois | Thermal Exfoliation of Fluorinated Graphite[END_REF][START_REF] Herraiz | Large-Scale Synthesis of Fluorinated Graphene by Rapid Thermal Exfoliation of Highly Fluorinated Graphite[END_REF][START_REF] Wu | Shock Exfoliation of Graphene Fluoride in Microwave[END_REF]. Unlike the CVD technique that has a tight process control and a relatively high cost, liquid-phase exfoliation is a simple, versatile means to produce ultrathin layers in large quantities [START_REF] Feng | Two-Dimensional Fluorinated Graphene: Synthesis, Structures, Properties and Applications[END_REF]. Herein, we show sonication-assisted exfoliation of (C2F)n in solvents yields few-layer platelets with an "Fdiamane"-like structure.

Material and Methods

Preparation of materials

Poly(dicarbon monofluoride) (C2F)n was synthesized using the conventional method. A monel boat containing Madagascar natural graphite was placed in a nickel reactor and heated to 350 °C under pure F2 gas flow for 8 h. (CF)n (GT1FF012) were purchased from ACS Material and used without further modification. Suspensions (10 mL) of (C2F)n (0.5 mg mL -1 ) in water or organic solvents were subject to sonication using a probe-sonicator (Branson Digital Sonifier 450; 20% amplitude, 1h; under ice bath and with a pulse effect of 10s on and 10s off to avoid overheating).

The suspensions in organic solvents were centrifuged (Centrifuge 5810, Eppendorf, 117 × g, 10 min) to remove non-dispersible species. The supernatants were separated affording stable suspensions.

Material characterization

XRD patterns were obtained using a PANalytical Xpert Multipurpose X-ray Diffraction System and analyzed using X'Pert HighScore Plus software. SEM images were recorded on a FEI Nova NanoSEM 450 FE-SEM. XPS spectra were acquired using an ESCALAB 250 Xi, Thermo Scientific. The 19 F solid-state NMR spectra were acquired on a Bruker AVANCE III 700 MHz spectrometer with a 16.4 tesla superconducting magnet, operating at a frequency of 658.88 MHz for the 19 F nucleus. Approximately 30 mg of sample was packed in a in 3.2 mm zirconia rotors fitted with Vespel ® caps and spun to 24 kHz at the magic angle in a double resonance 3.2 mm HX Low-gamma MAS probehead. The 19 F 90° pulse length of 2.5 μs was used with sufficient recycle delay to ensure full signal relaxation. The 19 F chemical shifts were referenced to CFCl3 at reference. The DMFIT software was used for deconvolution of the NMR spectra [START_REF] Massiot | Modelling One-and Two-Dimensional Solid-State NMR Spectra[END_REF]. AFM images were acquired on a Bruker Dimension ICON SPM operating in the tapping mode. A dispersion of (C2F)n in chloroform was drop-cast onto a freshly cleaved mica substrate and dried in air. AFM images were analyzed using the Gwyddion SPM software. TEM specimens were prepared by drop casting a (C2F)n dispersion onto a holey carbon-coated Cu grid (SPI Supplies, #2450-AB) and dried in air. TEM/STEM analysis was performed using a JEOL JEM-F200 Multi-Purpose FEG-S/TEM operating at an accelerating voltage of 80 kV. Image J and Pathfinder were used for processing the TEM and elemental mapping images. Raman spectra were recorded on a Renishaw InVia 2 Raman spectrometer at room temperature with a 532 nm laser excitation. For FTIR, pellets were prepared by compressing the well-mixed KBr and sample powders. FTIR was recorded on a Spectrum 100/Spotlight 400 spectrometer in the transmission mode.

Computational details

The DFT calculations for the (C2F) monolayer were performed using the Vienna Ab Initio Simulation Package (VASP) computational code, which implements the projector-augmented wave method (PAW), and the Perdew-Burke-Ernzerhof (PBE) form of the generalized gradient approximation that was taken for the exchange and correlation (XC) energy [START_REF] Kresse | Efficiency of ab-initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set[END_REF][START_REF] Kresse | Efficient Iterative Schemes for ab initio Total-Energy Calculations Using a Plane-Wave Basis Set[END_REF][START_REF] Perdew | Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF]. The PBEsol approach was also used, and the structural results are similar when compared to the PBE method (Table S1 in the Supporting Information). A large vacuum space between monolayers of 33.5 A ̊ is set to avoid interaction between the periodic images of the layers. In the ab initio calculations, only the two 2s and two 2p outermost valence electrons of C, and the two 2s and five 2p outermost valence electrons of F were explicitly taken into account, the rest of the electrons are considered frozen at the core. To ensure a convergence in the energy differences between phases of around 1 meV we used dense uniform 15x15x1 k-points Monkhorst-Pack grids to sample the Brillouin zone, and a cutoff in the kinetic energy of the plane-waves basis set of 600 eV. To obtain the unstrained configuration, the atomic positions and lattice vectors were fully relaxed (with residual forces less than 0.005 eV A ̊-1 ) and the optimized (relaxed) coordinates were then used for the electronic, elastic, vibrational properties, and core level binding energies.

The mechanical stability of the unstrained monolayer was checked in terms of elastic constants (Born stability criteria). In VASP calculations, the elastic constants for 2D materials can be obtained by normalizing the calculated constant with the lattice parameter c. By applying a set of suitable deformations to the monolayer and taking into account the reduced dimensionality, we obtained the in-plane elastic constants, where the Young's modulus (E) and Poisson's ratio (ν) can be obtained from the following relationships [START_REF] Vannucci | Anisotropic Properties of Monolayer 2D Materials: An Overview from the C2DB Database[END_REF]:

= ( - )/ =
To check the stability of 2D-(C2F) we have calculated the phonon dispersion curves in the whole BZ zone (a structure is locally stable only if the calculated frequencies of all its phonon modes are positive; otherwise, the presence of negative (imaginary) frequencies indicates an instability). To calculate the phonons we used the density functional perturbation theory (DFPT)

as implemented in the VASP code and a post processing Phonopy package to analyze the vibrational modes [START_REF] Togo | First Principles Phonon Calculations in Materials Science[END_REF].

The electron localization function (ELF) was calculated from the valence electronic density obtained within the PAW scheme implemented in VASP (where the core electrons are considered frozen and not explicitly considered). The visualization of the ELF isosurface was done with VESTA 3.5.5 [START_REF] Momma | VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data[END_REF]. We further performed a topological analysis of the all-electron density (the sum of the core and valence densities) calculated with VASP by means of the quantum theory of atoms in molecules (QTAIM), using the CRITIC2 code [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF][START_REF] Otero-De-La-Roza | Critic2: A Program for Real-Space Analysis of Quantum Chemical Interactions in Solids[END_REF]. For the calculation of the core level binding energies of the C1s and F1s core states with a half core-hole effects, we have applied the so-called final state method based in DFT [START_REF] Köhler | Density Functional Study of CO on Rh(111)[END_REF].

Results and Discussion

The synthesis of poly(dicarbon monofluoride) (C2F)n was first reported in 1979 [START_REF] Kita | Chemical Composition and Crystal Structure of Graphite Fluoride[END_REF] The surface composition and C-F bonding of (C2F)n was investigated using X-ray photoelectron spectroscopy (XPS). The survey spectrum (Fig. 2a) collected from (C2F)n shows that the composition is mainly C (60.2 at.%) and F (36.7 at.%). The difference in the fluorine-tocarbon (F/C) ratio (CF0.61 vs. CF0.5) could be owing to the presence of structural defects that allow carbon to be terminated with fluorine by forming CF2 and CF3 groups on the edges of the layers, and/or the presence of a small amount of (CF)n [START_REF] Kita | Chemical Composition and Crystal Structure of Graphite Fluoride[END_REF]. The XPS result for (CF)n (Fig. S1c)

shows a significantly higher F/C ratio of 1.0. Exfoliation of (C2F)n within various solvents (water, ethanol, toluene, chloroform, and DMF) was carried out using probe sonication, Fig. S2. (C2F)n can hardly be dispersed in water due to its hydrophobic nature, and most of the exfoliated material was floating on the water surface or sank to the bottom. Exfoliation within ethanol, chloroform, or DMF resulted in homogenous suspensions after removing any non-dispersible species using centrifugation. The Tyndall effect was observed for all the solvents (including water), suggesting the formation of colloidal dispersions of exfoliated (C2F)n. XPS spectra (Fig. S3) of the exfoliated (C2F)n show increased F/C ratios, 0.90, 0.85, 0.80, 0.82 and 0.90 for water, ethanol, toluene, chloroform and DMF, respectively. This could be owing to the exfoliation-induced fragmentation of (C2F)n that exposes more edges where CF2 and CF3 groups are present, and/or conversion of (C2F)n surface into (CF)n due to partial defluorination.

Raman spectroscopy and Fourier-transform infrared spectroscopy (FTIR) analyses were carried out for (CF)n, (C2F)n, and exfoliated (C2F)n, as shown in Fig. S1 andS4. For (CF)n, the Raman activity was quenched due to its wide band gap [START_REF] Wu | Shock Exfoliation of Graphene Fluoride in Microwave[END_REF], whereas for (C2F)n the Raman bands associated with carbon D (1324 cm -1 ) and G (1586 cm -1 ) bands were observed, which are close to the values measured in previous studies [START_REF] Gupta | Raman Scattering Study of Highly Fluorinated Graphite[END_REF]. In FTIR spectra, the strongest bands at 1192 cm -1 for (CF)n and 1222 cm -1 for (C2F)n can be assigned to the C-F stretching vibration of tertiary carbon atoms [START_REF] Rüdorff | IR-Spektrum des Kohlenstoffmonofluorids[END_REF]. A medium band at 1348 cm -1 was observed in the spectra, possibly corresponding to asymmetric stretching vibrations of peripheral CF2 groups [START_REF] Lagow | Some New Synthetic Approaches to Graphite-Fluorine Chemistry[END_REF]. The exfoliated (C2F)n show similar Raman and FTIR features to (C2F)n, which suggests that the "F-diamane" structure was retained.

Fig. 2c shows the 19 F solid state NMR spectra of (CF)n, (C2F)n, and exfoliated (C2F)n.

The presence of covalent C-F bonding in (C2F)n and (CF)n is confirmed by the position of the isotropic 19 F NMR peaks. The primary signal of the 19 F species was observed at the chemical shifts of -186 ppm and -187 ppm (vs CFCl3), for (CF)n and (C2F)n respectively. The (C2F)n and the exfoliated (C2F)n showed nearly identical 19 F spectra, indicating that the exfoliation process did not induce major chemical changes in the material. A relatively weak signal at -84 ppm was observed for (C2F)n, corresponding to terminal CF3 groups. An associated isotropic signal relating to CF2 groups at -120 ppm is also expected [START_REF] Zhang | Solid-State NMR and EPR Study of Fluorinated Carbon Nanofibers[END_REF][START_REF] Dubois | EPR and Solid-State NMR Studies of Poly(dicarbon monofluoride) (C2F)n[END_REF], but in the present case is overlapped by spinning sidebands. A resolved shoulder at -176 ppm was observed in the (C2F)n, while a similar signal albeit with a lower resolution is also evident in (CF)n. This signal is typical of (C2F)n layers [START_REF] Zhang | Solid-State NMR and EPR Study of Fluorinated Carbon Nanofibers[END_REF][START_REF] Dubois | EPR and Solid-State NMR Studies of Poly(dicarbon monofluoride) (C2F)n[END_REF], corresponding to the stacking of F-C(sp 3 )-C(sp 3 )-F. Note that in the extreme case, such a chemical shift is similar to fluorinated nanodiamonds, i.e. C(sp 3 )-C(sp 3 )-F, with F being located only on the surface with a low content (-164 ppm) [START_REF] Herraiz | A Multitechnique Study of Fluorinated Nanodiamonds for Low-Energy Neutron Physics Applications[END_REF]. Interestingly, from the deconvolution of the 19 Atomic force microscope (AFM) was used to examine the size and thickness of the exfoliated (C2F)n. Fig. 3a shows a typical image of the platelets from the chloroform dispersion.

The particle size in lateral dimension shows a roughly 100-fold decrease to 50-200 nm, which could result from exfoliation and/or centrifugation that removed larger flakes. AFM height analysis of the exfoliated (C2F)n (Fig. 3b) suggests a layered structure, with a single layer roughly 0.9 nm in thickness.

Fig. 4a shows a transmission electron microscopy (TEM) image of a typical exfoliated (C2F)n nanosheet. The corresponding selected area electron diffraction (SAED) pattern (Fig. 4b) has a six-fold symmetry, indicating a hexagonal symmetry for the exfoliated (C2F)n. Note that the intensities of the first-order diffraction peaks are higher than those of the second-order peaks.

This feature is consistent with the reported result for a suspended fluorinated bilayer graphene ("F-diamane"). For bilayer or few-layer graphene, the first-order diffraction peaks have lower intensities than the second-order peaks [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF]. A magnified TEM image (Fig. 4c), from the indicated area in Fig. 4a, shows the edge of the platelet, with a further zoomed-in image in the inset demonstrating the presence of 2-3 layers. Fig. 4d shows a plot of profile along the line in Fig. 4c inset, showing an interlayer distance of ~0.9 nm. The interlayer distance ranged from 0.8-0.9 nm along the lattice on the platelet edge. This value is consistent with the XRD and AFM results, and is in good agreement with the proposed model. The high-resolution TEM (HRTEM) image of an exfoliated (C2F)n and its corresponding fast Fourier transform (FFT) pattern (Fig. S6)

suggest a hexagonal crystalline order, which is consistent with the reported result [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF]. Elemental Fig. S7 shows the core-loss electron energy loss spectroscopy (EELS) spectrum of the exfoliated (C2F)n. The 1s-π* and 1s-σ* transition peaks can be seen in the C-K edge. Both peaks were observed even for highly fluorinated graphite samples [START_REF] Cavallari | Structural and Electronic Changes in Graphite Fluorides as A Function of Fluorination Rate: An XRS, PDF and DFT Study[END_REF][START_REF] Yoshida | Microscopic Characterization of the C-F Bonds in Fluorine-Graphite Intercalation Compounds[END_REF]. The presence of π* peak could arise from non-fluorinated graphitic species due to the less uniform fluorination of graphite flakes and/or partial defluorination during exfoliation. This is also evident in the XRD and XPS results discussed above. The relatively low intensity of π* peak compared to that of σ* peak indicates the disruption of the conjugated π-orbital network owing to fluorination.

The structure of (C2F)n appeared to be sensitive to electron beam irradiation, which is evident when comparing the TEM images (Fig. S8) taken from the same area before and after short-term exposure to the electron beam under TEM. Such sensitivity was also observed in a previous study where electron beam induced defluorination to bilayer graphene film [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF].

Nevertheless, we found the structure of the exfoliated (C2F)n stable after heating up to 300°C in air, as evident by XRD and 19 F NMR (Fig. S9).

DFT calculations were used to investigate the stability and C-F bonding of an "Fdiamane"-like structure, as well as its mechanical and electronic properties. The structure of an unstrained monolayer is shown in Fig. 5. The calculated bond lengths and cell parameters for the unstrained (C2F) monolayer are listed in Table 1.

Fig. 5 Crystal structure of a "F-diamane"-like (C2F) monolayer.

Table 1. Crystallographic description of the equilibrium (C2F) monolayer.

2D-(C2F) Space Group Atomic Position

Vo = 219. The local analysis of the electron density and its Laplacian at bond critical points (BCP) within the framework of the QTAIM method were used to investigate the covalency character of the chemical bonds in this structure. The C-C bonds are characterized by a higher electron density and larger negative Laplacian than the C-F bonds (see Table 2) indicating the typical covalent-shared interaction between the two C atoms, as seen from the ELF plots (see Fig. S11a).

As shown in Fig. S11b, the higher electronegativity of fluorine shifts the BCP (along the C(1)-F( 1)) path towards the C atom. A BCP was also found along the F(1)-F(1) path, which is characterized by a very small value for both the electron density and the (positive) Laplacian.

The C-F bonds with a smaller positive value of the Laplacian are understood to have partial ionic/covalent character.

A phonon dispersion curve of the 2D-(C2F) structure has been calculated using the DFTP method. As shown in Fig. 6, no imaginary frequency appears indicating that the monolayer is stable at ambient conditions. In order to interpret the binding energies of the core levels measured by XPS, we performed DFT calculations for the C1s and F1s core states in the final state approximation. The calculated values of the binding energies are given in Table 4. From these results we can see that the C atoms are bonded only to other C atoms (C1s-C) as the C-1s values are smaller than the case where they are attached to fluorine (C1s-F), a similar situation previously observed in a study on a C2F-sp 3 slab supercell [START_REF] Cavallari | Tuning C-F Bonding of Graphite Fluoride by Applying High Pressure: Experimental and Theoretical Study[END_REF]. The higher energy values of the core level obtained for the "F-diamane" monolayer can be related to the more covalent character of the C-F bonds than in the slab model (which is also supported by the calculated values of the Laplacian at the BCP).

Our DFT results are in fair agreement with the experimental values of the core level.

Conclusions

We have demonstrated a strategy for preparing ultrathin "F-diamane"-like structures by fluorinating graphite to poly(dicarbon monofluoride) (C2F)n, followed by "top-down" exfoliation of (C2F)n in a range of solvents. Different microscopic and spectroscopic analyses have been used to reveal the structure of the exfoliated (C2F)n nanosheets. Here we summarize the key experimental evidence for the presence of "F-diamane" structure: (i) The 19 F NMR result shows the two types of C-F bonds related to the stacking of F-C(sp 3 )-C(sp 3 )-F (-176 ppm) and F-C(sp 3 )-F (-187 ppm). Both C-F bonds are present regardless the number of layers, and their chemical shifts don't change with this number. This can be used as a direct evidence for the presence of interlayer C-C bonds in the exfoliated (C2F)n. (ii) The TEM images show a layered structure with an interlayer distance of 0.8-0.9 nm for the exfoliated (C2F)n, which is consistent with the XRD and AFM results, and agrees with the proposed (C2F)n structural model that is made of stacked layers of "F-diamane". (iii) The SAED and FFT patterns indicate a hexagonal symmetry for the few-layer (C2F)n, and the intensities of the first-order SAED peaks are higher than those of the second-order peaks. This is consistent with the reported TEM result for "F-diamane" [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF]. Here we also suggest that TEM must be used cautiously with specific settings such as low accelerating voltage and short exposure time, as the structure of (C2F)n was found to be sensitive to electron beam irradiation, and electron beam could induce defluorination [START_REF] Bakharev | Chemically Induced Transformation of Chemical Vapour Deposition Grown Bilayer Graphene into Fluorinated Single-Layer Diamond[END_REF]. The calculated topological parameters at the BCP suggested an increasing covalent character of the C-F bonds. This is related to higher energy values obtained for the core level, which is in good agreement with the experimental data. This method is suitable for preparing "F-diamane"-like materials in large quantities. Additionally, larger flakes of (C2F)n could be made by fluorination of HOPG or graphite with a larger flake size, which might allow a direct mechanical delamination of "Fdiamane" single layers, e.g. using a sticky tape.

  . The yield of (C2F)n depends on a number of factors, including the crystallinity of raw carbon, the fluorinating reagent, and the reaction temperature. Herein, (C2F)n was synthesized by heating graphite flakes to 350 °C under pure F2 gas flow (dynamic mode). SEM image shows the (C2F)n flakes range from 5-20 μm in size, Fig.1c. The X-ray diffraction (XRD) pattern, Fig.1d, suggests an interlayer distance of ~0.81 nm corresponding to the (001) diffraction line of (C2F)n.The peak at 42.1° is assigned to the (100) reflection (d = 0.21 nm). These values are in good agreement with the proposed structural model. The key distinguishing feature for (C2F)n compared to (CF)n is that the "stage-2" intercalated structure has a larger interlayer distance. The XRD pattern (Fig.S1b) for (CF)n suggests an interlayer distance of 0.64 nm. The relatively broad peaks for the (001) and (100) reflections suggests a high degree of structural disorder in (C2F)n and (CF)n. A peak at 26.6° observed for both (C2F)n and (CF)n corresponds to the (002) diffraction from residual non-fluorinated graphitic structure.

Fig. 1

 1 Fig. 1 The proposed (a) top-view and (b) side-view crystal structure of (C2F)n. (c) SEM image and (d) XRD pattern of (C2F)n flakes. The "*" symbols refer to the (002) diffraction from residual graphitic structure.

Fig. 2

 2 Fig. 2 XPS (a) survey and (b) C1s and F1s spectra for (C2F)n. (c) 19 F Solid state NMR of (CF)n, (C2F)n, and exfoliated (C2F)n.

  F spectra (Fig S5), the (C2F)n shows a relatively broad signal near -168 ppm suggesting the formation of such 19 F species.

Fig. 3

 3 Fig. 3 (a) AFM image of exfoliated (C2F)n platelets, with (b) the corresponding height analyses along the lines indicated in (a). A colloidal suspension containing exfoliated (C2F)n platelets is shown in the inset.

  maps of carbon and fluorine acquired in scanning transmission electron microscopy (STEM) mode indicate a uniform presence of carbon and fluorine elements in the exfoliated (C2F)n.

Fig. 4

 4 Fig. 4 (a) TEM image of an exfoliated (C2F)n platelet, with (b) the corresponding SAED pattern. (c) Magnified TEM image from the indicated area in (a), showing the edge of the platelet. The inset in (c) shows a further magnified image of the edge. (d) A plot of the profile along the line in (c) inset, indicating a ~0.9 nm interlayer distance. (e, f) Elemental maps of C and F from the exfoliated (C2F)n platelet.

5 Å 3 C2-C2 = 1 .C1

 531 554 A ̊ C2-C1 = 1.561 A ̊ C1-F1 = 1.378 A ̊ N°164 P-3m1 a = b = 2.559 A ̊; c = 38.711 A ̊ α = β = 90.00 γ = 120.00 1/3, 2/3, 0.431) The results for the elastic stiffness verify the Born's mechanical stability criteria. The calculated values of the in-plane Young's modulus and Poisson's ratio are 478.98 N m -1 and 0.10, respectively. Based on the electronic band structure calculation, we found a direct band gap at the Gamma point (4.02 eV, see Fig. S10) for the 2D-(C2F) structure.

Fig. 6

 6 Fig. 6 Phonon dispersion curves of 2D-(C2F) (a) and vibrational DOS (b).

Table 2 .

 2 Charge density ρ(r) and its Laplacian ∇ 2 ρ(r) calculated at the BCP, in atomic units.

	Chemical meaning	ρ (r)	∇ 2 ρ (r)
	C(2) -C(2)	0.2329	-0.5271
	C(2) -C(1)	0.2377	-0.5401
	C(1) -F(1)	0.2575	0.1259
	F (1) -F(1)	0.0148	0.0856

Table 3 .

 3 Calculated Raman-and Infrared-active phonon modes.

	Active modes	Raman (cm -1 )	Infrared (cm -1 )
	Eg	216.92, 421.63, 1247.23	-
	A1g	440.34, 1127.35, 1292.19	-
	Eu	-	329.66, 1176.01
	A2u	-	915.07, 1307.17

Table 4 .

 4 Core level energies obtained in the DFT calculations. The C1s and F1s energy values are referred to the electrostatic potential at vacuum.

	Core hole level site	Core level (eV)
	C1s (C)	289.1
	C1s (F)	290.5
	F1s (C)	692.7

Acknowledgment X.C. acknowledges the support from the Australian Research Council (DE180100294). S.R. acknowledges the support from Spanish MICINN Project No. PID2019-106383GB-C43. C.Z. is grateful for the award of Future Fellowship from the Australian Research Council (FT170100224). The authors thank Xiao Wang for drawing the crystal structure model in Fig. 1. The authors thank the Mark Wainwright Analytical Centre at the University of New South Wales.

Author Contribution

Declaration of competing interest

The authors declare that they have no known competing financial interests, personal relationships or organizations that could inappropriately influence the work reported in this paper.


Omar Gassara

Audrey Estublier
email: audrey.estublier@ifpen.fr

Bruno Garcia

Sonia Noirez

Adrian Cerepi

Corinne Loisy

Olivier Le Roux

Anelia Petit

Lena Rossi

Sean Kennedy

T Brichart

P Chiquet

L Luu

Van Lang

F André Duboin

J Gance

B Texier

B Lavielle

B Thomas

The Aquifer-CO2Leak project: Numerical modeling for the design of a CO2 injection experiment in the saturated zone of the Saint-Emilion (France) site

Keywords: CO 2, carbon capture and storage (CCS), saturated zone, porous media, experimental site, near-surface hydro-geological and hydro-geochemical

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The average global temperature on Earth has increased by 0.9 • Celsius since the beginning of the industrial era [START_REF] Thomas | Giss Surface Temperature Analysis (GISTEMP), NASA Goddard Institute for Space Studies[END_REF][START_REF] Hansen | Global surface temperature change[END_REF]. In particular, two-thirds of the warming has occurred at a roughly rate of 0.15-0.20 • C per decade since 1975. The phenomenon of global warming that is basically induced by the rise in global temperature has been the origin of numerous destructive effects on the planet: sea level has risen by 10 to 20 centimeters in certain area, the distribution of precipitation shifted throughout the 20th century, and extreme climate phenomena such as tsunamis or tropical storms became evermore frequent. The main contributor to these detrimental effects was identified as one of the main greenhouse gases, which is Carbon dioxide (CO 2 ). Indeed, in 2013, the atmospheric concentration of CO 2 exceed 400 parts per million (ppm) for the first time in recorded history [START_REF]Climate change: How do we know?[END_REF]. Plants of power, cement, steel or incineration are the main causes for this recent relentless rise in CO 2 concentrations in the atmosphere. However, these industries do not currently have alternative technologies to massively reduce their CO 2 emissions.

Carbon Capture and Storage (CCS) is an auspicious measure for reducing or preventing CO 2 from reaching the atmosphere. The idea consists in separating CO 2 from the emissions of stationary point sources such as power plants or factories and pumping the purified and concentrated fluid into suitable storage sites in the subsurface [START_REF]Total, Carbon capture and storage: the Lacq pilot. Project and injection period 2006-2013[END_REF][START_REF] Monne | The Lacq CCS Pilot, a First[END_REF][START_REF] Singh | A Review of Existing Carbon Capture and Storage Cluster Projects and Future Opportunities[END_REF][START_REF] Kapetaki | European Carbon Capture and Storage Project Network: Overview of the Status and Developments[END_REF][START_REF] Nataly | Energy related CO 2 emissions and the progress on CCS projects: A review[END_REF]. We talk about geological storage as CO 2 is trapped in reservoir rocks under the ground. The storage areas can be of different nature: deep saline aquifers, depleted petroleum fields or coal seams. Deep saline aquifers appear to be the storage reservoir with the highest capacity utilization (400 to 10,000 CO 2 Gton) [START_REF]Club-CO 2 , Carbon Capture, Use and Storage (CCUS): a solution for the future (2011) Dataset accessed 2018-10-31[END_REF]. Nevertheless, their structure and potential to permanently trap CO 2 remains very hard to assess; therefore a significant research effort must be undertaken to bring answers on their long-term potential for the prospects of CCS.

One of the main features of an aquifer site for CO 2 storage is to present a minimal probability of leakage (no natural fracturing, sealing properties of the caprock, well in good condition, etc.) [START_REF] Oldenburg | Coupled Vadose Zone and Atmospheric Surface-Layer Transport of Carbon Dioxide from Geologic Carbon Sequestration Sites[END_REF]. It has been proposed that a rate of 0.01 % per year of CO 2 leakage would be considered as the maximum tolerated risk for a CO 2 sequestration site [START_REF] White | Separation and capture of CO 2 from large stationary sources and sequestration in geological formations coalbeds and deep saline aquifers[END_REF]. Therefore, it is necessary to set up monitoring systems for CO 2 monitoring in order to detect these possible leaks, to assess their magnitude and their impacts and to predict possible corrective actions. Among the different monitoring techniques, the geochemical analysis that is the subject of this article, allows to monitor over time the fluids composition at soil surface and in wells [START_REF] Cohen | The CO 2 -Vadose Project: Experimental study and modelling of CO 2 induced leakage and tracers associated in the carbonate vadose zone[END_REF][START_REF] Gal | Study of the environmental variability of gaseous emanations over a CO 2 injection pilot: Application to the French Pyrenean foreland[END_REF][START_REF] Roux | The CO 2 -Vadose Project: Time-lapse geoelectrical monitoring during CO 2 diffusion in the carbonate vadose zone[END_REF][START_REF] Gilles | Monitoring pour le stockage du CO 2 et la production des hydrocarbures[END_REF][START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF][START_REF] Garcia | The CO 2 -Vadose Project: Numerical modeling to perform a geochemical monitoring methodology and baseline performance assessment for various geochemical variables (gas flux, gas composition, stable isotopes and noble gases) in the carbonate vadose zone[END_REF].

Monitoring CO 2 soil gas has been intensively studied in recent years (e.g. [START_REF] Gal | Study of the environmental variability of gaseous emanations over a CO 2 injection pilot: Application to the French Pyrenean foreland[END_REF], [START_REF] Beaubien | Monitoring of near-surface gas geochemistry at the Weyburn, Canada, CO 2 -EOR site 2001-2011[END_REF]). For example, Myers et al. [START_REF] Myers | Tracers-past, present and future applications in CO 2 geosequestration[END_REF] used gas tracers in particular as precursors of CO 2 leakage to identify the origin of the CO 2 detected. There are few studies on the interaction occurring between the geological storage formation and the gas phase. In their study [START_REF] Oldenburg | Coupled subsurface-surface layer gas transport and dispersion for geologic carbon sequestration seepage simulation, Tagungsbeitrag[END_REF], Oldenburg and Unger explained that depending on the type of release, the attenuation efficiency of the unsaturated zone could be drastically different. The evaluation of the natural flux and soil gas baseline concentration is of primary importance when studying the influence of natural processes on soil gas dynamics. Indeed, natural biological activity can induce variations of several percent in volume content of CO 2 in soil atmosphere ( [START_REF] Lewicki | Surface CO 2 leakage during two shallow subsurface CO 2 releases[END_REF], [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF], [START_REF] Schlomer | Baseline soil gas measurements as part of a monitoring concept above a projected CO 2 injection formation. A case study from Northern[END_REF]). The biological production of CO 2 is due to the organic matter degradation or roots respiration ( [START_REF] Schlomer | Baseline soil gas measurements as part of a monitoring concept above a projected CO 2 injection formation. A case study from Northern[END_REF]). Monitoring techniques for soil gas composition study have been extensively developed as a part of environmental and ecological applications. Different protocols for monitoring gas composition in soil have been developed for application as a routine tool ( [START_REF] Oldenburg | Coupled subsurface-surface layer gas transport and dispersion for geologic carbon sequestration seepage simulation, Tagungsbeitrag[END_REF], [START_REF] Schlomer | Baseline soil gas measurements as part of a monitoring concept above a projected CO 2 injection formation. A case study from Northern[END_REF]). Due to the natural fluctuation of biological activity induced by seasonal variation, the period of investigation to determine the baseline soil gas composition is of primary importance.

For example Gal et al. [START_REF] Gal | Study of the environmental variability of gaseous emanations over a CO 2 injection pilot: Application to the French Pyrenean foreland[END_REF] proposed a period of 18 months of monitoring for baseline evaluation as a reasonable timescale to evaluate the global natural cycle of CO 2 in soil. A baseline soil gas monitoring along several years could be required as well, since seasonal variation could differ year after year, especially in a context of climate change. In the experiment at the ZERT controlled release site, Cranfield experimental site, CO 2 gas was injected in shallow water and the released CO 2 was monitored at the soil surface with repeated grid measurements of soil CO 2 gas fluxes [START_REF] Lewicki | Surface CO 2 leakage during two shallow subsurface CO 2 releases[END_REF]. The authors concluded in their article to careful characterize the background CO 2 variability prior to CO 2 injection into the storage reservoir, to limit the total area of investigation by focus on features most susceptible to leakage (e.g., wells, faults) and to use various complementary CO 2 measurement techniques in a program of storage site monitoring.

Similarly, water-mineral chemical reactions can be a source or a sink of CO 2 . The carbon dioxide concentration is also particularly affected by soil temperature due to the strong dependence between temperature and chemical/biological reaction rates ( [START_REF] Barron-Gafford | The relative controls of temperature, soil moisture and plant functional group on soil CO 2 efflux at diel seasonal and annual scales[END_REF], [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF], [START_REF] Schlomer | Baseline soil gas measurements as part of a monitoring concept above a projected CO 2 injection formation. A case study from Northern[END_REF]). Concerning the liquid phase and the interface between the liquid phase and the gas phase (in the capillary zone in a natural context at subsurface location, between the saturated and the non-saturated zones), several techniques have also been developed (examples of sampling tubes, in-situ probes) but only few studies have put in evidence the influence of natural parameters variations (porosity, permeability, level of water table) on the CO 2 migration, its detection, and the degassing phenomenon that can occur.

Herein, this Aquifer-CO 2 Leak project is part of this research work necessary for the monitoring, the verification-characterization, and the geochemical detection of CO 2 leaks as well as the prediction of CO 2 behavior from reliable numerical models in near-surface environments, more particularly, in the saturated zone. It consists in experimenting CO 2 injection in not deep aquifer (phreatic) at Saint-Emilion in France on a site already well characterized from the hydrogeological and geochemical point of view. The aim is to get better understanding of CO 2 behavior in the saturated zone, to develop a know-how on the different aspects related to a CCS project (such as site characterization, injectivity control and leakage detection), and to test various geochemical and geophysical monitoring equipment and methods. In this article, we focus on the geochemical approach of this project that tests noble gases and a fluorescent molecule as tracers in the aquifer.

The main objective of this paper is to present the predictive modeling results of this small pilot-scale CO 2 injection experiment. The simulations have been carried out to design the pilot with the following purposes: (1) determining the volume of injected water enriched with dissolved CO2 and tracers, the injection rate and the amount of CO 2 and tracers to be injected in the saturated zone in such a way that their migrations can be followed with a sufficient precision, (2) predicting the carbonated water plume and the total time to return to the initial state of the domain, (3) determining the optimal location of the different geochemical sensors well for sufficient measurement coverage, and (4) defining as precise as possible a schedule for the measurements and/or sampling.

Site description and selection of the zone of interest

Saint-Emilion experimental site

The experimental site is an abandoned underground limestone in Saint-Emilion, in Nouvelle Aquitaine (France), exploited by "the rooms and pillars" method. It consists in an immense 70-hectare network of underground galleries that form a gigantic network of paths and developed on 3 levels. The carbonate formation is dated at Upper Oligocene (Stampien, 28-30 Ma) in Nouvelle Aquitaine basin and is composed of wackestone to grainstone facies characterizing an inner intertidal ramp platform. These carbonate facies are affected by a deep meteoric diagenesis process during 28-30 Ma. This geological emmersif event develops a good secondary porous network and a good porosity and permeability in the studied carbonate formation. The map of this underground quarry is shown in Figure 1 (a). The black lines represent the limits of each underground rooms or gallery. Pillars in gray are displayed on the map. The unsaturated zone, of about 18 m in depth, overlies the water table. The level of the water table Z wt is given by discontinuous blue lines and expressed in reference to the General Levelling of France (NGF). P5 and P6 denote respectively the reference well and the control well initially drilled in the site. The areas where drilling is possible (the height of the cavity is greater than 2 m) are marked in green and the areas where drilling is very complicated (the height is less than 2 m) in yellow. The underground water comes mainly from surface water infiltration due to effective precipitation. Some piezometric measurements were made near the reference well P5 that allowed us to determine the seasonal changes in the water table level and to estimate the hydraulic gradient at this area of the experimental site. More details about the geological properties can be found in [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF] [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF] and [START_REF] Garcia | The CO 2 -Vadose Project: Numerical modeling to perform a geochemical monitoring methodology and baseline performance assessment for various geochemical variables (gas flux, gas composition, stable isotopes and noble gases) in the carbonate vadose zone[END_REF] [START_REF] Garcia | The CO 2 -Vadose Project: Numerical modeling to perform a geochemical monitoring methodology and baseline performance assessment for various geochemical variables (gas flux, gas composition, stable isotopes and noble gases) in the carbonate vadose zone[END_REF]. All the gathered measurements are reported in Figure 1 (b) that depicts a cross-section of the part between the substratum and the gallery floor along the discontinuous black line shown in Figure 1 (a).

The domain is composed of 4 layered parts: (1) a layer of thickness that varies between 0 and 0.5 m of backfill; this backfill was added in order to ensure a flat gallery floor,

(2) a first layer of limestone with a maximum thickness of 1 m; the porosity in this layer varies between 0.35 and 0.45 and the permeability ranges from 10 to 14 D, (3) a second layer of limestone with a thickness of 3.5 m; the porosity is between 0.2 and 0.35 and the permeability varies from 6 to 10 D, and finally (4) an impermeable and nearly flat layer of clay named the sannoisian green marls; this layer is located at a total depth of 19.4 m (the soil altitude is 79.9 m NGF). Two extreme water table levels were measured depending on the hydrogeological recharge, a highest (January to February) and a lowest (June to July) water table. The maximum oscillation of the water table (zone of intermittent saturation) is 1.5 m. Furthermore, the measured hydraulic gradient i is of 5 % regardless the water table level.

Zone of interest

In order to define our simulation model, we have chosen a zone of interest on the site. The underlying idea is to cover the maximum of the green zone where the drilling is possible for both injection and observation. The selected zone is shown in red square in Figure 2 (a). Note that the dimensions of the domain are 20 m × 20 m and the model axes are oriented with the zone geometry. The Y-axis is parallel to the water flow direction in the aquifer (according to the hydraulic gradient), the X-axis is perpendicular to it and the Z-axis is set at the gallery surface. Furthermore, in order to estimate the water table level in this area, we extrapolated the measured hydraulic gradient (i = 5 %) to this zone and we kept the same flow direction in all areas. As shown in Figure 2 (a), the injection well is centered on the transverse (with respect to the flow direction) X-axis, that is its coordinates are (X = 10 m,Y = 14.8 m). Two

CooresFlow Software and modeling framework

General characteristics

In our publication, the numerical simulations allowing to realize the design of the pilot were carried out with CooresFlow ( [START_REF] Gallo | Longterm flow simulation of CO 2 storage in saline aquifers[END_REF], [START_REF] Estublier | Long-term fate of CO 2 in a saline aquifer: modeling issues[END_REF], [START_REF] Estublier | Simulation of a Potential CO 2 Storage in the West Paris Basin: Site Characterization and Assessment of the Long-Term Hydrodynamical and Geochemical Impacts Induced by the CO 2 Injection, Oil and Gas Science and Technology[END_REF]), IFPEN's in-house research simulator for reactive transport in porous media. It is developed in parternship with international industries and academics to model multiphase transport processes coupled with fluid-rock interactions, for various applications including CO 2 underground storage as well as enhanced recovery of hydrocarbon reservoirs. Numerous numerical simulators model non-isothermal multiphase multicomponent transport processes with or without fluid rock interaction ([27], [START_REF] Nordbotten | Uncertainties in practical simulation of CO 2 storage[END_REF]). Amoung them, CooresFlow aims to be a useful tool either for understanding complex processes occuring in porous media and for modeling them, or for simulating the fluid flow, the species migration and the rock reactivity for real sites with a geological model integrating both the structural complexity and numerous lithologies. Thus, it proposes particular functionalities: (1) modular physical models such as generic (e.g. EOS of Peng-Robinson) or user-defined (e.g. a new mathematical formula modeling a particular thermodynamical behavior) models, (2) unstructured meshes, (3) a friendly interface and (4) different complexity levels of coupling between the fluid flow model and the geochemical model associated to a parallel computing.

Mathematical formulation of a compositional monothermal single-phase flow model

The study concerns the simulation of water injection containing the following dissolved species (CO 2 , krypton and fluorescent molecule) in the saturated zone. To carry out this study, we used the following physical model: a compositional single-phase flow model that accounts for advection, dispersion and molecular diffusion within nonthermal context. It is worth noting that in this work, we will not consider the coupling between flow and reactive transport mechanisms. Indeed, the chemical reactions that may occur between solid and liquid have been neglected in a first time for the sake of simplicity. However, the rock-fluids interaction will be considered in the future works of this project and the effect of mineral dissolution and precipitation on a long term CO 2 migration will be studied and analyzed.

The compositional flow model used in this study is based on the following material balance equations:

∂C α j ∂t + div(J α j ) + Q α j = 0 ( 1
)
where α is the water phase, j the species index, C α j the molar concentration of the species j [mol.m -3 ], J α j the vector of molar flow rate of species

j [mol.s -1 .m -2 ], Q α j the source term (well) [mol.s -1 .m -3 ].
In addition, the following laws are used to describe fluid flows and transfers in the porous medium: Generalized Darcy's law for fluid flow:

V α = - K µ α (∇P α -ρ α g) (2)
Species transport by advection:

J α j, adv = C α j V α (3)
Fick's law [START_REF] Fick | On liquid diffusion[END_REF] for species diffusion/dispersion of species:

J α j, Diff = -D α j + Disp α • ∇C α j ( 4
)
Where P α is the pressure of the phase α [Pa], V α is the velocity of the α phase

[m.s -1 ], K is the medium permeability [m 2 ], µ α is the viscosity of the α phase [Pa.s],
ρ α is the density of the α phase [kg.m -3 ], g is the gravity constant [m.s -2 ], D α j is the diffusion tensor of a species j in the α phase [m 2 .s -1 ], Disp α is the dispersion tensor in the α phase [m 2 .s Concerning the initial conditions, the temperature T was set at 14 • C and the pressure P at the hydrostatic pressure as P = 1 bar at a depth of 1.25 m (at the higher situated east boundary), and the porous medium was supposed to be fully saturated with water carrying 15.5 mg/l of dissolved CO 2 . In addition, the water density ρ w and viscosity µ w were assumed to be constant at 1000 kg/m 3 and 1 cP, respectively, whatever the water is carrying CO 2 and krypton or not. Indeed, the impact of water density and viscosity differences in presence of dissolved gases on the flow is expected to be negligible in the studied case since at 1 bar and 14 • C the density of water enriched with CO 2 reaches only 1001 kg/m 3 and the viscosity is close to 1 cP.

The petro-physical properties of the domain were fixed at the minimum values measured in the layer 2 of limestone (see section 2.1, [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF]). The permeability k and the porosity φ were set constant as k = 6 D and φ = 0.2 ( [START_REF] Loisy | The CO 2 -Vadose Project: Dynamics of the natural CO 2 in a carbonate vadose zone[END_REF]) because of lack of heterogeneity measurements in the zone of interest. Table 1 summarizes the physicochemical properties of the different species dissolved in water that are considered in this simulation study. To complete the physical model, we have chosen at this stage of this study the following values for the field dispersivity by considering the scale 20m of the domain: the longitudinal dispersivity is as α L = 1 m and transverse dispersivity as α T = 0.33 m according to [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF] [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF].

A 8 cm diameter well was placed as given by Figures 2. The injection rate was estimated assuming that the injection was over the entire column at a constant volumetric flow rate and the injection velocity was close to the initial water velocity in the water table such that we limit overpressure around the well. Thus, the calculated injection

                
x H2O = 9.98 × 10 -1

x CO2 = 6.21 × 10 -4

x kr = 5.65 × 10 -5

x fluo = 1 × 10 -3 =⇒                  w H2O = 997, 220 (ppm)
w CO2 = 1, 515 (ppm)
w kr = 262 (ppm)

w fluo = 998 (ppm)
On one hand, by following these imposed values, the CO 2 and the krypton were injected at the maximum dissolved amount under the fixed temperature and pressure conditions. On the other hand, the injected amount of the fluorescent molecule was set at an arbitrarily chosen low value.

Results and discussions

In this section, we present the obtained simulation results and the discussions needed, by following this sequencing: first, a physical description of species concentra-tion dynamics in the cross-section YZ (Figure 2 (a)) is presented in both the injection and the post-injection stage. Then, we study the temporal evolution of the CO 2 mass fraction in the three observation wells and its plume at several time steps. Finally, we consider the concentration variation of the components across the XZ cross-section in order to better understand their migration and their spreading along the longitudinal and the transverse direction to the main water flow. Once the injection is stopped, the CO 2 concentrations spread rapidly in the direction of flow as the simulated profiles show smooth transitions with a maximum value that does not exceed 40 ppm after 3 days. Indeed, as time increases, the concentration profile decays, advects in the aquifer flow direction with high velocity, and spreads outward.

After 10 days, the injected CO 2 expands over the entire length of the section attesting for a great dispersion effect. This important concentration spreading is essentially due to dispersion, and more particularly to the longitudinal dispersivity, which remains the preponderant phenomenon in saturated ground-water flow as it dominates molecular diffusion.

We found that the same observations can be made for the krypton migration along the Y-Axis, as shown in Figures 3 (c) and (d). Similarly, the concentration of the krypton increases steeply during the injection stage to reach the maximum value w kr = 190 ppm, then it spreads rapidly and follows the same trend observed in CO 2 migration.

The only differences that have been found concern the maximum concentration values.

For both CO 2 and krypton, the maximum concentration corresponds to 73 % of the initial concentration that confirms the process of migration is mainly advection and dispersion and not diffusion. This observation applies also to the fluorescent molecule transport. For this, in the following of this work, we focus our analysis on the migration of a single component, which is CO 2 , since the same remarks remain valid for the various species. We are interested now in studying the temporal evolution of CO 2 concentration in the three observation wells (Obs1, Obs2 and Obs3). The obtained results are shown in Figure 4. Few hours after the injection stop, the injected CO 2 breaks through the first observation well (Obs1), then increases in time until reaching a maximum value of 37 ppm after 2 days, and finally, decreases rapidly to the initial value (15.5 ppm).

Concerning the second observation well (Obs2), the maximum value of CO 2 mass fraction that can be measured is 26 ppm and it can be recorded after 4 days. Located at 10 m from the injection point and at 4 m from the CO 2 plume center, the (Obs3) well is reached by CO 2 due to the transversal dispersion. The maximum value is 17 ppm at 9 days which is close to the initial value (15.5 ppm). The Figure 4 also shows that the system returns at its initial state after 20 days.

To get an idea about the dissolved CO 2 migration and plume temporal evolution, we present in Figures 5 (a To get more clarity about CO 2 migration, we study the dynamics of CO 2 mass fraction in the three cross-sections XZ of the observation wells. The corresponding results are reported in Figure 6. The maximum values that can be measured in each well and the plume extensions are still clear in this figure. Still some additional remarks can be made according to the CO 2 distribution. We found that the simulated profiles show a spreading Gaussian cloud which is centered on the main water flow direction and advances due to advection and dispersion as noted before. Indeed, as a function of time, the maximum value decreases between the three sections while the width of the spreading increases significantly. Note that the maximum value shown in Figure 6 (c) is 22 ppm which is different to the one recorded in Obs3 as the latter well is located at

Sensitivity study

The objective in this section is to provide information on the uncertainty associated with model simulations and to identify the most important parameters that have a strong influence on the model outputs. We begin this section by a discussion and a selection of the uncertain parameters in our study case which are due to lack of knowledge or even of measurements. Then, we focus on the one hand on the effect of the aquifer characteristics, which are the dip angle of the aquifer top-surface linked to the water table level, and on the other hand, on the effect of the petro-physical properties of the porous media (permeability, porosity, and dispersivity) on the CO 2 flow migration. In particular, the sensitivity results in terms of species concentration variations in space and time are presented and discussed in each simulation case.

Selection of the uncertain parameters

Until now, several parameters were estimated based on data gathered from surveys conducted on the site next to the zone of interest and from the literature or even by interpolating the available data. Indeed, in the primary model presented above, the permeability k and the porosity φ were fixed at 6 D and 0.2, respectively. In addition, we considered that the level of the water table is the highest measured one because we have scheduled to inject CO 2 in the aquifer in summer, period of high water table .   The longitudinal dispersivity (α L = 1 m) and the ratio α L /α T = 3, where α T denotes the transverse dispersivity, were also estimated from the literature [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF]. Finally, we note that we considered a dip angle of 5 % along the Y-axis and none along the X-axis (no available data and not expected to be high) in our model which remains a strong assumption that we must verify during the scheduled CO 2 injection. For the sake of clarity, we will refer to this case as the "reference" case.

Accordingly, we chose to conduct our sensitivity analysis by performing the following simulations by:

• adding a dip angle of 1 % along the X-axis, such that the groundwater flow will be changed and becomes two-directional flow, dips angles Y5%-X1%). In the reference case, the plume is elliptical and centered on the injection Y-axis (i.e. the main flow direction), whereas in the other case, the plume spreads over the horizontal plan and advances in the side of the low altitude of aquifer top-surface (from the upper-right to the lower-left boundaries). We note also that the width of the spreading is more important in the case of two dips angles (almost equal to 7 m versus 4 m in the reference case) as the one along the X-axis comes to enhance the effect of the transverse dispersion. Thus, 1 % dip angle along the X-axis changes significantly the plume shape as it is advected with the flow. Note that in this case, the third observation well becomes useless since the plume passes on the other side and in such a case we risk to observe nothing in this well. Or, if the dip angle was in the opposite direction, this third well becomes very important for the measurements. For this reason particular importance should be given to the well placement with regard to the value of this parameter which has to be estimated as the drilling progresses. The effect of the dip angle on the evolution of CO 2 concentration can be quantitatively evaluated via Figure 8. It shows the variation of the CO 2 mass fraction with time at the three observation wells according to each simulation case. As expected, the impact of the dip angle is clearer and more significant on the measurements in the third well, while in the other wells, the results are very similar to those obtained in the reference case. To finish with the effect of the aquifer characteristics, we consider the second parameter which is the water table level. Figure 8 shows that the model simulations are very sensitive to this parameter in virtually the three wells since the maximum values are reached in this case. Regarding the other parameters, dispersivity appears to be the parameter that has the most impact on CO 2 migration compared to permeability and especially when the longitudinal dispersivity is set at very low values (as given by the case α L = 0.3 m).

Let us now consider the sensitivity results: in the first well, the maximum value that can be measured is 51 ppm in the case of lowest water table or when α L = 0.3 m, and a minimum value of 28 ppm can be obtained in the case of k = 10 D or α L = 3 m.

However, the time of reaching these maximum values does not vary too much and it ranges over 1-3 days. Almost the same dependence can be noticed in the second well as the maximum value varies from 21 to 35 ppm and the time necessary for reaching it is after 3 to 5 days from the injection. For the third well, we found that the lowest water table has a significant impact on the temporal evolution of the CO 2 concentration; in this case, the maximum value is 19.5 ppm and it is reached after 9-10 days. This indicates a large spreading of the CO 2 plume in the transverse main direction of the water flow, as it will be confirmed hereinafter. This expected result is explained by the decrease of the transversal dispersivity flux along the vertical axis induced by the decrease of the aquifer thickness in case of low water table .   Figure 9 shows the variation of the CO 2 mass fraction across the sections XZ of the observation wells. These profiles are presented at the times ensuring the maximum concentration in each simulation case.

The spreading of the CO 2 concentrations is still clear in all cases and takes the form of Gaussian curve, but with varying widths at half maximum that range over 3-4 m in Obs1, 4-6 m in Obs2, and 4-12 m in Obs3. The maximum width is always obtained in case of simulation with lowest water table level. The dependence of the main model predictions on CO 2 concentration, that are the maximum CO 2 mass fraction that can be measured in each well (w max ), the time of reaching this peak (t max ) and the full width at half maximum of the Gaussian spreading (L c ), are summarized in Table 3. The method to calculate the full width at half maximum (L c ) is outlined in Figure 10. Note that the latter summary table is very useful in the field at the time of the experiment since it provides the times at which sampling and measurements can be performed without losing any precision and gives insights into CO 2 plume position and velocity at which it advances. The table also checks that in all the studied cases, CO2 and tracers can be detected by the measuring equipment since CO2 concentration is greater than ppm scale and krypton concentration is greater than ppb scale.

As a conclusion of this study, we advice to proceed with the following steps for placing the observation wells in the Saint-Emilion experimental site. During the well drilling phase, we should start by drilling three wells forming a rectangular triangle to measure the water table level at each well. The estimation of the dip angle along the X and Y-axis from these measurements will define the main direction of the aquifer flow. Then, other wells will be drilled to perform additional measurements to accurately monitor the CO 2 plume. They will be placed to cover a zone as a diamond-shape oriented in the main flow direction to optimize the site monitoring. Indeed, the species migrate both by advection and longitudinal and transversal dispersion. The diamond dimensions will be defined from the simulated species concentration taking into account the parameter uncertainty. Thus, a new fluid flow simulation using the measured data of the water table level and of the two dip angles should be performed to reduce the uncertainty of the diamond size.

Conclusions

The aim of this work was sizing both in time and space the planned experiment that consists in injecting dissolved CO 2 and tracer gases in water in the saturated zone of the Saint-Emilion experimental site. This objective was met by first determining the injection flow rate and the injected volume considering the site constraints. Thus, 200 l of solution containing the studied dissolved species was assumed to be injected at a volume rate of 2,520 ml/min. Then, the analysis of the simulation results allowed us to determine the size plume evolution with time, the date when the maximum of each species concentration is reached for three observation wells and the return time to the initial state of the system. The CO 2 plume size analysis showed that the average full width at half maximum is 3.5 m at the first observation well and 4.5 m at the second one, two days later due to transversal dispersion. We also observed that the maximum of CO 2 concentration is reached at about 2 days at 3 m and 4 days at 6 m from the injection well with a return at the initial state after 20 days.

Then, the sensitivity analysis of this work leads us to quantify the effect of the uncertain model parameters on the size plume and the date when the maximum of each species concentration is reached as well as to identify the sensitive parameters which have the strongest influence on the model predictions. The studied uncertain parameters were the dip angle along the X-axis, the longitudinal and transversal dispersivity, the water table level and the permeability and porosity parameters. The study concluded that the most sensitive parameters of this model are the dispersivity and the water table level. The permeability and porosity parameters and the dip angle along the X-axis have less effects. However, this study showed that this dip angle is a very important data that is necessary to measure to avoid drilling useless observation wells.

Finally, the results of this work come to verify that the proposed locations for injection and observation wells allow us to provide sufficient coverage of the measurements while considering the technical constraints for drilling on site and under some simplifying assumptions that have been studied and evaluated in this work. We also advice to start the well drilling phase by the measurement of the water table level at the three first wells and by the estimation of the dip angle along the X and Y-axis from these measures since these parameters are important for this site to know respectively the area dimensions for monitoring and the fluid flow direction. Other wells will be drilled to perform additional measurements to accurately monitor the CO 2 plume. They will be placed to cover a zone as a diamond-shape to optimize the site monitoring. The final position of these wells will be specified in the next papers. In total, eight wells are planned to be built for this site.

Limestone 1 :𝜙=0Figure 1 :

 11 Figure 1: (a) Map of Saint-Emilion site (underground view): the numbered pillars in gray are displayed on the map. Three different water table level (Z wt = 62, 62.5 and 63 m) were measured and shown in the figure by discontinuous blue lines leading to an hydraulic gradient i equal to 5 %. The flow direction is from the top left to the bottom right as given by the big blue arrow. P5 and P6 denote respectively the reference well and the control well initially drilled in the site. The areas where drilling is possible (the height of the cavity is greater than 2 m) are marked in green and the areas where drilling is very complicated (the height is less than

  observation wells are positioned parallel to the Y-axis and aligned with the injection well: the first one is at a distance of 3 m and the second at 6 m from the injection. A third observation well is placed at the following coordinates (X = 14 m,Y = 4.8 m), that is an offset of 4 m with respect to the flow direction. In order to get a clear idea about the studied zone disposition, we present in Figure 2 (b) a conceptual scheme of the underground at the YZ section. The aquifer shape was depicted by a top-angled trapezoid model geometry where the bottom is horizontal, and the left and the right corner are vertical. The dip angle of the aquifer top represents the hydraulic gradient i which was set at 5 %. Following this representation, the water table level Z wt in the domain varies between 62.75 and 63.75 m (NGF) leading to a maximum depth of 2.25 m with reference to the gallery floor. Note that the highest water table was chosen to depict the aquifer top. The model include an ambient formation water flow from the higher situated east boundary to the lower west boundary. Furthermore, the aquifer petro-physical properties are represented by the second layer of limestone (cf. Figure 1 (b)).

Figure 2 :

 2 Figure 2: Selection of the Zone of interest: (a) placement and dimensions of the selected zone on the experimental site of Saint-Emilion. The selected zone is shown in red square which covers the maximum of the green areas that offer an easier drilling for injection and observation. The position of the injection well is pointed by a black arrow and the three observation wells by red ones (Obs1, Obs2 and Obs3); (b) conceptual scheme of the YZ section of the selected zone that depicts the aquifer shape.

  CooresFlow is a THC (thermal hydro-chemistry) simulator that considers the following phenomena: (a) multiphase flow in porous media, with viscous and capillary forces, (b) transport of chemical components by advection, diffusion and dispersion, (c) transfers within the fluids or on the surface of the rock governed by local equilibrium, (d) chemical equilibrium or kinetic reactions to describe geochemical exchanges, (e) heat convection and thermal conduction, and (f) dynamic modification of the porosity and permeability of the porous medium over time.

- 1 ]

 1 characterized by α L the longitudinal dispersivity coefficient [m] and α T the transversal dispersivity coefficient [m].

4 .

 4 Figure 2 (b) shows the setup that was chosen for the three-dimensional simulation model. The simulation geometry was finely-meshed along the three axes X-Y-Z with constant space steps discretization: ∆X = 20 cm, ∆Y = 20 cm and ∆Z = 16 cm, leading to 100 × 100 × 20 cells. The cells size was carefully chosen to capture the physical phenomena detailed above as well as to limit the numerical dispersion. The domain is limited by boundary conditions: a hydrostatic pressure on both sides y = 0 and y = 20 m and a zero flux condition on the other faces, leading to an underground water velocity of 0.15 m/day (along the Y-axis). Concerning the initial conditions, the temperature T

  = 4.2 × 10 -5 m 3 /s = 2,520 ml/min. The injected water volume was estimated taking into account the porous volume of the pilot and the difficult conditions for moving experimental materials. Thus a total injected volume of 200 l was considered that represents 0.1 % of the pore volume PV (PV = 220 m 3 ) as the porosity is equal to 0.2. It leads to an injection time of 78 min. After the injection stage, we close the well and we simulate the migration of the dissolved species over 20 days. The imposed molar fraction x i and mass fraction w i (ppm) of each species i dissolved in the injected water are as follows:

Figures 3 (

 3 Figures 3 (a) and (b) show the dynamics of CO 2 mass fraction in the cross-section YZ of the injection well during both the injection and the post-injection stage. When the injection starts, the concentration of CO 2 increases by steps in the vicinity of the injection well until reaching a maximum value of w CO2 = 1, 100 ppm after a total injection time of 78 min. The latter maximum value is a little smaller than the imposed one due to the dispersion phenomenon taken into account in this study. At the end of the injection stage, the injected CO 2 is present at a maximum distance of 1.5 m from the injection well and even behind it (at the opposite of main water flow direction) due to a low overpressure induced by injection, but it does not reach the observation wells.

Figure 3 :

 3 Figure 3: Dynamics of CO 2 and krypton mass fractions along the section YZ of the injection well, during: (a) and (c) the injection stage, (b) and (d) the post-injection stage.

  ), (b) and (c) screenshots that show a top view of the CO 2 concentration distribution at times when the maximum of CO 2 is reached at each observation well during the post-injection stage. The effect of the longitudinal and transverse dispersivity is clear on CO 2 migration. At 2 days, the plume is of approximately elliptic shape with dimensions as follows: L L = 3 m and L T = 1 m, where L L denotes the longitudinal (main axis) spreading and L T the transverse (minor axis) one. After 4 days and when it reaches the second well, the CO 2 concentration spreads and the plume gains enormously in size as its dimensions become L L = 6 m and L T = 2 m. Finally, after 8 days from the injection, the CO 2 plume advances mainly along the Y-axis and covers a huge part of the domain (L L = 12 m and L T = 4 m) but with very low concentrations. Thus, the dissolved CO 2 plume just reaches the third observation well (Obs3) because it is located at 4 m from the plume center. It reaches the model boundary already after less than 10 days.

Figure 4 :

 4 Figure 4: Evolution of the CO 2 mass fraction in water versus time (in the observation wells).

Figure 5 :

 5 Figure 5: Zoom-in at the top view of the dissolved CO 2 plume at times when the maximum of CO 2 concentration is reached at each observation well during the post-injection stage: (a) after 2 days, (b) after 4 days and (c) after 8 days. The scales are adapted by figure.

Figure 6 :

 6 Figure 6: Dynamics of CO 2 mass fraction in the three cross-sections XZ of the observation wells: (a) Obs1 at 3 m, (b) Obs2 at 6 m and (c) Obs3 at 10 m from the injection well. Note that Obs3 is located at X = 14 m in figure (c), i.e. 4 m from the CO 2 plume center

Figure 7 :

 7 Figure 7: Zoom-in (top view) at the effect of 1 % dip angle along X-axis on the dissolved CO 2 plume comparing to the reference case at t = 10 days: (a) reference case with Y5% and (b) top aquifer with two dips angles Y5%-X1%.

Figure 8 :

 8 Figure 8: Evolution of the mass fraction of CO 2 versus time at the three observation wells: (a) at Obs1, (b) at Obs2 and (c) at Obs3. Note that the "reference" case corresponds to the simulation with k = 6 D, φ = 0.2, highest water table, α L = 1 m and single dip angle of 5 % along the Y-axis.

Figure 9 :

 9 Figure 9: Variations of the mass fraction of CO 2 in the sections XZ of the observation wells (at the time ensuring the maximum of CO 2 concentrations for each parameter, cf.Figure 8): (a) Obs1, (b) Obs2 and (c)

Figure 8 )

 8 : (a) Obs1, (b) Obs2 and (c)

Figure 10 :

 10 Figure 10: Method to calculate the full width at half maximum (L c ).

  single dip angle of 5 % along the Y-axis. (ppm) (m) (days) (ppm) (m) (days) (ppm) (m) (days) (ppm) (m) (days) (ppm) (m) (days) (ppm) (m)

  This study constitutes the first step in the modeling work of the Aquifer-CO 2 Leak project. Several points will be addressed within the future works after the injection experiment such as (1) the reactive transport modeling to consider the fluid rock interactions that occurs in carbonates of Saint Emilion, (2) the calibration of the fluid flow model by adjusting the uncertain parameters that was discussed herein and the calibration of the geochemical model and (3) the interpretation of the experimental measurements.The simulation work presented in this article and the scheduling of the future one after the experiment in Saint-Emilion show that simulation is essential at different steps of a project life whether it is a small-scale pilot as Saint-Emilion or a CO 2 storage at large scale. Simulation is necessary (1) to design the site for both the injection process (well position, rate, ...) and the monitoring network (observation wells, soil sensors to detect surface movement or gas leakage, ...) and (2) to predict the CO 2 plume migration and its impacts. At each step, the model is more complicated by integrating static and dynamic information measured in-situ and the uncertainties of the simulation results are reduced.

Table 1 :

 1 Physico-chemical properties of species dissolved in water at T = 14 • C and P = 1 bar.

		Molar mass Diffusion coefficient Solubility in water
		(g.mol -1 )	in Water (cm 2 .s -1 )	(KI)
	H 2 O	18.01	2.30×10 -5	-
	CO 2	44.01	1.41×10 -5	1,610
	krypton	83.80	1.51×10 -5	17,700
	fluorescent molecule 18.01	2.30×10 -5	-
	rate is Q inj

Table 3 :

 3 Summary of the sensitivity results of CO 2 and krypton migration to the uncertain parameters. The "reference" case refers to to the simulation conducted with k

	6
	=

m) in yellow; (b) Conceptual scheme of the aquifer near the reference well P5 and along the discontinuous black line shown in (a).

Acknowledgments

This research work was conducted within the Aquifer-CO 2 Leak project which is funded by the ADEME and the Région Nouvelle Aquitaine from France. The authors are grateful to IFPEN and research unit EA 4592 G&E ENSEGID-Bordeaux-INP for permission to publish this work and also all the partners of this project for useful discussions.

• considering the lowest level (in winter) of the water table Z wt : Z wt varies between 61.25 and 62.25 m (NGF) leading to a maximum depth of 3.75 m with reference to the gallery floor,

• using the maximum measured values of the permeability and the porosity in layer 2 of limestone (that is k = 10 D and φ = 0.35),

• considering the two limit values (maximum and minimum) of the longitudinal dispersivity found in the literature for the scale of our model [START_REF] Gelhar | A critical review of data on field-scale dispersion in aquifers[END_REF]:

Note that we chose these extreme values in order to figure out the possible variation range of the model outputs. In the following, single parameters were varied while the others of the "reference" case were kept constant in order to emphasize the effect of each parameter. Table 2 summarizes all the parameters values used in this sensitivity study.

Table 2: Parameters values selected for the sensitivity study. The reference values are marked with the special character † . Note that the Y-axis is parallel to the main water flow direction in the aquifer (according to the hydraulic gradient), the X-axis is perpendicular to it.

Properties

Uncertain

Sensitivity Results

In this section, we start with an analysis of the possible effect of 1 % dip angle along the X-axis on CO 2 plume. For this, we compare in Figure 7 the CO 2 plume at 10 days as obtained in each case (the reference case with Y5% and the case with two