Virginie Marry

Benjamin Rotenberg

Pierre Turq

Structure and dynamics of water at a clay surface from molecular dynamics simulation

Keywords: Water, clays, liquid-solid interface, hydrogen-bond, kinetics, molecular dynamics

We report a molecular dynamics study of the structure and dynamics of water at a clay surface. The negative charge of the surface and the presence of surface oxygen atoms perturbs water over two to three molecular layers, while the nature of the counterions (Na + or Cs + ) has only little effect. In the first molecular layer approximately half of the water molecules are H-bonded to the surface. We also analyze the H-bond network between surface water molecules. The diffusion of water molecules along the surface is slowed down compared to the bulk case. As far as the orientational order and dynamics of the water dipole are concerned, only the component normal to the clay surface is perturbed. We investigate the surface H-bonds formation and dissociation dynamics and its coupling to the release of molecules from the first molecular layer. We introduce a simple kinetic model in the spirit of Luzar and Chandler [Nature, 379, 55 (1996)] to allow for a comparison with bulk water dynamics. This model semi-quantitatively reproduces the molecular simulation results and suggests that H-bond formation is faster with the surface than in the bulk, while H-bond dissociation is slower.

INTRODUCTION

Clays are porous materials found in large quantities at the Earth's surface. Their mechanical (swelling), catalytic or retention (cationic exchange, anionic exclusion) properties explain their implication in a number of environmental and industrial processes. Elucidating the transport mechanisms of water and ions through clays is particularly important, since clays are considered as suitable materials for the geological storage of toxic and radioactive waste [START_REF] Andra | Évaluation de la faisabilité du stockage géologique en formation argileuse[END_REF][START_REF] Bradbury | Near Field Sorption Data Bases for Compacted MX-80 Bentonite for Performance Assessment of a High-Level Radioactive Waste Repository in Opalinus Clay Host Rock[END_REF]. Water at the surface of clays determines the hydrodynamic properties of flows through the porosity. This porosity is either saturated by an electrolyte solution (e.g. in soils) or unsaturated (e.g. in the caprock above a reservoir of gas such as CH 4 or CO 2 ) or saturated by both a hydrophobic and a water phase (e.g. in an oil reservoir). These examples pinpoint the need for a comprehensive understanding of the structure and dynamics of water at a clay surface.

Due to isomorphic substitutions these layered aluminosilicates bear a negative charge compensated by counterions. Upon drying of a clay suspension individual clay layers assemble to form anisotropic particles. The number of layers per particle depends on the nature of the clay and the history of the sample. Typical lateral extensions of a few ∼ 100 nm are to be compared to a thickness of ∼ 10 nm corresponding to a few individual layers. The space between layers gives rise to a nanoporosity called interlayer porosity containing counterions and water. Further assembly of particles into aggregates results in larger pores: Micropores between particles within micrometric aggregates and macropores between aggregates.

Most experimental studies of water at clay surfaces were conducted on clay at low hydration in order to avoid the signal of bulk water molecules. They thus mostly probed interlayer water. Its structure and dynamics has been investigated using X-ray diffraction combined with water adsorption experiments [START_REF] Cases | [END_REF][4][5][6], Nuclear Magnetic Resonance [7], Quasi-Eleastic Neutron Scattering [8][9][START_REF] Poinsignon | Proc. Int. Clay Conf[END_REF][START_REF] Swenson | [END_REF][12][13][14] or Infrared spectroscopy [15,16]. These experimental results suggest that interlayer water is more strongly hydrogen bonded than in the bulk and diffuse more slowly [17]. Surface specific techniques such as Sum Frequency Generation (SFG) and Second Harmonic Generation (SHG) have already proved extremely powerful in the study of water/air [18,19], water/membrane [20], water/metal [21] or water/mineral interfaces [22,23]. SHG has been employed to investigate the sorption of organic molecules at clay particle edges [24], but not water at basal surfaces.

