Alexandru Botan
email: alexandru.botan@upmc.fr

Benjamin Rotenberg
email: benjamin.rotenberg@upmc.fr

Virginie Marry

Pierre Turq

Benoit Noetinger

Alexandru Boţan

Benoît Noetinger

Carbon Dioxide in Montmorillonite Clay Hydrates: Thermodynamics, Structure, and Transport from Molecular Simulation

Keywords: carbon dioxide storage, clay minerals, swelling, diffusion, Monte Carlo, molecular dynamics

come

Introduction

Clay minerals are of great practical importance in soil science, groundwater hydrology, natural gas and petroleum reservoir engineering, storage of carbon dioxide or toxic chemical and radioactive waste. The role of a clay formation in gas or oil reservoirs is to trap the buoyant fluid in a lower porous formation (often water-filled carbonates). In this context, the most important properties are the low hydraulic permeability of the clay caprock and its ability to retain mobile species. Up to now, the main focus in the literature on CO 2 storage has been on the macrosopic two-phase flow in the host formation, assuming that the integrity of the caprock will be preserved so as to prevent the CO 2 plume from leaking toward the surface. The physics of this transport is very rich, as the introduction of CO 2 perturbs the equilibrium between carbonate rocks and dissolved carbonate ions, which may result in the rock dissolution, pore opening, carbonate reprecipitation, etc. 1,2 The interaction between the CO 2 plume and the clay caprock is also a crucial issue. However, studies on the microscopic scale are necessary to understand the complex interplay between pore size, fluid composition in the pores, and transport of mobile species. Indeed, interactions on the molecular level between water, CO 2 and the clay surface determine how the pore size and composition change when the buoyant bulk water/CO 2 mixture reaches the clay cap-rock. A scenario that one wants to avoid is the dehydration of clay particles induced by the presence of CO 2 , leading to their shrinkage. The resulting fractures would then increase the permeability of the caprock, thus reducing its ability to retain the fluid in the reservoir. 3 These properties depend strongly not only on temperature and pressure, mainly controlled by the burial depth, but also on the clay mineral composition, as was recently discussed by González Sánchez et al. 4 For example, the microscopic diffusion coefficient of water, as measured by Quasi-Elastic Neutron scattering, is similar to that of bulk water in compacted kaolinite and pyrophillite (uncharged clays), but considerably smaller in charged clays such as smectites and illites. While all clay minerals are layered aluminosilicates, their structure depends on their charge (neutral or negative), as the latter also implies the presence of compensating counterions between layers. The presence and nature of counterions is a key factor in determining how water and other molecules might penetrate into the interlayer space, a process generally accompanied by an increase in the interlayer distance and thus known as clay swelling. [5][6][7] Ab-initio simulations have shown that hydration of clay surfaces is a function of the structure of the crystal lattice. [8][9][10] Numerous classical molecular simulations and neutron scattering experiments have found strong dependence of swelling and hydration of clay minerals on the interlayer cation size and charge. [11][12][13][14][15][16][17][18][19][20][21] Such interlayer nanopores are absent in uncharged, non-swelling clays, and mobiles species such as water are present only in larger pores between clay particles (stacks of clay layers).

Recently, much attention has been paid to the study of the adsorption of different molecules and ions onto hydrated clay surfaces. These investigations are, in general, experimentally difficult and may be supplemented by using molecular simulation. Molecular simulations have been used

to study the atomic-scale interactions between ions and clay surfaces, e.g. in the case of uranium sorption 22,23 or cation exchange. 24,25 Adsorption is also possible in the case of methane [26][27][28][29] and carbon dioxide, 30 for which stable hydrate complexes were observed experimentally in the presence of smectites under pressures as low as 10 bar and at temperatures as high as 300 K. Titiloye and Skipper [31][32][33] used molecular simulation to investigate the structure and dynamics for given compositions of methane-water mixtures as interlayer species under pressure and temperature corresponding to reservoir conditions. They observed that an increase in pressure (in the reservoir) or in size of the cations causes a larger swelling and an increase in the self-diffusion coefficient for the methane molecules. Yang and Zhang also reported molecular dynamics simulations of the structure and diffusion coefficient of carbon dioxide in dry, uncharged clay-like slit pores. 34 Cole et al. investigated the influence of pressure on the properties of dry supercritical CO 2 in muscovite nanopores, using molecular simulations and small-angle neutron scattering. 35 However, the above-mentioned simulation studies were carried out at fixed number of molecules, pressure and temperature (NPT ensemble) and do not provide information on the amount and composition of the interlayer water-gas mixture, which vary with burial depth and pore size. 36,37 Since it is the chemical potential of the adsorbing molecules (water, methane, carbon dioxide) that is fixed by the reservoir, previous studies of the thermodynamics of clay swelling as a function of relative humidity have been performed using grand-canonical simulations, i.e. in the µ H 2 O V T ensemble. 14,38 The mixture composition in clays differs from that of the bulk mixture under the same conditions because the phase diagram is modified by the clay/fluid interactions. In order to address this issue, Odriozola et al performed simulations in the µ H 2 O µ C 2 H 6 PT ensemble to study the behaviour of ethane in hydrated montmorillonite interlayers. 39 This ensemble allows a direct measurement of the water and ethane content in the clay in equilibrium with given reservoir conditions. The ethane concentration observed in the simulations was very low and did not affect the swelling. Since the solubility of carbon dioxide in water is higher than that of ethane and even methane, 30,40 one could expect a priori a stronger effect of CO 2 on clay swelling.

In this paper we report a molecular simulation study of carbon dioxide in hydrated interlayers of Na-montmorillonite for two reservoir conditions, with the same temperature T =348 K but different pressures, namely P=25 and 125 bars. The equilibrium interlayer distance and composition are first obtained from Monte Carlo simulations in the µ H 2 O µ CO 2 V T ensemble, where the chemical potentials for water and carbon dioxide correspond to bulk phase coexistence in the reservoir. Molecular Dynamics simulations are then used to investigate the interlayer structure and dynamics.

Computational methods

General strategy

In order to simulate the behaviour of a clay caprock in contact with a CO 2 reservoir, we considered a Na-montmorillonite clay in equilibrium with H 2 O-CO 2 mixtures under two sets of P/T conditions close to relevant ones for geological storage, and under which two bulk phases coexist: H 2 O-rich liquid and CO 2 -rich gas at P=25 bar and T =348.15 K on the the one hand, H 2 O-rich liquid and CO 2 -rich supercritical fluid at P=125 bar and T =348.15 K, on the other hand. To assess the influence of CO 2 on the thermodynamics, structure and dynamics of hydrated clay we also considered Na-montmorillonite clay in equilibrium with pure water at P=125 bar and 348.15 K.

The grand canonical Monte Carlo method (µV T ensemble) provides an ideal computational tool for investigating the adsorption of molecules in porous materials. [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] In the present work, the clay interlayer exchanges both water and carbon dioxide molecules with reservoirs which set their chemical potentials. The latter are fixed by the coexistence between the CO 2 plume and the liquid H 2 O which fills the clay pores (liquid-vapour equilibrium). In the case of the clay phase the pressure has only normal component, therefore the external pressure is also considered as pressure applied normal to the clay surface P app . The equilibrium states of the system, characterized by the interlayer distance z and their composition, correspond to the minimum of the swelling free energy (F): 11

∆F = -L x L y ∫ z z 0 [P(z ′ ) -P app ]dz ′ , ( 1
)
where L x and L y are the dimensions of the simulation box along x and y axes, respectively. P(z ′ ) is the pressure as a function of the interlayer distance z ′ .

GCMC simulations require the knowledge of the chemical potentials µ H 2 O and µ H 2 O in the reservoirs. Numerous experimental and numerical works provides the data for pressure, temperature and composition of mixtures under reservoir conditions, but no data is available for the chemical potential. However, the chemical potential can be estimated from these data in µV T ensemble by running a series of simulations at fixed µ H 2 O and µ CO 2 and choosing the ones which give the right density and composition of the mixture for both coexisting bulk phases.

The equilibrium states are subsequently simulated using Molecular Dynamics to investigate the interlayer structure and dynamics.

Model and methods

The model of clay used in the calculations is the sodium-saturated Wyoming-type montmorillonite with unit formula [START_REF] Skipper | [END_REF]

U(r i j ) = U LJ +U C = 4ε i j [ ( σ i j r i j ) 12 - ( σ i j r i j ) 6 ] + q i q j 4πε 0 r i j ,
where r i j is the distance between sites i and j of different molecules, q i is the partial charge of the site, σ i j and ε i j are LJ parameters deduced from the conventional Lorentz-Berthelot combining rules. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] These parameters and charges can be found in references. [44][START_REF] Berendsen | Interaction models for water in relation to protein hydration[END_REF][START_REF] Harris | [END_REF]48 Monte Carlo simulation in the grand canonical ensemble (µ H 2 O µ CO 2 V T ) were performed using configurational bias techniques (CBMC). 14 The main idea is to divide the trial molecule insertion into two insertion steps: first of the central atom (O for water and C for carbon dioxide), then of other atoms (H or O). We generate k 1 trial configurations for the first atom and select one with since the k th one corresponds to the actual position which must be included in the set of trial configurations. The acceptance probability for the insertion step reads:

probability p 1 = exp(-β ∆u LJ 1 )/W 1 , where β = 1/(k B T ),
P N→N+1 acc = min ( 1, V Λ 3 (N + 1) W 1 k 1 W 2 k 2 exp[β (µ -∆U C )]
) , and for the removal:

P N→N-1 acc = min ( 1, Λ 3 N V k 1 W 1 k 2 W 2 exp[β (∆U C -µ)]
) .

Here Λ is the de Broglie wavelength, V the volume of the simulation box, N the number of molecules, µ the chemical potential, and ∆U C the change in electrostatic energy induced by the molecule insertion (removal). This algorithm was applied for H 2 O and CO 2 molecules with k

H 2 O 1 = k CO 2 1 = 50, k H 2 O 2 = 25 and k CO 2 2
= 50. The acceptance probability for the translation and rotation steps is based on the standard Metropolis criterion. 49 The simulations to estimate the chemical potentials from the density and composition of the mixture were performed for both coexisting bulk phases separately in a cubic simulation box. The and values for z were taken in the range 24 to 35 Å, corresponding to interlayer distances from 12 to 17.5 Å. The spherical cutoff radius is equal to the half of the smallest box side. Long-range electrostatic interactions were computed using Ewald summation. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] The system was equilibrated for 2 • 10 7 steps and 2 • 10 7 -5 • 10 7 additional steps were used for the calculation of the interlayer composition and pressure.

The pressure is calculated from the expression:

P = k B T < N tot > V - 1 L x L y ⟨( δU δ z )⟩ . ( 2
)
Here N tot is the total number of interlayer molecules, δU the change in configurational energy induced by an infinitesimal change δ z in the box size in the z direction (normal to clay layers). The pressure was computed every 100 steps with ghost volume changes δ z taken randomly within the range -0.01:0.01 Å. In the case of bulk phases one changes the volume along all directions. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] Molecular dynamics simulations were carried out in the canonical ensemble using the DLPOLY software package. [START_REF] Smith | The DLPOLY 2 user manual[END_REF] Initial configurations were taken from Monte Carlo simulations at the equilibrium distances,with the relative horizontal position of the clay surfaces fixed to the most probable one (as determined from the Monte-Carlo simulations). The equations of motion were integrated using a 1 fs time step and the SHAKE algorithm. [START_REF] Ryckaert | [END_REF] The temperature was fixed using a Nosé-Hoover thermostat. Configurations were sampled every 0.1 ps, and the production runs were 3.0 ns. Selfdiffusion coefficients along the clay surfaces were calculated using the two-dimensional Einstein relation:

D ∥ = lim t→∞ ⟨∆x(t) 2 + ∆y(t) 2 ⟩ 4t , ( 3
)
where < ∆x(t) 2 > is the mean-square displacement of the particle in the x direction.

Results and discussion

Thermodynamics Monte Carlo simulation in the µ H 2 O µ CO 2 V T ensemble requires the knowledge of the chemical potentials µ H 2 O and µ CO 2 , which are fixed by the temperature T and pressure P in the reservoir.

Thus, we first determined the chemical potentials of water and carbon dioxide so as to obtain in grand-canonical µ H 2 O µ CO 2 V T density and composition of both bulk phases close to experimental data 52,[START_REF] Gillespie | Vapor-Liquid and Liquid-Liquid Equilibria: Water-Methane, Water-Carbon Dioxide, Water-Hydrogen Sulfide, Water-nPentane, Water-Methane-nPentane[END_REF] and molecular simulations in the Gibbs ensemble. 47 1, the composition of both phases (CO 2 -poor and CO 2 -rich) in the reservoir, as well as the density of the CO 2 -rich phase (gas vs. supercritical fluid) are not the same at 25 and 125 bar.

Figure 2 shows how the pressure normal to the clay layers varies with the interlayer distance.

Each data point, indicated with a symbol, is taken from 1 to 5 independent simulations, while the lines are guides to the eye. A quantitative difference is observed in the first pressure minimum.

This results in a larger free energy barrier for the transition from monolayer to bilayer state for the mixture, as shown in Figure 3. We calculated the swelling free energy per unit area for each chosen set of chemical potentials using Eq. ( 1). For each reservoir condition, there are only two that there is no difference (within the error bars) between pure water and the mixture for bilayer states, while the monolayer spacing is slightly larger for pure water. The bilayer state for all curves is favored with a free energy much lower than that of the monolayer.

The CO 2 interlayer mole fraction for mixtures at 25 and 125 bar are presented in Figure 4. The horizontal lines represent the CO 2 mole fraction in the bulk H 2 O-rich phase. In contrast to the reservoir (bulk), there is only one stable phase in the interlayer, with a CO 2 mole fraction larger in the clay phase than in the bulk one. We checked that these results do not depend on the initial composition of the simulated system. From a more quantitative point of view, the monolayer for pure water in the reservoir consists of 6.2 water molecules per sodium ion; for the mixture at This observation persists despite very long simulation times and the use of temperature annealing, and is not a simulation artefact. A similar feature was found for an ethane-water mixture in clay. 39 Nevertheless the decrease in CO 2 mole fraction from the monolayer to the bilayer mainly originates from the increase in water molecules.

These equilibrium states determined by grand-canonical Monte-Carlo were then simulated using Molecular Dynamics to investigate the interlayer structure and dynamics. For these simulations, the relative horizontal position of the clay surfaces was fixed to the most probable configuration as determined during the MC simulations and will be discussed bellow. to the surface, despite the fact that contrary to the monolayer case there would be enough space to accommodate a perpendicular orientation.

Structure

In order to clarify the situation, we plot in Figure 6 the x-y distributions for O CO 2 and C atoms.

The oxygen atoms tend to arrange near the center of the clay hexagonal cavities, whereas C atoms have a small tendency to stay close to O clay atoms. In the bilayer state, one oxygen of the CO 2 molecule enters in the hexagonal cavities, as can be inferred from the location in the x-y plane of the oxygen atoms that are closest to the surface (|z| ∈ [2.0:4.5] Å see Figure 5). The projected C-O CO 2 distances (in the x-y plane) are different for monolayer and bilayer states, which is in agreement with the density profiles: In the former case the distance is larger and corresponds to CO 2 parallel to the surface, while for the bilayer state CO 2 molecule orientations are distributed rather uniformly around the center of the cavity. These observations are in agreement with MD simulations for carbon dioxide in dehydrated uncharged clays. 34 This means that the charge of the clay surface and the associated presence of counterions and water molecules does not seem to play an important role in this case. Furthermore, horizontal movements of the clay layers, performed using MC simulation, show that in the monolayer case CO 2 molecules fully determine the position of each layer with respect to the other, a feature observed with Sr 2+ counterions. 13 The preferred position is shown in Figure 6, where the two clay layers are represented by blue and red-yellow hexagons, respectively. Such an arrangement can be explained by the fact, that both oxygen atoms of the same CO 2 molecule tend to arrange near the center of the hexagonal cavity at the same time (one per surface) as illustrated on Figure 6. Note that we have not found any preferred positions for bilayer state. This means that CO 2 molecules near one surface do not influence the position of the second one, and is consistent with the entrance of the CO 2 molecule in the hexagonal cavity.

Transport

The transport properties of the interlayer fluid can be characterized by the self-diffusion coefficients of the species parallel to the clay surface (D ∥ ). The trajectories obtained from molecular dynamics simulation were analyzed and D ∥ were calculated as the longtime limiting slope of mean square displacements (MSD) versus time (see Eq. ( 3)). Table 2 reports the values of the diffusion coefficients calculated from the MSD between 0.5 and 1.0 ns. The standard deviation was calculated over four independent simulations. The experimental diffusion coefficients of H 2 O

given in Table 2 were obtained by collaborators in synthetic hectorite at T =347 K by neutron spin echo. [START_REF] Dubois | Articles in progress[END_REF] Although hectorite clay is different from montmorillonite, it has the same swelling properties and gives similar diffusion coefficients at ambient temperature. 55 All data, given in the table for comparison, correspond to P=1 bar. Diffusion coefficients generally decrease with an increase in pressure. 56 This means that the ratio of our simulations over the corresponding experimental data for bulk phase, which is already equal to ∼1.4, should probably be higher if we compared with experimental data at higher pressures (see Table 3). This discrepancy can be traced back to the inaccuracy of the SPC model for describing the dynamical properties of water. 57 Nevertheless, we can estimate the influence of CO 2 on the diffusion of interlayer species by analyzing the ratio D ∥ /D bulk , where D bulk are the values of D simulated in the bulk solution. This allows to compare with the previous studies of H 2 O-clay systems. 58 In the absence of CO 2 , we find that D ∥ /D bulk for water and counterions are smaller than in Ref. 58, both in the mono-and bilayer cases, as expected for the higher pressure investigated in the present work.

