Sébastien Le Crom
email: lecrom.sebastien@gmail.com

Christophe Tournassat

Jean-Charles Robinet

Virginie Marry
email: virginie.marry@sorbonne-universite.fr

Influence of Water Saturation Level on Electrical Double Layer Properties in a Clay Mineral Mesopore: a Molecular Dynamics Study

Molecular dynamics simulations were conducted to evaluate the effect of water desaturation on ions and water density distribution and diffusion properties in a clay mineral mesopore. Simulations were carried out using the Polarizable Ion Model (PIM) force field. Desaturation of the pore was responsible for significant changes in localization and mobility of water molecules and ions, especially in most unsaturated conditions, for which charge overcompensation by cations adsorbed on the surface through inner-sphere complexes were observed. While cations and anions mobility decreased monotonously with increasing desaturation, pore-averaged diffusion coefficient of water molecules exhibited a complex evolution with an increased diffusion coefficient value at low desaturation level, and then a drop of this value at high desaturation level.

Water diffusion was influenced by two antagonistic processes, i.e. acceleration of water molecules diffusion at water/air interface; and slowdown at water / clay mineral interface. These observations can be put in perspective with differences in macroscopic diffusion properties among cationic, anionic and water tracers observed in experiments conducted on desaturated clay mineral rock samples.

Introduction

Clay minerals have a wide range of industrial applications from catalytic activities to energy and ecological engineering. In particular, clay-rich geological formations are considered as potential host rocks for radioactive waste repositories, and as cap-rocks for CO 2 sequestration in geological settings. [START_REF] Tournassat | Natural and Engineered Clay Barriers[END_REF][START_REF] Bourg | Sealing Shales versus Brittle Shales: A Sharp Threshold in the Material Properties and Energy Technology Uses of Fine-Grained Sedimentary Rocks[END_REF][START_REF] Neuzil | Can Shale Safely Host U.S. Nuclear Waste? Eos[END_REF] Clay mineral particles are constituted of stacks of aluminosilicate layers. In clay minerals such as smectites and illite, layers are negatively charged because of heterovalent isomorphic substitutions taking place in their structure. This negative charge is compensated by cations located next to the external surfaces of the particles as well as in their interlayer spaces. [START_REF] Bergaya | Handbook of Clay Science. General Introduction: Clays, Clay Minerals[END_REF][START_REF] Bergaya | Handbook of Clay Science. General Introduction: Clays, Clay Minerals[END_REF] At the external surfaces that are in contact with pore aqueous solution, a structuring of the aqueous phase forms, which is commonly called electrical double layer (EDL), and which differs from the bulk solvent phase. In EDL, anions and cations are not distributed uniformly because of charge balance requirements, but also because of specific Van der Waals interactions between the mobile species (water and ions) and the solid. [START_REF] Sposito | The surface chemistry of natural particles[END_REF] The EDL is usually conceptually subdivided in, at least, two regions: the compact region next to the surface, and the diffuse layer (DL) region. [START_REF] Sposito | The surface chemistry of natural particles[END_REF][START_REF] Lee | Hydrated cation speciation at the muscovite (001)-water interface[END_REF][START_REF] Lee | Monovalent ion adsorption at the muscovite (001) solution interface: relationships among Ion coverage and speciation, interfacial water structure, and substrate relaxation[END_REF][START_REF] Marry | Structure and dynamics of water at a clay surface from Molecular Dynamics simulation[END_REF][START_REF] Rotenberg | Water and ions in clays: Unraveling the interlayer/micropore exchange using Molecular Dynamics[END_REF][START_REF] Bourg | Molecular Dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl-CaCl2) solutions[END_REF][START_REF] Tournassat | Comparison of Molecular Dynamics simulations with triple layer and modified Gouy-Chapman models in a 0.1 M NaCl -montmorillonite system[END_REF] A reliable description of the EDL properties is necessary to model transfers of water and aqueous solutes in clayey media because of the strong coupling between surface electrostatic properties and fluxes of solvent and ions. [START_REF] Tournassat | Reactive transport modeling of coupled processes in nanoporous media[END_REF] In the context of underground radioactive waste storage or CO 2 sequestration, unsaturated transient conditions are predicted to take place during both operation and after-closure phases. Modeling fluid and solute transport in the presence of unsaturated conditions is thus an important aspect to make reliable predictions of the evolution of the above mentioned systems with time. [START_REF] Churakov | Mobility of Na and Cs on Montmorillonite Surface under Partially Saturated Conditions[END_REF][START_REF] Gimmi | Water retention and diffusion in unsaturated clays: Connecting atomistic and pore scale simulations[END_REF] As water and ions distribution and dynamics in the EDL cannot be measured directly, [START_REF] Bourg | Stern layer structure and energetics at mica-water interfaces[END_REF] atomic simulations such as molecular dynamics simulations (MD) are widely used for the purpose of characterizing ions distributions and dynamics in nano-and mesopores. In the present study, we investigated the evolution of the electrical double layer properties as a function of increasing desaturation, using molecular dynamics simulations, with an emphasis on water and ion organization and mobility at mineral-water and water-air interfaces. However, the reliability of simulation predictions depends on the ability of the force field, i.e. the set of parameters and equations used to calculate interactions between atoms, to correctly describe the physics of these systems. In a previous study, [START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF] we showed that simulations run with a polarizable force field resulted in diffuse layer properties predictions that were significantly different from those calculated with non polarizable force fields. In particular, the use of a polarizable force field led to an increased order of water molecules and ion organization at the vicinity of clay mineral surfaces compared to results obtained with non polarizable force fields. For this reason, we conducted our simulations with a polarizable force field.