Computer simulations have provided microscopic insights into the interlayer properties of clays. They allowed for the interpretation of experimental data related to structural [6,[25][26][27][28][29][30][START_REF] Sposito | Proc. Nat. Acad. Sci[END_REF][START_REF] Park | [END_REF], thermodynamical [33][34][35][36][37][38][39][40][41][42] and dynamical [13,14,[43][44][45][46][47] properties of clay interlayers. The present study is devoted to the external surfaces corresponding to the interface between a solution saturating clay micropores and the clay particles. More precisely, a simplified picture of clay microporosities distinguishes two types of external surfaces. Basal surfaces are found along the clay layers and differ from interlayer surfaces by the absence of another layer facing them. Lateral surfaces are very different since they consist of broken edges of layers and allow for the exchange of water and ions between the micropore and the interlayer. We studied the lateral surfaces of Montmorillonite clay particles in a previous paper [48] and consider here only basal surfaces. Due to the anisotropy of Montmorillonite particles the latter form the major part of the external surfaces. Previous studies on basal surfaces of Pyrophillite, Laponite, Montmorillonite and Beidellite clays focussed on the properties of ions in their vicinity, in particular their possible sorption at the mineral surface [49][50][51]. The structure and diffusion of water at the surface of other layered minerals such as Brucite, Gibbsite, Hydrotalcite, Muscovite, Talc [52][53][54][55] or Tobermorite [56] have also been investigated.

The purpose of the present paper is to study with molecular simulations the influence of a clay surface on the structure and dynamics of water, with a particular attention to the H-bond network formed (a) between water molecules and the surface and (b) between surface water molecules. After presenting the simulated system and the simulation methods, we first describe the structural properties of water. We then explore the diffusion and the orientational dynamics, the kinetics of formation and dissociation of H-bonds with the surface, and finally the coupling of these events with the release of water molecules from the first molecular layer near the surface.

SYSTEM AND METHODS

System description

Each clay layer consists of one sheet of octahedral aluminum oxide between two sheets of tetrahedral silicon oxides. A fraction of Al 3+ and Si 4+ are substituted by cations of lower valence, Mg 2+ and Al 3+ respectively. The Montmorillonite studied in this article contains only octahedral substitutions. Its unit cell formula is : C 0.75 Si 8 (Al 3.25 Mg 0.75 )O 20 (OH) 4 . The negative charge is thus 0.75 e per unit cell. It is located in the middle of the clay layer and is compensated by 0.75 monovalent counterions C + = Na + or Cs + . The atomic structure of the clay layer was taken from X-ray diffraction measurements [START_REF] Brindley | Crystal Structures of Clay Minerals and their X-ray Identification[END_REF][START_REF] Maegdefrau | [END_REF]. The layer thickness is 6.54 Å. The horizontal dimensions of the simulation box are 35.88 × 41.44 Å2 corresponding to 4×8 unit cells. The simulation box contains one layer cut into two halves as illustrated in Fig. 1. Periodic boundary conditions in all dimensions are used. We thus simulate a periodic stacking of infinite clay layers, as was done e.g. by de Carvalho and Skipper [49]. Including an additional layer and a corresponding interlayer to account for the presence of the "bulk" clay particle [50] does not lead to significant changes in the distribution of water and ions near the surface.

The space between the layers is filled with water, counterions and salt. Figure 1 shows the case of a pore containing 1600 water molecules, 24 Na + counterions and 28 NaCl pairs as added salt. This corresponds to a concen-tration in the middle of the pore of ∼ 1 M. The box size in the direction normal to the surfaces is 42.9 Å (see below). The nature of the counterion has a crucial effect on the interlayer properties. In order to investigate the influence of the counterion on the structure and dynamics of water at basal surfaces, we also simulated a system with Cs + counterions and CsCl as added salt. We will show in the next sections that the distance between surfaces is large enough to ensure that both surfaces are independent of each other for the considered salt concentration.

FIG. 1: The simulation box contains one clay layer of 4 × 8 unit cells. The layer is divided in two halves and periodic boundary conditions in all direction are used. The negative charge of the mineral is compensated by 24 Na + counterions. The system also includes 1600 water molecules and 24 NaCl pairs. Green : Al and Mg; yellow : Si; red : O; white : H; blue : Na; pink : Cl; silver : Water molecules.

We use the SPC/E water model [59], known to accurately reproduce structural and dynamical properties of water, in particular its dielectric constant and diffusion coefficient [60]. It was also successfully used to describe the hydrogen bond dynamics in bulk water [61]. This point is particularly important for our study, since we want to assess the effect of clay surfaces on the H-bond structure and dynamics. In this model the molecule is rigid with an O-H distance of 1.0 Å and an H-O-H angle of 109.47 • . Clay layers are considered as rigid molecules.

The force field used to describe the interactions between clay atoms, water atoms and ions consists of pairwise interactions. The interaction between two atoms contains a direct electrostatic (Coulomb) contribution and a Lennard-Jones potential to account for short-range repulsion and long-range dispersion forces :

V ij = q i q j 4π 0 r ij + 4ε ij σ ij r ij 12 - σ ij r ij 6 (1)
where q i and q j are the charges of atoms i and j, and r ij the distance between them. The Lennard-Jones parameters σ ij and ε ij are computed from individual σ ii and ε ii using Lorentz-Berthelot mixing rules. We use the individual parameters for ions in SPC/E water of Dang et al. which were shown to correctly describe the dynamics of ions in bulk water [62]. The parameters for clay atoms are the ones of Smith [30]. Within this force field the Lennard-Jones parameters for surface oxygen atoms are identical to those of water oxygen atoms. The combination of all these parameters, summarized in table I, has already proved to describe reasonably well a number of experimental results on the thermodynamical, structural and dynamical [13,14,30,[44][45][46][47]

Monte-Carlo simulations and equilibrium distance

The equilibrium distance between surfaces for the above water and salt content is first determined by Monte-Carlo simulation in the isothermal-isobaric (N, σ zz , T ) ensemble [START_REF] Frenkel | Understanding Molecular Sim-ulations, From Algorithms to Applications[END_REF] with fixed number of water molecules N , stress normal to clay layers σ zz = 10 5 Pa and temperature T = 298 K. In order to ensure consistency with the periodic boundary conditions, rotation of the clay layers are excluded. Only translation of clay layers parallel to each other (leaving the intersurface distance unchanged) or in the direction perpendicular to the surface (changing this distance and hence the volume of the simulation box) are allowed. The equilibrium vertical dimension of the box is 42.90 ± 0.05 Å for the Na + counterion and NaCl case. The distance is then fixed to this value in subsequent molecular dynamics simulations. The same value was taken for Cs + with CsCl as added salt. Considering the difference in density between the two salt solutions, the error made in doing so is smaller than 1 Å.

Molecular Dynamics simulations

Molecular dynamics simulations are performed in the canonical (N, V, T ) ensemble using the DLPOLY simulation package [START_REF] Smith | Daresbury Laboratory[END_REF]. As starting configurations we use the outcome of the above described Monte-Carlo simulations. The equation of motion for the rigid SPC/E are integrated using the SHAKE algorithm [START_REF] Ryckaert | [END_REF] and we use a time step of 1 fs. For the determination of the density profiles averages are taken over 2.5 ns. Water dynamical properties are computed from 260 ps simulations sampled every 50 fs. Results are systematically compared to that in bulk solutions at the same 1 M salt concentration. For these simulations we use 5 NaCl or CsCl pairs in 256

H 2 O.
Given the geometry of the system we expect dynamical properties to vary with the distance from the clay surface. The pore is thus divided into layers, defined in the next section, in which each quantity of interest is evaluated. Here we explain how the diffusion properties are computed in each layer. The presence of the surface breaks the symmetry of the diffusion tensor which needs to be described by one component normal to the surface (D ⊥ ) and two equal components parallel to the surface (D ). The component parallel to the surface D is evaluated in each layer using the method of Liu et al. [66]. This method considers the mean-square displacement ∆y 2 (t) = [y(t) -y(0)] 2 along y (resp. x) of molecules remaining in the i-th layer between 0 and t as :

D (z i ) = lim t→∞ [∆x 2 (t) + ∆y 2 (t)] S i C (t) 4t S i C (t) , (2)
where z i corresponds to the center of the layer and S i C (t) = 1 if the molecules remains continuously in the ith layer between 0 and t, while S i C (t) = 0 otherwise. Thus S i C (0) is the fraction of molecules in the i-th layer and S i C (t) / S i C (0) is the survival probability of molecules in that layer. S i C (t) decreases faster as the width of the layer decreases. Hence the statistical accuracy on the meansquare displacement (MSD) can become rather poor after a few picoseconds if too thin layers are considered. However times required for the MSD in (2) to reach a linear regime can be of this order of magnitude. We found that a good trade-off could be obtained by choosing layers wider than 3 Å (see the next section) and performing a linear regression of the MSD between 8 and 16 ps. The linear behavior of the MSD is illustrated in Fig. 2.

A relatively simple evaluation of the normal component D ⊥ can also be obtained in regions where the density is constant from the autocorrelation of the functions [66]

Ψ i n (t) = sin nπ z(t) -z i min z i max -z i min , (3)
where z i min and z i max define the i-th layer of width L = z i max -z i min and n is an integer. When the density is constant the autocorrelation

g n (t) = Ψ i n (t)Ψ i n (0)S i C (t) / S i C (0) decays as exp(-n 2 π 2 L 2 D ⊥ t) so that D ⊥ is given by D ⊥ (z i ) = - L nπ 2 lim t→0 1 t ln Ψ i n (t)Ψ i n (0)S i C (t) S i C (0) . (4)
For pure diffusion, the decay is exponential at all times.

Because numerical data from simulations can be noisy as

FIG. 2: Mean-square displacement along the surface [∆x 2 (t) + ∆y 2 (t)] S i C (t) / S i C (t)
of water oxygen atoms in the surface (solid line) and central (dashed line) layers. The regression for the diffusion coefficient is done between 8 and 16 ps, where the MSD grows linearly with time. The inset shows that at long times (t > 60 ps) the MSD for the surface layer deviates from a linear regime. This is due to the faster decay of the survival probability S i C (t) / S i C (0) in that layer, leading to larger statistical inaccuracy at long times.

g n (t) becomes very small (i.e. larger n and longer times), we use the limit t → 0 to determine the slope. However this "short time" regime can be relatively long for n = 1. We will show in the next section that the water density is constant in the central layer, so that the above procedure can be applied. Figure 3 shows that ln g n (t) is indeed linear for n = 1 and 2 in the case of water in the central layer. In addition, the slope for n = 2 is four times larger than for n = 1. This allows (a) to conclude that the motion of water molecules in the direction perpendicular to the surface is diffusive in the central layer and (b) to determine the value of D ⊥ from the slope. The data for n = 2 are noisier than for n = 1 because the values of g 2 are much smaller than g 1 at long times.

The calculation of various other correlation functions in each layer is also based on survival probabilities in the layer. The details of their definitions will be given when needed in the following sections. Unless otherwise stated, the error estimates are computed using the block averaging method [START_REF] Frenkel | Understanding Molecular Sim-ulations, From Algorithms to Applications[END_REF].

STRUCTURAL PROPERTIES

Density profiles

The negative charge of the surface results in an excess of cations in its vicinity and a depletion layer for the anions. This well known effect is usually described in terms of double layer, with a fixed "Stern" layer of counterions at the surface partly screening the electrostatic interactions and a more diffuse layer containing both counterand coions. This screening is often described within the framework of Poisson-Boltzmann theory [67]. Molecular

(t) = Ψ i n (t)Ψ i n (0)S i C (t) / S i C (0)
in the case of water in the central layer with Na + ions. g1 and g2 decay exponentially and the decay rate for n = 2 is four times larger than for n = 1 (results are reported divided by n 2 ). Thus the motion of water in the direction perpendicular to the surface is diffusive in the central layer and the value of D ⊥ can be computed from (4).

simulations allow to go beyond this continuous solvent representation and account for the complex interplay between ions, solvent molecules and the surface.

The density profiles for cations, anions and water oxygen and hydrogen atoms as a function of the distance normal to the surface are displayed in Fig. 4. The vertical dashed-dotted lines indicate the location of the surface oxygen atoms. The double-layer picture is seen to be reasonably well suited, as suggested by the large excess of counterions near the surface and the depletion of coions. However there is a difference between Na + and Cs + ions: The maximum in the latter case is located closer to the surface than in the former. This is due to the fact that Cs + ions form mainly inner-sphere complexes (ISC) with the clay surface, in which surface oxygen atoms belong to the first coordination shell of the cation, whereas Na + ions mainly form outer-sphere complexes (OSC) in which only water oxygen atoms belong to the first coordination shell. This fact has been observed experimentally as well as in molecular simulation of clay interlayers [68]. Both types of complexes are illustrated by typical snapshots from our simulations in Fig. 4.

The distributions of water oxygen (O W ) and hydrogen atoms (H W ) are almost independent of the nature of the counterions, because most water molecules do not belong to the solvation shell of cations. Only the results for the Na + case are reported in the bottom part of Fig. 4. They reveal a strong ordering of water molecules near the surface over two to three water layers. In particular hydrogen atoms can come closer to the surface than oxygen atoms. We finally note that with the high ionic strength used for our simulations water and ionic profiles rapidly reach a constant bulk value. This confirms that we are indeed studying independent (i.e. non-interacting) clay surfaces. We now turn to the lateral distribution of O W and H W atoms for molecules located in the first molecular layer, defined by a distance between O W and the plane containing the surface oxygen atoms smaller than 4 Å. This cutoff distance corresponds to the location of the minimum between the two first oxygen density maxima. A shorter cut-off distance of 2.3 Å was used for hydrogen atoms to include only the closest atoms from the surface, since a second peak in the hydrogen density is found within the first molecular layer (see the bottom part of Fig. 4). Figure 5 reports the results for one clay unit cell, with darker regions corresponding to higher densities. Red and yellow circles indicate surface oxygen and silicon atoms, respectively. The latter are located below oxygens as can be seen in the snapshots of Figs. 4 and7. Hydrogen atoms have a clear propensity to reside over surface oxygens and avoid being above cavity centers and silicon atoms. Oxygen atoms exhibit a smoother density profile, with a preference for sites near (but not above) surface oxygen atoms and in the cavity center.

Water orientation and H-bond network

The previous observations on the position of water oxygen and hydrogen atoms, both vertically (H atoms can go closer to the surface than O) and laterally (H are located mainly above surface oxygen atoms) suggest that water molecules in the first layer form hydrogen bonds with the surface oxygen atoms. We now further investigate this point by studying the orientation of water molecules in the first molecular layer near the surface. To this end we first define a hydrogen bond with a surface oxygen O S by the same criterium as the one used in the study of bulk water [61,69,70]: The distance O S O W between the two oxygens must be smaller than 3.5 Å, and simultaneously the angle between the O W H W and O W O S vectors must be smaller than 30 • . Though simple this definition was shown to be a reasonable choice [71].

Table II summarizes the proportions of surface water molecules (i.e. with an O atom within 4 Å of the surface) as a function of the number of H-bonds with the surface (O S ), with water molecules in the surface layer (O 1 W ) and water molecules in the second layer (O 2 W ). Approximately 50% of surface molecules form one H-bond with surface oxygen atoms (type A molecules). Almost 45% donate one H-bond to at least one surface molecule without forming a H-bond with the surface (type B). Type B molecules can be further divided into molecules that are H-bond donors to two surface molecules (14% of all surface molecules, type B 2 ), to only one surface molecule (12%, type B 1 ) and to one surface molecule and one water molecule in the second water layer (11%, type B 1 ). Other possible configurations each contribute to less than 2% of all surface molecules.

The organization of surface molecules can be described in terms of their orientation with respect to the surface, defined by two angles illustrated in Fig. 6. The first one θ ∈ [0, π] is the angle between the dipole µ of the water molecule and the normal to the surface (z axis). The average orientation in the first layer is characterized by cos θ = µ z /µ = -0.40 ± 0.01 with Na + counterions and -0.325 ± 0.005 with Cs + counterions. The negative sign indicates that the dipole of water molecules are on average pointing toward the surface. This is consistent with the fact that H atoms can approach closer than O atoms. Note however that θ is smaller than the maximum angle θ = π that would be obtained for a dipole normal to the surface with two H atoms closer to the surface than the water O atom. In a clay interlayer with Na + and two molecular layers (i.e. one per clay surface) the average orientation is cos θ = -0.345 ± 0.009. The other relevant angle is β ∈ [0, π/2] between the normal n to the molecular plane and the normal to the surface. The distribution of both angles are reported in Fig. 6 for molecules H-bonded to the surface (type A) and in Fig. 7 for other surface molecules (type B).

Figure 6 clearly shows that the molecular plane of surface H-bonded molecules is almost perpendicular to the surface, since the maximum probability is found around β ∼ 90 • . Nevertheless, thermal fluctuations around the vertical position without destroying the H-bond are possible. The most probable value of cos θ corresponds to θ ∼ 140 • , consistent with the fact that the water oxygen atom is not located above the H-bonding hydrogen, but slightly shifted (see Fig. 5). For the SPC/E water model which has an angle between OH bonds of 109.5 • , a vertical OH bond would indeed result in a maximum around 180 -109.5/2 ≈ 125.25 • . Typical configurations for molecules of type A (H-bonded to the surface) can be seen in the snapshots of Fig. 7.

The angular distribution is more complex in the non surface H-bonded case (type B, Fig. 7). For these molecules the maximum probability is obtained along the circle cos 2 θ + cos 2 β = 1 corresponding to situations where the plane containing both µ and n is orthogonal to the clay surface (in that case θ = 90 • ± β). Along this circle the most probable orientations are found between θ ∼ 105 and 140 • . All molecules H-bonded to two surface molecules (type B 2 ) are oriented as such, as can be seen in the bottom part of Fig. 7. This configuration is similar to H-bonds in bulk water as illustrated on the bottom snapshot of the same figure. Surface molecules of type B 1 , forming only one H-bond with surface molecules, can rotate about their donating OH bond. This results in a much more diffuse distribution centered around an ellipse (dotted line in the top part of Fig. 7). Let us conclude this section on structural properties by remarking that the preferential orientation of the dipole toward the surface is not surprising: The presence of surface oxygen atoms leads to hydrogen bonds with the surface, and the overall negative charge of the clay generates an electric field oriented toward the surface along which the dipole tends to align. In turn, this ordering of the first molecular layer of water affects the structure of the second layer by favoring H-bond donation from water molecules in the second layer to oxygen of surface molecules. Nevertheless the structure of this second layer is already closer to the bulk one.

DYNAMICAL PROPERTIES Diffusion

We start our study of the dynamics of water near a clay surface by investigating its diffusion properties. In order to evaluate how the components D ⊥ and D evolve with the distance to the surface we first need to define layers across the pore from the density profiles. The surface layers have already been defined in the previous section as corresponding to water oxygen atoms less than 4 Å apart from the plane containing the surface oxygen atoms. The cut-off corresponds to the location of the density minimum between the two first maxima. The second layer is comprised between 4 and 7.5 Å from the surface, i.e. between the above minimum of the density profile and the next one. A third layer can be defined by including the water molecules between 7.5 and 13 Å from the surface, corresponding to the part of the density profile that still hasn't reached the bulk value. Finally, the central layer corresponds to all other molecules, at distances greater than 13 Å from both surfaces, where the density profile has reached the constant value of the bulk liquid. These four regions are delimited by vertical dotted lines in Fig. 8.

D and D ⊥ are then evaluated in each layer as explained in the simulation methods section. Results for D are reported in the top part of Fig. 8 for both Na + and Cs + counterions. In both cases diffusion is slowed down near the surface compared to bulk water. The diffusion coefficient value in the central layer is comparable to the bulk result for the SPC/E model (2.4±0.1 10 -9 m 2 s -1 ) [72] and the experimental one (2.3±0.1 10 -9 m 2 s -1 ) [START_REF] Krynicki | Faraday Discuss[END_REF]. The slight deviation comes from the presence of salt at a relatively high concentration.

A striking difference between both types of counterions is observed in the second layer where D is smaller in the Na + case, even smaller than the value in the first layer for the same ion. Such a difference probably arises from the different hydration properties of the cations. The cations have a similar diffusion coefficient in the first layer (0.9±0.1 10 -9 m 2 s -1 ). Thus both tend to slow down the molecules belonging to their hydration shell (partly located in the second water layer). However Cs + ions form mostly inner-sphere complexes with the surface and therefore perturb less the water molecules in the second water layer than Na + ions. The latter are mainly located between the first and second water layers and form outer-sphere complexes with the surface. In addition, it is known that water molecules in the first hydration shell of Na + ions are more strongly bound than in the Cs + case and have a longer residence time. This time has been computed by Koneshan [62] with the same water and ions model that we have used: It is on the order of 20 ps for Na + , approximately twice longer than for Cs + .

Comparison of the results obtained for D ⊥ in the central layer (the only region where the density is constant) with different eigenfunctions Ψ n (n = 1, 2) confirms that Water diffuses more slowly on the surface than in the bulk. This slowing down is particularly pronounced in the Na + case for molecules belonging to the second layer.

the motion is diffusive in the central regions (see Fig. 3). We found values of 2.05 ± 0.05 10 -9 m 2 s -1 with Na + counterions and 2.20 ± 0.05 10 -9 m 2 s -1 with Cs + . These values are slightly smaller than the bulk value and than the D values. This indicates that although the bulk density has been reached the diffusion normal to the surface is to some extent still perturbed by the confinement by the surfaces. Using a larger simulation box size should allow to recover the bulk value in the central layer.

Water rotation

The previous section was devoted to the translational dynamics of water molecules. We now turn to the rotational dynamics. It can be quantified by the following correlation functions of the components of the water dipole µ:

C i ⊥ (t) = µ z (0)µ z (t) S i I (t) µ z (0) 2 S i I (t) (5a) C i (t) = [µ x (0)µ x (t) + µ y (0)µ y (t)] S i I (t) [µ x (0) 2 + µ y (0) 2 ] S i I (t) (5b)
where S i I (t) = 1 if the molecule is present in the i-th layer both at time 0 and time t and S i I (t) = 0 otherwise. At variance with S i C (t) introduced for the computation of the diffusion coefficient, such a definition allows for molecules to leave and come back into the i-th layer between 0 and t and characterizes an intermittent survival probability. As we explained earlier the continuous survival probability S i C (t) / S i C (0) decreases very fast and results in poor statistics for times longer than 20 ps. We thus consider S i I (t) to define the correlation functions ( 5), although the orientation of a molecule is influenced by its stay in neighboring layers. Results for Na + and Cs + are very similar and only the former are reported in Fig. 9, which also indicates the results for a bulk 1 M NaCl solution (in that case C ⊥ = C ) and the survival probability P i (t) = S i I (t) / S i I (0) in each layer. 5). The survival probabilities P i (t) = S i I (t) / S i I (0) are also indicated (dashed-dotted lines). The rotation of molecules about the normal to the surface (C ) is not hindered by the surface (C ∼ C b ), at variance with the rotation in the other direction (C ⊥ ). The plateau at times longer than 10 ps in the surface layer suggests that molecules oscillate around a particular orientation. In the central layer all components behave as in the bulk salt solution.

The orientational dynamics in the central layer (d) is the same than in a bulk salt solution. As water molecules get closer to the surface the normal component C ⊥ (t) decays more slowly, while the tangential component C (t) is hardly affected. The characteristic decay time of C ⊥ increases by a factor of two when moving from the bulk to the surface layer (a). Furthermore, C ⊥ does not decay to 0 but reaches a plateau value, at least on the 60 ps timescale. The limit values are C ⊥ (∞) ≈ 0.5 for Na + and 0.4 for Cs + , which compare well with the ratio of equilibrium averages µz 2 µ 2 z = 0.49 ± 0.06 for Na + and 0.36±0.02 for Cs + . This suggests that the dipole of water molecules with a long survival time in the surface layer oscillates around a particular orientation. Good candidates for such molecules are the one forming H-bonds with the surface (type A).

To further investigate this hypothesis we studied the distribution of dipole orientations as a function of first exit time from the surface layer. The normalized distributions are reported in Fig. 10. Although very broad for short exit times it narrows progressively toward a distribution similar to that of type A molecules, also indicated on this figure. These observations indicate that molecules that are not H-bonded to the surface leave the surface layer more rapidly than surface H-bonded molecules.

H-bonds formation and rupture

The fact that the surface H-bonded molecules have a longer residence time in the surface layer calls for further study of the H-bonds formation and dissociation kinetics. The last two sections complete the description of the surface water dynamics by adressing this issue. The H-bond kinetics in bulk water has been investigated in details by various approaches [61,69,70,74,75]. Some of them have also been used to investigate the water dynamics in the coordination shell of ions [76] or at the surface of micelles and proteins [75]. Apart from the one of Laage and Hynes [61,77] that counts switching events from a Hbond with a given molecule to one with another molecule, all approaches consider individual H-bonds. We thus define a H-bond characteristic function h(t) = 1 if a given water hydrogen atom is H-bonded to the surface at time t and 0 otherwise. In the spirit of previous works [74][75][76]78] we now introduce two time-correlation functions

A IS (t) = h(0)h(t)S I (t) h(0)S I (0) and (6a)
A CS (t) = h(0)h(t)S C (t) h(0)S C (0) , (6b)
where the two functions S I and S C have already been described and refer here to the surface layer only. They characterize the intermittent and continuous presence in the surface layer, respectively. S C allows to decouple the H-bond formation and dissociation reactions from the displacement into and out of the surface layer, whereas S I contains information on the exchange of molecules between the surface layer and the rest of the system. The notation A is associated with molecules forming an Hbond with the surface (see previous section). Both A IS and A CS give probabilities that the H is bonded to the surface at time t knowing that it was at time 0. They thus characterize the so-called intermittent lifetime of Hbonds, because they are allowed to break and reform between 0 and t. Other choices are possible: Counting only bonds that were never broken between 0 and t leads to a continuous H-bond lifetime [74,76,78,79]. Figure 11 reports A CS (t) and A IS (t) in the Na + counterion case. The probability that a molecule has never left the surface between 0 and t vanishes at long times. Hence A CS goes to zero at infinity. If the ratio between the volume of water and the area of the surface is large, the probability that a molecule comes back to the surface also vanishes at long times and A IS goes to zero at long times. Although at short times the two correlation functions coincide A CS decays faster than A IS . This is only because molecules leaving the surface between 0 and t contribute to the latter but not to the former. The decays of both function is not monoexponential, even in the simpler case of A CS . This is partly due to libration and inter-oxygen vibration at the sub-picosecond timescale, but also to more subtle effects which we now analyze along the lines introduced by Luzar and Chandler for bulk water [69]. 6) and (7). The inset illustrates the quantitative agreement between the simulation results for AIS and BIS and Eq. ( 8) with k1 = 0.14 ps -1 and k-1 = 0.06 ps -1 for times longer than 2 ps (shorter times are not shown).