The results given in Table 2 indicate that the presence of carbon dioxide molecules reduces the mobility of interlayer H 2 O and Na + . To reach this conclusion, one must separate the effect of confinement by the clay layers, which can be deduced from the simulations in the absence of CO 2 . This confinement results in a decrease of the diffusion coefficient, relative to the bulk, by a factor ∼ 2 in the bilayer and ∼ 10 in the monolayer cases, respectively. These values are consistent with the experimental data at 347 K (see Table 2). In the presence of CO 2 , this reduction is even more dramatic, with a factor of ∼3 for the bilayer, and ∼25 for the monolayer. Thus, the diffusion of H 2 O and Na + in the presence of CO 2 for the monolayer is about 2.5 times smaller than in the absence of CO 2 and only 1.5 times for the bilayer. Since the CO 2 content is higher in the former case (see Figure 4), we conclude that the more CO 2 molecules in the interlayer, the smaller the diffusion coefficient of all species. Such a behaviour is of course not unexpected, as crowding by bulky molecules such as CO 2 hinders the diffusion of other species. Finally, it is worth noting that the smallest D ∥ /D bulk ratios are obtained for the CO 2 molecules themselves. This can be explained by the apparently stronger interactions of these molecules with the surface, reflected in the molecular profiles discussed in section Structure. Such a behaviour is reminiscent of the observed difference between isovalent cations, for which it has been observed that larger ions (e.g. Cs + compared to Na + ) tend to remain closer to the surface than smaller ones, and that their diffusion in the interlayer is more reduced compared to the bulk. 59

Conclusion

We used grand-canonical Monte Carlo and Molecular Dynamics simulations to study the Namontmorillonite hydrates in equilibrium with carbon dioxide rich reservoirs. From the µV T study, we determined H 2 O and CO 2 contents as a function of the interlayer distances. These results show that the hydrated clay system is capable of adsorbing CO 2 molecules. The swelling free energy curves indicate that for the considered conditions neither swelling nor shrinkage should be observed. For monohydrated systems containing carbon dioxide, CO 2 molecules lie in the center of interlayers, parallel to the clay surface, and they stabilize a particular position of the clay surfaces relative to each other. For the bilayer, which contains more water molecules but no additional carbon dioxide, the linear O=C=O molecules are located close to the clay surface with one O atom above the center of hexagonal cavities and slightly tilted relative to the surface. The presence of interlayer CO 2 inhibits the diffusion of all the mobile species.

In the future, we will investigate how much the results depend on the force field and describe the system, using e.g. the CLAYFF force field. 60,61 Another interesting aspect is the effect of the reservoir composition, taking into account the presence of salt and dissolved HCO - 3 and CO

  Na + 0.75 [Si 8 ](Al 3.25 Mg 0.75 )O 20 (OH) 4 .The layered structure of the clay is simulated by two clay layers, using periodic boundary conditions, as done previously e.g. in Refs. 11,14,43,44. One of the clay layers is situated in the centre of the simulation box, while the other is divided in two half-layers located at the top and bottom of the box, as shown in Figure1. Na + cations, water and carbon dioxide molecules are distributed in the interlayer spaces. Each clay layer consists of 8×4 clay unit cells, totalling 1280 atoms, and is treated as a rigid body. The water and carbon dioxide molecules are simulated using the SPC[START_REF] Berendsen | Interaction models for water in relation to protein hydration[END_REF] and EPM2 46 models, respectively. The combination of these two models was shown to correctly predict the phase behaviour of H 2 O-CO 2 mixtures in the range of pressures and temperatures considered in the present work.47 Pairwise additive Lennard-Jones and Coulomb potentials are used to model interactions between particles:

1 j=1

 1 ∆u LJ 1 is the change in the LJ contribution to the energy produced by the first atom insertion, W 1 = ∑ k exp(-β ∆u LJ 1 ) is the Rosenbluth factor. Then we generate k 2 trial configurations for the other atoms. An orientation of the molecule is selected with probability p 2 calculated similarly to p 1 , but with ∆u LJ 2 the change in the LJ contribution to the energy for all atoms of the molecule, except the central one. The same procedure applies for the molecule's removal, though in this case we generate only k-1 trial configuration,

  box size was taken to correspond to 400-500 molecules of H 2 O and CO 2 for each considered density. An approximate set of chemical potentials was first obtained by trying 50 (µ H 2 O ; µ CO 2 ) pairs in the range [-49:-42]×[-40:-34] kJ/mol. More precise values were then determined by dichotomy (see next section for a discussion of sensitivity). In the case of the clay phase the simulation box corresponds to a periodically replicated parallelepiped with angles adjustable on-the-fly, which allows us to shift clay layers in the horizontal directions without destroying the structure. This is necessary since the relative arrangement of the clay surfaces has been shown to depend on the interlayer distance and content.13,48 The x-y dimensions of the simulation box are 41.44×35.88 Å 2 ,

  free energy minima points: The first one represents a monolayer state (a layer one molecule thick), the second a bilayer (two molecular layers); all other intermediate states are unstable. They are found for the mixture under 25 bar at 12.3 (monolayer) and 15.4 Å(bilayer), for the mixture under 125 bar at 12.3 and 15.5 Å, and for pure water at 125 bar at 12.4 and 15.3 Å. It should be noted

Figure 5

 5 Figure 5 displays the density profiles for water oxygen (O H 2 O ) and hydrogen (H H 2 O ), carbon dioxide oxygen (O CO 2 ) and carbon (C) and sodium (Na) atoms perpendicular to clay surface. The data for the mixtures at 25 and 125 bar are very similar and we report only the results at 125 bar. The profiles obtained in the mono-and bilayer cases are compared to the clay system with pure water only. As can be seen from the figure, CO 2 molecules hardly influence the distribution of other atoms. The distance between the C peaks and the surface are almost the same in the mono-and bilayer states. This indicates that C have a tendency to remain close to the clay surface, and is

  The radial distribution functions (r.d.f.) for theNa-O H 2 O , O clay -C and O H 2 O -C pairs for both interlayer spacings at P=125 bar are shown in Figure 7. Changes in the Na-O H 2 O r.d.f. due to the presence of CO 2 are very small. The first peak for O H 2 O -C is situated at 3.92 Åin both the clay interlayer and in the bulk phase. The nearest-neighbour coordination number was calculated by integrating over the first peak of the correlation function. C atoms are surrounded by approximately the same number of oxygen atoms for both interlayer hydration state (15.9 for mono-and 16.7 for bilayer, respectively). Interestingly, however, the nature of these atoms differ: There are 4.9 O H 2 O and 11 O clay around each carbon atom in the monolayer case, and 11 O H 2 O and 5.7 O clay in the bilayer case.

64 Figure 1 :

 641 Figure 1: Schematic representation of the simulation box.

Figure 2 :

 2 Figure 2: Pressure as a function of the interlayer distance at T = 348.15 K. The black solid, red dashed and blue dashed-dotted lines are interpolations of the data points for reservoir conditions corresponding to the mixture at 25, 125 bar and pure water at 125 bar, respectively.

Figure 3 :

 3 Figure 3: Clay swelling free energy per unit area plotted as a function of the interlayer distance.

Figure 4 :

 4 Figure 4: CO 2 mole fraction plotted as a function of interlayer distance for Na + -montmorillonite at T =348.15 K. The horizontal lines indicate the CO 2 mole fraction in the bulk H 2 O-rich phase. The grey stripes visualize the location of the stable states.

Figure 5 :

 5 Figure 5: Density profiles of the C, O, H and Na atoms for the mixture (first column) and pure water (second column) at P=125 bar, T =348.15 K. The first and second rows correspond the monolayer and bilayer states, respectively. The vertical dashed lines represent the clay surface oxygen atoms.

Figure 6 :

 6 Figure 6: Distribution of O CO 2 (left) and C(right) atoms parallel to the clay surface. Light regimes correspond to high density. The first and second rows correspond the monolayer and bilayer states, respectively. In the monolayer case, a typical configuration of the CO 2 molecule is indicated. In the bilayer case, results for O CO 2 correspond only to the atom closest to the surface (|z| ∈[2.0:4.5] Å see Figure 5). Results are given for one unit cell with dimensions a × b = 5.18 × 8.97 Å 2 . Oxygen surface atoms are indicated for one clay surface by red circles, silicon atoms by yellow circles. The location of the second surface is indicated by the blue hexagons.

Figure 7 : 21 Figure 1 .

 7211 Figure 7: Radial distribution functions. Results are compared to a bulk ionic solution for Na-H 2 O, and to the bulk CO 2 -poor liquid for O H 2 O -C.

  app = 25 bar mixture, P app = 125 bar pure water, P app = 125 bar16 17

Figure 2 .

 2 Figure 2.

Figure 6 .

 6 Figure 6.

Figure 7 .

 7 Figure 7.

  Both quantities are most sensitive to changes in the chemical potential of CO 2 . Depending on the bulk phase (H 2 O-rich or CO 2 -rich), notable changes are observed either in the density or in the mole fraction. In the liquid H 2 O phase a change of 0.25 kJ/mol in µ CO 2 leads to changes in mole fraction up to 20 % (in its absolute value) and in density of only 1 %. For the CO 2 -rich phase such variations of µ CO 2 result in significant density changes (up to 30 %) and negligible changes in the mole fraction (1 %). Changing µ H 2 O by 0.25 kJ/mol leads to small changes in density or mole fraction (less 3 %). The retained values for the chemical potential are shown in Table1. The chemical potential for pure water at T =348.15 K is found with the same approach to be -46.95 kJ/mol. These values were then used for three different sets of GCMC simulations in Na-montmorillonite clay at 348.15 K, corresponding to three different reservoir conditions: 1) pure water at 125 bar, which serves as a reference (before injection of CO 2 ), 2) CO 2 -H 2 O mixture at P app =25 bar, and 3) CO 2 -H 2 O mixture at 125 bar. As can be seen in Table

  25 bar, 5.2 H 2 O and 0.6 CO 2 ; and for the mixture at 125 bar, 4.9 H 2 O and 0.8 CO 2 . For the bilayer under the same reservoir conditions we find 11.7 H 2 O per Na + , 11.6 H 2 O and 0.3 CO 2 , and 11.4 H 2 O and 0.6 CO 2 , respectively. These data show that the total interlayer density increases with pressure. It should be noted that the number of CO 2 molecules is larger for the monolayer.