Materials and Methods

Systems

All unsaturated pore simulation boxes were generated on the basis of simulations previously carried out in water saturated condition and in the presence of 0.1 mol•L -1 NaCl, for which full information is available in our previous study. [START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF] Simulation boxes had all the same dimensions, i.e. 20.72 × 35.88 × 54.70 Å3 in the x, y and z directions respectively. Therefore, each system was composed of an identical 4×4 montmorillonite layer unit with 0.75 charges per unit cell, 12 sodium cations that compensated the layer charge in addition of two pairs of Na + and Cl -ions and n water molecules with n ranging from 72 (most unsaturated system) to 1122 (water saturated system) (Figure S1). Desaturation was achieved by removing water molecules from the center of the pore. Ions present initially in the removed aqueous phase were moved back into the aqueous phases so that there was the same number of Na + and Cl -ions in front of each surface bordering the pore. Because only water molecules were removed during the desaturation process, NaCl concentration in the aqueous phase increased proportionally with desaturation (Table 1). The three most unsaturated systems had 6, 12 and 18 water molecules per charge-compensating Na + cations, which was commensurate with the number of water molecules per cations in monohydrated, bihydrated and trihydrated sodium montmorillonite interlayer pores respectively. [START_REF] Tesson | Classical polarizable force field to study hydrated charged claysand zeolites[END_REF][START_REF] Ferrage | Classical polarizable force field to study dry charged clays and zeolites[END_REF] The number of gas molecules in our system can be evaluated to be less than one in the simulation boxes for 300 K and 1 bar conditions (ideal gas law). Consequently, no oxygen or nitrogen molecules were added in the gas phase, thus air phase was considered to be equivalent to vacuum in this study.

Table 1: Systems description with the number of water molecules per clay mineral layer unit, the total number of water molecules, the percentage of water molecules with respect to the saturated system and the concentration of NaCl in the aqueous phase.

N H 2 O / unit cell N H 2 O N H 2 O / N H 2 O , sat (%) c N aCl (mol•L -1 )

Molecular Dynamics Force Field

