Water and Ion Dynamics in Confined Media: A Multi-Scale Study of the Clay/Water Interface

: This review details a large panel of experimental studies (Inelastic Neutron Scattering, Quasi-Elastic Neutron Scattering, Nuclear Magnetic Resonance relaxometry, Pulsed-Gradient Spin-Echo attenuation, Nuclear Magnetic Resonance Imaging, macroscopic diffusion experiments) used recently to probe, over a large distribution of characteristic times (from pico-second up to days), the dynamical properties of water molecules and neutralizing cations diffusing within clay/water interfacial media. The purpose of this review is not to describe these various experimental methods in detail but, rather, to investigate the speciﬁc dynamical information obtained by each of them concerning these clay/water interfacial media. In addition, this review also illustrates the various numerical methods (quantum Density Functional Theory, classical Molecular Dynamics, Brownian Dynamics, macroscopic differential equations) used to interpret these various experimental data by analyzing the corresponding multi-scale dynamical processes. The purpose of this multi-scale study is to perform a bottom-up analysis of the dynamical properties of conﬁned ions and water molecules, by using complementary experimental and numerical studies covering a broad range of diffusion times (between pico-seconds up to days) and corresponding diffusion lengths (between Angstroms and centimeters). In the context of such a bottom-up approach, the numerical modeling of the dynamical properties of the diffusing probes is based on experimental or numerical investigations performed on a smaller scale, thus avoiding the use of empirical or ﬁtted parameters.