Table 1 :

 1 2- Estimated chemical potentials in

	µ H 2 O µ CO 2 V T simulation (bulk phase) for coexistence curve of H 2 O-CO 2 mixtures, based on and experimental data. 52,53 The statistical uncertainties in the last digits are given in parentheses.	T, P, µ H 2 O , µ CO 2 , ρ H 2 O-rich , x H 2 O-rich CO 2 ρ CO 2 rich , (kJ/mol) (kJ/mol) (g/cm 3 (K) (bar) ) mole fraction (g/cm 3 ) x CO 2 -rich CO 2 mole fraction	348.15 25.33 ---0.0055 -0.9894	348.15 25.33 ---0.0054 -0.9818	348.15 126.66 -----0.9915	348.15 25 --0.944(8) 0.0034(13) 0.0401(5) 0.9777(70)	348.15 125 --0.948(10) 0.0140(24) 0.313(18) 0.9940(8)	348.15 24.23(7) a -46.85 -39.25 0.934(3) 0.004(1) 0.039(1) 0.982(1)	348.15 122.35(98) a -46.85 -35.35 0.940(5) 0.014(1) 0.298(3) 0.994(1)
	+Gibbs ensemble simulation 47		Experiment 52,53			Simulation 47	NPT+Gibbs simulation	µ H 2 O µ CO 2 V T	(this work)
	NPT

a

The pressure was calculated in the CO 2 -rich phase.

Table 2 :

 2 Diffusion coefficient data obtained in this work for the clay phase, T=348.15 K and some of those presented by Malikova et al58 and collaborators.54 The standard deviations in the last digits are given in parentheses.

		bulk
		D ∥ /D		0.02							0.19
	2 CO s -1
		D ∥ /10 -9 m 2	0.24 (4)	0.17 (2)						1.6 (2)	1.5 (1)
		bulk
		D ∥ /D		0.03	0.08						0.30	0.46
	Na +	s -1
		D ∥ /10 -9 m 2	0.19 (2)	0.13 (2)	0.36 (3)		0.41(5)	0.70(8)		1.3 (1)	1.4 (1)	2.1 (1)		1.3(1)	3.7(5)
		bulk
		D ∥ /D		0.04	0.12	0.13	0.18	0.19			0.39	0.48	0.35	0.56	0.60
	2 O s -1 H D ∥ /10 -9 m 2	0.51 (2)	0.32 (2)	0.99 (2)	0.75	0.83(10)	1.6(1)		3.5 (1)	3.2 (1)	3.9 (1)	2.1	2.0	4.2
		Monolayer	CO 2 -H 2 O-clay, P=25 bar, this work	CO 2 -H 2 O-clay, P=125 bar, this work	H 2 O-clay, P=125 bar, this work	H 2 O-clay, T=347 K, experiment a 54	H 2 O-clay, T=330 K, simulation a 58	H 2 O-clay, T=375 K, simulation a 58	Bilayer	CO 2 -H 2 O-clay, P=25 bar, this work	CO 2 -H 2 O-clay, P=125 bar, this work	H 2 O-clay, P=125 bar, this work	H 2 O-clay, T=347 K, experiment a 54	H 2 O-clay, T=320 K, simulation a 58	H 2 O-clay, T=360 K, simulation a 58

a P=1 bar.

Table 3 :

 3 Diffusion coefficient data obtained in this work for the bulk phase, T=348.15 K and some of those presented by Easteal et al,62 Thomas and Adams.63 The standard deviations in the last digits are given in parentheses.

		s -1
	CO 2	D/10 -9 m 2
		s -1
		m 2
	Na +	D/10 -9		4.6 (4)			3.5
		s -1
	H 2 O	D/10 -9 m 2	8.2 (1)	8.1 (1)	5.646	6.582
			CO 2 -H 2 O, P=125 bar, this work Na + in water, P=125 bar, this work	62 H 2 O, T=343 K, experiment a	62 H 2 O, T=353 K, experiment a	Na + in water ab	CO 2 -H 2 O, experiment a

ions. We are also investigating large (interparticle) pores to clarify the crossover between the surface dominated regime in the interlayer, where only the water-rich phase was found to be stable to the bulk regime where two phases coexist.

Acknowledgement

We thank the Institut Français du Pétrole for permission to publish these results. We are grateful to Prof. J.-P. Hansen, Prof. J.-F. Dufrêche and Dr. C. Nieto-Draghi for valuable discussions.

Graphical TOC Entry


Sami Tazi

Benjamin Rotenberg
email: benjamin.rotenberg@upmc.fr

Mathieu Salanne

Michiel Sprik

Marialore Sulpizi

Absolute acidity of clay edge sites from ab-initio simulations

I. INTRODUCTION

The acido-basic properties of clay minerals play a key role in their surface chemistry.