The explicit consideration of polarizability can drastically change the predictions of aqueous phase properties at mineral-water [START_REF] Tesson | Classical polarizable force field to study hydrated charged claysand zeolites[END_REF][START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF]47 and at air-water interfaces. 48-51 Consequently, we used the Polarizable Ion Model (PIM) force field. This force field was initially developed for the modeling of electrolyte solutions, 52 and was extended to various clay and zeolite minerals. [START_REF] Tesson | Classical polarizable force field to study hydrated charged claysand zeolites[END_REF]47,53,54 In PIM, each atom is treated as a polarizable charged ball. PIM was combined with the four-site Dang-Chang water model, [START_REF] Dang | Molecular dynamics study of water clusters, liquid, and liquid-vapor interface of water with many-body potentials[END_REF] which uses point dipoles to handle water molecule polarizability. Because Cl -were forced to approach the clay mineral surface with the increase of desaturation, it was necessary to add interaction parameters for Cl - and clay mineral layer atoms in PIM (Table S1). PIM parameters fitting procedure followed the method of Madden et al. [START_REF] Aguado | Ewald summation of electrostatic multipole interactions up to the quadrupolar level[END_REF][START_REF] Madden | From first-principles to material properties[END_REF] It consisted in an iterative process on several independent and characteristic configurations, with variation of the force field parameters (using Minuit code [START_REF] James | Minuit -a system for function minimization and analysis of the parameter errors and correlations[END_REF] ) to minimize differences between dipole and force values calculated by DFT and by classical MD force field. These differences were measured by χ 2 µ and χ 2 F for dipoles and forces error functions respectively (Eq S1 and Eq S2 in supporting information). The values of these error functions (0.109 and 0.572 for χ 2 µ and χ 2 F respectively) were of the same order of magnitude as those obtained during previous fits of the PIM force field parameters. [START_REF] Tesson | Classical polarizable force field to study hydrated charged claysand zeolites[END_REF]47,53,54 Simulation details

Simulations were run with the CP2K 59 package. Periodic boundary conditions were used in all directions. The cutoff radius for short-range electrostatic and repulsion-dispersion interactions was set to 10 Å. The long-range electrostatic interactions were computed using the Ewald method 60 with an accuracy of 10 -6 , an α parameter of 0.35 and 11, 19 and 29 k vectors in the x, y and z directions respectively. The SHAKE algorithm 61 was used to keep the water molecules rigid.

Unsaturated systems were first equilibrated for 3 ns in the N V T ensemble at 300 K using a Nosé-Hoover thermostat, gradually increasing the timestep from 0.1 fs to 0.5 fs and the thermostat constant from 10 fs to 1000 fs. The total simulation time of the production phases for the unsaturated pores was set to 40 ns in the N V T ensemble at 300 K with a thermostat constant of 1000 fs and with a timestep of 0.5 ns. We reduced the timestep from 1 fs for the saturated system to 0.5 fs for the unsaturated system in order to avoid numerical errors generated by species acceleration at water-air interface. Indeed, we observed that a timestep of 1 fs was responsible for water molecules to approach too closely from each other, thus causing the overall energy of the system to be unstable, and stopping the simulation.

Data post-processing

Radial distribution functions, density distributions profiles perpendicular to and density distributions maps parallel to the montmorillonite surface were averaged on the whole trajectories. Diffusion coefficients were calculated on subsets of 5 ns trajectories, a procedure which allowed to calculate standard errors. Diffusion coefficients parallel to the basal surfaces (xy plane), D || , were calculated as a function of the distance to the surface from Mean Square Displacements (MSD) and survival probability chosen as the uninterrupted residence in the layer (methodology available in Liu et al. [START_REF] Liu | On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid vapor interface of water[END_REF] and Marry et al. [START_REF] Marry | Structure and dynamics of water at a clay surface from Molecular Dynamics simulation[END_REF] ). The aqueous phase volume was gradually reduced with the increase in water desaturation. Thereby, D || of aqueous species were calculated only if the statistics, i.e. the number of the given species atoms in a slice parallel to the surface, was sufficient. For similar statistical reasons, the evaluation of the Cl -diffusion coefficient as a function of their distance from the surface was not carried out because of their small amount (2 atoms) in the systems.Because of the presence of periodic boundary conditions, hydrodynamic effects result in computation of diffusion coefficients that vary as a function of simulation box sizes for otherwise equivalent bulk systems. [START_REF] Tazi | Diffusion coefficient and shear viscosity of rigid water models[END_REF] Simonnin et al. [START_REF] Simonnin | Diffusion under confinement: hydrodynamic finite-size effects in simulation[END_REF] proposed an equation to correct the computed value of the diffusion coefficients for confined systems between two solid planes simulated in an orthorhombic box with a square base dimension L and with a pore width H. However, to our knowledge, no correction has been proposed in the case of an aqueous phase confined between a solid and a gas phase. In addition, in the saturated system, this correction would account for only 5 % of the computed value for water molecules. Therefore, no correction was made to the diffusion coefficients presented in this study.