We first note that once its H-bond with the surface is broken a water molecule is most likely still in the surface layer. This is obvious if the bond breaks after a reorientation event, which is the most probable dissoci-ation pathway [61]. However it is also holds if the dissociation is caused by a displacement of the molecule along the surface or even a moderate displacement away from the surface. This is because the surface layer corresponds to a larger distance (4 Å) than the hydrogen-bond (3.5 Å). Therefore a bond dissociation event is more likely to transform a type A molecule into a type B molecule than to lead to the direct removal of the molecule from the surface layer. A similar discussion had been introduced in the bulk water case by Luzar and Chandler who considered H-bonded and non H-bonded waters in the first shell around a given molecule [69]. Following their approach, we introduce the probability for a water H atom to be in the first layer (i.e. with the O atom of the same molecule less than 4 Å away from the surface) without forming a H-bond with the surface at time t, given that it was forming such a bond at time 0 :

B CS (t) = h(0)[1 -h(t)]S C (t) h(0)S C (0) , (7)
with the same function S C described above. A similar definition can be made using S I . Figure 11 also reports the evolution of B IS and B CS . Assuming that the rupture and formation of H-bonds with the surface can only lead to or come from molecules in the first layer we can write

dA dt = -k 1 A + k -1 B , (8)
where first order kinetics for both events has also been assumed. Since the H-bond formation and dissociation events do not depend on whether the molecules has left the surface layer in the past, the same rates k 1 and k -1 should be obtained if one considers A IS and B IS on the one hand, or A CS and B CS on the other hand. The inset of Fig. 11 shows that Eq. ( 8) holds for A IS and B IS with k 1 = 0.14 ± 0.01 ps -1 and k -1 = 0.06 ± 0.01 ps -1 , at least for times longer than 2 ps which are the only ones reported. The uncertainty on all rate constants are estimated visually by varying their values until the linear relationship defining them (between dA dt , A and B) is not satisfied anymore. The provided uncertainty is thus only indicative. The simple two-state model does not account for libration and inter-oxygen vibrations at very short-time scales. Thus deviation from Eq. ( 8) at short time is expected and was also observed in the bulk water case [69]. Equation ( 8) is also satisfied for A CS and B CS with k 1 = 0.15 ± 0.01 ps -1 and k -1 = 0.07 ± 0.01 ps -1 , very close to the above values, as illustrated in the inset of Fig 12. These results further support the hypothesis of first order kinetics for both surface H-bond formation and rupture. We finally note that the H-bond dynamics in the Cs + counterion case is very similar with k 1 = 0.14 ± 0.01 ps -1 and k -1 = 0.05 ± 0.01 ps -1 in the intermittent case (IS) and k 1 = 0.16 ± 0.01 ps -1 and k -1 = 0.06 ± 0.01 ps -1 in the continuous case (CS).