The increased acidity of water confined in smectite clay interlayers [START_REF] Bergaya | Handbook of clay science[END_REF], for example, contributes to its catalytic properties. The protonation state of the edges of clay layers has an impact on both the stability of colloidal clay suspensions [START_REF] Tombácz | Colloidal behavior of aqueous montmorillonite suspensions: the specific role of pH in the presence of indifferent electrolytes[END_REF] and on the sorption of ions [Baeyens andBradbury 1997, Bradbury andBaeyens 1997]. This is particularly important since the retention of contaminants by argillaceous rocks is pivotal in the current concepts for the geological disposal of long-lived high-level nuclear waste in several countries such as Switzerland, France and Belgium [ANDRA 2005, Bradbury and[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]. Clays are lamellar aluminosilicate minerals. While the negative charge of their basal surface builds up because of isomorphic substitutions in the mineral lattice, that of edge sites (which arise from the finiteness of the sheets) is due to the deprotonation or protonation of surface groups. The protonation state is then controlled by the acidity of the various groups and the pH of the solution.

Several difficulties arise when trying to address the question of acidity from an experimental point of view. The most problematic one for the assignment of individual pK a 's to the different surface groups is that macroscopic measurements such as titration experiments cannot distinguish between them. Moreover, Bourg et al. have shown recently that the analysis of such experiments must be performed with great care, in particular to correctly take into account the surface charge at the onset of titration [START_REF] Bourg | Modeling the acid-base surface chemistry of montmorillonite[END_REF]]. In the wait for the use of site-specific techniques (e.g. spectroscopic ones), the most promising approach has been to interpret experiments within the framework of a conceptual model, relying on assumptions (e.g. on the reactions to consider or how to decouple electrostatic from chemical contributions) and associated parameters (equilibrium constants for the various reactions) [Baeyens andBradbury 1997, Bradbury andBaeyens 1997]. The predictions of the model are then compared to the limited and global experimental data. In the best case, a favourable comparison could still occur by a compensation of errors between the different assumptions and parameters underlying the prediction [START_REF] Bourg | Modeling the acid-base surface chemistry of montmorillonite[END_REF]]. In order to limit this risk, Tournassat et al. have proposed to constrain a large number of parameters from structure-based models [START_REF] Tournassat | The titration of clay minerals ii. structure-based model and implications for clay reactivity[END_REF]]. More precisely, the density of each site can be predicted from the structure of the mineral layer and the corresponding acidity constants estimated using semi-empirical bond-valence based models such as MUSIC [START_REF] Bickmore | Bond-valence methods for pKa prediction. ii. bond-valence, electrostatic, molecular geometry, and solvation effects[END_REF], 2004[START_REF] Brown | Recent developments in the methods and applications of the bond valence model[END_REF][START_REF] Hiemstra | Intrinsic proton affinity of reactive surface groups of metal (hydr)oxides: The bond valence principle[END_REF][START_REF] Tournassat | The titration of clay minerals ii. structure-based model and implications for clay reactivity[END_REF]]. Nevertheless, there has been to date no direct way to probe the charge density on the edge and to measure the deprotonation free energies of the different sites. Another indirect strategy to determine the overall protonation state of clay edges has been proposed by Zhao et al. who found using Atomic Force Microscopy that the force between a silica tip and a surface consisting of muscovite edges turned from repulsive at pH∼10 to attractive at pH∼6, from which they concluded that these edge surfaces had a point of zero charge of pH 7-8. [START_REF] Zhao | Probing surface charge potentials of clay basal planes and edges by direct force measurements[END_REF] In this respect atomistic simulations have emerged as a powerful tool to address the microscopic properties of clays. While force-field based, classical molecular simulations have for example provided insights on bulk clays (including interlayer water and ion dynamics [START_REF] Ferrage | Hydration properties and interlayer organization of water and ions in synthetic na-smectite with tetrahedral layer charge. part 2. toward a precise coupling between molecular simulations and diffraction data[END_REF][START_REF] Marry | Microscopic simulations of interlayer structure and dynamics in bihydrated heteroionic montmorillonites[END_REF][START_REF] Marry | Microscopic simulation for structure and dynamics of water and counterions in a monohydrated montmorillonite[END_REF][START_REF] Rinnert | Hydration of a synthetic clay with tetrahedral charges: A multidisciplinary experimental and numerical study[END_REF], swelling [START_REF] Boek | Monte Carlo molecular modelling studies of hydrated Li-, Na-and K-smectites: Understanding the role of potassium as a clay swelling inhibitor[END_REF][START_REF] Tambach | Molecular simulations of swelling clay minerals[END_REF]] and ion exchange [Rotenberg et al. 2009, Teppen and[START_REF] Teppen | Hydration energy determines isovalent cation exchange selectivity by clay minerals[END_REF]), or on ion sorption [START_REF] Greathouse | Molecular dynamics simulation of uranyl(vi) adsorption equilibria onto an external montmorillonite surface[END_REF]] and water at basal surfaces [START_REF] Marry | Structure and dynamics of water at a clay surface from molecular dynamics simulation[END_REF][START_REF] Rotenberg | Molecular explanation for why talc surfaces can be both hydrophilic and hydrophobic[END_REF][START_REF] Wang | Effects of substrate structure and composition on the structure, dynamics, and energetics of water at mineral surfaces: A molecular dynamics modeling study[END_REF][START_REF] Brown | Recent developments in the methods and applications of the bond valence model[END_REF], their application to clay edges have been so far limited [START_REF] Rotenberg | Water and ions in clays : Unraveling the interlayer/micropore exchange using molecular dynamics[END_REF][START_REF] Rotenberg | Molecular simulation of aqueous solutions at clay surfaces[END_REF]. One of the main reasons is that the description of interfaces often requires to address polarization effects and, in the case of clay edges, reactivity. Electronic structure calculations provide a direct way to include such effects describing both the solid and the liquid in a consistent way. In the last years Density functional Theory (DFT)-based molecular dynamics simulations have provided a first principle understanding of bulk clays such as the properties of interlayer ions [Boek andSprik 2003, Suter et al. 2008] and their vibrational properties [START_REF] Larentzos | An ab initio and classical molecular dynamics investigation of the structural and vibrational properties of talc and pyrophyllite[END_REF], or of water at basal surfaces [START_REF] Bridgeman | Ab-initio total energy study of uncharged 2:1 clays and their interaction with water[END_REF][START_REF] Tunega | Theoretical study of adsorption sites on the (001) surfaces of 1:1 clay minerals[END_REF], 2004]. More recently, the acidity of interlayer water [Churakov andKosakowski 2010, Liu et al. 2011] and the structure and reactivity of clay edges [START_REF] Bickmore | Ab initio determination of edge surface structures for dioctahedral 2:1 phyllosilicates: implications for acid-base for reactivity[END_REF][START_REF] Kremleva | Uranyl adsorption at (010) edge surfaces of kaolinite: A density functional study[END_REF][START_REF] Liu | Surface complexes of acetate on edge surfaces of 2:1 type phyllosilicate: Insights from density functional theory calculation[END_REF], 2012] have also been investigated. In particular, Churakov has studied water confined between edges of pyrophyllite. He evidenced transient proton exchange between surface groups mediated by surface water molecules [START_REF] Churakov | Structure and dynamics of the water films confined between edges of pyrophyllite: a first principle study[END_REF]] and proposed a first attempt to provide an acidity scale of the surface groups from their deprotonation energy in vacuum [START_REF] Churakov | Ab initio study of sorption on pyrophyllite: structure and acidity of the edge sites[END_REF]]. In this context, it has also been suggested that ab-initio simulations could provide refinements of the bond-valence based models of acidity by introducing the relaxation of bond lengths [START_REF] Bickmore | Ab initio determination of edge surface structures for dioctahedral 2:1 phyllosilicates: implications for acid-base for reactivity[END_REF]] (whose importance was already suggested by [START_REF] White | Analysis and implications of the edge structure of dioctaedral phyllosilicates[END_REF]) or the solvation structure (coordination and bond lengths) [START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF]] determined from DFT simulations.

One possible approach to compute the acidity of surface sites is to compute the potential of mean force along a suitable reaction coordinate, as was done by Leung et al. for silica [Leung et al. 2009] and goethite [START_REF] Leung | Predicting the acidity constant of a goethite hydroxyl group from first principles[END_REF]. DFT-based ab-initio simulations can also be combined with thermodynamic integration to determine the free energy of proton transfer from a given site to a water molecule, which provides a direct measure of the acidity of this site. This strategy has been successfully employed to calculate the pK a of simple aqueous compounds [START_REF] Cheng | Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics[END_REF], Costanzo et al. 2011[START_REF] Mangold | Absolute pK(a) values and solvation structure of amino acids from density functional based molecular dynamics simulation[END_REF][START_REF] Sulpizi | Acidity constants from dft-based molecular dynamics simulations[END_REF], 2008] as well as the surface acidity and point of zero charge of simple oxides in contact with water, namely titania [START_REF] Cheng | Acidity of the aqueous rutile TiO2(110) surface from density functional theory based molecular dynamics[END_REF], quartz [START_REF] Sulpizi | The silica-water interface: how the silanols determine the surface acidity and modulate the water properties[END_REF] and alumina [START_REF] Gaigeot | Oxide/water interfaces: how the surface chemistry modifies interfacial water properties[END_REF]]. In the present work, we apply this strategy to provide the first direct calculation of the pK a of the different groups which are present on the (010) face of pyrophyllite. We first introduce the simulated system and computational method used to determine the deprotonation free energy. We then compute the pK a of the three edge sites present on this surface, namely SiOH, AlOH and AlOH 2 . We further analyze the changes in structure and solvation of each site upon deprotonation and relate them to the reorganization free energies. Finally, we compare our results to previous estimates for the same edges and for similar sites on other mineral surfaces, and discuss the predictions of the MUSIC model.

II. METHODS

A. Acidity constant calculations

The acidity constants of the different reactive groups present on the clay edge are computed using the reversible proton insertion/deletion method [START_REF] Cheng | Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics[END_REF], Costanzo et al. 2011[START_REF] Mangold | Absolute pK(a) values and solvation structure of amino acids from density functional based molecular dynamics simulation[END_REF][START_REF] Sulpizi | Acidity constants from dft-based molecular dynamics simulations[END_REF], 2008]. pK a 's are estimated from the reaction free energies for the transfer of a proton from the surface to a water molecule in the bulk part of the sample. For example, if we consider the silanol site as a reference, this reads:

-SiOH + H 2 O aq → -SiO -+ H 3 O + aq (1)
The free energy of this transfer process is computed using thermodynamic integration. The charge of the acidic proton of a surface OH group is gradually switched off transforming it into a neutral dummy ("ghost") hydrogen atom which is tethered to the oxygen but otherwise does not interact with the rest of the system. Simultaneously a similar dummy hydrogen attached to a water molecule is charged up, creating a hydronium H 3 O + in the Eigen form. Following the approach discussed in detail in Ref. [START_REF] Costanzo | The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode[END_REF]], the reaction free energy for the deprotonation process is calculated according to the following integral [START_REF] King | Investigation of the free energy functions for electron transfer reactions[END_REF]:

∆A = 1 0 dη ∆E η ( 2
)
where ∆E is the vertical energy gap, defined as the potential energy difference between products P (-SiO -+ H 3 O + aq in Eq. 1) and reactants R (-SiOH + H 2 O aq ) for instantaneous configurations of a molecular dynamics trajectory. The subscript η indicates that the averages are evaluated over the restrained mapping Hamiltonian

H η = (1 -η) H R + η H P ,
where η is a coupling parameter which is gradually increased from 0 (-SiOH + H 2 O aq ) to 1 (-SiO -+ H 3 O + aq ). This Hamiltonian also contains a harmonic restraining potential V restr , independent of η, keeping the dummy atoms close to the equilibrium position of the H + nucleus in the protonated site and hydronium, as explained in the next section. The Simpson rule (three point approximation) is used to evaluate the integral in Eq. 2:

∆A TP = 1 6 ( ∆E 0 + ∆E 1 ) + 2 3 ∆E 0.5 (3)
This requires the generation of three trajectories corresponding to values of η=0, 0.5, 1.

Simulation length for each value of η are reported in Table I. The three point formula is often a good compromise between computational cost and accuracy of the free energy change [START_REF] Hummer | Calculation of free energy differences from computer simulations of initial and final states[END_REF]Szabo 1996, Sulpizi andSprik 2008]. However, if the curvature of ∆E η is large, more integration points may be required. According to the derivation in Ref. [START_REF] Costanzo | The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode[END_REF]] the expression for the pK a , in the low temperature limit to the quantum vibrational partition function of the acid proton, reads:

2.30k B T pK a = k B T ln c o Λ 3 H + + 1 0 dη ∆E η (4)
where c o = 1 mol dm -3 is the unit molar concentration, and Λ H + is the thermal wavelength of the proton. The logarithm of the product c o Λ 3 H + accounts for the liberation entropy of the proton and is responsible for a correction of -3.2 pK a units to the thermodynamic integral.

Finite system size effects are of course a point of concern in these calculations [START_REF] Costanzo | The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode[END_REF]]. In the present case, since the system is overall neutral, we do not expect issues related to the presence of net charges (and compensating background) as also discussed in Ref. [START_REF] Churakov | Ab initio study of sorption on pyrophyllite: structure and acidity of the edge sites[END_REF]]. The dimensions of the system considered here (see below) are similar to previous studies of other mineral surfaces for which finite site effects could be estimated as being smaller than the uncertainty on the final pK a value.

When considering the intrinsic acidity of different sites on the same surface, one must ensure that less acidic sites will not be reprotonated by proton transfer from more acidic ones.