Additional simulations with a non polarizable force field

We carried out simulations on unsaturated systems with the popular non polarizable force field ClayFF. [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF] It is not in the scope of this article to conduct a comprehensive comparison of force fields. However some of the results obtained with ClayFF reinforced the analyses made on the systems under investigation with PIM. Full details on ClayFF simulations are available in the SI.

Results and discussion

Interfacial water structuration Distribution profiles of oxygen (OW) and hydrogen (HW) atoms in water molecules was subdivided in four layers (or regions) in the direction normal to the clay mineral surface (Figure 1). In water saturated conditions, the three first layers exhibited features that were characteristic of water structuring (density peaks), and the fourth layer, in the center of the pore, was a bulk-like region in which OW and HW density profiles were almost flat at a value equal to bulk water density.

In desaturated systems, a water / gas interface appeared, with water covering the surface of the clay mineral and gas (vacuum) staying in the middle of the pore, in agreement with the observation that water remains at the surface and hydrates counterions of hydrophilic charged clay minerals. [START_REF] Marry | Structure and dynamics of water at a clay surface from Molecular Dynamics simulation[END_REF][START_REF] Churakov | Mobility of Na and Cs on Montmorillonite Surface under Partially Saturated Conditions[END_REF][START_REF] Simonnin | Mineral and ion specific effects at clay water interfaces: Structure, diffusion, and hydrodynamics[END_REF] As the system was desaturated, the interface moved closer to the clay mineral surface, and was characterized by a progressive decrease of OW and HW density from bulk density down to zero over a distance of about 5 Å, hereafter described as a transition zone. Distributions of water density in the water saturated part of the system were identical for all simulations down to 27 H 2 O per unit cell. OW distributions computed for these systems exhibit a peak doublet in the first layer that is indicative of two water populations with different orientations relative to the surface and to Na + cations. [START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF] Below this threshold value, further desaturation influenced water distribution profiles at water-clay interface. In the most unsaturated case, the decrease in intensity of OW first density peak indicated that the clay mineral surface was partially dehydrated (Figure 1a). A similar evolution of water structure at the interface as a function of desaturation was also observed in previous simulations using a non polarizable force field. [START_REF] Churakov | Mobility of Na and Cs on Montmorillonite Surface under Partially Saturated Conditions[END_REF] In less extreme conditions, the density peak doublet 5.75 and 6.6 Å merged into a wide single peak (Figure 1a). A peak doublet reappeared in the most unsaturated case, with maxima located at distances closer to the surface (5.65 Åand 6.35 Å) than in water saturated conditions. Concomitantly, the two HW density peaks closest to the surface drifted away about 0.1 Å from the surface as the surface was dehydrated, suggesting a slight change in water structure at clay mineral surface.

