Keywords: Swelling clay, Water adsorption, Atomistic modeling, Grand canonical, Na-montmorillonite

Swelling clays are known to exhibit crystalline swelling at the layer scale (∼ nm): X-Ray Diffraction (XRD) identifies discrete jumps in inter-layer spacing upon change in relative humidity (RH) corresponding to the thickness of a water layer. Continuum theories are not adapted to describe systems at the scale of a few water molecules, and much effort has been dedicated to understand and predict crystalline swelling from molecular simulation. Grand canonical (GC) molecular simulations combined with thermodynamic stability analysis is one of the most appropriate techniques to predict swelling in function of RH from atomistic simulation. Yet, usual GC Monte Carlo is quite inefficient, and improved algorithms are needed to reach reasonable accuracy. In this work, we report results of configurational bias GC molecular simulations of Namontmorillonite, one of the most studied swelling clays. This dataset covers a wide range of relative fugacities above and below vapor saturation, and a wide range of spacings, encompassing all hydration states. It makes it possible to derive a complete hydration diagram in function of external pressure and relative fugacity. Such a diagram extends literature results to undersaturated conditions and can be used to anticipate the crystalline swelling upon any complex loading. Of particular interest is the free swelling upon change in RH, for which our predictions compare well with experimental XRD results available in the literature. The constitutive behavior of the confined water has been little investigated so far, and existing studies were limited to the undrained stiffness (at constant water content). Since Gibbs-Duhem equation does not hold for confined water, proper characterization requires 3 independent moduli, in contrast with bulk water which is characterized by the bulk modulus only. Here, we provide estimates of the 3 moduli. The constitutive behavior of confined water significantly differs from that of bulk water, but the bulk behavior is recovered in the capillary domain at spacing larger than 19 Å.

Introduction