Similarly, for sites whose pK a is outside the water domain (0-14), one must prevent proton transfer from water to the deprotonated sites, to really sample the products of the surface deprotonation reaction. This can be achieved in simulations by introducting additional restraining potentials on O-H bonds that are not considered for the deprotonation. Since these constraints are present in both the initial and final states, they do not contribute to the free energy difference (as for the tethering of the dummy hydrogens) in Eq. 2. Essentially, the pK a value reflects the absolute acidity of a given site, conditional to the fact that the rest of the surface remains unchanged, i.e. as if this site was isolated on the surface. Once these values have been determined, one can draw a predominance diagram to predict the actual state of the surface as a function of pH.

B. System and simulation details

We simulate a finite neutral pyrophyllite particle in contact with a slab of bulk liquid water along its (010) face. The simulation box contains 312 atoms, namely 12 aluminum, 24 silicium, 80 oxygen (of which 8 form SiO edge sites, 8 AlO edge sites and 8 OH groups), 28 hydrogen atoms (of which 8 form SiOH edge sites 4 AlOH edge sites, 8 AlOH 2 edge sites and the 8 last the OH groups), 56 water molecules. The unit cell is orthorhombic with dimensions L x =10.36 Å, L y =9.35 Å and L z =31.0 Å. Periodic boundary conditions allow to simulate an infinite stack along y of sheets which are infinite along the x direction and finite along the z one (See Fig. 1).

The initial configuration for the DFT-based molecular dynamics is prepared using a classical molecular dynamics simulation in which the clay particle is fixed and the water equilibrated in the NVT ensemble at 300 K. The SPC/E model [START_REF] Berendsen | The missing term in effective pair potentials[END_REF] for water and the CLAYFF force field [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF] for clay are used and the equations of motion for the rigid water molecules are integrated using the SHAKE algorithm [START_REF] Ryckaert | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF]. DFT-based simulations are performed with the PBE functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], analytic Goedecker-Teter-Hutter [START_REF] Goedecker | Separable dual-space gaussian pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic separable dual-space gaussian pseudopotentials from h to rn[END_REF] ]. For the DFT calculations, the latter uses an algorithm based on the hybrid Gaussian and plane wave method [START_REF] Lippert | A hybrid gaussian and plane wave density functional scheme[END_REF]]. In order to analyze the electronic charge redistribution upon deprotonation, we have computed the density-derived atomic partial charges (DDPAC) according to Ref. [START_REF] Blöchl | Electrostatic decoupling of periodic images of plane-wave expanded densities and derived atomic point charges[END_REF]. A charge calculation is performed for configurations along the molecular dynamics trajectory, and average values are reported.

The harmonic restraining potential V restr used for the thermodynamic integration to maintain the dummy proton close to the equilibrium position reads:

V restr = n d i=1 1 2 k d,i (d i -d 0,i ) 2 + nα j=1 1 2 k α,j (α j -α 0,j ) 2 (5)
where the first sum corresponds to restraining the distance d from the nearby oxygen atom to the equilibrium value d 0 in the protonated form and the second to the angles involving the dummy proton. The equilibrium values are determined from a prior simulation without restraints, while the restraining parameters have been chosen according to the prescriptions of Ref. [START_REF] Sulpizi | Acidity constants from vertical energy gaps: density functional theory based molecular dynamics implementation[END_REF]. Details on the number of terms and on the corresponding parameters are reported for each site in Table II. The choice of a particular water molecule to form the hydronium is not relevant provided that this molecule is beyond the second solvation shell of the acidic surface proton under investigation.

III. RESULTS AND DISCUSSION

A. Structure of the interface

The structure of the (010) face of pyrophyllite interface with a water film has already been investigated in detail by Churakov [START_REF] Churakov | Structure and dynamics of the water films confined between edges of pyrophyllite: a first principle study[END_REF]], and we provide here the main features to serve as a reference for the discussion of the deprotonated sites. Moreover, we have used different simulation conditions, in particular a larger clay fragment and a wider water region, and this preliminary analysis allows to investigate the possible differences. The water density is equal to the bulk density except for the adsorbed layers, and its structure differs slightly from that of liquid water because of confinement. Nevertheless, it remains liquid-like. As can be seen in Fig. 1 to surface groups. Within small quantitative variations, all these results confirm the findings of Churakov [START_REF] Churakov | Structure and dynamics of the water films confined between edges of pyrophyllite: a first principle study[END_REF]]. Overall, the coordination number of the oxygen atom on these groups is appproximately 3 for AlOH 2 (2 H and 1 Al), 4 for AlOH (3 H and 1 Al) and 3 for SiOH (2 H and 1 Si).

B. Acidity of surface groups

We now turn to the calculation of the pK a of the edge sites using the thermodynamic integration method explained above, during which the acid proton in each OH group on the clay edge is gradually removed and transferred into the aqueous solution. The convergence of the three vertical energy gaps corresponding to η = 0, 0.5 and 1 is illustrated for the SiOH site in Fig. 3. The pK a of the different surface groups are given in Table IV. We find that the most acidic site is SiOH, with a pK a of 6.8 ± 0.4, closely followed by AlOH 2

(pK a = 7.6 ± 1.3). Finally the AlOH group is the less keen to loose a proton, with a pK a of 22.1 ± 1.0. In order to understand the relative acidity of these sites, it is instructive to analyze the contributions of the vertical energy gaps ∆E AH and the reorganization free energies λ, to the overall reaction free energies ∆A, which are summarized in Table IV. The vertical energy gap measures how difficult it is to deprotonate the site without changing the structure (bond lengths, solvation) around it, while λ = ∆E AH -∆A quantifies the stabilizing effect of structural relaxation (including solvation) [START_REF] Hummer | Calculation of free energy differences from computer simulations of initial and final states[END_REF].

For each reaction, all these quantities include an identical contribution of the protonation of a water molecule to form a hydronium, so that differences between them really reflect what happens on each site.

The deprotonation of SiOH is stabilized by the arrival of two water molecules donating H-bonds, increasing the coordination number of the oxygen atom from approximately 3

(1 Si, 1 bonded H and 1 accepted H-bond) to 4 (1 Si and 3 accepted H-bond), without change in the distance of the first neighbours R O•••H (see Table III). This allows an efficient stabilization of the deprotonated site (large λ) which contributes favorably to the reaction.

Other possible contributions to the stabilization include the charge redistribution and the relaxation of the Si-O bond length. As can be seen in Table V, there is no significant change in the partial charge of the Si and O upon deprotonation, while to the other hydroxyl present on the same Al before deprotonation. Nevertheless, there is an asymmetry between these two sites from the point of view of solvation: Whereas AlOH on the neutral surface is coordinated by two surface water molecules, there is only one around the deprotonated AlOH 2 site. This arises from a crowding effect, as no more than three molecules can approach the ligands of the Al: This steric constraint implies that only one molecule can stabilize the deprotonated site, since two others already donate H-bonds to the neighbouring OH group. Overall, the balance between these two effects is hard to predict and the reaction free energy is only 0.05 eV larger from that for SiOH, leading to a moderate pK a difference.

In order to evaluate the integral (2), the computationally more efficient end point approximation is not sufficiently accurate because of the non-linear response of the solvent and resort to the three-point formula (3) is necessary. As already commented in our recent work [START_REF] Cheng | Acidity of the aqueous rutile TiO2(110) surface from density functional theory based molecular dynamics[END_REF][START_REF] Cheng | Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics[END_REF], Costanzo et al. 2011, Gaigeot et al. 2012[START_REF] Mangold | Absolute pK(a) values and solvation structure of amino acids from density functional based molecular dynamics simulation[END_REF][START_REF] Sulpizi | Acidity constants from dft-based molecular dynamics simulations[END_REF], 2008, Sulpizi et al. 2012] the deprotonation event involves a highly non-linear response of the solvent reorganization. An easy to monitor indicator of deviations from linearity is the dependence of the variance of the energy gap fluctuations on protonation state. Comparing the standard deviation of the gap distribution in the protonated and deprotonated trajectory, denoted by σ AH and σ A-respectively in Table IV, we notice that σ A-is larger than σ AH . An additional observation is that if σ AH is the same for all the systems investigated, σ A-, however, exhibits some variations and we therefore can not count on cancellation of corrections to linear response when comparing the different systems.

C. Comparison with previous estimates

As mentioned previously, direct comparison with experimental data is not possible, since such measurements for individual sites on a given clay edge are to date not available. However we can compare our data with previous estimates obtained with empirical methods based on bond-valence models, as will be discussed below, and previous DFT-based estimates of the deprotonation enthalpies ∆H vacuum dep for the same sites in vacuum by [START_REF] Churakov | Ab initio study of sorption on pyrophyllite: structure and acidity of the edge sites[END_REF]]. In the latter study, only one value, for AlOH 2 , was directly calculated from the relaxed deprotonated surface while the others were estimated from an empirical relation involving the partial charge of the oxygen atom (also determined from DFT calculations).

These estimates concluded that ∆H vacuum dep was most negative for AlOH 2 , followed by SiOH and AlOH. We find that the vertical energy gaps on the solvated surface are consistent with this order, i.e. that it is easier to deprotonate AlOH 2 , followed by SiOH and AlOH, for configurations of the interfacial water corresponding to the protonated sites. However, the relaxation of the surface water molecules also contributes to the reaction enthalpy (and so does the formation of a solvated hydronium ion, which is common to the three deproto-nation reactions). Overall, we measure reaction enthalpies ∆H solv dep = E A --E AH , with the energies corresponding to averages of the deprotonated (resp. protonated) state along the deprotonated (resp. protonated) trajectory, for the reaction between weak or very weak acids with water, which are smaller than the corresponding statistical error of ≈ 0.05 eV.