This change in water structure was confirmed by the proportion of OW located above the middle of the ditrigonal cavity that was higher in the most unsaturated case (Figure S2). In water saturated conditions, chloride ions were repelled from the clay mineral surfaces and located in the middle of the pore in agreement with the expected repulsive electrostatic interactions at negatively charged surfaces. On the opposite, sodium ions were attracted to the surface, and formed concentration peaks (Figure 2), the intensity of which decreased as the distance from the surface increased. [START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF] In the presence of unsaturated conditions, the water/gas interface pushed ions towards the surface. As a result, Cl -concentration distribution became more and more compact, and its center of mass shifted towards the surface with increasing desaturation levels. In the most dehydrated state, Cl -concentration distribution had the form of a sharp peak at about 7 Å from the surface. Anions remained closer to the water / gas interface than cations. Two hypothesis can be formulated to explain this observation: 1) anions are repelled from the negatively charged clay mineral surface by electrostatic forces; and/or 2) while small hard ions like Na + are repelled by the water / air interface because of the repulsive interaction with their image charge, [START_REF] Onsager | The Surface Tension of Debye-Hückel Electrolytes[END_REF] more polarizable ions like Cl -experience an attractive force towards the water / air interface due to van der Waals interactions 48 (see snapshot Figure S3). Na + ions were pushed towards the clay mineral surface with increasing desaturation level, leading to the growth of the Na + concentration peaks that was the closest to the clay mineral surface, but without much consequences on the position of Na + concentration peaks for water saturation conditions above 13.5 H 2 O per unit cell. At higher levels of desaturation, Na + ions came closer to the surface and screened the electrostatic repulsion between Cl -and the clay surface (Figure 2). New peaks arose, with cations standing between the negatively charged surface and the anions. In these new concentration peaks, cations were partially dehydrated, and they were coordinated directly to clay mineral surface atoms (see Table S2 and Table S3 for Na + first shell coordination composition), which is an indication of the formation of inner sphere complexes. These inner sphere cations were also observed in the absence of anions. [START_REF] Churakov | Mobility of Na and Cs on Montmorillonite Surface under Partially Saturated Conditions[END_REF] However, the presence of chloride anions increased the proportion of inner-sphere cations, and favored their appearance at lower desaturation levels. In the most unsaturated case, the peak the closest to the surface corresponded to cations located above the ditrigonal cavities and the second peak corresponded to cations above surface oxygen atoms (Figures 3 a andb). These two preferred positions had already been observed in our previous work carried out in water saturated conditions with the ClayFF force field. [START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF] Because the Na + -surface distance was the same, we infer that the Na + preferred positions with respect to the surface does not depend on the force field, despite the deformation of the surface cavities observed with PIM. [START_REF] Tesson | Classical polarizable force field to study hydrated charged claysand zeolites[END_REF][START_REF] Crom | Influence of Polarizability on the Prediction of the Electrical Double Layer Structure in a Clay Mesopore: A Molecular Dynamics Study[END_REF]53,54 Both cations and anions density peaks in region 1 (Figure 2) were located at a same distance from the surface (≈ 7 Å). Ions present in this region exhibited preferential locations on the surface (Figures 3 c andd): Cl -ions were located above Si atoms, while Na + ions were located near surface oxygen atoms of the ditrigonal cavity, leaving the center of the cavities to water molecules and Na + inner sphere complexes. Radial distribution functions g Cl-N a and g Cl-HW (Figure S4 and Figure S5) evidenced a decrease in the number of water molecules in Cl -hydration sphere, while the number of nearest Na + neighbors (located at a distance of 2.8 Å from Cl -) increased to reach an average coordination number greater than 1 in the most unsaturated system (Table S4 and Table S5). Consequently, Cl -ions formed ion pairs with Na + ions. In order to better visualize the relationships between the different species in very unsaturated systems (with 4.5 H 2 O per unit cell system or 13.5 H 2 O per unit cell), distribution profiles were plotted for all species on the same graph (Figure S6).

Ions distributions at interface

Above mentioned changes in ions distributions have a direct influence on the electrostatic properties within the pore. Surface charge compensation (CC) was calculated as follows:

CC = 100 L x L y Q z 0 (C Na + (u) -C Cl -(u)) du (1)
where L x and L y are the dimensions of the box in x and y directions, Q is the charge of clay mineral surface in elementary charge unit (which corresponds to the total charge of the clay mineral sheet divided by two), C Na + (z) and C Cl -(z) are the concentrations of cations and anions as a function of z in number of particles per unit volume. CC is zero at the clay mineral surface and 100 % in the middle of the pore because the system is overall electrically neutral and symmetric.