The fact that Eq. ( 8) is satisfied for both counting procedures (intermittent or continuous) and both counterions with approximately the same rate constants in all cases allows us to conclude that the surface H-bond formation and rupture can reasonably be described as first order kinetic processes. An additional confirmation of consistency is found with the equilibrium ratio between (a) surface H atoms (belonging to water molecules with an O less than 4 Å apart from the surface) forming an Hbond with surface O atoms and (b) surface H atoms not bonded to the surface. This ratio hS CS / (1 -h)S CS is 0.42 ± 0.01 and 0.39 ± 0.01 for Na + and Cs + , respectively. The function S SC is evaluated at time 0 only and ensures that we count only molecules in the surface layer. The same result is obtained with S IS , since both functions coincide at t = 0. This population ratio should be equal to k -1 /k 1 . Using the above values of the rates one finds for Na + k -1 /k 1 = 0.4 ± 0.1 and 0.5 ± 0.1 with the intermittent and continuous counting procedures, respectively. For Cs + one finds 0.4±0.1 (IS) and 0.4±0.1 (CS). These results compare well with the equilibrium population ratios and give further credit to the numerically determined rate constants.

In order to understand the influence of the surface on the H-bond dynamics we also performed a similar analysis for H-bonds between water molecules in a bulk 1 M NaCl solution. The bulk analog of the surface layer is the first coordination shell of the water molecule [69]. We chose the same inter-oxygen distance cut-off of 4 Å to define the B state. This is slightly larger than the 3.5 Å used by Luzar et al. which corresponds to the first coordination shell of a water molecule. However our choice is necessary for a meaningful comparison with the surface case. It also has the advantage of being not too close to the cut-off criteria defining the H-bond, thereby reducing the effect of inter-oxygen vibration on the rupture short-time kinetics. The bulk rates computed from A IS and B IS are k b 1 = 0.26 ± 0.01 ps -1 and k b -1 = 0.03 ± 0.01 ps -1 for H-bond dissociation and formation, respectively. For a bulk 1 M CsCl solution the rates are very similar: k b 1 = 0.28 ± 0.01 ps -1 and k b -1 = 0.04 ± 0.01 ps -1 . These results indicate that H-bonds with the surfaces form faster and dissociate more slowly than H-bonds with water molecules in the bulk. This constitutes one of the main conclusions of the present study. Both effects might be attributed to the average electric field pointing toward the negatively charged clay surface. Indeed, this field favors orientations of the water molecules with a dipole pointing in the direction of the surface, as was observed in the previous sections. Furthermore the orientation of the molecule fluctuates less under an electric field. Since reorientation is the main cause of H-bond formation and dissociation we should expect a faster formation of Hbond with the surface (favored by a dipole in the direction of the surface) and a slower dissociation (hindered by smaller fluctuations). However this hypothesis is difficult to test directly. A possibility would be to reproduce the same study for an uncharged clay such as pyrophyllite. Note that too high an electric field would align the dipole with the normal to the surface. In such a configuration no H-bond with the surface can be formed.