For that purpose, some of the dynamical studies presented here were performed with synthetic clays (laponite and saponite) in order to avoid impurities within the clay network. In the case of natural clays (kaolinite, beidellite, and vermiculite), the samples were purified according to standard procedures [60] and size separated by centrifugation. Under such conditions, clay/water interfaces are well characterized systems and ideal to investigate the influence of charged interfaces on the thermodynamical and dynamical properties of confined fluids.
Finally, numerical modeling is required to relate the organization of these clay/water interfaces and the dynamical information extracted from the various experimental data obtained by INS, QENS, PGSE, MRI, and macroscopic diffusion experiments (see Figure 1). In that context, an atomic description of the clay/water interface is required for modeling the structural and dynamical properties of these interfacial systems probed in the crystalline swelling regime. For that purpose, quantum calculations are used to directly determine the dynamical properties of confined ions and water molecules by using ab initio DFT quantum methods [61,62]. Quantum calculations are also useful to evaluate atomic properties of the clay network, such as atomic partial charges, by using semi-empirical [63] or ab initio [64] quantum calculations. These parameters are further exploited in classical statistical treatments of the clay/water interfaces by using classical Monte Carlo (MC) simulations and Molecular Dynamics (MD) to investigate the structural and dynam-ical properties of the confined ions and water molecules within the crystalline swelling regime [2,12,17,28,40,41,46,[63][64][65][66]. By contrast, in the so-called osmotic regime, one may easily neglect the atomic structure of the clay network and the solvent molecules by using classical Brownian Dynamics (BD) in order to determine the effective mobility [67] of the neutralizing counterions in the presence of structureless charged surfaces [68]. In the same context, numerical integration of differential equations can also be used to study the mobility of diffusing probes at the macroscopic scale (Refs. [54,69], and references therein). In the present review, these techniques are described following a path of increasing probed time-scales, together with the associated numerical methods, in order to highlight the methodology used to get information regarding the multi-scale properties of clay/water interfaces. The purpose of that study is to exploit such a multi-scale analysis of the dynamical properties of ions and water molecules diffusing within clay sediments in order to develop a bottom-up approach of such complex interfacial dynamical processes. In such analysis, we avoid to use fitted parameters derived from empirical models by performing realistic numerical modeling exploiting structural or dynamical properties, such as average residence time or local mobility, extracted from preliminary studies performed on a smaller scale.

Inelastic Neutron Scattering
The intensity of the INS spectra is directly proportional to the so-called Generalized Density of States, labeled GDOS(ω), that is the density of states of a system taking into account the scattering cross section of the various atoms: where M j and σ j are, respectively, the atomic mass and the cross section of the atom labeled j. An important point is that the scattering cross section of Hydrogen is almost 2 order of magnitude that of the other atoms of the clay (Si, Al, O), so the GDOS(ω) gives mainly information on the vibrational modes of the H atoms. For the hydrated systems, the simulated spectra were calculated with Density Functional Theory Molecular Dynamics simulations since the confined water molecules remain liquids and, thus, disordered compared to monocrystalline clay network. After solving classical Newton equations, the GDOS(ω) was approximated by evaluating the velocity auto-correlation function of the atoms: where v j (t) is the instantaneous velocity of atom j. Before comparing the predictions obtained from the ab initio MD simulations with experimental data, the contributions from the dry clay network is subtracted [18]. The resulting spectra are displayed in Figure 3 for two hydration states corresponding, respectively, to the formation of one and two hydration layers [70]. In both longitudinal and perpendicular directions, the agreements between predicted g H−water (ω) and experimental spectra are semi-quantitative with a fair matching of the characteristic frequencies (see Figure 3), partially validating the use of the above-mentioned approximation in the treatment of such local motions of the confined water molecules.
(a) (b) Figure 3. Comparison between the experimental GDOS water (ω) and simulated g H−water (ω), see Equation (2), for INS spectra of hydrated saponite sediments along the directions parallel (a) and perpendicular (b) to the clay layers (see text). Reprinted with permission from Ref. [18]. Copyright (2017) American Chemical Society.

Quasi-Elastic Neutron Scattering
Quasi-Elastic Neutron Scattering (QENS) experiments are performed [12,28] to probe the diffusive motions of the confined water molecules and their orientational dependence regarding the clay surface. These measurements are performed with the same saponite synthetic clay at various hydration states and neutralized by sodium or calcium counterions. Since the dynamical range probed by these QENS measurement extend up to a few hundred pico-seconds (see Figure 1), the molecular modeling of the clay/water interface is performed by using the classical CLAYYFF force field [64], with a minor correction of its Lennard/Jones parameters introduced to improve the agreement between the structure of the clay/water interface obtained by X-Ray and neutron scattering experiments and simulated data [2,57,70]. Grand Canonical Monte Carlo (GCMC) simulations [71] are first performed to determine the number of confined water molecules as a function of the experimentally fixed water chemical potential [12,28] and the layer-to-layer distance previously determined by X-Ray Scattering (XRS) [70]. The mobility of the confined water molecules and neutralizing sodium or calcium counterions is then modeled by performing numerical simulations of Molecular Dynamics (MD) [72] (see Figure 4), starting from an equilibrium configuration previously determined by GCMC simulations. The QENS spectra displayed in Figure 5 are directly evaluated from the Fourier Transform of the Intermediate Scattering Function (ISF): convoluted by the spectral resolution of INS experimental set-up [28]. As displayed in Figure 5, an excellent agreement is obtained in all cases between the experimental spectra and the simulated numerical data. As a consequence, one may be quite confident in the validity of the water mobility determined by these numerical simulations on a timescale slightly shorter than nano-second. The water self-diffusion coefficients displayed in Table 1 are obtained by integrating the velocity auto-correlation function [73] of the water molecules, allowing to determine the local water mobility in the direction parallel and perpendicular to the clay surface [74]. In Figure 5, the wave numbers q probed by these QENS experiments vary between 0.3 and 1.8 Å −1 , thus restricting the corresponding diffusion length between 21 and 3.5 Å. Furthermore, the same QENS measurements probe a limited dynamical range characterized by a spectral window of 0.25 meV with a resolution of 0.0025 meV (see Figure 5), corresponding to maximum diffusion times of 132 ps. By taking into account the order of magnitude of the apparent mobility of the confined water displayed in Table 1, this maximum diffusion period restricts the order of magnitude of the probed diffusion length to typically a few Angstroms. Furthermore, as displayed in Figure 4, the time-scale probed by the QENS measurements totally matches the dynamical range probed by MD numerical simulations. As a consequence, the agreement between these QENS experimental data and the numerical results obtained by MD simulations depends not only on the quality of the numerical model used in that study but also on the matching between the time-scales and diffusion-ranges probed by both experimental and numerical methods. As a consequence, the self-diffusion coefficients reported here in the direction perpendicular to the clay surface quantifies only the local mobility of the confined water molecules and does not correspond to macroscopic displacements totally impeached by the confining clay platelets (see Figure 6a). By contrast, the water mobility of the water molecules quantified by the same QENS measurements in the direction parallel to the clay surface can be used to determine the order of magnitude of the residence time of the same water molecules (labeled τ C in Figure 6a) confined in the interlamellar space between two clay platelets. As displayed in Table 1, confinement under low electrostatic coupling (i.e., for 0.8 electron per unit cell) simply reduces the water mobility by a factor ten by comparison with bulk water. By contrast, at higher electrostatic coupling (i.e., for 1.4 electron per unit cell), the water mobility may be reduced by two orders of magnitude. These numerical results illustrate how geometric and energetic constraints can significantly modify the dynamics of confined water molecules.  Table 1. Water dynamical parameters, D // and D ⊥ , respectively, the self-diffusion coefficient defined in the direction parallel and perpendicular to the clay surface, deduced from MD simulations for all samples investigated at 300 K. The data are from Table 3 in the Ref. [28] with permission from Ref. [28]. Copyright (2012) American Chemical Society. Na + -saponite 0.8 e/uc-1W 2.4 ± 0.20 Na + -saponite 0.8 e/uc-2W 7.2 ± 0.50 3.0 ± 1.0 Na + -saponite 1.4 e/uc-1W 0.6 ± 0.05 Na + -saponite 1.4 e/uc-2W

Nuclear Magnetic Resonance Relaxometry
In the last decade, Nuclear Magnetic Resonance (NMR) relaxometry appeared as a powerful tool [75,76] to investigate dynamics of fluids under confinement by measuring the frequency variation of the NMR relaxation rates. Usual NMR spectrometers allow to measure easily the relaxation of the nuclear magnetization in the direction parallel and perpendicular to the static magnetic field, leading, respectively, to the longitudinal and transverse relaxation rates. For bulk liquids, these two independent measurements lead to the same relaxation rates. By contrast, the transverse relaxation rates of confined fluids largely exceed their longitudinal counterpart due to their large frequency variations. The purpose of NMR relaxometry is to probe that frequency variation of the relaxation rates in order to extract dynamical information on the molecular motions (rotation and diffusion) inducing the fluctuations of the magnetic coupling monitoring the NMR relaxation processes. For I = 1/2 spin nuclei [35,75,76], field-cycling NMR spectrometers are perfectly appropriate to perform such investigations since the order of magnitude of the probed relaxation rates is generally compatible with the time required to switch the magnetic field of these spectrometers (typically, milli-seconds). As an example, field cycling NMR relaxometry measurements were only performed with dilute dispersions of synthetic clays [37]. By contrast, the same field cycling spectrometers become inefficient for quadrupolar nuclei (spin I > 1/2) because of the large enhancement of their relaxation rates under confinement. In that context, spin-locking relaxation measurement [77][78][79][80] is a powerful tool since it may be performed with usual spectrometer and is simply limited by the short dead time required to generate a sequence of NMR pulses (typically, micro-seconds).
Furthermore, a complete understanding of the NMR relaxation processes requires a quantum analysis of the spin. For I = 1/2 spin nuclei, this is easily performed by using the four Pauli matrices [75] to describe the three independent states of the magnetization (M x , M y and M z ) completed by the unit operator, in order to obtain a complete basis set. For quadrupolar nuclei (spin I > 1/2), this basis set becomes incomplete and new operators must be added [81] (five for I = 1 spin nuclei like 2 H, twelve for I = 3/2 spin nuclei like 7 Li or 23 Na and up to sixty for I = 7/2 spin nuclei like 133 Cs). These supplementary operators describe new quantum states, also called coherences, accessible to the nuclear magnetization by using specific pulse sequences [82][83][84]. By measuring the corresponding multi-quanta relaxation rates [85], it becomes possible to separately identify the contributions from the different relaxation mechanisms of confined quadrupolar nuclei [10,41,43,44,51,80,86]. Finally, by performing selective multi-quanta spin-locking relaxation measurements [85], one can extract dynamical information on confined fluids from the frequency variations of the relaxation rates of these additional coherences.
In the case of heavy water, the electrostatic field gradient felt by the deuterium atom is oriented along the −→ OD director. As a consequence, the NMR relaxation rate of these deuterium atoms may be determined by the correlation time describing the reorientation of that −→ OD director by reference with the static magnetic field. In bulk water, that correlation time is about 3.6 pico-seconds, leading to a NMR relaxation rate of 2.4 s −1 . By applying the same analysis to the NMR relaxation rate of confined water molecules, the measured relaxation rate (typically, 10 4 s −1 ; see Figure 7) requires reorientation correlation time of the −→ OD director with an order of magnitude of 15 nano-seconds, i.e., three order of magnitude larger than the value measured by QENS (see Michot et al. [12]). That interpretation is, thus, invalid since it neglects the influence of confinement on the specific orientation of the confined water molecules (see Figure 6d), leading to the splitting of the deuterium resonance line (see Figure 6b,c). As a consequence, the decorrelation of the electric field gradient felt by the confined water molecules is incomplete until they remain confined inside the interlamellar space of the same clay sediment (see Figure 6a). By contrast, after desorption from one oriented micro-domain, the residual quadrupolar coupling vanishes and the long-time decrease of the correlation function monitoring the water quadrupolar coupling is then driven by the water residence time τ C , as displayed in Figure 6a.   Figure 6. (a) Schematic view of the multi-scale organization of the hydrated clay-sediment induced by the coexistence of various micro-domains resulting from the stacking of hydrated clay platelets; (b) variation of the 2 H NMR spectra of the confined water molecules as a function of the sediment orientation in the static magnetic field; (c) variation of the residual quadrupolar splitting (see Equation (5)), extracted from the spectra displayed in Figure 6b; (d) snapshot illustrating one equilibrium configuration obtained by Grand Canonical Monte Carlo simulations of the water molecules and neutralizing Na + counterions confined between fragments of saponite platelets. Reprinted with permission from Ref. [80]. Copyright (2014) American Chemical Society.
In that context, a first set of 2 H multi-quanta spin-locking relaxometry experiments were performed to determine the average residence time of water molecules confined within the interlamellar space of clay platelets forming aggregates (Figure 6a) of dense sediments of Na + -beidellite clay [80]. As illustrated in Figure 6b, the apparent splitting ν Q of the 2 H resonance line varies as a function of the sediment orientation β LF into the static magnetic field, according to the expected relationship ( Figure 6c): where β LF is the angle between the direction of the static magnetic field and the normal n to the surface of the clay film. This residual splitting of the 2 H NMR resonance line results from the specific orientation of the confined water molecules as illustrated by an equilibrium configuration of the confined ions and water molecules obtained by Grand Canonical Monte Carlo simulations (Figure 6d). By measuring the specific relaxation rates of the so-called T 10 (s), T 11 (a, s), T 20 , and T 22 (a, s) coherences [82] of this I = 1 spin nucleus as a function of the sediment orientation into the static magnetic field (Figure 7), it becomes possible to quantify the relative contributions of the quadrupolar and heteronuclear dipolar NMR relaxation mechanisms [80]. Three specific angles (0 • , 30 • , and 90 • ) are then selected to determine the time evolution of four independent coherences (labeled T 11 (s), T 21 (a), T 21 (s), and T 22 (a) [82][83][84] under spin-locking conditions by applying an irradiation field of calibrated irradiation power ω 1 . For I = 1/2 spin nuclei, such spin-locking relaxometry measurements probe only a single angular velocity (2 ω 1 ) for each irradiation power [77][78][79]. By contrast, in the presence of a residual quadrupolar coupling (ω Q ), I = 1 spin nuclei probe three different angular velocities: for a single irradiation power ω 1 [82]. Some raw data are displayed in Figure 8 in addition to their Fourier transforms illustrating the specific angular velocities (k 1 , k 2 , and k 3 ) probed by these spin-locking relaxometry measurements. As also illustrated by Equation (6), the three apparent velocities (k 1 , k 2 and k 3 ) vary as a function of the sediment orientation which modulates the residual quadrupolar splitting (ω Q ); see Equation (5). As a consequence, it becomes possible to investigate a broad dynamical range by using a limited number of irradiation powers ω 1 (see Figure 9). Figure 9 exhibits a transition between a plateau, at low angular velocities, and a power-law decrease, at high angular velocities [80]. By using the water self-diffusion coefficient of the confined water molecules previously determined by QENS measurements under equivalent conditions (see Section 2.2), numerical simulations of Brownian Dynamics are performed to determine the average residence time of water molecules confined within the interlamellar space between two beidellite platelets. That residence time fully corresponds [43] to the critical time τ C evaluated from the measured angular velocity ω C (τ C = 1/ω C ; see Figure 9). This result illustrates how multi-quanta spin-locking NMR relaxometry of quadrupolar nuclei largely extends up to ∼ 10 −4 s, the time-scale previously probed by QENS measurements.  The same spin-locking relaxometry measurements were performed by 7 Li NMR spectroscopy to study the mobility of lithium counterions within aqueous dispersions of synthetic laponite clays [10]. Dense dispersions are prepared by uni-axial oedometric compression; more details are given in the Supplementary Materials. Above the Isotropic/Nematic phase transition [5], the 7 Li NMR spectra exhibit a splitting of the resonance line (see Figure 10) that perfectly cancels at the so-called magic angle; see Equation (5). This behavior results from the occurrence of a single nematic phase, with a director parallel to the compression axis. Since 7 Li is a I=3/2 spin nucleus, six different angular velocities may now be probed by using 7 Li NMR relaxometry with a single irradiation power ω 1 [84] (see Figure 11): As illustrated in Figure 12, a broad range of angular velocities is then investigated by using a limited number of angular velocities. Figure 12 also exhibits the dispersion curve of the NMR relaxation rates. Here, we also detect a transition between a low-frequency plateau and a high-frequency continuous decrease. Numerical simulations of Brownian Dynamics (BD) are then required to interpret these results. In order to reproduce ionic condensation on the clay surface, a mean force potential [73] is evaluated from preliminary Monte Carlo simulations of the distribution (see Figure 13) of thousands Lithium counterions neutralizing the negative charge of an isolated laponite platelet. By using the weak overlap approximation [68], BD simulations are further performed for modeling lithium diffusion within a collection of partially oriented laponite platelets. As illustrated in Figure 14, the decorrelation of the memory function, induced by the fluctuations of the quadrupolar relaxation mechanism [41], evolves according to a t −1 power law for diffusion time longer than 10 2 ns. This behavior results from the strong electrostatic coupling between the highly charged clay platelet and the lithium counterions trapped in the so-called condensation layer (see Figure 14). This phenomenon is responsible for the low frequency logarithmic decrease of the 7 Li NMR relaxation rates displayed in Figure 12. As a consequence, two different confinements may be detected by NMR relaxometry within such solid/liquid interfacial systems, induced either by geometrical constraints, as in the case of water molecules adsorbed within the interlamellar space of dense sediments, or by energetic coupling, as in the case of ionic condensation within dilute clay dispersions.

Pulsed-Gradient Spin-Echo Attenuation
Attenuation of the NMR signal measured after an echo pulse-sequence encoded by using gradients of the static magnetic field leads to dynamical information on the mobility of the NMR probes in the direction parallel to the applied field gradient [87]. These attenuation measurements directly probe the Intermediate Scattering Function; see Equation (1) in the Supplementary Materials for wave-vectors q defined by the relationship: where γ is the gyromagnetic ratio of the NMR nucleus (2.6752 × 10 8 rad/s for 1 H and 1.037×10 8 rad/s for 7 Li), and δ is the duration of the applied field gradient and G its strength; see Figure S2 in the Supplementary Materials for more details.
Because of the maximum strength of the applied field gradient (typically, G max = 1.6 T/m) and its duration (typically, δ = 500 µs), the maximum resolution defined as resolution max = 1/q max is not better than 15 µm, thus corresponding to macroscopic mobility measurements. Since the magnetic field gradient may be applied along any direction within the sample, 1 H PGSE NMR measurements were first performed to detect the Isotropic/Nematic phase transition of laponite aqueous dispersions [11,50,88] by measuring the anisotropy of the water self-diffusion tensor. These measurements were possible thanks to the great anisotropy of the laponite platelets (diameter∼300 Å, thickness∼10 Å). As illustrated in Figure 15, the Isotropic/Nematic transition of these laponite suspensions does not vary as a function of the neutralizing monovalent counterions and is well reproduced by numerical simulations of Brownian Dynamics mimicking water diffusion within a collection of partially oriented hard disks [50]. The same PGSE measurements were performed by 7 Li NMR spectroscopy within aqueous dispersions of laponite platelets neutralized by lithium counterions. By contrast with the continuous decrease of the measured water mobility normalized by the bulk water mobility, Figure 16 exhibits a strong initial decrease of the lithium mobility by reference to its value measured within simple LiCl aqueous solutions. This phenomenon is the fingerprint of the ionic condensation of the lithium counterions at the surface of the highly charged laponite platelets. Numerical simulations of Brownian Dynamics may reproduce that phenomenon by including an effective mean force potential as described previously (see Figure 13). While these BD simulations agree perfectly with experimental data at low clay concentrations (see Figure 16), large differences appear at higher concentrations because of the limitation of the weak overlap approximation [68] implied in these BD simulations [50].
Because of the time scale required to build these magnetic field gradients (typically, 60 µs), such PGSE measurements become inefficient to measure the mobility of water molecules or counterions confined within the interlamellar space of clay aggregates due to the corresponding enhancement of their transverse relaxation rates. A last illustration of the use of 1 H PGSE concerns the analysis of anisotropy of water diffusion in compacted clay samples made of micrometer-sized kaolinite particles [51,52]; see the Supplementary Materials for more details. Since kaolinite is unswelling clay, these measurements are appropriate to detect the macroscopic mobility of the water molecules diffusing around the clay particles, leading to the so-called tortuosity factor [89,90]. In addition, this mobility is also strongly impacted by the preferential orientation of lamellar particles in the system (see Figures 18 and 19) [52]. These properties can be quantified based on X-Ray Scattering technique on bi-dimensional detector plate (2D-XRS) (see Figure 17) and the analysis of intensity distribution along the diffraction ring of the 001 reflection [52]. This measurement allows extracting the order parameter P 2 through the second order Legendre polynomial function [91], a parameter ranging from 0 for isotropic organization of particle to 1 if particles are perfectly oriented. The influence of preferred orientation for a constant porosity value can reveal a reduction of water mobility by a factor 2 when increasing anisotropy in particle orientation (see Figure 18). Such anisotropy in water mobility is well highlighted using 1 H PGSE-NMR technique and can be quantitatively reproduced by DB simulations using virtual porous media mimicking the size, shape, and orientation of the particles in the porous media (see Figure 18). By accounting for the self-diffusion of water probes in the interparticle volume and elastic collisions with particle interfaces, the simulated data are found to well reproduce experimental measurements (see Figures 18 and 19), thus allowing to decouple the contribution of porosity and particle orientation on the overall water diffusion process.  Dynamics are used to simulate the water self-diffusion tensor within virtual porous media, shown as inserts, mimicking these clay samples (see text). Reprinted from Ref. [52]. Copyright (2020), with permission from Elsevier. Figure 19. Influence of the preferred orientation of particles P 2 on the water self-diffusion tensor measured (red) and simulated (blue) within a collection of clay samples with the same porosity. The z-direction is perpendicular to the preferred alignment of clay platelets. Reprinted from Ref. [52]. Copyright (2020), with permission from Elsevier.

Magnetic Resonance Imaging
Magnetic Resonance Imaging (MRI) experiments were also performed by using 1 H NMR spectroscopy to obtain long time-scale information on the macroscopic mobility of water molecules within clay porous media made of kaolinite particles [51]; see the Supplementary Materials for more details. The experimental conditions are selected to probe a total length of 25 mm. The porosity of the dry kaolinite sample (φ dry = 0.52) is evaluated by helium pycnometry [51]. The clay sample is initially saturated with bulk water, leading to a total length (noted L 1 ) of 12.3 mm (see Figure 20a). The MRI experiment starts just after the addition of a limited amount of heavy water (length L 2 = 7 mm) on the top of the hydrated clay sediment, leading to a total length (L = 19.3 mm) fitting the experimental window. As illustrated in Figure 20a, the 1 H water concentration profiles exhibit some attenuation of the 1 H NMR signal at the bottom of the clay sample because a significant section of the clay sample (5.3 mm) is located outside the optimal detection area of the detection coil. We also detect alterations of the water concentration profile near the clay/heavy water interfacial area, further limiting the useful length of the hydrated sample to 5 mm. By integrating the 1 H NMR signal over the corresponding limited area (see Figure 21), we obtain information of the long time-scale, τ ≈ (40 ± 2) × 10 3 s, implied in this exchange between bulk and heavy water. Numerical simulation are performed [51] to fully validate that analysis (see Figures 20b and 21) by using the porosity of the dry clay sample and the water mobility previously measured by 1 H PGSE attenuation measurements (see Section 2.4 above). That result clearly illustrates how MRI experiments may be used as complementary tool to probe the mobility of confined molecules over a broad macroscopic time-scale (from minutes up to day).   Figure 21. Direct comparison of the time evolutions of the measured and simulated numbers of confined water molecules. Experimental and numerical data are fitted by a biexponential law: where m i are constants. Reprinted with permission from Ref. [51]. Copyright (2018) American Chemical Society.

Macroscopic Diffusion Experiments
Out-Diffusion experiments were performed to determine the time-scale characterizing the exchange between Ca 2+ , initially located in the interlayer of vermiculite disks and neutralizing its charge, and Na + or Sr 2+ , present in aqueous solutions [53,54] (see Figure 22 and the Supplementary Materials for more details). Firstly, the resulting exchange rates between these counterions were well predicted by BD simulations using self-diffusion coefficients for interlayer and aqueous bulk cations, calculated from preliminary MD simulations and the geometry of both aqueous reservoir and vermiculite disks (see Figure 22a) in the case of the Ca 2+ /Na + exchange [53]. The results of this study showed that the time-evolution of the concentration gradient of interlayer Ca 2+ was the main driving force interpreting exchange rates measured at all aqueous reservoir salinity investigated, while limitation at the solid/liquid interface, due to diffusion in the aqueous reservoir, can acts as a limiting factor to interpret exchange rates measured at low salinity conditions. As displayed in Figure 22b, the same experimental data were also well reproduced by a classical finitevolume model based on continuous differential equations, constrained by a self-diffusion coefficient of the slowest interlayer cation. That approach uses also selectivity coefficients obtained independently in batch and describing the ionic exchange at thermodynamic equilibrium (see Ref. [54] for more details). Using the two numerical approaches (BD and classical finite-volume model based on continuous differential equations), data were always based on self-diffusion coefficients obtained from preliminary numerical simulations of Molecular Dynamics (MD) without the use of any fitted parameter [53,54]. As displayed in Figure 22, the agreement between experimental and numerical data validates the use of self-diffusion coefficients obtained from MD simulations to predict the dynamic of the exchange between interlayer cations and aqueous ones over several days (up to 20 days) and then the multi-scale approaches proposed. As a second illustration of long-time diffusion process, Through-Diffusion experiments are commonly used to measure the flux of a tracer, traversing a compacted sample and to extract an effective diffusion coefficient. Such effective diffusion coefficient is directly linked to the pore diffusion parameter commonly derived from NMR measurements probing only the porosity accessible for the investigated probe. A large number of experiments were performed in literature to investigate the influence of porosity on the diffusion of water and ions in porous media made of clayey particles [69,[92][93][94]. Water diffusion coefficients obtained using Through-Diffusion experiments and PGSE-NMR analyses performed with compacted kaolinite samples, characterized by the same porosity and preferential particle orientation, revealed similar results (see Figure 23), despite the different time scale probed by the two methods (i.e., ms for PGSE-NMR versus several days for Through-Diffusion [52]). This finding provides interesting information regarding representative time and length scales for the investigation of water dynamics in clay porous media. Moreover, this consistency allows taking advantage of both techniques for the study of clay/water interfaces. On the one hand PGSE-NMR cannot be used when sample contains paramagnetic elements, whereas Through-Diffusion is insensitive to the chemical composition of the sample. On the other hand, PGSE-NMR provides the full tensor of water dynamics whereas Through-Diffusion experiments requires to be repeated when changing the direction of diffusion probed. For the analysis of clay porous media with both interparticle and interlayer porosity (dual porosity media), Through-Diffusion was successfully applied on vermiculite to investigate the role played by the different types of porosities and preferential orientation of particle on water dynamics [55,94]. Among the large set of different numerical methods to interpret the experimental data [69,95,96], Asaad et al. [55] used BD simulations on representative virtual porous media to mimic the diffusion of water probes in both interparticle and interlayer volume, using self-diffusion coefficients extracted from MD simulations. The good agreement between experiments and simulations (see Figure 23) allowed highlighting the marginal role of interlayer volume on the flux of water tracer traversing samples, and the similar role of preferential orientation of particle when dealing with single or dual porosity media. Figure 23. Influence of the particle orientation and interlayer porosity on the apparent mobility of water molecules (i.e., effective diffusion coefficient) diffusing through clay porous media in the direction perpendicular to the preferred orientation of the clay particles. Data are normalized to the self-diffusion coefficient of water in bulk solution. The data obtained for kaolinite (single porosity of ∼0.26) by 1 H NMR PGSE attenuation (empty yellow triangle) and TD experiments (full yellow triangle) are taken from Asaad et al. [55] and Tertre et al. [94], respectively. TD data obtained for vermiculite (dual porosity) for a similar interparticle porosity are shown as full red circles [55]. Experimental data are compared with results from BD simulations (empty squares) taking into account interlayer diffusion (orange symbols) or only the mobility in the interparticle volume (black symbols) [55]. Reprinted from Ref. [55]. Copyright (2021), with permission from Elsevier.

Conclusions
This review focuses on dynamical properties of the water molecules and neutralizing counterions confined within clay particles. By exploiting a large panel of complementary experimental studies, it becomes possible to investigate a large time scale varying between pico-second and days, implying local motions, as well as macroscopic displacements, within the clay aggregates. In addition to these experimental studies, we also performed multi-scale modeling of these clay/water interfacial systems in order to better extract dynamical information from the raw experimental data and further validate the conclusions of our analysis. This allowed deriving a complete bottom-up analysis of the dynamical properties of ions and water molecules confined within dense clay sediments. In that context, the interpretation of such complex dynamical phenomena is based on numerical modeling exploiting realistic models and dynamical information (residence time, local mobility) extracted from preliminary measurements performed on a smaller scale. Such complementary approaches should be useful for the analysis of other solid/liquid interfacial systems, such as cementeous and zeolitic materials, porous silicates, membranes and micelles, and aqueous dispersions of macromolecules and biopolymers.
Supplementary Materials: This Supplementary Material details the materials and the experimental methods presented in this review. The following are available at https://www.mdpi.com/article/10 .3390/colloids5020034/s1. Figure S1: (a) X-ray diffraction measurements of the swelling of montmorillonite neutralized either by Na + (blue) or by Ca 2+ (red) counterions and dispersed respectively in NaCl or CaCl 2 aqueous solutions. The data are issue from Ref. [1]. (b) Structures of kaolinite and swelling clays (montmorillonite, laponite, saponite, beidellite and vermiculite). Figure S2: Schematic view of the pulse sequence used to perform NMR Pulsed-Gradient Spin-Echo attenuation measurements in order to probe the influence of clay on the water and Li + ions mobility. The red and blue arrows illustrate the time domains corresponding to the transverse (T 2 ) and longitudinal (T 1 )