Clays are some of the most abundant minerals in the Earth crust and the swelling of clays upon hydration is a well-known issue for a variety of applications 1,[START_REF] Bergaya | Handbook of clay science[END_REF] . At low water content, when submitted to water vapor, it has long been evidenced by X-Ray Diffraction that the swelling of clay is 'crystalline', in the sense that it occurs by sudden changes in the basal spacing at the clay layers scale (~nm) [START_REF] Mooney | Adsorption of Water Vapor by Montmorillonite. II. Effect of Exchangeable Ions and Lattice Swelling as Measured by X-Ray Diffraction[END_REF][START_REF] Norrish | Crystalline Swelling of Montmorillonite: Manner of Swelling of Montmorillonite[END_REF] . The c-axis spacing can be accurately measured by XRD while controlling the water activity either through the relative humidity (RH) of the vapor [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites; 2, The Li+ , Na+ , K+ , Rb+ and Cs+ -exchanged forms[END_REF] or through the ionic strength of a liquid solution 6,7 . The spacing is found to increase step-by-step: for instance, the spacing of Na montmorillonite (Na-mmt) increases from ~9.5 nm (dry state) to ~12.5 nm (water content around 0.05-0.15 g water per g clay), to ~15.5 nm (0.2-0.3 g water / g clay), and to ~18.5 nm (>0.4 g water / g clay). Interestingly, the differences in basal spacing (~3 Å) corresponds to the size of a water molecule so that one can readily identify those hydration states as the formation of water layers inbetween clay mineral layers. For this reason, the hydration states are referred to as 'xW' states with 'x' the number of water layers: 0W is the dry state (~9.5 nm), 1W the one water layer state (~12.5 nm) etc.. Therefore, crystalline swelling appears to be the consequence of hydration transitions between those xW states upon a change of water activity.

The RH of hydration transition strongly depends on the nature of the clay, in particular the layer charge and its localization in the mineral layers (surface vs. center), and on the nature of the counter-ion in the interlayer space. It is most often inhomogeneous within a sample since clay mineral compositions are never perfectly homogeneous [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites; 2, The Li+ , Na+ , K+ , Rb+ and Cs+ -exchanged forms[END_REF] . Detailed analysis of XRD results suggests the coexistence of different hydration states at the same RH 8,9 , and significant hysteresis between adsorption and desorption can be observed [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites; 2, The Li+ , Na+ , K+ , Rb+ and Cs+ -exchanged forms[END_REF]7 . While XRD provides detailed information about the deformation at the clay layer scale, measuring the complementary thermodynamic quantity, namely stress, is much more challenging. Surface force apparatus are limited to easily isolated surfaces (e.g., mica) 10 , but is not easily applicable in general. More macroscopic measurements of stress and stiffness are often limited to the domain of osmotic swelling (much larger water content) 11 , and suffer from the difficulty to relate macroscopic measurements to nanometric properties 12,13 . Molecular simulation offers an interesting alternative to experiments to address the elementary mechanisms of swelling at the nanoscale. The very first studies of hydrated clays dates back to the 1990's 14,15 , and, since then, molecular simulation has been widely used to investigate all sorts of physical properties at the layer scale from structure analysis to electrolyte diffusion, vibrational properties and mechanics [START_REF] Greathouse | Developments in Clay Science[END_REF] . Regarding the hydro-mechanical couplings at the origin of crystalline swelling, early works by Molecular Dynamics (MD) were able to recover the experimental relationship between water content and basal spacing with plateaus in basal spacing indicative of a crystalline swelling [START_REF] Boek | Molecular Modeling of Clay Hydration: A Study of Hysteresis Loops in the Swelling Curves of Sodium Montmorillonites[END_REF][START_REF] Boek | Monte Carlo Molecular Modeling Studies of Hydrated Li-, Na-, and K-Smectites: Understanding the Role of Potassium as a Clay Swelling Inhibitor[END_REF] . Investigation of the energetics of water adsorption shows oscillations with minima at the basal spacing of the plateaus which further confirms the crystalline nature of swelling [START_REF] Smith | Molecular Computer Simulations of the Swelling Properties and Interlayer Structure of Cesium Montmorillonite[END_REF] . Yet, usual MD simulates closed systems, i.e., at fixed water content, and therefore cannot provide information on the propensity of a clay to swell in function of the RH, which is the key information for practical application. To this end, one has to consider grand canonical (GC) simulations, that is, simulations in which the water content can fluctuate and one controls the chemical potential instead (i.e., the RH). GC simulations are quite inefficient due to extremely low acceptation rates of particle insertion and deletion, and it appears necessary to adopt peculiar GC algorithms to reach reasonable accuracy (e.g., configurational bias Monte Carlo). For this reason, GC studies of clay swelling are much less frequent than MD studies. First GC studies were not accurate enough to identify the crystalline swelling [START_REF] Delville | Structure of liquids at a solid interface: an application to the swelling of clay by water[END_REF][START_REF] Delville | Adsorption of vapor at a solid interface: a molecular model of clay wetting[END_REF] , but later studies did observe oscillations in the solvation pressure isotherms as a function of basal spacing which confirms the existence of the discrete xW hydration states [START_REF] Karaborni | The Swelling of Clays: Molecular Simulations of the Hydration of Montmorillonite[END_REF] . Early GC studies were limited to a single basal spacing or a single water chemical potential, usually corresponding to liquid water in standard conditions, and thus were not able to estimate how basal spacing evolves with RH [START_REF] Delville | Structure of liquids at a solid interface: an application to the swelling of clay by water[END_REF][START_REF] Delville | Adsorption of vapor at a solid interface: a molecular model of clay wetting[END_REF][START_REF] Karaborni | The Swelling of Clays: Molecular Simulations of the Hydration of Montmorillonite[END_REF][START_REF] Shroll | Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling[END_REF][START_REF] Chávez-Páez | Monte Carlo simulations of Wyoming sodium montmorillonite hydrates[END_REF][START_REF] Hensen | Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation[END_REF][START_REF] Smith | Molecular simulations of hydration and swelling in clay minerals[END_REF][START_REF] Whitley | Free energy, energy, and entropy of swelling in Cs-, Na-, and Sr-montmorillonite clays[END_REF] . Combining GC and isobaric simulations, Hensen et al. [START_REF] Hensen | Why Clays Swell[END_REF] first related clay swelling and RH from molecular simulation. Large domains of meta-stability appear in the swelling curves, indicative of the swelling hysteresis. Such combination of GC and isobaric simulations has been used to investigate the role of the counter-ion, and of the nature of the clay [START_REF] Teich-Mcgoldrick | Swelling Properties of Montmorillonite and Beidellite Clay Minerals from Molecular Simulation: Comparison of Temperature, Interlayer Cation, and Charge Location Effects[END_REF] . Hydration transitions are peculiar types of phase transition and disentangling stable and metastable states requires proper thermodynamic analysis by integrating the solvation pressure isotherm with respect to basal spacing to derive the so-called swelling free energy [START_REF] Shroll | Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling[END_REF] . Integration of the solvation pressure provides the thermodynamic potential minimum at equilibrium, so that stability and meta-stability can be identified as global or local minima, respectively. Note that such thermodynamic integration was also applied to MD results [START_REF] Young | Simulations of Clay Mineral Swelling and Hydration: Dependence upon Interlayer Ion Size and Charge[END_REF] (fixed water content) and the corresponding thermodynamic potential is the Gibbs free energy which characterizes stability at fixed water content, instead of fixed RH.

The advantage of constant volume GC simulation with thermodynamic analysis 23 over constant pressure GC simulation [START_REF] Hensen | Why Clays Swell[END_REF] is that one can predict stable and metastable states at any external pressure from a single solvation pressure isotherm. This route was first applied to GC results of hydrated clays by Shroll et al. [START_REF] Shroll | Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling[END_REF] to predict how stable basal spacing evolves with the external pressure (dehydration of clays induced by a mechanical loading). Following the same approach, Whitley et al. [START_REF] Whitley | Free energy, energy, and entropy of swelling in Cs-, Na-, and Sr-montmorillonite clays[END_REF] addressed the effect of the nature of counter-ion in the interlayer space and Smith et al. [START_REF] Smith | Molecular simulations of hydration and swelling in clay minerals[END_REF] investigated the effect of layer charge, but these studies remain limited to a single chemical potential (liquid water at ambient conditions).

Predicting how stable and metastable basal spacings evolves with RH requires many GC simulations varying both chemical potential and basal spacing. Such complete picture was first proposed by Tambach et al. [START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF] , focusing on the 1W-2W transition over 5 different RH.

This approach provides a systematic methodology to estimate the crystalline swelling quantitatively from molecular simulation, and address issues of practical interest such as the role of counter-ion, layer charge and nature of clay [START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF][START_REF] Yang | Layer Charge Effects on Adsorption and Diffusion of Water and Ions in Interlayers and on External Surfaces of Montmorillonite[END_REF][START_REF] Li | Molecular Simulation Study of Montmorillonite in Contact with Water[END_REF] or the hysteresis 34 . The same route can be used to estimate the complete xW hydration phase diagrams in function of chemical potential, external pressure and temperature 35 . Recent investigations have addressed the case of inhomogeneous hydration states, the importance of solid flexibility, and the importance of force field parameterization [36][37][38][39][40] .

While swelling and hydration transitions have been the focus of many molecular simulation studies, the mechanical stiffness of the hydrated clays have attracted much less attention.

Existing studies investigated only the undrained stiffness ('undrained' means that water content is hold constant, by opposition to 'drained' when water is free to flow out of the system to maintain a constant water chemical potential). Results show a strong anisotropy with differences of one to two orders of magnitudes between the direction dominated by the water layer stiffness and that dominated by the mineral layer stiffness 13,41 . Experimental measurements appears consistent with the 'weak' (water dominated) directions 13 , which suggests that those directions are of particular importance with respect to the macroscopic behavior.

The undrained stiffness is not enough to fully characterize the mechanical behavior of an hydrated phase. As we shall explain in section 'Discussion/Confined water properties', the adsorbed water no more satisfies the Gibbs-Duhem equation, and, as a consequence, proper mechanical characterization involves three independent moduli (including the usual rigidity with respect to deformation), whereas bulk water is characterized by its bulk modulus only 42 . Recently, we proposed a first rough estimation of all three moduli 43 based on previous simulation results 35 . These estimations are limited to water activities above 1 corresponding to saturated conditions typical of geomechanics, and, to the best of our knowledge, no estimation exists in the crystalline swelling domain.

In this paper, we perform GC simulation of Na-mmt varying both basal spacing and water activity at 300K, and we follow the thermodynamic route [START_REF] Shroll | Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling[END_REF] to determine the domains of stability and meta-stability of each hydration phase. To the best of our knowledge, we report in this paper the largest dataset of GC results of Na-mmt covering 21 different water activities both below (17) and above (4) saturation with c-axis spacing between 9.25 Å and 23.25 Å. From these results, we derive the stable and metastable swellings and water contents as function of RH in unjacketed conditions ('unjacketed' is the term used in mechanics to refer to a porous medium immersed in a fluid so that the mechanical stress in the solid reduces to the bulk pressure of the fluid, what is commonly known as a free swelling experiment). The predicted free swelling agrees well with existing XRD experimental studies. The large dataset of GC results makes it possible to draw a complete phase diagram of the hydration states, both below and above saturation, thus extending that proposed by Honorio et al. 35 (above saturation only). Moreover, we provide the first detailed estimates of the three moduli characterizing the mechanics of the 1W and 2W hydration states both below and above saturation. Nevertheless, current accuracy of our result does not allow reliable estimations for the 3W state and for RH below 0.2 (cavitation of the bulk phase). The raw results of our GC simulation (pressure and water content isotherms) are made available in a public repository (https://doi.org/10.5281/zenodo.4633697).

In section 'Methods', we provide all the details of the GC molecular simulation we perform. In section 'Results', we present the direct results of the molecular simulation (pressure and water content). In section 'Discussion/Hydration transitions', we perform the thermodynamic analysis to determine the domains of (meta)stability of the different hydrated states. The unjacketed swelling and water content are derived, as well as the hydration phase diagram. In section 'Discussion/Confined water properties', we provide the estimates of the three moduli characterizing the behavior of the hydrated states.

Methods

The most common clays studied by molecular simulations to investigate swelling are Arizona

montmorillonite (M x [Si 8 ][Al 4-x Mg x ]O 20 (OH) 4 ), Wyoming montmorillonite (M x+y [Si 8-y Al y ] [Al 4-x Mg x ]O 20 (OH) 4 ) and beidellite (M y [Si 8-y Al y ][Al 4 ]O 20 (OH) 4
), which are dioctahedral smectites of low charge (x + y between 0.5 and 1.25, and usually 0.75 or 1 in most studies). Beidelite contains tetrahedral substitutions only (close to the surface of the mineral layer), whereas Arizona montmorillonite contains octahedral substitutions only (at the center of the mineral layer), and Wyoming montmorillonite is intermediary, most often with 2/3 of octahedral substitutions and 1/3 of tetrahedral substitutions. In this work, we simulate Na-Wyoming montmorillonite, which corresponds to the intermediary situation, with x = 0.5 and y = 0.25 corresponding to an average layer charge of x + y = 0.75, and with a common counter-ion (Na). This model clay, of chemical formula Na 0.75 [Si 7.75 Al 0.25 ][Al 3.5 Mg 0.5 ]O 20 (OH) 4 , has been studied by many to investigate clay swelling 13,15,[START_REF] Boek | Molecular Modeling of Clay Hydration: A Study of Hysteresis Loops in the Swelling Curves of Sodium Montmorillonites[END_REF][START_REF] Boek | Monte Carlo Molecular Modeling Studies of Hydrated Li-, Na-, and K-Smectites: Understanding the Role of Potassium as a Clay Swelling Inhibitor[END_REF][START_REF] Chávez-Páez | Monte Carlo simulations of Wyoming sodium montmorillonite hydrates[END_REF][START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF]34 For the simulated systems, the unit cell is replicated 4 and 2 times along the a and b axis, respectively; so that the lengths of the simulated systems in the x and y directions are 20.64 Å and 17.93 Å, respectively. For basal spacings below 17 Å, the unit cell was replicated twice along the c-axis, i.e., two layers are simulated, whereas for larger basal spacings a single unit cell was considered in the c-axis. Considering 2 layers at small basal spacings is critical to avoid size effects [START_REF] Chávez-Páez | Monte Carlo simulations of Wyoming sodium montmorillonite hydrates[END_REF] . When two layers are simulated, we enforce that the two basal spacings are strictly equal by recentering each layer at its initial center of mass at each simulation step. Note that doing so we disregard inhomogeneous hydration states, in which neighboring clay layers can coexist with different basal spacings 36 . Heterogeneity within a layer is also suggested, and would be an essential aspect of the mechanisms of transition between hydration states 37,38 . Inhomogeneous hydration states are likely, at least because layer charge distribution is never perfectly homogeneous. Imposing periodic boundary conditions and maintaining equal basal spacings are therefore simplifications with respect to reality.

To simulate the hydration of the Na-mmt we use the Clay-FF force field 44 , which provides an accurate description of the electrostatic interactions and is therefore well-adapted to model the swelling (dominated by non-bonded interactions) [START_REF] Greathouse | Developments in Clay Science[END_REF] . This potential allows full atomic flexibility of the mineral layers and therefore an unambiguous computation of the stress (no frozen zone that would act as a dead volume). Accounting for the solid flexibility is obviously essential to address the deformation of the minerals when hydration is inhomogeneous 37,38 , but it appears also important for homogeneous cases when the minerals are flat 40 , as is the case in the present work. The original clayFF uses the SPC model to simulate water.

Here, we use the rigid SPC/E model instead 45 , which is more consistent with experiments regarding the mechanics of liquid water (pressure vs. density). Water SPC/E is known to underestimate notably the saturation vapor pressure/fugacity 46 (998.70 Pa vs. 3 536.8 Pa).

As we shall see, the magnitudes of the pressure isotherms in this work are in the range of hundreds of MPa, so that this inaccuracy in vapor pressure appears indistinguishable in the isotherms, and has almost no impact on the stability analysis. In contrast, it is essential to predict correctly the liquid pressure which can reach hundreds of MPa and prevails at large basal spacings. Moreover, an accurate liquid pressure is important in our estimation of the moduli describing the adsorbed water mechanics (see section 'Discussion/Confined water properties'). We apply a cutoff at 8.5 Å for the Lennard-Jones interactions. The Ewald summation is used to compute the long range electrostatic interactions, with a cutoff of 8.5 Å for the direct part, and a target relative error in forces of 10 -7 . All simulations are performed with LAMMPS software 47 (https://lammps.sandia.gov/), which is modified to include the configurational bias Monte Carlo algorithm (see hereafter).

As mentioned in the introduction, usual GC Monte Carlo simulations of water in clays suffer from a very low acceptation rate of the insertion/deletion trials, and it is necessary to use peculiar GC methods to reach reasonable computational efficiency. In this work, GC simulations are performed following a configurational bias scheme inspired from that proposed by Hensen et al. [START_REF] Hensen | Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation[END_REF] (for the MCY water model) and reused in many studies [START_REF] Hensen | Why Clays Swell[END_REF][START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF]34 (MCY and TIP4P water models). Configurational bias Monte Carlo consists in biasing the insertion/deletion trials in favor of more probable configurations, while correcting for the bias in the acceptation probability of the Monte Carlo step 48 . The scheme proposed by Hensen et al. [START_REF] Hensen | Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation[END_REF] consists in biasing first the position of the oxygen atom of water, and then the orientation (positions of the hydrogen), based on the non-electrostatic (short range) part of the MCY potential. And the electrostatic part is considered only when calculating the acceptation probability correcting the bias. Doing so, the computationally expensive part of the potential (electrostatics) is computed only once for the acceptation probability. This biasing scheme follows a similar logic as the double cut-off scheme of Vlugt et al. 49 but the split of the potential is based on the nature of the interaction rather than on a cut-off. In both cases, the bias takes into account the short range repulsion of the potential, which turns out to be the dominant factor for efficient GC trials. We reuse the scheme of Hensen et al. [START_REF] Hensen | Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation[END_REF] to perform GC simulations of the SPC/E water model. Since the hydrogens of water SPC/E include only electrostatic interactions, the application of the scheme for the orientational part reduces to choosing a random orientation (no bias), and only the trial position of the oxygen is biased. The scheme is the following:

Insertion trials 1. k O trial positions for the oxygen atoms are generated randomly: r O 1 , ..., r O k O . The non- electrostatic (Lennard-Jones) part U O LJ r O i
of the potential energy of these oxygens with the existing atoms is computed, and one of them, r O k , is chosen randomly with probability

p trial r O k = exp -βU O LJ r O k W O , with W O = i exp -βU O LJ r O i the Rosenbluth factor and β = (k B T ) -1 .
2. The positions of the hydrogens r H 1 and r H 2 are obtained by choosing randomly an orientation.

The insertion attempt is accepted with probability

p acc (N → N + 1) = min 1, V N +1 W O k O βf exp -βU H 2 O elect r O k , r H 1 , r H 2 , where U H 2 O elect r O k , r H 1 , r H 2
is the electrostatic interaction of the trial water molecule with the existing atoms, and f is the fugacity related to the chemical potential µ according to µ = kT ln f q kin (T ) , with q kin (T ) the kinetic part of the molecular partition function at temperature T 48 .

Deletion trials

1. One of the existing oxygens of water molecules (r O 1 ) is selected for trial deletion and k O -1 additional trial oxygen positions are generated randomly: r O 2 , ..., r O k O . The Rosenbluth factor for the non-electrostatic part of the potential is computed:

W O = i exp -βU O LJ r O i .
2. The positions of the hydrogens r H 1 and r H 2 are obtained by choosing randomly an orientation.

The selected water molecule is deleted with probability

p acc (N → N -1) = min 1, N V k O W O 1 βf exp βU H 2 O elect r O 1 , r H 1 , r H 2 .
(Note that there is a typo in the paper of Hensen et al. [START_REF] Hensen | Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation[END_REF] : N + 1 at the numerator should be replaced by N 48 ) One can check that this scheme ensures the super-detailed balance 48 , which guarantees that the molecular configurations are sampled according to Boltzmann statistics. We implemented this scheme in LAMMPS. The validity of the implementation was checked against usual GC Monte Carlo for bulk water (Figure 1, bottom). The computational efficiency of the GC simulations is improved by a factor of 5 for the simulation of bulk water SPC/E. The performance is even greater for GC simulations of hydrated clays, with improvement sometimes exceeding a factor of 10. Since each of our GC simulations of hydrated clays requires a few weeks of calculation (single processor), the use of the configurational bias scheme appears critical to our work. In practice, hybrid GCMC/MD simulations are performed in which 50 configurational bias GCMC insertion / deletion trials are performed every 100 MD steps. No GCMC translations or rotations are performed, the displacements of the atoms are handled by the MD. Each MC step considers k O = 100 possible positions of oxygen. The MD part of the simulation considers a velocity-Verlet integration with a timestep of 1 fs. A Langevin thermostat is applied to impose a constant temperature of 300K, with a damping parameter of 100 fs. The SHAKE algorithm is used to impose the rigidity of the water molecules. The systems are equilibrated over 1.5 × 10 6 configurational biais GCMC steps (150 × 10 6 oxygen positions), and the thermodynamic averages are produced over 7 to 20 × 10 6 configurational biais GCMC steps (700 to 2000 × 10 6 oxygen positions) depending on the computational costs of the systems (more expensive systems are computed over less steps). The accuracy of the computed thermodynamic averages are estimated by the block averaging techniques 48 , with blocks of 10 5 configurational biais GCMC steps (10 7 oxygen positions).

Before addressing the hydration of clay by GC molecular simulation, it is necessary to establish the correspondence between chemical potential µ and the bulk properties, namely vapor and liquid pressures (P v and P l ) and densities (ρ v and ρ l ) for the water SPC/E model.

Of particular interest are the properties of the superheated liquid, between saturation and cavitation. Superheated water is relevant for hydrated clays, because one expects the molecular simulation results to converge to capillary water at large basal spacings; and capillary water is superheated liquid water which is stabilized by the Laplace pressure across the menisci separating the liquid and vapor phase in unsaturated porous media. We display in Figure 1 (top) the results of a series of canonical simulations of bulk water SPC/E at 300K. The predicted density-pressure curve reproduces reasonably well the equation of state of Wagner and Pruß 50 . This equation of state is the reference one for the IAPWS (International Association for the Properties of Water and Steam). It is accurately calibrated on experimental data above saturation, but it remains theoretical in the superheated domain due to the lack of experimental data [START_REF] Shmulovich | Experimental superheating of water and aqueous solutions[END_REF] . Note that, in the simulation results, cavitation is observed at a density of 875 kg/m 3 . Systems at smaller densities are inhomogeneous with a vapor bubble in the liquid. Comparison with the IAPWS equation of state applies to the homogeneous liquid only, i.e., at densities higher than 875 kg/m 3 . The relationship between chemical potential and pressure/density can be obtained directly from GC simulation or indirectly from canonical simulations followed by a thermodynamic integration of the Gibbs-Duhem equation (dµ = dP/ρ). The results are displayed in Figure 1 (bottom). Note that we do not report our results in terms of the chemical potential µ, but in terms of the fugacity f /f sat relative to the fugacity at saturation. The relative fugacity coincides with the water activity (i.e., RH in unsaturated conditions) which is more easily interpreted than chemical potential, and fugacity is the natural control parameter in GC Monte Carlo (see acceptation probabilities above). The results of Figure 1 (bottom) exhibit good correspondence between the two approaches (direct GC simulation, and indirect canonical simulation with Gibbs-Duhem integration). Note that the thermodynamic integration requires to determine first the saturation point which is obtained as the intersection between the liquid and vapor pressure curves: f sat = 998.70Pa. From this series of preliminary calculations, we can establish the list of fugacities that we will consider later on for the simulation of clay hydration. The list of fugacities and associated bulk properties are listed in Table 1. The selected fugacities correspond to prescribed RH below saturation and to prescribed bulk liquid pressures above saturation. Note that below a relative fugacity of 20%, the superheated liquid is unstable, so liquid properties cannot be reported. Note also that the water vapor behaves almost as an ideal gas at the temperature considered (300K).

Results

In GC simulations, one simulates systems at fixed fugacity (chemical potential) and layer spacing. One can fully characterize the hydro-mechanical couplings by computing the conjugated thermodynamic quantities, namely the pressure normal to the layer and the water content. We report in this section those two quantities for a variety of fugacities and spacings. The raw results are also made available in a public data repository (https://doi.org/10.5281/zenodo.4633697). Note that we report our results in terms of the z-axis spacing (d) and not of the c-axis spacing (c). The two spacings are related according to d = c sin (β), with β = 100.46 • , and they differ slightly (sin (β) ≈ 0.983). While the basal spacing is the relevant spacing with respect to crystallography, the z-axis spacing is the relevant spacing from the mechanical point of view, since it is the spacing conjugated to the normal pressure P zz . Note also that we introduce the 'inter-layer' spacing h as the spacing in-between the mineral layers accessible to water. While the z-axis and c-axis spac- (bottom) Fugacity-pressure and fugacity-density curves obtained either by direct GC simulations (configurational bias or usual Mont Carlo) or by canonical simulation and thermodynamic integration of the Gibbs-Duhem equation. Saturation corresponds to the intersection between the liquid and vapor pressures.

Table 1: Fugacity (f ) conditions considered for the study of water SPC/E adsorption in Na-mmt at 300K, and properties of the bulk vapor and liquid phases at those fugacities: vapor and liquid pressures (P v and P l ) and densities (ρ v and ρ l ) and liquid bulk modulus (K l ). The relative fugacity f /f sat with respect to the fugacity at saturation (f sat ) quantifies the water activity which coincides with the concept of relative humidity in the undersaturated domain (f ≤ f sat ). Note that the vapor phase properties are almost that of an ideal gas. MPa). A Savitzky-Golay filter is applied to smoothen the scatter of the results (lines). In what follows, we systematically use the filtered curves (stability analysis, hydration diagram, confined water properties). For sake of readability, the results of Figure 2 cover only about one relative fugacity out of two. We display in Figure 3 the filtered curves at all the relative fugacities considered, with a specific highlight on the convergence to the bulk liquid pressure at large spacings.

f (Pa) f /fsat (%) Pv (Pa) ρv (g/m 3 ) P l (MPa) ρ l (kg/m 3 ) K l (MPa) 0 0 0 0 - - - 0.
Consistently with previous GC studies of hydrated clay, we observe that the pressure isotherms oscillate with the spacing, the decreasing branches corresponding to the xW (meta)-stable hydration states. While the 0W (dry) and 1W hydration states are observed at all fugacities, the 2W state appears only for relative fugacities f /f sat larger than 1%. As for the 3W hydration state, none of the GC studies mentioned in the introduction addressed it.

Not only it requires exploring large spacings, but it also requires a high accuracy in pressure because the associated oscillation in the isotherm does not exceed a few tens of MPa. Here, the 3W state is clearly visible at spacings around 18Å and for relative fugacities of 100% or above. The fact that the 3W state appears between 90% and 100% RH is consistent with experiment (2W-3W transition observed only above 90% of RH 6,7 ). At spacings larger than about 19 Å, no more variation in the isotherms are visible, and the corresponding normal pressures are almost equal to the bulk water pressures. More precisely, at relative fugacities above bulk cavitation (30% or higher), we recover the bulk liquid pressure (dotted lines in Fig. 3 bottom), and at relative fugacities below bulk cavitation (20% or smaller) we recover the bulk vapor pressure (i.e., almost 0 at the scale of the figure). This is consistent with the usual theory of capillary water: for an hydrophilic porous medium, one expects the smallest pore to be filled by the superheated liquid in unsaturated conditions. Applying Kelvin's law and assuming perfect wetting, the pore diameter below which capillary water is expected at a relative fugacity of 30% is about 9Å. Taking into account the thickness of the adsorbed layers, this is well above the largest inter-layer spacing considered in this study. Note however that recent studies suggest that the surface of montmorillonite can turn hydrophobic if the surface is devoid of counter-ions [START_REF] Szczerba | Structure and Dynamics of Water-Smectite Interfaces: Hydrogen Bonding and the Origin of the Sharp O-D W /O-H W Infrared Band From Molecular Simulations[END_REF][START_REF] Szczerba | Intrinsic hydrophobicity of smectite basal surfaces quantitatively probed by molecular dynamics simulations[END_REF] , and, in such cases, no capillary regime would be expected. In the present simulations, counter-ions are present at the surface, so that a full surface wetting it is expected along with a capillary regime. Our results, however, are not accurate enough to identify an osmotic effect (double layer theory) that would induce a deviation from the capillary pressure [START_REF] Israelachvili | Intermolecular & Surface Forces[END_REF] . Osmotic effects are typically in the range of a few MPas and reaching an accuracy of 1 MPa requires peculiar molecular simulation set up that do not use GC steps (see [START_REF] Yang | Swelling pressure of montmorillonite with multiple water layers at elevated temperatures and water pressures: A molecular dynamics study[END_REF] for instance, that addresses the 3W hydration and osmotic effects by mean of such set up). Water density ρ In Figure 4, we display the water density in the pore as a function of spacing and relative fugacity. Again, the points are the raw results, and the lines are obtained by application of the Savitzky-Golay filter. The typical accuracy of the results is ±3kg/m 3 (block averaging estimate). Results near the bulk cavitation (f /f sat = 20%) are typically the least accurate with a typical error reaching ±10kg/m 3 . In Figure 4, only about one relative fugacity out of two is displayed for readability. In Figure 5, the filtered curves are displayed for all the considered relative fugacities, and the domains corresponding to the different hydration states are highlighted with colored areas. We provide a special highlight on the domain of liquid-like densities in Figure 5 (bottom) which shows that the confined water density does converge with respect to the spacing d, i.e., a peak in the density curve. This is indeed what is observed here. These oscillations of density also emerge from the structuration of the confined water in layers, which has been widely investigated in the literature [START_REF] Greathouse | Developments in Clay Science[END_REF] . One can note that the peak of the 1W state exhibits an inhomogeneous evolution with relative fugacity. For basal spacings below 12 Å the water molecules are more strongly adsorbed than for basal spacings above 12 Å, so that desorption occurs at lower fugacities. As for the pressure isotherms, a clear transition appears for relative fugacities between 20% and 30%, that corresponds to the cavitation of the superheated bulk liquid. Above cavitation, the system converges to a liquid-like form at large basal spacings (capillary domain), which is consistent with the capillary theory and Kelvin's law. Below cavitation, the superheated bulk liquid is no more stable and the system converges to a vapor-like form at large basal spacings. Yet, the densities are still much higher than that of a vapor because there remains a surface adsorption.

(kg/m 3 ) f/f sat (undersaturated) f/f sat (supersaturated) 30% 50% 70% 100% 204% 407% 0% 0.1% 1% 10%

Discussion

Hydration transitions

In this section, we perform a stability analysis of the hydrated clay, in order to predict the crystalline swelling and adsorption isotherm in unjacketed conditions. This analysis follows the same route as that used by Schroll et al. [START_REF] Shroll | Molecular dynamics simulations in the grand canonical ensemble: Application to clay mineral swelling[END_REF] and Tambach et al. [START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF] . Since we investigated a broader range of basal spacings and relative fugacities, we are able to predict all hydration transitions from 0W to capillary. Moreover this methodology is not limited to unjacketed conditions, but any external pressure can be considered, which makes it possible to draw a complete hydration diagram describing both the stable and meta-stable domains of the xW states as a function of relative fugacity and external pressure. While unjacketed conditions are the preferred conditions for laboratory tests, such hydration diagrams are relevant to the (geo-)mechanical applications in which clay is not necessarily free to swell, but is confined by surrounding soils and rocks.

Let us first recall the basics of the stability analysis. In GC simulations, the hydrated clay is under thermal and osmotic equilibrium with its environment so that one controls the temperature T (300K in this work), the chemical potential µ of water (or equivalently the relative fugacity), and the spacing (d). The thermodynamic equilibrium of this system corresponds to the minimization of an hybrid potential λ which is the combination of the Helmholtz free energy for the solid mineral and of the Grand potential for the water:

λ = min N H 2 O F A -µ N H 2 O A
with F the total Helmholtz free energy. The energy balance of this potential upon changes in chemical potential and spacing is given by dλ| T = -P zz d(d) -

N H 2 O
A dµ. This energy balance leads to the Maxwell relation we discussed in the previous section ( ∂Pzz

∂µ d,T = ∂(N H 2 O /A) ∂d µ,T
). The thermodynamic potential as a function of the spacing is readily obtained by integration of the pressure isotherm of Fig. 3:

λ (d) -λ (d 0 ) = - d d 0 P zz d(d).
Let us now consider the hydrated clay in mechanical equilibrium with its environment, so that one imposes a pressure P ext instead of the spacing d. Such conditions are more realistic with respect to field application, since one never controls a spacing at the nanometer scale. In that case, thermodynamic equilibrium corresponds to the minimization of the so-called swelling free energy ω which is related to λ according to a the Legendre transform ω = min d (λ + P ext d). While the stable spacing corresponds to the global minimum of ω (d) = λ (d) + P ext d, local minima correspond to meta-stable states. Therefore, analyzing the pressure isotherms of Fig. 3, one can determine all stable and meta-stable hydration states at any prescribed external pressure P ext . From such stability analysis, one concludes that only convex portion of the potential λ (d) can be (meta-)stable, i.e., only decreasing branches of the pressure isotherm can be observed under pressure control: the xW hydration states.

Of particular interest is the situation in which the clay is free to swell at equilibrium with the water vapor, the usual conditions in most adsorption experiments. This corresponds to unjacketed conditions in which the external pressure is equal to the bulk vapor pressure: 6 (top left), we display how the swelling free energy ω (d) = λ (d)+P ext d profiles evolve with relative fugacity. One can readily see the different minima corresponding to the stable and meta-stable hydration states, and it appears that the global minimum evolves toward higher hydrations states as relative fugacity increases, i.e., the system exhibits crystalline swelling. We estimated the fugacity of the hydration transitions by interpolating between successive free energy profiles. Similarly, we also estimated the limits of metastability (spinodal). The estimated stable swelling and the domains of meta-stability of each hydration states are displayed in Figure 6 (top right). The 0W-1W stable transition is found to occur at a low relative fugacity (9%), while the 1W state remains metastable upon desorption down to very low relative fugacity (0.018%) and the 0W state always remains metastable upon adsorption. As for the 1W-2W transition, stable change occurs at an to the capillary domain is expected as one reaches saturation. We also display in Figure 6 (top right) XRD experimental measurements of the spacing and molecular simulation results reported in the literature. Experimental equilibration is often long and one may follow in part the metastable branches so that reported hydration transitions are expected somewhere inbetween the stable path and the spinodal points (limits of meta-stability). XRD experiments on Na-mmt usually report 0W-1W transitions at relative fugacities (= RH) below 20%, 1W-2W transition around 60%-70% and 2W-3W transition above 90% 3,5,8,9,57-64 . This is well predicted from our GC molecular simulations. Few molecular simulation works do report a swelling curve, and we could find only two works for Na Wyoming montmorillonite [START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF][START_REF] Li | Molecular Simulation Study of Montmorillonite in Contact with Water[END_REF] , none of which considers the exact same force field. These results compare well with ours in terms of spacing plateaus. [START_REF] Tambach | Molecular Simulations of Swelling Clay Minerals[END_REF] tend to underestimate the RH of the 0W-1W transition, while [START_REF] Li | Molecular Simulation Study of Montmorillonite in Contact with Water[END_REF] tend to overestimate the RH of the 1W-2W transition. More generally, this raises the question of the most appropriate force field to predict swelling. In this work, we consider clay-FF with water SPC/E, which appears to predict reasonably well the swelling curve.

P ext = P v . In Figure
Recent studies have suggested to modify clay-FF (increase the steric repulsion with bridging oxygen) in order to improve the predicted X-ray and neutron diffraction spectra and reduce the stability of the 0W state 38,39,52 . Regarding swelling, this is expected to increase the basal spacing of the dry state, and to reduce the RH of the 0W-1W transition (0W less stable). It does not seem that these changes would provide any significant improvement in the comparison with experimental swelling. Accordingly, from a swelling perspective, the interest of using the modified clay-FF is non obvious.

From the stability analysis, one can also derive the water content in unjacketed conditions.

We display in Figure 6 (bottom left) the water content as a function of relative fugacity (adsorption isotherm). Crossing the swelling and adsorption isotherms, one can also derive the spacing as a function of water content (Fig. 6 bottom right). The adsorption isotherm appears very similar to the swelling isotherm with a step-like structure. Accordingly, the swelling and the water content appear to be strongly correlated. Actual measurements of the water content are generally not as much step-like as the spacing measurements. The difficulty comes from the fact that one measures the overall water content and discriminating between inter-layer water, surface adsorption outside clay layers, and capillary water, is not obvious [START_REF] Bérend | Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites; 2, The Li+ , Na+ , K+ , Rb+ and Cs+ -exchanged forms[END_REF] . In this respect, XRD measurements of spacing are more relevant to confront our results than adsorption isotherm because it addresses specifically the clay layer scale.

Regarding the water content -swelling curve, it has been investigated accurately by many in the molecular simulation literature since it can be estimated from MD at fixed water content (does not require GC simulations) [START_REF] Greathouse | Developments in Clay Science[END_REF] . Our results are well consistent with the literature, with a first plateau corresponding to the 1W state around 0.05-0.15 g water / g clay, a second plateau for the 2W state around 0.2-0.3 g water / g clay and then a spacing almost linear with the water content [START_REF] Mooney | Adsorption of Water Vapor by Montmorillonite. II. Effect of Exchangeable Ions and Lattice Swelling as Measured by X-Ray Diffraction[END_REF][START_REF] Norrish | Crystalline Swelling of Montmorillonite: Manner of Swelling of Montmorillonite[END_REF][START_REF] Teich-Mcgoldrick | Swelling Properties of Montmorillonite and Beidellite Clay Minerals from Molecular Simulation: Comparison of Temperature, Interlayer Cation, and Charge Location Effects[END_REF][START_REF] Calvet | Hydratation de la montmorillonite et diffusion des cations compensateurs[END_REF][START_REF] Fu | Changes in the Properties of a Montmorillonite-Water System during the Adsorption and Desorption of Water: Hysteresis1[END_REF] . Here, we can also distinguish the stable, meta-stable and unstable portions of the curve, which cannot be done by MD. A majority of water contents are actually unstable, and a large portion of the plateaus are actually metastable. Also, while the potential hysteresis is large in the swelling and adsorption isotherms, very little hysteresis is visible in the water content -spacing curve.

Repeating the stability analysis at arbitrary external pressures, we can determine a complete hydration diagram of the Na-mmt (Fig. 7). The diagram displays the stable transitions between hydration states (plain lines) and the domains of meta-stability of each state (colored areas) delimited by the spinodals (dashed lines). The dotted lines represent the hydration transition between two hydration states when both states are metastable only. The domain of stability of the bulk water is restricted to a line (black plain line) because pressure and fugacity are not independent properties for a bulk fluid in isothermal conditions (Gibbs- The stable transitions are given by the plain lines. The domains of meta-stability are displayed by the colored areas delimited by the spinodals (dashed lines). The dotted lines correspond to the transitions between two hydrations states in a domain where both hydration states are metastable. Note that the domain of stability of the bulk fluid is limited to the black plain line, since fugacity and pressure are not independent properties for a bulk fluid (Gibbs-Duhem relation). The black '+' line corresponds to the superheated bulk liquid (metastable phase of the bulk water), and is also the hydration domain of the capillary water (again restricted to a line). The 2W-3W transition is not displayed here because it is almost superimposed with the bulk liquid isotherm (less than 10 MPa above). The inset (top) provides the swelling free energy barriers between two hydration states at the point of hydration transition.

Duhem relation: dP = ρdµ = ρkT df f ). Likewise, the domain of the superheated liquid (metastable bulk water) is restricted to a line (black '+' line) and is also the domain of the capillary hydration state. The hydration diagram can be used to predict the evolution of the hydration state upon any external loading in fugacity and/or pressure, both regarding stable hydration and regarding the potential hysteresis. For instance, the unjacketed adsorption and desorption of Figure 6 correspond to a loading that follows exactly the bulk vapor line.

Other loadings of interest are the drying under restricted tension that leads to desiccation cracking, or the consolidation-induced dehydration.

In practice, in geological conditions, the reservoir fluid is not pure water and is likely to contain other fluids. Of particular interest is the case of carbon dioxide (CO 2 ) for carbon storage applications. Such fluid mixture would be biphasic with one phase mostly composed of liquid water and dissolved CO 2 and the other phase mostly composed of gaseous or supercritical CO 2 with dissolved water. In clays, the water-rich phase is expected to be the wetting phase and therefore is expected to fill the smallest pores, which we study here.

The solubility of CO 2 in liquid water reaches 3 CO 2 molecules for 100 water molecules at pressures of 100 MPa 65 . These CO 2 may have some influence on the predicted hydration phase diagram, but we think this effect should be somehow limited owing to the small CO 2 content.

If the determination of the stable and metastable states is clear, to what extent metastable states are actually maintained in practice is a complex issue. It depends on the magnitude of the energy barrier that must be overcome to perform the hydration transition, in comparison with the thermal agitation energy. The mechanisms of minimum energy barriers can have many origins (e.g., heterogeneous bubble nucleation for a superheated liquid). In the particular case of clay hydration, a possible mechanism is the deformation of the clay mineral to accommodate a progressive hydration transition over an entire clay particle, like the propagation of a wave [START_REF] Honorio | Flexibility of nanolayers and stacks: implications in the nanostructuration of clays[END_REF] . In this respect, the magnitude of the energy barrier for an entire clay particle is expected to scale with the length of a clay layer, and not with its area.

Interestingly, thermal agitation of the clay layer ∼ kT √ N also scales with the length, which suggests that the energy barrier per unit length of the wave is a good indicator of the propensity to meta-stability. Balancing the energy of deformation of a clay layer and the energy barrier between hydration states, it appears that a wave-like hydration transition is expected to occur over a few nanometers only [START_REF] Honorio | Flexibility of nanolayers and stacks: implications in the nanostructuration of clays[END_REF] . Therefore, a relevant indicator of meta-stability is the swelling free energy barrier per nm 2 , compared to the thermal agitation energy (kT ). We estimated the energy barrier that appears in the swelling free energy profiles (Figure 6 and 0.5 • kT /nm 2 . Accordingly, the energy barriers are of the order of or higher than the thermal agitation energy in all cases, and meta-stability appears possible for all hydration transitions. This is consistent with the hysteresis observed in experiments (e.g., 7 for the 2W-3W transition). The energy barriers reported are the maximum at the external pressures of hydration transitions. In contrast, the energy barriers fall to zero at the spinodal (decrease almost linear with external pressure from the hydration transition to the spinodal, see for instance 35 Fig. 8). So meta-stability is not expected to hold until the spinodal. The large magnitude of the energy barrier for the 0W-1W transition suggests that most of the meta-stability domain can be observed, whereas the moderate magnitude for the 1W-2W and 2W-3W transitions suggests that only part of the meta-stability can be observed.

Confined water properties

While hydration transitions have attracted much attention in the literature, the mechanics of the hydration states have been little investigated. In particular, while the mechanical behavior of bulk water can be described by a single modulus, namely the bulk modulus

K b = ρ ∂P ∂ρ T
(inverse of compressibility), the mechanics of the confined water is more complex because the Gibbs-Duhem equation does not hold anymore for the adsorbed water 42 . The Gibbs-Duhem equation dP | T = ρdµ, which arises from the extensivity of the behavior of the fluid with respect to volume, entropy and number of particles, imposes a one-to-one relationship between pressure and chemical potential at fixed temperature. In particular, it imposes ∂P ∂V µ,T = 0, which is clearly not true for the water confined in clays since the pressure isotherms vary with the spacing (Fig. 2). Moreover, bulk water exhibits a one-to-one relationship between pressure and density so that the bulk modulus characterizes both the rigidity with respect to volume and to number of particles (K b = N ∂P ∂N V,T = -V ∂P ∂V N,T ). This does not hold anymore for the confined water and one must distinguish the two rigidities.

As we rederive it hereafter, a complete mechanical description of confined water requires 3 independent moduli, and, in this section, we provide estimates of those 3 moduli for the 1W and 2W hydration states. Estimates for the 3W state were not accurate enough and are not reported here.

Let us consider a closed volume of bulk water in isothermal conditions, for which Gibbs-Duhem equation holds. The thermodynamic state of the water depends on a single variable, the density, and the constitutive behavior is fully described by a single moduli, the bulk In what follows, we report estimates only for relative fugacities higher than 30% (which still covers a majority of situations for practical applications).

∂N

∂V P b ,T for the confined water. This quantity also corresponds to the coupling term in the constitutive behavior (2):

K N V K N N N V = ∂P ∂µ V,T = ∂N
∂V µ,T . The superheated liquid water is considered as the reference bulk state. One recovers a value of 1 (limit of the superheated liquid) at large spacings corresponding to the capillary domain. are averages over several relative fugacities. We distinguish the domains 30%-100% and 100%-407%, and error bars are provided for the standard deviation. The Maxwell relation can be accurately verified when the superheated liquid is taken as the reference bulk state (Fig. 8). Considering the vapor as the reference bulk state yields similar results, but with a higher sensibility to RH since the vapor density ρ b evolves significantly. Since ∂P ∂P b V,T can be estimated only on average over a range of RH, considering the superheated liquid as reference bulk state appears more appropriate.

In Figure 8, not only the Maxwell relation is well verified, but one can also see that the derivatives converge to 1 in the capillary domain, which confirms that one recovers the case of the superheated bulk water. Deviations from 1, characteristic of the non-validity of Gibbs-Duhem equation, are significant for spacings below 16 Å (1W and 2W hydration states). In what follows, we consider the derivation with respect to spacing ( ∂N ∂V µ,T ) rather than the derivation with respect to chemical potential ( ∂P ∂µ V,T

) to estimate the coupling coefficient

K N V K N N .
We report in Figure 9 the estimated values of the K IJ for the 1W and 2W hydration states. The drained rigidity

K d V V = -V ∂P ∂V µ,T , coupling coefficient K N V K N N = V N ∂N
∂V µ,T , and inverse chemical rigidity 1

K N N = V N 2 ∂N ∂µ V ,T
are derived from the pressure and water content isotherms. The undrained stiffness K V V = -V ∂P ∂V N,T = K d V V + (K N V ) 2 K N N is derived from the three other K IJ . In the limit of a bulk fluid verifying Gibbs-Duhem equation, K d V V = 0 and all the other moduli are equal to the bulk modulus K V V = K N V = K N N = K b . Here, we clearly observe K d V V = 0 and K V V , K N V , and K N N independent. Deviation from the bulk behavior is more pronounced for the 1W state than for the 2W state. The magnitudes K IJ are also independent of the spacing and we report here the averages over all the capillary domain. As expected, the behavior in this domain compares very well with that of the bulk liquid: K V V = K N V = K N N = K b and K d V V = 0. Osmotic effects, as predicted from the double layer theory, are not excluded but are beyond the accuracy of our GC simulations.

Conclusion

In this paper, we report results of grand canonical molecular simulations of a model Wyoming Na-montmorillonite at 300K, for basal spacings covering all the hydration states and for various chemical potentials both below and above bulk water vapor saturation. This dataset is large enough to derive a full hydration diagram of the Na-montmorillonite in function of external pressure and water chemical potential. This diagram provides the domains of (meta-)stability of the different hydration states and can be used to predict the evolution of the hydration for any arbitrary evolution of pressure and chemical potential. The case of the free swelling/shrinkage upon change in relative humidity (unjacketed response) is considered in detail. Stability analysis provides the stable swelling curve, the domain of hysteresis and the energy barriers to overcome meta-stability. These results are consistent with the existing literature on Na-montmorillonite (both experiments and simulation works).

The molecular simulation results are also used to estimate the moduli characterizing the constitutive behavior of confined water. While the constitutive behavior of bulk water is fully characterized by a single modulus (bulk modulus), the characterization of confined water requires 3 independent moduli, because the Gibbs-Duhem equation does not hold anymore for confined water. The estimates of these moduli are in the range of a few GPa and can reach up to twice the bulk modulus of the bulk liquid. Interestingly, the rigidity of coupling K N V = N ∂P ∂N V,T = -N ∂µ ∂V N,T exhibits values quite similar to the bulk modulus of the bulk liquid, whereas the other moduli vary significantly. At spacings larger than 19Å and fugacities above bulk cavitation (f /f sat ≥ 30%), we recover the behavior of the superheated bulk liquid: the pressure and density isotherms converge to the bulk pressure and density, and all the constitutive moduli converge to the bulk modulus. This is consistent with the capillary theory. Osmotic effects (double layer theory) are not excluded but are beyond the accuracy of the grand canonical molecular simulations.

Figure 1 :

 1 Figure 1: Bulk properties of water SPC/E at 300K. (top) Density-pressure curve obtained by canonical simulations compared to the equation of state of Wagner and Pruß 50 (IAPWS-95).(bottom) Fugacity-pressure and fugacity-density curves obtained either by direct GC simulations (configurational bias or usual Mont Carlo) or by canonical simulation and thermodynamic integration of the Gibbs-Duhem equation. Saturation corresponds to the intersection between the liquid and vapor pressures.

Figure 2 :

 2 Figure 2: (top) Pressure P zz normal to the Na-mmt layer as a function of the z-axis spacing d at a temperature of 300K. The points are the results from the GC simulation, and the lines are Savitzky-Golay filter of the raw results. The estimated error is displayed for half of the curves (typical error ±9 MPa). Only about one relative fugacities out of two is displayed here, see Figure 3 for all fugacities. (bottom) Snapshots of sampled atomic configurations for the different stability domains.

Figure 4 :

 4 Figure 4: Confined water density in Na-mmt as a function of the z-axis spacing d at a temperature of 300K. The points are the results from the GC simulation (typical error ±3 kg/m 3 ), and the lines are Savitzky-Golay filter of the raw results. Only about one relative fugacities out of two is displayed here, see Figure 5 for all fugacities.

Figure 7 :

 7 Figure7: Hydration diagram of Na-mmt at 300K: (top) overview over the complete range of pressures and relative fugacities, (bottom) zoom on the 2W, 3W and capillary domains. The stable transitions are given by the plain lines. The domains of meta-stability are displayed by the colored areas delimited by the spinodals (dashed lines). The dotted lines correspond to the transitions between two hydrations states in a domain where both hydration states are metastable. Note that the domain of stability of the bulk fluid is limited to the black plain line, since fugacity and pressure are not independent properties for a bulk fluid (Gibbs-Duhem relation). The black '+' line corresponds to the superheated bulk liquid (metastable phase of the bulk water), and is also the hydration domain of the capillary water (again restricted to a line). The 2W-3W transition is not displayed here because it is almost superimposed with the bulk liquid isotherm (less than 10 MPa above). The inset (top) provides the swelling free energy barriers between two hydration states at the point of hydration transition.

  top left) for the different hydration transition (inset in Fig. 7 top). The energy barrier strongly depends on the type of hydration transition, and moderately depends on the relative fugacity. Typical values are about 0.07 N/m, 0.008 N/m, and 0.002 N/m for the 0W-1W, 1W-2W, and 2W-3W transitions, respectively. Reduced to a nm 2 and compared to the thermal agitation energy, these barriers correspond to 17.6 • kT /nm 2 , 2 • kT /nm 2

  modulus: dP = ρdµ = K b dρ ρ . Let us now consider a volume of confined water in isothermal conditions, for which Gibbs-Duhem equation does not hold. The complete description of the thermodynamic state requires distinguishing the number of molecule N and the volume V as independent state parameters. The first derivatives of the Helmholtz free energy F with respect to N and V define the thermodynamic quantities conjugated to the two control variables: the pressure P = -∂F ∂V N,T and the chemical potential µ = ∂F ∂N V,T . The complete description of the constitutive behavior of the confined water requires 3 independent moduli to evaluate the changes of P and µ upon changes of N and V : the rigidity with respect to volume K V V = -V ∂P ∂V N,T , the 'chemical' rigidity K N N = N 2 V ∂µ ∂N V ,T , and the rigidity of coupling K N V = N ∂P ∂N V,T = -N ∂µ ∂V N,T (note that we take advantage of a Maxwell relation in the range of several hundreds of GPa, whereas the confined water rigidity is in the range of a few GPa 41 . Since the K IJ are first derivatives of the pressure and water content with respect to spacing and chemical potential, the quality of our estimations depends on how much GC results are resolved in spacing and chemical potential. It appears that our results are well resolved in spacing (every 0.25 Å), but less in chemical potential, especially in the domain of low chemical potentials (Tab. 1). As a consequence, we could obtain reasonably accurate estimates of the K IJ only for relative fugacities above 30% (above bulk cavitation).

Figure 8 :

 8 Figure 8: Verification of the Maxwell relation ∂P ∂µ V,T = ∂N ∂V µ,T ⇔ ∂P ∂P b V,T = 1 ρ b

First, we investigate

  the validity of the thermodynamic representation of confined water by checking the Maxwell relation ∂P ∂µ V,T = ∂N ∂V µ,T . The quantity involved in this relation also corresponds to the coupling coefficient in the constitutive equations (2): K N V µ,T . This Maxwell relation establishes a relation between the pressure and water content isotherms and can be used to derive one isotherm from the other 67 . Deviations from this relation would point to an incomplete thermodynamic description, e.g., the need for more state variables. In practice, we prefer rewriting the Maxwell relation as follows: ∂P ∂P b V,T = 1 ρ b ∂N ∂V P b ,T , where we replaced the chemical potential µ by the bulk pressure P b and bulk density ρ b thanks to the Gibbs-Duhem equation (subsript b refers to the bulk properties at the same µ and T ). Since the GC results are more resolved in spacing than in relative fugacity, the estimations of 1 ρ b ∂N ∂V P b ,T are much more precise than that of ∂P ∂P b V,T . Accordingly, the values reported for ∂P ∂P b V,T

Figure 9 :

 9 Figure 9: Estimates of the K IJ moduli characterizing the constitutive behavior of confined water.

( 6 )

 6 Slade, P. G.; Quirck, J. P.; Norrish, K. Crystalline Swelling of Smectite Samples in Concentrated NaCl Solutions in Relation to Layer Charge. Clays and Clay Minerals 1991, 39, 234-238. (7) Laird, D. A.; Shang, C.; Thompson, M. L. Hysteresis in Crystalline Swelling of Smectites. Journal of Colloid and Interface Science 1995, 171, 240-245. (8) Ferrage, E.; Lanson, B.; Sakharov, B. A.; Drits, V. A. Investigation of smectite hydration properties by modeling experimental X-ray diffraction patterns: Part I. Montmorillonite hydration properties. American Mineralogist 2005, 90, 1358-1374. (9) Ferrage, E.; Lanson, B.; Michot, L. J.; Robert, J.-L. Hydration Properties and Interlayer Organization of Water and Ions in Synthetic Na-Smectite with Tetrahedral Layer Charge. Part 1. Results from X-ray Diffraction Profile Modeling. The Journal of Physical Chemistry C 2010, 114, 4515-4526. (10) Israelachvili, J. N.; Pashley, R. M. Molecular layering of water at surfaces and origin of repulsive hydration forces. Nature 1983, 306, 249-250. (11) Yan, L.; Low, P. F.; Roth, C. B. Swelling Pressure of Montmorillonite Layers versus H-O-H Bending Frequency of the Interlayer Water. Clays and Clay Minerals 1996, 44, 749-756. (12) Ortega, J. A.; Ulm, F.-J.; Abousleiman, Y. The effect of the nanogranular nature of shale on their poroelastic behavior. Acta Geotechnica 2007, 2, 155-182. (13) Ebrahimi, D.; Pellenq, R. J.; Whittle, A. J. Nanoscale Elastic Properties of Montmorillonite upon Water Adsorption. Langmuir 2012, 28, 16855-16863. (14) Delville, A. Modeling the clay-water interface. Langmuir 1991, 7, 547-555. (15) Skipper, N. T.; Refson, K.; McConnell, J. D. C. Computer simulation of interlayer water in 2:1 clays. The Journal of Chemical Physics 1991, 94, 7434-7445.

  ings are unambiguously defined, the inter-layer spacing results from an arbitrary partition of the z-axis spacing between the mineral layer and the electrolyte. Here, we consider a thickness of the mineral layer equal to 9.36 Å, which corresponds to the z-axis spacing in the 0W state (dry) in the relaxed configuration (P zz = 0 MPa, T = 300K). Accordingly, we define the inter-layer spacing as h = d -9.36Å. Other definitions are possible (e.g., based on Van der Walls radii, spacing at which adsorption starts, etc.) and would lead to inter-layer spacing that would differ by ±0.2Å. The inter-layer spacing is used in this work to estimate the confined water density ρ = M H 2 O N Ah , with N the number of mole of water, M H 2 O the molar mass and A the area of the mineral layer. It is also involved in the calculation of the moduli characterizing the mechanics of confined water (see section 'Discussion/Confined water properties').

	00998699	0.001	0.0100726	7.25755E-05		-		-	-
	0.0998699 0.998699 9.98699 19.9740 49.9349 99.8699 199.740	0.01 0.1 1 2 5 10 20 500 0 500 Pressure P zz (MPa) 1000	0.100726 1.00726 10.0726 20.1452 50.3631 100.726 201.621	7.25755E-04 7.25755E-03 7.25755E-02 0.145151 0.362878 0.725755 1.45375 f/f sat (undersaturated) 30% 50% 70% 100% 0% 0.1% 1% 10%	-------f/f sat (supersaturated) 204% 407%	-------	-------
	299.610	30 1000	302.573		2.18243		-157.570		892.454	878.864
	399.480 499.349	40 50	10	403.549 504.567 12	14	2.91153 3.64155 16 Spacing d ( ) 2W	18	-121.409 -92.6449 20	22	921.579 940.079	1337.36 1554.25
	599.219	60		605.679 1W		4.37334		-68.7367	953.954 capillary	1709.61
	699.089	70	0W	706.938		5.10763		-48.2569 3W		965.044	1834.48
	798.959	80		808.381		5.84457		-30.3293		974.256	1940.04
	898.829	90		909.998		6.58329		-14.3769		982.115	2031.91
	998.699	100	1011.71		7.32145	1.01171E-03	988.952	2113.48
	1432.95	143.482	-		-		50		1011.17	2390.70
	2040.70	204.336	-		-		100		1031.40	2661.60
	2887.11	289.087	-		-		150		1050.03	2928.75
	4060.63	406.592	-		-		200		1067.33	3193.58

  Complete series of pressure isotherms (top) and zoom on the domain of large spacings (bottom). The decreasing branches of the oscillations correspond to the xW hydration states. For spacings larger than 19 Å, we refer to the capillary domain since no more oscillations can be distinguished and the normal pressure appears almost equal to the bulk pressure (dotted lines): superheated liquid above cavitation (f /f sat ≥ 30%), vapor below cavitation (f /f sat ≤ 20%).

						f/f sat (undersaturated) 5% 10% 20% 30% 40% 50% 0% 0.001% 0.01% 0.1% 1% 2%	60% 70% 80% 90% 100%	f/f sat (supersaturated) 143% 204% 289% 407%
			0W		1W		2W	3W	capillary
	Pressure P zz (MPa)	100 0 100 200 300 400 500					f/f sat (undersaturated) f/f sat (supersaturated) 5% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 143% 204% 289% 407% 0% 0.001% 0.01% 0.1% 1% 2%
		200
		300	0W	10	12 1W	14	16 Spacing d ( ) 2W	18 3W	20	22 capillary
	Figure 3:

  evolves toward higher hydration states as relative fugacity increases which corresponds to the crystalline swelling (black squares). The local minima are meta-stable states. The stable evolution and domains of meta-stability of spacing (top right) and water content (bottom left) are derived from the stability analysis. Combining both isotherms, one obtains the curve relating water content and spacing which exhibits little hysteresis (bottom right).

	intermediate relative fugacity (67%), the 2W state remains meta-stable upon desorption
	down to 36%, and the 1W state always remains metastable upon adsorption. Finally the
	2W-3W transition occurs above 90% of relative fugacity: the stable transition is expected
	at 97.7%, the 3W states remains meta-stable upon desorption down to 90.4%, and the
	2W state is always metastable upon adsorption. Finally, a transition from the 3W state
			0W	1W		2W	3W	capillary	f/f sat (undersaturated) 30% 40% 50% 60% 70% 80% 90% 100% stable path 0% 0.001% 0.01% 0.1% 1% 2% 5% 10% 20%	Spacing d ( )	14 16 18 20 22	metastable limit 0W 1W 3W 2W stable path 0W 1W 2W 3W	1W Cases et al. (1992) 3W Tambach et al. (2004) 2W experiments Mooney et al. (1952) Fripiat et al. (1965) Keren and Shainberg (1975) Watanabe and Sato (1988) Li et al. (2019) simulations Fu et al. (1990) Laird et al. (1995) Dios Cancela et al. (1997) 2W Ferrage et al. (2005)
									f/f sat (supersaturated) 143% 204% 289% 407%		12 10	1W	0W
														0	20	40 Relative fugacity f/f sat (%) 60	80	100
	Water content m H 2 O /m solid (kg/kg)	0 0.0 0.1 0.2 0.3 0.4 0.5	1W stable path 0W 1W 2W 3W	20	40 Relative fugacity f/f sat (%) 60 0W 2W metastable limit 0W 1W 2W 3W	80	1W	3W	100 2W	Spacing d ( )	0.0 10 12 14 16 18 20 22	0.1 stable path 0W 1W 2W 3W	0.2 Water content m H 2 O /m solid (kg/kg) 0.3 0.4 metastable limit 0W 1W 2W 3W	0.5

Figure 6: (top left) Evolution of the swelling free energy profiles ω (d) = λ (d) + P ext d with the relative fugacity at 300K in unjacketed conditions (external pressure equal to the bulk vapor pressureP ext = P v ). Thermodynamic equilibrium (minimum of ω)

Acknowledgement

Funding from the ALLUVIUM and ARPENTONS projects are gratefully acknowledged.

The ALLUVIUM project is part of the FUTURE consortium and is funded in the framework of the Programme Investissement d'Avenir supervised by the French National Research Agency (ANR). The ARPENTONS project is part of the NEEDS/MIPOR initiative (CNRS, ANDRA, CEA, EDF, IRSN, BRGM).

of F ). The constitutive equations describing the behavior upon changes of N and V are:

The precise definitions of the K IJ were chosen to have units of pressure, like the usual bulk modulus for bulk water. As a consequence, when Gibbs-Duhem equation is valid, all the K IJ reduce to a single value K IJ = K and one recovers the usual behavior (dP = ρdµ = K dρ ρ ). The constitutive equations (1) are not adapted to the GC conditions of the molecular simulation (control of µ instead of N ). Rearranging equations (1) we obtain:

with

∂V µ,T the rigidity with respect to volume in drained conditions (µ fixed). The moduli that appear in the constitutive equations (2) can be estimated directly from the GC molecular simulation results by computing the derivatives of the pressure and adsorption isotherms (Figs. 2 and4). Note, however, that the volume V of the confined water is somehow ambiguous since it is based on how the interlayer spacing h is defined (see discussion at the beginning of section 'Results'). Here, we consider h = d-9.36Å where 9.36Å corresponds to the thickness of the mineral layer estimated as the z-axis spacing in the 0W state (dry) in the relaxed configuration (P zz = 0MPa, T = 300K). Other choices could lead to mineral thicknesses that can differ by ±0.2Å. This corresponds to an accuracy of the interlayer spacing of ±7%, ±3.5%, and < ±2% for the 1W, 2W, and capillary states, respectively. Accordingly, except for K N V which does not depend on the interlayer spacing, the values reported in this section for the different moduli are relevant up to this maximum accuracy. Note also that our analysis also assumes the solid much more rigid than the interlayer fluid. This assumption is well verified since estimated values of the solid rigidities are of the moduli are in the range of a few GPa, i.e., in the range of the bulk modulus of liquid water. Yet, the estimated values vary significantly within the 1W and 2W domains from 0 to about twice the bulk modulus of liquid water. All moduli are increasing with fugacity and no particular change can be observed at the transition from undersaturation (f /f sat ≤ 100%) to supersaturation (f /f sat ≥ 100%). Interestingly, the rigidity of coupling

∂V N,T exhibits much less variability than the other moduli, and its average value over the 1W and 2W domains, respectively, compares well with the bulk modulus of liquid water (Fig. 10). There is, though, a surprising 'bump' in the estimate of K N V at the highest fugacity (f /f sat = 407%) in the 1W domain. We believe this is an irrelevant artifact due to the inaccuracy of our estimation of derivatives at the boundary of our domain of study. To satisfy thermodynamic stability in osmotic conditions (P and µ fixed), the K IJ moduli cannot take any value but must verify K V V ≥ 0, K N N ≥ 0, and

We observe in Figure 9 that these conditions are well satisfied in the 1W and 2W domains. The most restrictive condition is K d V V ≥ 0, and one readily observes that this condition is no more satisfied outside the 1W and 2W domains.

The average K N V over the 1W and 2W domains is also quite similar to the bulk modulus. This is not the case of the other K IJ of the 1W and 2W water.

We provide in Figure 10 the estimated K IJ moduli in the capillary domain. Since the pressure and water density are independent of the spacing in this domain (Figs. 2 and4), the