It thus seems that the agreement between the deprotonation enthalpy in vacuum and the deprotonation free energy is somewhat fortuitous. In a later study of water confined between pyrophyllite edges, Churakov reported transient proton exchange between SiOH or AlOH 2 on the one hand and AlOH on the other hand, via surface water molecules [START_REF] Churakov | Structure and dynamics of the water films confined between edges of pyrophyllite: a first principle study[END_REF]].

While we did not observe such rare events during the course of our simulations, possibly due to subtle changes in the organization of water molecule in a twice larger water slab in our case, this suggests that the acidity of these sites is comparable and indeed larger than that of AlOH.

Most previous acidity estimates for edge sites of clays follow from bond-valence based empirical models such as MUSIC [START_REF] Hiemstra | Intrinsic proton affinity of reactive surface groups of metal (hydr)oxides: The bond valence principle[END_REF]]. One of the major appeals of such models is their ability to obtain estimates of pK a for oxides or hydroxides in a very simple way. The main assumptions underlying the MUSIC model are (i) that the pK a of an acid can be estimated from the structure of the corresponding base (ii) an empirical linear relation between the pK a and the valence of this oxygen and its neighbours, fitted on experimental pK a 's of acids in homogenous solution, which implicitly includes information on the protonated state (iii) estimates for the coordination number and contribution of each neighbour to the valence of the oxygen. The application of this method from solutes to oxide surfaces, in particular clay edges, has already been considered. Tournassat et al. have applied this strategy using the original method, which makes use of a priori values for the coordination numbers and valence contribution of bonded H atoms, received H-bonds and metals [START_REF] Tournassat | The titration of clay minerals ii. structure-based model and implications for clay reactivity[END_REF]]. The latter depends on the metal-oxygen distance and Bickmore et al. have demonstrated that surface relaxation after deprotonation (in vacuum) may contribute to the change in valence, hence in pK a [START_REF] Bickmore | Ab initio determination of edge surface structures for dioctahedral 2:1 phyllosilicates: implications for acid-base for reactivity[END_REF][START_REF] Bickmore | Bond-valence methods for pKa prediction: critical reanalysis and a new approach1[END_REF]. Moreover, it was subsequently shown that the MUSIC model could be improved by incorporating structural information (bond length and coordination numbers) from DFT calculations [START_REF] Bickmore | Bond-valence methods for pKa prediction. ii. bond-valence, electrostatic, molecular geometry, and solvation effects[END_REF]], eventually leading to the revised MUSIC model [START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF].

In Table VII we compare different pK a estimates from the MUSIC model using different values for the coordination numbers and valence units for the neighbouring H (attached to the O or received from H-bonds) and metal (see Appendix for details of the MUSIC model).

In the first version, we report the results of Tournassat et al. who used the a priori recommendations of the original MUSIC model for coordination numbers and valences [START_REF] Tournassat | The titration of clay minerals ii. structure-based model and implications for clay reactivity[END_REF]]. In model II, coordination numbers from the DFT simulations are used. Finally, in the revised version (model III), both coordination numbers and valences are determined from the DFT simulations [START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF]]. The estimates of the original model are in remarkable agreement with that of the present work: The value for AlOH is almost the same as ours, out of the water domain; that for SiOH and AlOH 2 are on the alkaline side, with a pK a for the latter approximately one unit larger than for the former, although both are larger than ours. Interestingly, for the modelling of titration curves, which was the main purpose of Tournassat et al. , a smaller value than the one predicted by MUSIC was used for SiOH (8.2 instead of 9.1) to obtain a better agreement with the experiments. This suggests that our prediction of a more acidic SiOH group might reflect more accurately the experimental situation.

As can be seen in Table VII, the actual coordination numbers are not very different from the a priori estimates, so that the prediction of model II are still reasonable (although shifted in the wrong direction for SiOH and AlOH 2 ). Introducting the valences from DFT, computed as explained in Ref. [START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF]] only worsens all predictions, yielding in particular very acidic SiOH and AlOH sites. From this comparison we conclude that the original method provides the closest agreement with the prediction of thermodynamic integration, but this agreement is fortuitous, since including information on the correct structure and valence, which are the physical ingredients underlying this model, does not improve the result. This suggest that including such information on the deprotonated state only is not sufficient to capture the free energy of reaction, which accounts for the changes between the initial and final states. As discussed above, the reaction free energies for all three sites include a large contribution of the reorganization free energies and the structure of the deprotonated state might be very different from that of the protonated one, especially in the AlOH case for which deprotonation is followed by the departure of the neighbouring OH 2 group (this prediction for the intrinsic acidity of this site assumes that no proton transfer between surface groups occurs, while this would happen in practice because the other groups are more acidic). This particular case confirms that such empirical models, while very useful thanks to their simplicity, should be used with caution when applied to surfaces. The proximity of the various sites results in correlations between their solvation shell, especially when several sites share the same metallic ion (Al in the present case).

D. Comparison with other surfaces

It is interesting to compare our findings with previous calculations on simple oxides. For example, for the (0001) surface of quartz some of us have found 2 types of silanols, "outof-plane" and "in-plane", with different acidity constants, namely 5.6 and 8.5 [START_REF] Sulpizi | The silica-water interface: how the silanols determine the surface acidity and modulate the water properties[END_REF]]. Here the silanols are preferentially out-of-plane. As can bee seen in Table III and Fig. 2, they donate a H-bond to water (n H•••Ow = 0.92) and accept H-bonds mainly from water H (0.80 H-bonds) and occasionally from nearby surface groups (0.26 H-bond from AlOH 2 ). Such strong H-bonds donated by silanols to water have also been observed on the (0001) quartz surface [START_REF] Gaigeot | Oxide/water interfaces: how the surface chemistry modifies interfacial water properties[END_REF], Sulpizi et al. 2012]. The higher pK a of the silanol groups on the (010) edges of pyrophyllite with respect to the out-of-plane silanols on the quartz (0001) surface could be due to the nearby surface groups, and in particular the less strong H-bond accepted from AlOH 2 .

We can also compare the results for the aluminol acidity with previous calculations on the (0001) surface of alumina. The main difference in this case is in the coordination of the oxygen with the metal. Indeed in the present case, the OH group in AlOH is singly coordinated by aluminum, while on the (0001) surface of alumina the OH group of the Al 2 OH aluminol is doubly coordinated by Al. This difference also has a quite strong impact on the solvation structure of the the respective OH groups. In the case of Al 2 OH aluminols, the OH donates a not so strong H-bond to water and accepts a quite strong one from another water molecule, while on the (010) edge of pyrophyllite AlOH groups do not donate any H-bond to water and only accept one. Overall, the stronger interaction of this surface hydroxyl group with water in alumina correlates with a smaller pK a of 16.6 [START_REF] Gaigeot | Oxide/water interfaces: how the surface chemistry modifies interfacial water properties[END_REF]], compared to 22.1 in the present case.

IV. CONCLUSION

We have determined the acidity of the SiOH, AlOH and AlOH 2 groups present on the (010) edge of pyrophyllite, using DFT-based molecular dynamics simulations combined with thermodynamic integration. We have found that the most acidic group is SiOH, closely followed by AlOH 2 and that AlOH should not deprotonate in water. We have shown that the deprotonation of each site is followed by moderate changes in bond lengths and charge distributions, but sometimes more consequent changes in solvation. The solvation of surface sites and hence their acidity is strongly influenced by the proximity of other sites, in particular for AlOH and AlOH 2 which share the same Al. Introducing the DFT-derived bond lengths and valence units in the MUSIC worsens its predictions, which are surprisingly (but fortuituously, since it relies on incorrect assumptions in this case) in reasonable agreement with our estimates. The pK a values suggest that around pH=7 the sorption of ions on this face should occur only on the SiO -and deprotonated AlOH 2 which, considering the presence of the other hydroxyl group, can be seen as Al(OH) - 2 , sites. Several natural extensions of this work include the study of additional sites, e.g. SiOH + 2 , AlOH + 2 or protonated bridging SiOAl sites, and different faces of pyrophyllite, in particular the (110) one, and of substituted clays such as montmorillonite. Indeed, the presence of substitutions in the tetrahedral (Si 4+ → Al 3+ ) and octahedral (Al 3+ → Mg 2+ ) layers influences the affinity for water [START_REF] Liu | Atomic-scale structures of interfaces between phyllosilicate edges and water[END_REF]] and the acidity of surface groups. Another direct application could be the use of the pK a values to predict titration curves using e.g. the analytical approach of Bourg et al. [START_REF] Bourg | Modeling the acid-base surface chemistry of montmorillonite[END_REF] or the simulation strategy of Delhorme et al. [START_REF] Delhorme | Acid-base properties of 2:1 clays. i. modeling the role of electrostatics[END_REF]. This would provide an indirect way to confront our results to experimental data. A direct test of our predictions could eventually come in the near future thanks to the rapid development of surface sensitive experimental probing the solid-liquid interface, such as Sum Frequency Generation (SFG) spectroscopy [START_REF]Liquid interfaces probed by second-harmonic and sum-frequency spectroscopy[END_REF][START_REF] Shen | Surface properties probed by 2nd harmonic and sum-frequency generation[END_REF]].