Another possible origin of the slowing down of H-bond dissociation at an interface is the decrease of the number of adjacent water molecules, as observed by Liu et al. at the water/air interface [80]. Breaking a bond requires to find a partner molecule to form a new bond [61]. Thus a smaller number of neighbors should indeed decrease the H-bond dissociation rate k 1 . The water/air interface is a particularly good illustration of the coupling between the formation and dissociation event on the one hand, and pair diffusion of the partner molecules on the other hand, a mechanism first evidence in bulk water by Luzar and Chandler [69,70]. Neglecting this distinction leads in that case to the (wrong) conclusion that the H-bonds have a shorter lifetime, because the diffusion at the interface is faster than in the bulk [80]. We now examine in more detail the release of water molecules from the surface layer and its coupling to the H-bond kinetics.

H-bonds kinetics and release from the surface layer Equation ( 8) is not sufficient to describe the evolution of A(t) and B(t). To obtain a closed set of equations we now consider the evolution of type B molecules. This is easier in the case where molecules leaving the surface layer are not counted anymore, i.e. for A CS and B CS . The non-bonded H atoms can either form a H-bond with the surface as already described or leave the surface layer. Assuming once more a first order kinetics with rate constant k 0 for this last process we have

dB CS dt = k 1 A CS -(k -1 + k 0 )B CS . (9)
When considering molecules reentering the surface layer (IS case) one should also include an additional source term. This would result in an open set of equations. The solution of Eqs. ( 8) and ( 9) with initial conditions A CS (0) = 1 and B CS (0) = 0 reads:

A CS (t) = s + -k 1 √ ∆ e -s-t - s --k 1 √ ∆ e -s+t B CS (t) = k 1 √ ∆ e -s-t -e -s+t , (10a)
where

s ± = k 1 + k -1 + k 0 ± √ ∆ 2 ∆ = (k 1 + k -1 + k 0 ) 2 -4k 1 k 0 ≥ 0 . (10b)
Using the values of k 1 and k -1 determined above solely from dACS dt we find that the simulation results are semiquantitatively reproduced with a rate k 0 ≈ 0.09 ps for the release from the surface layer. Simulation results are compared to the two-state model in Fig. 12. Considering the number of approximations underlying this simple kinetic model, such an agreement is very satisfactory.

At long times the solutions of Eqs. ( 8) and ( 9) satisfy lim

t→∞ A CS /B CS = ( √ ∆ + k 0 + k -1 -k 1 )/2k 1 .
In the case of slow water release from the surface layer compared to H-bond formation and dissociation (k 0 k -1 , k 1 ) this expression reduces to ratio k -1 /k 1 . This not the case with the values determined above which lead to a limit (A CS /B CS ) ∞ ≈ 0.72. The simulation results for the ratio A CS /B CS also reach a constant value of 0.70±0.05. This further supports the validity of this simple model and illustrates the robustness of the numerically determined rate constants.