In the simulations carried out with the three most desaturated conditions, an overcompensation of the surface charge took place in the region z = 7 Å -10 Å (Figure 2c). In this region, the surface and the slice of fluid closest to the surface formed a positively charged wall, a feature that cannot be captured by simple continuum theories such as the Gouy-Chapman model. Continuum-scale modeling of unsaturated clay mineral systems requires thus the consideration of more elaborated models. Note that this phenomenon of overcom-pensation was also observed in simulations run with the non-polarizable force field ClayFF, for which overcompensation was more pronounced with increasing NaCl salt concentrations (Figure S7 and Table S6). with first steps of increasing desaturation (Figure 4 and Table S7). As desaturation proceeded further, water diffusion decreased down to 1.67 × 10 -9 m 2 • s -1 in the most unsaturated state.

This two-step evolution was driven by the relative influences of two competitive effects: (1) water molecules at the water/air interface diffused more quickly than in the bulk; 66 the proportion of water molecules located at the interface in the 4 th layer increased (yellow points on Figure 4) with desaturation, and thus the average diffusion coefficient increased too; on the contrary (2) the interactions of water molecules with cations and with the clay mineral surface slowed water diffusion down; in the 1 st water layer closest to the surface, water molecules had a diffusion coefficient of 3.5 × 10 -10 m 2 • s -1 , except for the most unsaturated states for which the diffusion coefficient in the first layer increased to 1.5 × 10 -9 m 2 • s -1 because of the presence of the water/air interface located in the same layer in these conditions. This two-step evolution of pore-averaged water diffusion coefficients as a function of desaturation had already been observed by Churakov with a non polarizable force-field. [START_REF] Churakov | Mobility of Na and Cs on Montmorillonite Surface under Partially Saturated Conditions[END_REF] Ionic species mobility ) whole pore 0 inner-sphere 1 outer-sphere 2 outer-sphere 3 outer-sphere pore center An increase in cation diffusion coefficient was observed near the water/air interface, similarly to observations made for water molecules, albeit less significant (Figure 5, Table S8). The proportion of cations near the water/air interface was much lower than that close to the clay mineral surface. As a result, the slowdown near the clay mineral surface counterbalanced more efficiently the acceleration near the water/air interface, and the average diffusion coefficient decreased with increasing desaturation levels. In the most unsaturated system, the diffusion coefficient was equal to 2.2 × 10 -10 m 2 • s -1 , i.e. about seven times lower that in bulk water (D N a + = 1.5 × 10 -9 m 2 • s -1 ) [START_REF] Tazi | Diffusion coefficient and shear viscosity of rigid water models[END_REF] and about five times lower than in saturated conditions. In fact, in the conditions corresponding to 4.5 and 9 H 2 O per unit cell, cations located in region 0 (inner sphere cations) were immobile, and contributed strongly to the observed slowdown of cations diffusion.

