Bridging molecular and continuous descriptions: the case of dynamics in clays

The theory of transport in porous media such as clays depends on the level of description. On the macroscopic scale, hydrodynamics equations are used. These continuous descriptions are convenient to model the fluid motion in a confined system. Nevertheless, they are valid only if the pores of the material are much larger than the molecular size of the components of the system. Another approach consists in using molecular descriptions. These two methods which correspond to different levels of description are complementary. The link between them can be clarified by using a coarse-graining procedure where the microscopic laws are averaged over fast variables to get the long time macroscopic laws. We present such an approach in the case of clays. Firstly, we detail the various levels of description and the relations among them, by emphasizing the validity domain of the hydrodynamic equations. Secondly, we focus on the case of dehydrated clays where hydrodynamics is not relevant. We show that it is possible to derive a simple model for the motion of the cesium ion based on the difference on time scale between the solvent and the solute particles.


INTRODUCTION
Clays are multi-scale materials. The description of the dynamics of such complex systems is to be made at different scales. On the macroscopic scale, hydrodynamics equations are very convenient to model the fluid motion and the diffusion of the various species through this porous medium. Most of the time, the Navier-Stokes equation is coupled to the Poisson equation and the Fick's law in order to compute the macroscopic properties of the system as a function of the geometry. Yet the link among these macroscopic laws and the microscopic properties of the material has to be clarified. They are valid only if the pores of the material are very large, but the exact condition has to be determined. In addition, macroscopic parameters (zeta potential, effective viscosity, diffusion coefficients, adsorption coefficients) are required.
Further approaches are based on molecular descriptions. With the development of computer science, it is possible to perform molecular dynamics simulations in the case of such a confined system. The precise motion of the particles is obtained at short time and short length scales. In that case the difficulty arises from the fact that the boundary conditions of the simulation do not correspond to the exact experimental system since long-scale porosities are ignored and the link with the experiments has to be clarified.
Thus the two methods are complementary since they allow to describe the complex structure of the An Acad Bras Cienc (2010) 82 (1) material. At short scales molecular dynamics can be used whereas macroscopic laws are required at large scales. The link between the descriptions still requires further developments. Rigorous strategy can be obtained if one uses coarse-graining procedures where the microscopic laws are averaged over fast variables to get the long time macroscopic laws. We present such approaches in the case of clays.
First, we will focus on hydrated clays. We will briefly introduce the various levels of description that can be used in order to describe montmorillonite clays. Comparison between molecular dynamics and macroscopic laws allows the determination of the validity domain of commonly used theories, such as the Poisson-Boltzmann equation. In the case of hydrodynamics, the comparisons show that the macroscopic laws are valid if slip boundary conditions are taken into account.
Second, the case of ionic diffusion for dehydrated clays is analyzed. In that case, the classical equations (e.g. Nernst-Planck model) are not valid any more because of the size of the pore. Continuous solvent models based on bulk equation are not valid because of the layering of the solvent, which depends on its molecular structure. We have introduced a coarse-graining procedure, which generalizes the latter to the case of confined systems. We take into account the molecular structure of the water to compute (i) the interactions between the particles and (ii) their mobility. These preliminary calculations show that it is possible to describe the motion of the ions in a continuous solvent model, even in the case of dehydrated clays, provided that the mass of the ion is much larger than that of the solvent (cesium ion e.g.).