The present study opens the way to the accurate simulation of the sorption of ions onto edge sites of clay minerals. Not only does it allow to predict the protonation state of the surface, which is a key factor. It also provides insights for the development of accurate force fields for classical molecular dynamics simulations. For example, the fact that after deprotonation of AlOH 2 both hydroxyl groups are identical suggest that polarizable ion models, which have been successfully developped for the description of molten and crystalline oxides and hydroxides [START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF]] might be well suited to also describe edge sites of aluminosilicate minerals. Resort to classical molecular dynamics will then allow to sample configurations over sufficiently long time scales to investigate the sorption/desorption processes.

Site

  pseudopotentials and DZVP-MOLOPT [VandeVondele and Hutter 2007] basis sets with a cut-off of 280 Ry. Only the Γ point of the Brillouin zone was considered and simulations were performed at 350 K, which is the usual condition to simulate water under ambient conditions with the PBE functional. Both the force field and DFT based molecular dynamics simulations were performed with the CP2K package [CP2K developers group, http://cp2k.org

  , the clay plane is tilted with respect to the directions of the simulation box. This tilt allows the bottom and top surfaces of the clay sheet to be placed relative to each other (taking into account the periodic boundary conditions) in accordance with the geometry of the bulk clay unit cell. This orientation of the clay is not frozen and fluctuates during the simulations, as a result of the relatively small size of the clay particle considered here. However the solid-liquid interface always remains perpendicular to the z direction and the two surfaces are equivalent.Three different deprotonable groups are present on the (010) edge, namely SiOH, AlOH 2 , and AlOH as can bee seen on the insert of Fig.1. Their solvation can be analyzed in terms of the coordination numbers of the O and H surface atoms by water (and other surface atoms) and the distance to the corresponding neighbours. The results are summarized in

  FIG. 1: Simulated system: The (010) pyrophyllite surface faces a slab of liquid water. The insert illustrates the three types of deprotonable sites present on this surface: SiOH, AlOH and AlOH 2 . The latter share the same aluminum atom. Si atoms are represented in yellow, Al in green, O in red and H in white. The simulation box is also indicated, and periodic boundary conditions in all directions are used.

FIG. 3 :

 3 FIG. 2: Schematic view of the hydrogen bond network at the neutral surface. Average numbers of hydrogen bonds for each type of bond are indicated.

  The MUSIC method provides pK a estimates based on the bond valence of the deprotonated site, using the coordination numbers (h.b.) and valence (val., s) of donor (d) and acceptor (a) hydrogens and of the metallic ion (M e -O). In the original version (model I), coordination numbers and valences are estimated a priori, as in Ref.[START_REF] Tournassat | The titration of clay minerals ii. structure-based model and implications for clay reactivity[END_REF]]. In model II, coordination numbers from the DFT simulations are used. In the revised version[START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF]] (model III), both coordination numbers and valences are determined from the DFT simulations. Details of the MUSIC calculations can be found in the Appendix. Results from TableIVare also reported for convenience.

FIG. 4 :

 4 FIG.4: When the AlOH site deprotonates, the adjacent OH 2 group departs from the aluminum, which remains 5-fold coordinated, and donates a H-bond (blue line) to the AlO -group. Two additional water molecules (not shown) donate H-bonds to this oxygen before and after deprotonation.

Table III

 III

and the hydrogen bond patterns are depicted in Fig.

2

. In particular, silanol groups (SiOH) donate strong H-bonds and accept weaker ones, while AlOH 2 groups only donate H-bonds and AlOH groups only accept them. Moreover, some H-bonds are formed between surface groups, e.g. between AlOH 2 and SiOH. Meanwhile, AlOH groups do not exchange H-bonds with other surface groups, accept ∼1.8 H-bonds from water and donate none. Finally, AlOH 2 groups do not accept any H-bonds, donate ∼0.65 to water and ∼0.25

  Table VI indicates that the shortening (from 1.64±0.03 to 1.60±0.03 Å) of the Si-O bond is very limited.The deprotonation of AlOH involves a larger vertical gap and a smaller reorganization free energy, compared to SiOH. This means that this site is both intrinsically more difficult to deprotonate and less stabilized by subsequent reorganization of the surface. When AlOH deprotonates, the coordination number of the oxygen hardly changes, however the nature of the neighbours changes, from 1 Al, 1 bonded H and 2 accepted H-bonds to 1 Al and 3 accepted H-bonds. One also observes a shortening of the distance to the donating water molecules (R O•••H decreases from 1.65 to 1.60 Å) and of the Al-O bond (from 1.86 to 1.73 Å) as well as a moderate charge redistribution (with an increase in the partial charge of O from -0.9 e to -1.1 e). The most spectacular change, however, is the departure of the OH 2 group adjacent to AlOH as a water molecule, with a resulting decrease in the coordination of Al from 6 to 5, as can be seen in Fig.4. Such a phenomenon was recently investigated on the neutral surface by Liu et al. , who found a free energy cost of 3 kcal/mol (≈ 5 k B T ) for the departure of the water molecule (hence a relatively low probability for this event on the neutral surface)[START_REF] Liu | Atomic-scale structures of interfaces between phyllosilicate edges and water[END_REF]]. It is this departure which opens the way to the third donating water molecule to stabilize the deprotonated site. Nevertheless, all these structural reorganizations are not sufficient to counterbalance the effect of an otherwise larger vertical gap, leading to a larger pK a .

Finally, for AlOH 2 the two terms play in opposite directions. While the vertical gap ∆E AH is smaller than for SiOH, possibly thanks to the presence of a remaining H to stabilize the charge on the deprotonated site, the reorganization is less efficient (smaller λ). In that case, the coordination number of oxygen decreases slightly from 3 (1 Al and 2 bonded H) to approximately 2.75 (1 Al, 1 bonded H and 0.75 accepted H-bond). The change is limited due to steric constraints, as no additional molecule can approach this site. The final situation corresponds to a transfer of the proton from AlOH 2 to a water molecule that was receiving a H-bond from this site. This molecule then transferred one of its initial H atoms to another molecule, according to the Grotthus mechanism, leading to a hydronium ion further away from the surface and a deprotonated site stabilized by accepting a H-bond from the surface water whose dipole is now oriented toward the surface. There is no change in the partial charges and a decrease in the Al-O bond length (from 2.05 to 1.85 Å). The resulting deprotonated group is in fact identical (in terms of geometry and partial charges)

TABLE I :

 I Simulated time (in ps) for each value of the coupling parameter η and each site.

		Simulation time (ps)
		η = 0.0 η = 0.5 η = 1.0
	SiOH	17.0	22.1	18.1
	AlOH	8.6	8.7	11.4
	AlOH 2 19.1	17.6	17.7
	Site	n d d 0 k d n α	α 0	k α
	AlOH 1 0.98 0.2 1	116 0.2
	SiOH 1 1.01 0.2 1	118 0.2
	AlOH 2 1 1.03 0.2 2 109/121 0.2
	H 3 O + 3 1.00 0.2 2	111 0.2

TABLE II :

 II List of the parameters employed in the harmonic potential Eq.

TABLE IV :

 IV Site ∆E AH σ AH ∆E A -σ A -∆E 0.5 σ 0.5 ∆A pK a λ SiOH 7.57 1.03 -6.42 1.50 0.61 0.60 0.60 6.8±0.4 6.97 AlOH 7.97 1.02 -6.95 2.05 2.02 0.58 1.52 22.1±1.0 6.45 AlOH 2 7.22 1.05 -6.15 2.24 0.71 0.49 0.65 7.6±1.3 6.57 Average vertical deprotonation gap ∆E for the protonated (AH), deprotonated (A -

Acknowledgments

The authors thank Rodolphe Vuilleumier for useful discussions and acknowledge financial support from the Agence Nationale de la Recherche under grant ANR-09-SYSC-012. This work was performed using HPC resources from GENCI-CINES/IDRIS (Grants c2010086491, c2011086491 and c2012086491). This work was carried out under the HPC-EUROPA2 project (project number: 228398) with the support of the European Commission Capacities Area -Research Infrastructures Initiative. S. T. is grateful to Alexey Cheptsov for technical assistance during his visit.

Appendix: Acidity estimates with the MUSIC model

The MUSIC model allows to estimate the pK a of a site from the coordination number of the oxygen by hydrogen atoms (m H on the site, n H from other molecules), their valencies (s H for the former, 1 -s H for the latter) and the metal valencies s M e i as:

In this expression, A = 19.8 is a constant whose value is determined from regressions of experimental data for aqueous acids, V = -2 is the oxygen valency, while the metal valencies are determined from the corresponding metal-oxygen distances r M eO i as exp((r 0i -r M eO i )/b) with b = 0.37 Å and r 0i reference bond lengths calculated from bulk crystal structures. In this work, we compare our results from thermodynamic integration with three variants of the MUSIC model where the coordination numbers and valencies are determined according to different prescriptions (see text). The coordination numbers and valencies reported in Tab. VII from DFT calculations are computed using a cut-off of r c = 2.375 Å as [START_REF] Machesky | Surface protonation at the rutile (110) interface: Explicit incorporation of solvation structure within the refined music model framework[END_REF]]:

with n the coordination number and