The evolution of Cl -diffusion coefficient with saturation was similar to that observed for cations (Figure 6, Table S9), but it exhibited a larger amplitude of variation. Indeed, average Cl -diffusion coefficient value decreased by a factor of about eight compared to the bulk water value (D Cl -= 1.9 × 10 -9 m 2 • s -1 ) [START_REF] Tazi | Diffusion coefficient and shear viscosity of rigid water models[END_REF] and six compared to the water saturated conditions. This evolution may seem surprising since anions, unlike cations, are a priori not strongly attracted to negatively charged surfaces. Nevertheless, formation of ion pairs in such spatially constrained conditions may be ultimately responsible for the loss of Cl - mobility. Diffusion experiments have evidenced that anions diffusion is more strongly reduced than cation diffusion in the presence of water unsaturated conditions. Numerous processes were put forward to explain this observation, including differences in diffusion pathway disruptions among cations anions and neutral molecules. [START_REF] Savoye | New experimental approach for studying diffusion through an intact and unsaturated medium: a case study with Callovo-Oxfordian argillite[END_REF][START_REF] Savoye | Mobility of cesium through the callovo-oxfordian claystones under partially saturated conditions[END_REF][START_REF] Savoye | Experimental study on diffusion of tritiated water and anions under variable water-saturation and clay mineral content: comparison with the Callovo-Oxfordian claystones[END_REF][START_REF] Savoye | Effect of Water Saturation on the Diffusion/Adsorption of 22 Na and Cesium onto the Callovo-Oxfordian Claystones[END_REF][START_REF] Savoye | How mobile is tritiated water through unsaturated cement based materials? New insights from two complementary approaches[END_REF] Our study emphasizes that mechanisms taking place at the surface, such as formation and immobilization of ion pairs, may also be, partly, responsible for the observed differences between anions and cations diffusion properties evolution with desaturation. As for water diffusion, our study corroborates the fact that a limited desaturation is not responsible for a decrease of water molecules selfdiffusion coefficient. On the contrary, a limited desaturation increases its value, an effect that is antagonistic to the effect of diffusion pathway disruption on the macroscopic effective diffusion coefficient value. The relative influence of these two competing effects was probed by Gimmi et al., [START_REF] Gimmi | Water retention and diffusion in unsaturated clays: Connecting atomistic and pore scale simulations[END_REF] who concluded that the decline in effective water diffusion due to the increase in tortuosity brought about by the addition of gas in the clay matrix is not compensated by the the increase of water mobility at the vapor-liquid interface. / unit cell and salt concentration).

Conclusion

• Clay charge compensation distribution profiles in the direction normal to the clay mineral surface as function of pore saturation for systems simulated with the non polarizable ClayFF force field.

• D || of water as a function of pore saturation.

• D || of sodium ions as a function of pore saturation.

• D || of chloride ions as a function of pore saturation.

Figure 1 :

 1 Figure 1: Density distribution profiles normal to clay mineral surface as a function of pore saturation: (a) water oxygen atoms, and (b) water hydrogen atoms. The clay mineral oxygen distributions are shown in brown, and the center of mass of the surface oxygen is represented by a brown dotted line at z = 3.27 Å. The surface water layers are delimited by the dashed black lines. Insets in the two panels focus on the layers located between 2.5 Å and 8 Å and help to visualize the distributions of atoms densities next to the clay mineral surface.

Figure 2 :

 2 Figure 2: Density distribution profiles normal to clay mineral surface as a function of pore saturation: (a) Na + , (b) Cl -, and (c) surface charge compensation in the EDL. Clay mineral layer oxygen distributions are shown in brown, and the center of mass of surface oxygens is represented by a brown dotted line at z = 3.27 Å. The different peaks observed on Na + distributions are delimited by dashed black lines. Insets in panels (a) and (b) focus on the region located between 2.5 Å and 10 Å and help to visualize the distributions of ions next to the clay mineral surface. Inset in panel (c) focuses on the region located between 6 and 12 Å and helps to visualize the charge overcompensation for the three most unsaturated pores.

Figure 3 :Figure 4 :

 34 Figure 3: Atomic density maps averaged over a unit cell of the clay mineral surface of (a) Na + in the first peak of inner-sphere cations (region 0 inner-sphere ), (b) Na + in the second peak of inner-sphere cations (region 0 inner-sphere ), (c) Na + in the first peak of outer-sphere cations (region 1 outer-sphere ) and (d) Cl -in the region 1 outer-sphere . Si (yellow), O (red) and H (white) clay mineral atoms were drawn to show the position of the ditrigonal surface cavities.

Figure 5 :

 5 Figure 5: D || of sodium ions as a function of pore saturation. The definitions of the different regions are given on the ions distribution profiles (Figure 2). The coefficient values and associated errors are averaged values of two symmetric layers.

Figure 6 :

 6 Figure 6: D || of chloride ions as function of the pore saturation.