FIG. 12: Evolution of the population of surface H-bonded molecules (A) and surface molecules not H-bonded to the surface (B) for Na + counterions. Only molecules that never leave the surface layer between 0 and t are considered in ACS(t) and BCS(t) defined by Eqs. ( 6) and (7). Simulation results (symbols) are compared to the predictions of the two-state model described in the text (lines) with rate values: k1 = 0.15 ps -1 , k-1 = 0.07 ps -1 and k0 = 0.09 ps -1 . The semi-logarithmic scale shows the exponential decay of ACS and BCS at long times. The inset illustrates the quantitative agreement between the simulation results and Eq. ( 8) for times longer than 2 ps (shorter times are not shown).

The numerical values of k 1 , k -1 and k 0 are all of the same order of magnitude. Therefore all processes are coupled. The benefit of distinguishing the intermittent and continuous countings is now more evident, since this allows to decouple the formation and rupture of H-bonds from the release of molecules from the surface layer. The characteristic release time 1/k 0 ≈ 11 ps is also the time where A IS and A CS start to diverge, as can be seen in Fig. 11. In addition, the hypothesis that only water molecules that are not H-bonded to the surface may leave the surface layer is validated by the agreement with the simulation results. This in turn explains the results of Fig. 10, namely that molecules staying longer in the surface layer have an orientational distribution similar to that of surface H-bonded molecules: Surface H-bonded molecules first need to break their H-bond before leaving the surface. Thus they tend to stay longer in the surface layer. The time it takes to reach the distribution of surface H-bonded molecules in Fig. 10 is also of the order of 1/k 0 ≈ 11 ps. This constitute an additional proof of consistency of this interpretation and of the numerical value of the release rate k 0 .

A final remark should be made about the kinetic model developed in this section. Although it reproduces semiquantitatively the molecular dynamics results, one should not forget that it is too simple. In particular we showed in the previous section that surface molecules that are not H-bonded to the surface (type B) could be H-bonded to zero, one or two surface molecules. This distinction is not included in the kinetic model. The good predictions of this model suggest that there is a fast equilibrium between all these forms compared to the surface H-bond formation and dissociation. The first order kinetics assumption for all processes is also a strong one. It is only valid for "not too short" times, as was mentioned earlier. With these caveats in mind, we still argue that this two-state model captures the dynamics of most relevant processes involving water molecules in the first layer near the clay surface. It also allows for a comparison of the surface and bulk behaviors, which was the main purpose of the present study.

CONCLUSION

We have studied the structure and dynamics of water molecules near a clay surface using molecular simulations. We showed that the both structure and dynamics are perturbed over two to three molecular layers only and that the nature of the counterion does not have a large effect on them. We demonstrated that in the first molecular layer approximately half of the water molecules are H-bonded to the surface and we analyzed the H-bond network between surface water molecules. The diffusion of water molecules along the surface is slowed down compared to the bulk case. As far as the orientational order and dynamics of the water dipole are concerned, only the component normal to the clay surface is perturbed. This effect can be understood by the presence of the electric field in the direction of the negatively charged surface and the presence of surface oxygen atoms leading to H-bonds between water molecules and the surface. Finally we investigated the surface H-bonds formation and dissociation dynamics and its coupling to the release of molecules from the first molecular layer. A simple kinetic model was introduced to allow for a comparison with bulk water dynamics. It semi-quantitatively reproduces the molecular simulation results and suggests that H-bond formation is faster with the surface than in the bulk, while H-bond dissociation is slower. Both effects can again be traced back to the electric field in the direction of the surface, which promotes orientations favorable to H-bonding. The slower H-bond dissociation can also be explained by a decrease in the number of adjacent water molecules to form a new bond [80].

Although obtained for a specific type of surface these conclusions are likely to apply to other types of mineral surfaces. Our observations are indeed related to the overall negative charge of the surface and to the presence of surface oxygen atoms, two features encountered with other common minerals. Results obtained here for Mont-morillonite clay should thus be transferable to other clay minerals. We expect the observed effects of the surface to be more pronounced as its charge density increases, e.g. with micas. In addition, the localization of the substitutions within the clay layers could play a role in the H-bond structure and dynamics. We expect that moving the charge defects from octahedral to tetrahedral sheets (as in Beidellite), i.e. closer to water molecules, will enhance the observed effect of the surface. It would be of particular interest to investigate positively charged surfaces such as Layered Double Hydroxides, also known as "cationic clays" [START_REF] Duan | Layered Double Hydroxides (Structure and Bonding)[END_REF]. In that case the average electric field and the hydrogen covered surfaces should both favor orientations where the water molecules accept H-bonds from surface hydroxyl groups. This would completely modify the structure and dynamics of the surface water molecules.