LEVELS OF DESCRIPTIONS
Hydrated smectite clays such as montmorillonite are lamellar mineral crystals composed of charged layers separated by an aqueous solution. They exhibit special features towards hydration and ion fixation (Karaborni et al. 1996). Clay lamellae form thin platelet-shaped particles of diameter close to several hundreds of Ångströms. However this lamellar geometry is valid only at small length scales. At larger scales, the structure is more complex and it leads to multi-porosities. Thus the de-scription of such systems requires a multi-scale strategy. Many theoretical studies at different levels of description have been performed. They are summarized in Figure 1.
Ab-initio molecular dynamics (Boek and Sprik 2003, Tunega et al. 2004, Benco et al. 2001) provides information about the electronic degrees of freedom, but because of its computational cost it is restricted to the smallest time and length scales. Classical Monte-Carlo (Skipper et al. 1995) or Molecular Dynamics simulations (Malikova et al. 2004, Arab et al. 2004, Greathouse et al. 2000 are able to describe larger systems. For example, the mechanism of crystalline swelling at low hydration is well reproduced by these techniques where the various atoms and ions are considered explicitly. Nevertheless, very large systems (e.g. for high hydration, or when different platelets have to be taken into account) cannot be treated by this technique. The use of alternative methods such as Brownian dynamics simulations ) is inescapable. In that case, only a part of the degrees of freedom (e.g. the motion of the (free) ions) is considered. At large scales, continuous methods (e.g Unfortunately, the link among these methods is at first sight far from obvious. Yet this issue is very important since every method requires the knowledge of the parameters of the model it is based on. These parameters could be obtained -at least in principle -from the simulations at the lower level of description. This self-consistent approach is not commonly used for clays, though. Most of the numerical studies obtain their parameter from experiments or from semi-phenomenological considerations. Only a few attempts have been made towards that direction for clays , Porion et al. 2003. The general method, which bridges the gap among the various models, is the elimination of fast variables (Van Kampen 1985). Every level of description can be obtained from the preceding one by averaging over An Acad Bras Cienc (2010) 82 (1)
Pore Averaged Laws Similarly, in the case of Brownian dynamics, the motion of the solute particles (i.e. ions) is obtained by considering that the only relevant variables are the degrees of freedom of the solute particles. This Brownian description which treats the solvent as a continuum is valid only if the solvent molecules are much faster than the solute particles.
The elimination of the fast variables is the key for the derivation of the various models of Figure 1, except for the last one. Indeed, the experimental macroscopic law of porous media (such as Darcy's law) are theoretically obtained differently, by averaging over the geometry of the material, for example thanks to the homogenization procedure of Moyne and Murad Murad 2002, 2003).

HYDRODYNAMIC LIMIT
It should be noted that the hydrodynamic equations themselves correspond to the elimination of fast vari-ables. Thus the laws of non-equilibrium thermodynamics (de Groot and Mazur 1983) such as the Navier-Stokes equation for the convection and the Fick's law for the diffusive motion of the solute particles, are valid in the hydrodynamic limit (Ailawadi et al. 1971, Zwanzig 2001 where the hydrodynamic modes become slow. The latter corresponds to the limit of large distances r or small wave-vectors k, as can be understood from the following argument.
Non-equilibrium thermodynamics predicts the evolution of particle and momentum distributions in space. For example, in the case of diffusion, it gives the evolution of the density of particles in space for the various components of the systems. The microscopic density of a species α reads: (1) where the ith particle of type α has position r i (t) at time t. The Fourier transform is: In the hydrodynamic limit k → 0, ρ α (k, t) becomes time-independent because it tends towards N α , the total number of particles of type α. For k = 0, ρ α (k, t) = N α depends on time. But its evolution gets slower as k → 0. Thus for k → 0, the density of particle becomes a slow variable. A similar argument can be made for the average velocity (or momentum) of the particles: because of the linear momentum conservation law. Thus in the hydrodynamic limit (small wave-vectors or large distances), hydrodynamic laws are valid because the hydrodynamic quantities become slow, on the grounds of the conservation laws.
An Acad Bras Cienc (2010) 82 (1) The accuracy of hydrodynamic descriptions can be checked by comparing their predictions to the results obtained from molecular dynamics Turq 2003, Dufrêche et al. 2005b). Some of the results are given in Figure 2 for Na montmorillonite clays with an interlayer spacing equal to 47.5 Å. In that case, macroscopic hydrodynamic models generally calculate counterion distributions with the help of the Poisson-Boltzmann equation. Even for that relatively small interlayer spacing, this hydrodynamic model appears to be valid. The main differences are the oscillations given by the microscopic model, which are characteristic of the discrete nature of the solvent molecules. The electroosmotic flow can be calculated from molecular dynamics as well. In that case simple hydrodynamic theory fails to reproduce the microscopic results. In fact, it seems that the problem comes from the stick boundary conditions. One can show that microscopically, slip boundary conditions have to be taken into account. More precisely, the parallel component of the hydrodynamic fluid velocity at the interface is not zero. It is actually given by: where z is the direction orthogonal to the interface and δ is the Navier slip length that quantifies slipping. For montmorillonite, δ is equal to a few Ångströms, which is very small, but it is enough to modify the electroosmotic flow. The reason comes from the fact that electro-osmosis is created by surface charges, so that it strongly depends on the boundary conditions. Conversely, hydrodynamic motions created by bulk forces (such as a pressure gradient in the case of a Poiseuille flow or of a streaming potential experiment) depend less on the slip length, as it has been recently confirmed by Joly et al. (Joly et al. 2006).

DEHYDRATED CLAYS
At very low interlayer spacings, hydrodynamic models inside the pores do not hold anymore. For example, in the important case of monohydrated clays, there is only one layer of water molecules between two clay sheets, so there is no chance that hydrodynamic equations work. Yet it is not impossible to derive a continuous solvent model for the dynamics of the ions in that case. Indeed, the hydrodynamic quantities are not slow anymore, because the distance between the layers is not large enough, but if the ions are much heavier than the water molecules, their dynamics is much slower than the dynamics of the solvent. Thus, the position and the velocity of the solute particles are slow variables and a coarse-graining procedure can be performed. The motion of the ions is given by a Brownian model. The diffusion of heavy ions in clays at low hydration is particularly important because of the practical applications.
In particular, it corresponds to the case of nuclear waste storage.

BROWNIAN MODEL
Most theories of transport of electrolyte solutions, which generalized the pioneering work of Onsager (Onsager and Fuoss 1932) are actually based on Brownian description. At first sight it may be surprising that Brownian theories can be applied for simple ions that are not quite heavier than the solvent molecules. However, in practice, for bulk solutions, the model is satisfactory even for small ions as Li + , probably because of the phenomenon of hydration. A shell of solvent molecules becomes strongly attached to each ion so that it is probably not too bad an approximation to consider a given ion together with its solvation shell as a single heavy entity. Different Brownian models can be implemented for the dynamics of ions in an aqueous solution. They are summarized in Figure 3. At the microscopic scale, discrete solvent descriptions consider all the different particles, the solvent and the ions. The dynamics (Hansen and McDonald 1986) is driven by Newton's law (which corresponds to the Liouville's equation in phase space). If solute particles can be considered to be much heavier than the solvent molecules, Brownian models (based on continuum solvent descriptions) are valid. There are mainly two classes of Brownian description (Résibois 1968). The first one corresponds to the Langevin approach. The relevant observables are the velocity (or the momentum) and the position of the ions. The evolution in phase space is given by a Fokker-Planck equation, as it can be shown by a Kramers Moyal expansion (Risken 1996, Van Kampen 1992. The second one is the Smoluchowski approach. It is valid for time scales longer than the relaxation time τ v = 1/γ of the velocities of    (Altenberger and Friedman 1983, Sung and Friedman 1983, Dufrêche et al. 2005a). The reason is that the characteristic time of ions in bulk solution is the Debye time of the electrolyte (Turq et al. 1992). This time is typically of a few nanoseconds, it is much longer than the relaxation time of the velocities (closer to the picosecond) so that the Smoluchowski approach is valid.
In the case of monohydrated clays the situation is completely different because of the presence of the mineral layers. The characteristic time of the motion of the ion is not the Debye time, because the ion diffusion is actually driven by the surface-ions interactions (averaged over the solvent configuration). Molecular dynamics simulations (Malikova et al. 2004) show that the dynamics follows a hopping mechanism (site-site diffusion), with a characteristic time close to the picosecond, which is of the order of magnitude of τ v . Thus, contrary to bulk solutions, the Smoluchowski's equation is not valid. The correct Brownian model is based on the Langevin level of description.
There is another price to pay compared to the bulk case. The site-site diffusion is coupled to the hydration/dehydration of the water. Thus, it is not possible anymore to consider the shell of water molecules as a part of a Brownian particle, because the solvation changes during the diffusion process. As a matter of fact, Brownian theory can only be applied only to real heavy ions. The case of small ones cannot be studied from Brownian models, because the time scale of the hydration process is close to the time scale of diffusion.

LANGEVIN DYNAMICS SIMULATION
We applied a Langevin description to the case of a clay with Cs + counterions and a water monolayer. The separation of timescales between the slow solute (ion) and fast solvent (water) dynamics is such that one can An Acad Bras Cienc (2010) 82 (1) average over the solvent degrees of freedom. Within the Langevin approach, the solute/solvent interaction is modelled by two parameters. The dynamical properties depend on a single friction coefficient γ (in s −1 ): The friction force on an ion is proportional to its velocity (F solv = −mγ v). The important point is that γ is different from the bulk values because the solvent configuration is completely modified by the geometry in the monolayer. In particular, the orientational configuration of the water molecules is modified by the sheets so that the dielectric part of the friction (Hubbard andOnsager 1977, Bagchi 1998), which is an important contribution to the bulk friction, is completely different. The second parameter is the external (free energy) potential V e f f which models the solute/surface interactions.
The Langevin equation which gives the evolution of the velocity v is where F is the external force, i.e. the opposite of the gradient of V e f f . (t) is a random force. It can be modelled by a white noise whose correlation function is related to γ (Risken 1996). We solved this equation numerically using a first order algorithm. Such a numerical solution method is called the Langevin dynamics simulation method.
The Langevin equation is equivalent to the Fokker-Planck (FP) equation which gives the evolution of the probability density function f (x, v, t) of finding an ion at a given position with a given velocity. We also solved numerically the FP equation by a Lattice Fokker-Planck method (Moroni et al. 2006) which we used in this particular context of ionic mobility in compacted clays . The two methods (Langevin dynamics and Lattice Fokker-Planck simulations) were found to give exactly the same results. An important issue is the value of the McMillan potential V e f f (McMillan and Mayer 1945), averaged over the configurations of the solvent. It is exactly defined as the N particles distribution function of the solute when it is infinitely dilute. For the sake of sim-plicity, we performed preliminary calculations where V e f f is given by the low density limit: with ρ(r) is the average concentration of counterions obtained by molecular dynamics. This simplified model neglects the correlation among the ions.
V e f f , which was assumed to be a 2D potential is represented in Figure 4. This potential clearly indicates that we have a site-site diffusion for the dynamics of the cesium ions, the depth of the potential being a few times the thermal energy k B T . Langevin dynamics gives exactly the same ion distribution (not represented). It is a direct consequence of the fact that the canonical distribution function is a stationary solution of the Fokker-Planck equation.

DIFFUSION OF IONS
We calculated the velocity correlation function Z (t) (Hansen and McDonald 1986) which characterizes the diffusion of the ions in the system for different values of the friction γ . The results are presented in Figure 5.
At large friction γ = 10 ps −1 , Z (t) can be divided into two parts. First at short times a simple relaxation corresponds to the friction of the solvent. Then there An Acad Bras Cienc (2010) 82 (1) is a long time negative tail which corresponds to the jumps from site to site. If the friction decreases and approaches γ = 1 ps −1 , oscillations appear in the velocity correlation function. They correspond to the oscillations of the ions in the layer sites. They are typical of the Langevin regime for which the velocity relaxation is not faster than the characteristic time of the interactions. The two times are close to 1 ps. The results obtained from Langevin dynamics can be compared to the Z (t) calculated from molecular dynamics. Of course, the velocity correlation function is quantitatively different because of the simplicity of this preliminary model: the friction is isotropic and the ion correlations are neglected. However the curves are qualitatively similar. Discrete solvent simulations recover small oscillations. The value of the friction can be obtained at small time scales. γ seems to be between 2 and 5 ps −1 . It corresponds to a domain where Z (t) starts to exhibit very small velocity oscillations.

CONCLUSION
We briefly describe the different approaches used to model the dynamics of clays. Every level of description is valid if the variables it models (relevant variables) are slow variables. In the case of hydrodynamics, the macroscopic laws are valid as long as the typical size of the system is larger than a few nanometers, if proper boundary conditions are taken into account. For monohydrated systems where there is just one layer of solvent molecules it is still possible to use a continuous solvent model if the ion is heavy. Then the dynamics is obtained from a Brownian model. Preliminary results indicate that this method may be applied for Cs + . Palavras-chave: argila, modelagem multiescala, dinâmica molecular, hidrodinâmica, movimento Browniano.