•

  Desaturation drastically modified the localization and the diffusion of mobile species (water and ions) in most unsaturated conditions, corresponding to a monolayer (4.5 H 2 O per unit cell), a bilayer (9 H 2 O per unit cell) or a trilayer of water (13.5 H 2 O per unit cell). In particular, cations engaged in immobile inner sphere complexes located in ditrigonal cavities, resulting in a local overcompensation of the surface charge. In parallel, while water diffusion coefficient remained relatively unaffected by the extent of desaturation, ions mobility was greatly reduced in the most unsaturated simulation cases. The complexity of mutual interactions between water anion-cation-clay mineral surface and water-air interface pleads for the development of advanced continuum models to take into account desaturation effect on transport properties in clay-rich media. Description of the systems simulated with the non polarizable ClayFF force field (N H 2 O

•

  Mean square displacement (MSD) of water molecules as a function of time and as a function of distance from clay mineral surface (system: pore with 36 H 2 O / unit cell). Molecular Dynamics modeling study. Geochimica et Cosmochimica Acta 2006, 70, 562-582. (31) Greathouse, J. A.; Cygan, R. T. Molecular Dynamics simulation of uranyl(vi) adsorption equilibria onto an external montmorillonite surface. Physical Chemistry Chemical Physics 2005, 7, 3580-3586. (32) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J.; Hou, X. Molecular modeling of the structure and energetics of hydrotalcite hydration. Chemistry of Materials 2001, 13, 145-150. (33) Kalinichev, A. G.; Kirkpatrick, R. J. Molecular dynamics modeling of chloride binding to the surfaces of calcium hydroxide, hydrated calcium aluminate, and calcium silicate phases. Chemistry of Materials 2002, 14, 3539-3549. (34) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J. Molecular modeling of water structure in nano-pores between brucite (001) surfaces. Geochimica et Cosmochimica Acta 2004, 68, 3351-3365. (35) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J.; Cygan, R. T. Structure, energetics, and dynamics of water adsorbed on the muscovite (001) surface: a Molecular Dynamics simulation. The Journal of Physical Chemistry B 2005, 109, 15893-15905. (36) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J. Asymmetric hydrogen bonding and orientational ordering of water at hydrophobic and hydrophilic surfaces: A comparison of water/vapor, water/talc, and water/mica interfaces. The Journal of Physical Chemistry C 2009, 113, 11077-11085. (37) Kirkpatrick, R. J.; Kalinichev, A. G.; Hou, X.; Struble, L. Experimental and Molecular Dynamics modeling studies of interlayer swelling: water incorporation in kanemite and ASR gel. Materials and Structures 2005, 38, 449-458.

Acknowledgement

This work is part of the Ph.D. project of S.L.C., which was funded by the French National Agency for Radioactive Waste Management (Andra). This study benefited from HPC ressources granted by GENCI (TGCC and Cines, Grant No. A0060410107) and Sorbonne Université (MeSu). V.M. thanks the CNRS project NEEDS-DARIUS. Additional fundings (C.T.) from the EC Horizon 2020 project EURAD under Grant Agreement 847593 (WP FUTURE), and from a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" Programme LabEx VOLTAIRE, 10-LABX-0100 are acknowledged.

Supporting Information Available

• Snapshots of the studied systems (saturated and unsaturated).

• Minimization functions (forces and dipoles) used in the PIM force field fitting procedure.

• Cl --Clay mineral atoms interaction parameters in PIM force field.

• Water oxygen atoms (OW) density maps in the first layer for systems 45 H 2 O / unit cell and 4.5 H 2 O / unit cell.

• Snapshot of half a pore (36 H 2 O / unit cell) with focus on ionic distributions.

• Position and number of OW in the first peak of the radial distribution functions of OW water atoms around Na + as function of pore saturation.

• Na + first shell coordination composition.

• Radial distribution functions and integrals of Na + around Cl -as function of pore saturation.

• Position and number of Na + in the first peak of the radial distribution functions of Na + around Cl -as function of pore saturation.

• Radial distribution functions and integrals of HW around Cl -as function of pore saturation.

• Position and number of HW in the first peak of the radial distribution functions of HW around Cl -as function of pore saturation.

• Non polarizable force field simulations: simulation details description.

• Water atoms and ions distribution profiles normal to clay mineral for the 4.