FIG. 3 :

 3 FIG. 3: Logarithm of the correlation function gn(t) = Ψ i n (t)Ψ i n (0)S i C (t) / S i C (0)in the case of water in the central layer with Na + ions. g1 and g2 decay exponentially and the decay rate for n = 2 is four times larger than for n = 1 (results are reported divided by n 2 ). Thus the motion of water in the direction perpendicular to the surface is diffusive in the central layer and the value of D ⊥ can be computed from (4).

FIG. 4 :

 4 FIG. 4: Distribution of ions (top) and water (bottom). The values for Cl -and OW have been scaled by a factor of two for clarity. The vertical dashed-dotted lines indicate the location of surface oxygen atoms. The water oxygen (OW) and hydrogen (HW) distributions are almost independent of the counterion and only results for Na + are shown. They indicate a strong ordering of water at the clay surface, over two to three molecular layers. Cs + cations form mainly inner-sphere complexes with the surface (top-right), while Na + form preferably outer-sphere complexes (bottom-right).

FIG. 5 :

 5 FIG. 5: Distribution of water oxygen (OW) and hydrogen (HW) atoms over the clay surface. Only molecules belonging to the first layer (distance smaller than 4 Å for OW, 2.3 Å for HW) are considered. Darker regions correspond to higher probabilities. Results are given for one unit cell with dimensions a × b = 5.18 × 8.97 Å2 . Oxygen surface atoms are indicated by red circles, silicon atoms (below oxygens, see Figs.4 and 7) by yellow circles.

  FIG. 5: Distribution of water oxygen (OW) and hydrogen (HW) atoms over the clay surface. Only molecules belonging to the first layer (distance smaller than 4 Å for OW, 2.3 Å for HW) are considered. Darker regions correspond to higher probabilities. Results are given for one unit cell with dimensions a × b = 5.18 × 8.97 Å2 . Oxygen surface atoms are indicated by red circles, silicon atoms (below oxygens, see Figs.4 and 7) by yellow circles.

  (i.e. with an O atom less than 4 Å away from the surface) as a function of the number of H-bonds donated to the surface (OS), to water molecules in the surface layer (O 1 W ) and to water molecules in the second layer (O 2 W ). Only the configurations contributing to more than 2% are shown. Results only slightly depend on the nature of the counterions (Na + or Cs + ).

FIG. 6 :

 6 FIG.6: Orientational distribution for surface H-bonded molecules (type A). The definition of both angles is illustrated on the right part of the figure. θ is the angle between the dipole moment µ of the molecule and the normal to the surface (z axis). β is the angle between the normal n to the molecular plane and the normal to the surface. The maximum probability corresponds to almost vertical OH bonds, with a water oxygen not exactly above the surface oxygen (see Fig5and the snapshots in Fig.7).

FIG. 8 :

 8 FIG. 8: Diffusion coefficient of water molecules along the surface (D ) as a function of the distance to the surface for Na + (•) and Cs + ( ) counterions. The vertical dotted lines indicate the boundaries between the layers (labelled a to d), and their relation to the water oxygen density profiles (bottom).Water diffuses more slowly on the surface than in the bulk. This slowing down is particularly pronounced in the Na + case for molecules belonging to the second layer.

FIG. 9 :

 9 FIG. 9: Correlation functions for the water dipole orientation in a bulk 1 M NaCl solution (C b , solid line) and its components along (dotted line) and normal (dashed lines) to the surface for water molecules in the first (a), second (b), third (c) and central (d) layers defined in Fig. 8. The dipole correlation functions are defined by Eqs. (5). The survival probabilities P i (t) = S i I (t) / S i I (0) are also indicated (dashed-dotted lines). The rotation of molecules about the normal to the surface (C ) is not hindered by the surface (C ∼ C b ), at variance with the rotation in the other direction (C ⊥ ). The plateau at times longer than 10 ps in the surface layer suggests that molecules oscillate around a particular orientation. In the central layer all components behave as in the bulk salt solution.

FIG. 10 :

 10 FIG. 10: Angular distribution of the dipole for surface water molecules as a function of first exit time from the layer. All molecules can leave the surface layer after a short time. Only the ones with an OH bond almost normal to the surface remain in this layer for longer times. The distribution of molecules leaving at long times progressively narrows toward the equilibrium distribution of surface H-bonded molecules (type A, dashed lines).

FIG. 11 :

 11 FIG. 11: Evolution of the population of surface H-bonded molecules (A) and surface molecules not H-bonded to the surface (B) for Na + counterions. The correlation functions AIS, ACS, BIS and BCS are defined by Eqs. (6) and(7). The inset illustrates the quantitative agreement between the simulation results for AIS and BIS and Eq. (8) with k1 = 0.14 ps -1 and k-1 = 0.06 ps -1 for times longer than 2 ps (shorter times are not shown).

TABLE I :

 I properties of bulk clays at low hydration (interlayer). The long-range electrostatic interactions are computed via three dimensional Ewald summation. Charges qi and atomic Lennard-Jones parameters σii and εii. The pair Lennard-Jones parameters are obtained from the Lorentz-Berthelot mixing rules.

	Molecule	Atom	qi (e) σii ( Å) εii (kJ.mol -1 )
		Al	3.0	0.0	0.0
		Mg (Subst. Al) 2.0	0.0	0.0
	Clay	Si O (Surface)	1.2 -1.0 3.166 1.84	13.18 0.650
		O (O h )	-1.424 3.166	0.650
		O (T d )	-0.8 3.166	0.650
		H	0.424	0.0	0.0
	Water	O H	-0.848 3.166 0.424 0.0	0.650 0.0
		Na	1.0	2.587	0.418
	Ions	Cs	1.0	3.883	0.418
		Cl	-1.0 4.401	0.418

TABLE II :

 II Distribution of water molecules in the surface layer

Acknowledgments

Fruitful discussions with Damien Laage (Ecole Normale Supérieure, Paris) and Daan Frenkel (AMOLF, Amsterdam) are gratefully acknowledged. B.R. acknowledges financial support from the Agence Nationale pour la Gestion des Déchets Radioactifs (ANDRA, France).