B Rotenberg

J.-F Dufrêche

P Turq

Frequency-dependent dielectric permittivity of salt-free charged lamellar systems

We present a new model to analyze dielectric spectroscopy measurements on charged lamellar systems, with the following improvements with respect to the hitherto available models: (i) it does not rely on the hypothesis of local electroneutrality, and allows to treat the salt-free case; (ii) the chemical exchange governing the partition between free and bound ions is properly taken into account; (iii) a fully analytical solution is provided. The variation of the frequency-dependent dielectric permittivity with both thermodynamic and kinetic characteristics of the free-bound ion equilibrium is presented. In particular, the relative weights of both relaxation modes (exchange and transport), and their characteristic frequencies are discussed. This study opens the way to the analysis of systems for which the usual models are irrelevant, such as salt-free clay gels or membranes.

Introduction

The dielectric properties of heterogeneous materials consisting of inclusions of a conducting medium in an insulating matrix have attracted the attention of scientists for more than a century [START_REF] Maxwell | Treatise on Electricity and Magnetism[END_REF]. Early theories, known as Maxwell-Wagner [START_REF] Wagner | [END_REF] and later Maxwell-Wagner-Sillars [START_REF] Sillars | Proc. Inst. Electr. Eng[END_REF] theories were based solely on macroscopic electrodynamics. Later on, considerable effort has been devoted to developing theoretical models based on a microscopic description of the charge carriers in the conducting inclusions, and the complex permittivity could be related with some success to the motion of charge carriers. In particular, its frequency dependence was derived from the dynamics of charge relaxation, either directly in the conducting medium, often an electrolyte consisting of anions and cations dispersed in a solvent [START_REF] Debye | [END_REF]5], or from the relaxation of charges accumulated at its boundaries [6,7,8,9,10,11,12,13,14,15]. The relaxation of the ionic cloud in an electrolyte was first described by Debye, and its peculiar frequency dependence is now known as the Debye-Falkenhagen effect [START_REF] Debye | [END_REF]. More recently, our understanding of the frequency dependent conductivity of bulk electrolytes has been improved by the use of Mode Coupling Theory [5]. In solid ionic conductors too, the concepts of ionic atmosphere and of its relaxation seem to be appropriate to interpret the frequency-dependent conductivity [16]. Charge accumulation at boundaries in presence of an applied AC electric field was largely discussed in various contexts such as materials between (blocking or partially non-blocking) electrodes [6,7,8], saline inclusions in rocks [9] or lamellar systems [10,11,12], particularly biological membranes.

Though often described in terms of equivalent electrical circuits, the frequency dependent conductivity of electrolytes, colloidal suspensions or emulsions, can be efficiently interpreted on a more physical basis at the Smoluchowski level by assuming that the transport of the charged species is governed by the Nernst-Planck equation [13,14]. The pioneering work of Trukhan [15] in the context of conducting inclusions in an insulating material, both in planar and spherical geometries, opened the way to the study of lamellar systems [10,11,12]. The resolution of the problem is in the previous cases greatly simplified by assuming that electro-neutrality is locally satisfied, in the absence of the applied field. Space charges then arise only near the electrodes, under the influence of the field. This picture is only justified when the Debye screening length, characterizing the range of electrical interactions, is much smaller than the system size. In the electrochemical context, this is achieved by the presence of a "supporting electrolyte". Even in presence of an electrolyte, recent studies have shown that under an alternating field, electro-neutrality could be perturbed over distances much larger than the Debye length if the diffusion coefficients of anions and cations are different, in dramatic contrast with the results for a DC field [17,18]. The authors were able to provide a fully analytical solution for insulating spheres suspended in an electrolyte [17], and numerical results for charged particles [18]. Two important conclusions are: (i) AC fields can induce violations of electro-neutrality over large distances, resulting in substantial contributions to the dielectric behavior of the suspension, even for systems that are neutral at equilibrium; (ii) the effect can be up to five orders of magnitude larger for charged particles, i.e. when a space charge is present at equilibrium.

From the theoretical point of view, all investigations of the polarization (either at electrodes or of lamellar systems) have assumed so far that the system under consideration was locally electro-neutral. In many systems however, when no "supporting electrolyte" is added, charges arising from the dissociation of surface sites, either acidic or cation compensated sites (such as oxygen atoms of mica or clay sheets, or the hydrophilic head of an ionic surfactant) can be found in the vicinity of the charged surfaces, mobile in the solvent which thereby becomes a region of space charge, and an internal electric field arises. Electro-neutrality is in that case only guaranteed by the presence of the (oppositely) charged surfaces. The equilibrium profile of the counter-ion distribution depends on the charge density of the surface, and is determined by the compensation of diffusion (due to concentration gradients) and migration (due to the internal electric field): it leads to the well-known Poisson-Boltzmann solution [19,20]. The dynamical response to an applied field of this charge distribution was apparently never studied, though it is necessary if one wants to investigate the permittivity of such systems, which contain no added salt. A possible reason for this absence of consideration might be found in the difficulty of the task. Indeed, the coupling of this distribution with an electric field does include a contribution which is absent for locally neutral system, namely that of the internal electric field "at rest".

On the experimental side, many systems belonging to this category (no added salt) have received over the last decades a growing interest. In particular, the dynamics of cations in clays was recognized as a crucial point because these materials are considered as possible (engineered and geological) barriers for the long term storage of nuclear wastes disposal [START_REF] Andra | Tome 1: Contexte et objet[END_REF]. Among other studies, dielectric spectroscopy has proved very promising.

However, there is a need for a theoretical background to interpret the experimental data and extract microscopic information from these macroscopic measurements. The hitherto available models are restricted to the case of added salt and are therefore not relevant for this purpose.

Another important feature of the ionic dynamics, which concerns all the above-mentioned cases, is the possibility of adsorption of the ions at the surfaces. This adsorption is generally ion-specific (not only charge-dependent), and both the chemical equilibrium at the surfaces can influence the dynamics of the charge relaxation in the liquid phase. Indeed, in bulk electrolytes, it has been shown that chemical association between ions could influence the charge dynamics [22], while relaxation of the ionic cloud do influence diffusion limited reactions [23]. Similarly bound (to a surface) and free ions don't have the same mobility (the latter being zero for bound ions), so that the exchange between both forms affects the overall dynamics of the mobile ions. The influence of such an exchange on dynamics has been reported in the case of water in clays [24], and especially on dielectric properties for exchangeable polar molecules in general [25] and for water bound to biological molecules [26] in particular. The exchange of cations is something different, but a parallel can be drawn with these studies, since it also affects the dynamics, and hence the dielectric properties. In micellar solutions, the exchange between cations bound to the surface (immobile) with others diffusing around the micelles was shown to lead to a complex dielectric response [27,28]. In turn, one can ask if it is possible to extract the thermodynamic (distribution constant between bound and free ions) and dynamic (kinetic constants of the adsorption and desorption reaction) characteristics of the exchange reaction from the dielectric permittivity of the whole system, which is a probe of the global charge dynamics.

The scope of the present paper is to show that despite the complexity of the problem, a (relatively) simple analytical solution can be found in the case where charged species in the liquid phase stem only from dissociation of the surfaces, and when a chemical equilibrium dictates the partition between free and bound ions. The solution is presented and the influence of structural and kinetic parameters on the complex permittivity is discussed. The possibility of extracting microscopic information from its experimental determination is also considered. The paper is organized as follows. In section 1, the microscopic model is presented. In section 2, the equilibrium concentration profile and internal electric field are recalled. Section 3 presents the fully analytical solution for the complex permittivity of the system; the resulting expression, particularly the influence of the various parameters, is discussed in section 4. Finally, section 5 summarizes the major findings of the present work.

Microscopic model

We only consider the ionic motion normal to the charged surfaces. Within the meanfield treatment adopted here, there is no concentration gradient parallel to the surfaces, and therefore no diffusive flux (nor internal electric field) along the planes (E = 0). We restrict ourselves to the homo-ionic case (i.e. only one type of counter-ions), with counterions characterized by their charge z and their diffusion coefficient D. We also restrict ourselves to the case of no added salt: the only ions between the charged surfaces come from the dissociation of surface sites. Electro-neutrality is guaranteed by the presence of the charged surfaces.

The ionic distribution is described in terms of a two-state model: ions can be either adsorbed on the surfaces (Stern layer), or mobile in the interlayer space. Such a distribution is depicted in figure 1. This model accounts for the possible formation of inner-sphere complexes at the surface, where some solvent (water) molecules of the ion's coordination shell are replaced by atoms of the charged surface. Such a distribution is characterized by the concentration profile c(x, t) (in m -3 ), and the adsorption on both surfaces Γ ± (t) (in m -2 ). Γ + (resp. Γ -) refers to the surface at x = +L/2 (resp. x = -L/2), where the origin of the x axis normal to the surfaces has been taken in the middle of the interlayer space (L is the interlayer distance).

A particular ion located near a charged surface (at x = ±L/2), can switch between the two states following first order kinetics, i.e. the number of mobile ions adsorbing per unit surface and per unit time is k 1 c(±L/2, t) (k 1 in m.s -1 ), whereas the number of adsorbed ions desorbing per unit surface and per unit time is k 2 Γ ± (t) (k 2 in s -1 ). The prime is here to underline that k 1 and k 2 don't have the same dimensionality. At equilibrium, we simply have k 1 c 0 (±L/2) = k 2 Γ 0 . The relaxation towards chemical equilibrium has a time constant:

τ χ = 1 k 1 /L + k 2 (1)
One can define the distribution constant K d as the ratio between the number of bound ions per unit volume vs. mobile ions per unit volume.

K d = bound mobile (2)
With this definition, K d is not the equilibrium constant of the exchange reaction. We will discuss later the relation between K d and the kinetic constants. The bare charge of the surface is -σ 0 (without loss of generality we consider that the counter-ions are cations, and thus σ 0 > 0), and their effective charge -σ = -σ 0 + zΓ 0 = -σ 0 /(1 + K d ). K d will be shown to have a large influence on the dynamics of the charge distribution in the presence of an applied electric field, and hence on the polarization of the interlayer.

The solvent is treated as a continuum of dielectric constant S , and all dielectric relaxation mechanisms associated with water, such as dipolar rotation of water molecules are therefore outside the scope of the present model. Their contribution to the dielectric spectrum should be added to the one that we are calculating here, and the dielectric constant of the whole interlayer is hence larger than that of the pure solvent (see section 4). This assumption makes sense if the relaxation time of the solvent (typical time for a solvent molecule reorientation) is much shorter than the relaxation time(s) of the ions.

Indeed, in that case the two processes correspond to different time scales and can be considered as decoupled. Such an assumption should be verified a posteriori, from the experimental measurements. The value of S is debatable, since for confined water (small L), there is no clear reason for it to be equal to its bulk value. Nevertheless, there has been some evidence that a mesoscopic model such as Poisson-Boltzmann involving the bulk value of S can provide reasonable results, at least as far as the equilibrium profile is concerned, for distances as short as 15 Å [19].

Equilibrium concentration profile

The charge of the surfaces attracts the counter-ions, thereby enhancing their local concentration. This density inhomogenity gives rise to an internal electric field and a diffusive flux. At equilibrium, diffusive and migration (electric) fluxes cancel out. The total flux j(x) is then identically zero:

j(x) = -D dc dx + zD k B T c(x)E(x) = 0 ( 3
)
where E is the electrical field, deriving from the electrical potential V (E = -dV /dx).

k B is the Boltzmann constant, and T the temperature. Furthermore, the local Gauss theorem reads:

dE dx = - d 2 V dx 2 = zc(x) 0 s (4)
Inserting ( 4) in (3), and using the electro-neutrality condition (the total number of ions compensates the total charge of the surfaces), one finds for the concentration profile the Poisson-Boltzmann solution [19,20]:

c 0 (x) = 1 2πL B α 2 cos 2 (αx) (5)
with

2πL B σ z = α tan αL 2 (6)
where L B = z 2 /4π 0 s k B T is the Bjerrum length , and for the internal electric field:

E 0 (x) = z 0 s α 2πL B tan(αx) (7)
These quantities are both symmetric with respect to the x = 0 plane, which means that c 0 (x) (scalar) is even while that E 0 (x) (component of a vector) is odd. The constant α can be interpreted as an inverse screening length, generalizing the Debye inverse screening length κ to the case of a confined electrolyte. The dimensionless parameter αL hence determines the "softness" of the ionic profile (see equation ( 5)). Equation ( 6) implies that for a given bare charge density σ 0 , αL implicitly depends on K d via the effective charge σ. Indeed, it can be rewritten as:

αL 2 tan αL 2 = L 4λ = L 4 σ 0 z 0 S k B T × 1 1 + K d (8)
where one has introduced the (K d dependent) Gouy length λ, defined as the distance from the charged plane at which the electrostatic energy of an ion is equal to the thermal energy (neglecting the influence of the other ions). The ratio L/λ completely defines the parameter αL, and hence the equilibrium profile. 1. It is clear from figure 2b that when L/λ is small, the concentration profile is very flat, while for large values of L/λ it is sharply increasing in the vicinity of the charged planes. Such different ionic distributions will respond very differently to an applied electric field, as will be shown in section 4.

For the same lamellar system (same bare charge σ 0 , same interlayer spacing L, same counter-ion charge z), the specific interactions between the ions and the surface (which could be different for e.g. Na + and Cs + ), as described here by the distribution constant K d , tune the effective charge, and hence the magnitude of the electrical (non-specific) interactions, which in turn controls the whole shape of the density profile.

One can express the concentration at the charged surface c 0 (L/2) as a function of the total amount of mobile ions in the interlayer, by integrating (5). Since at equilibrium we furthermore have k 1 c 0 (L/2) = k 2 Γ 0 , we can use these two relations to link K d to the kinetic constants, with the result:

2k 1 Lk 2 = K d sin αL αL = K d 2πL B L 2 C0 (αL) 2 + (πL B L 2 C0 ) 2 (9)
Where one must not forget that αL implicitly depends on K d . C0 denotes the average concentration in the interlayer; using (6), we have:

C0 = α tan(αL/2) πL B L = 2σ Lz (10)
In the limit of weak electrostatic interactions (αL → 0), either because the bare charge density is small, or because the adsorption is strong (K d 1) -resulting in a screening of the interactions-or finally because the interlayer spacing is very small (see (8)), the concentration profile becomes very "flat", and almost equal to its average value C0 . In this case one simply has α 2 = 2πL B C0 , as expected in the Debye-Hückel limit. Moreover, in this limit equation ( 9) simplifies to

K d = 2k 1 /k 2 L.
The limitations of the Poisson-Boltzmann equation [29] arise from the fact that it neglects ionic correlations beyond a mean-field treatment of the electrostatic interactions.

In particular, the ionic size, which is a fundamental parameter for the description of ion transport in bulk electrolyte solutions [30], is neglected. For low to moderate charge lamellar systems such as clays [31,32], the hard sphere repulsion between the ions adds small oscillations to the Poisson-Boltzmann concentration profile, as can be shown from continuous solvent Monte-Carlo simulations. Yet the latter do not correspond to the oscillations seen in discrete solvent molecular dynamics simulations [19]. In that case the ionic concentration profile still oscillates around the PB curve but the period of the oscillations is different: they are in fact reminiscent of the layering of the solvent. Such an effect could be approximately taken into account by adding an oscillating term to the external potential experienced by the ions. The magnitude of such oscillations is not large, though. Therefore, the Poisson-Boltzmann treatment is probably not too bad an approximation, at least for low to moderate charge surfaces and large interlayer distances.

Steady-state polarization -Permittivity

In dielectric spectroscopy measurements, an alternating electric field of small magnitude is applied to the sample, at a given frequency ω. Ions move under the applied field, but are stopped by the clay sheets. This is strictly analogous to the situation of an electrolyte between two blocking electrodes, excepted that the only ions present in our case are the counter-ions compensating the charge of the surfaces. This motion leads to an excess of ions at the surface (and a lack at the other surface), so that the local chemical equilibrium is no longer satisfied.

After a certain time, a steady state is reached, and all physical quantities evolve with an angular frequency ω, possibly with a phase lag with respect to the applied field E app .

In a typical experiment, the frequencies range from a few Hz to a few GHz (10 9 Hz), so that the propagation time of the electromagnetic interactions L/c 10 -16 s (for L of a few hundreds of nm) is negligible compared to the period of the oscillations of the electric field. For higher frequencies, which are experimentally accessible, one starts to probe the solvent relaxation, in which we are not interested here. The approximation of quasi-stationary state is therefore perfectly justified for computing the electric field.

The presence of the electric field breaks the symmetry with respect to the interlayer mid-plane, and an electrical polarization appears. That is precisely this polarization, which also evolves with the frequency ω, that is measured. For sufficiently weak applied fields, P (ω) is proportional to E app , which allows to define the dielectric constant r by

P (ω) = 0 r E app (ω) (11)
It is commonly split into real and imaginary parts as r = r -i r (where i 2 = -1). The rate of absorption of energy at a given frequency is proportional to r .

We now calculate the polarization as a function of frequency, in the linear response regime. The excess quantities (above c 0 , E 0 and Γ 0 ), are denoted by c 1 , E 1 and δΓ. In complex notation, for a quantity f (x, t) we have:

f (ω) = f 0 (x) + f 1 (x)e iωt (12)
the reference value f 0 being zero for P , E app , δΓ and j, while it is given by ( 5) and ( 7) for c and E; for Γ ± it is simply Γ 0 . One can show that (when the steady state is reached) c 1 is an odd function of x, while E 1 and j are even, and δΓ + = -δΓ -. Furthermore, c 1 , E 1 , δΓ, j and P are all proportional to E app (first order terms). Higher order terms have to be left aside in calculating the permittivity (linear response).

The polarization is the first moment of the charge distribution. In computing the latter, one must not forget the contribution of the charge excess at the surfaces:

P = z δΓ + -δΓ - 2 + z L L/2 -L/2 xc 1 (x)dx (13)
Using ( 4) and the expression of the internal electric field at the boundaries, E 1 (±L/2) = ∓zδΓ ± / 0 s , one simply finds after integrating by parts:

P = - 0 s L L/2 -L/2 E 1 (x)dx (14)
The flux of counter-ions up to first order in E app is (the zeroth order term vanishes):

j(x) = -D dc 1 dx + zD k B T [c 0 (x)E 1 (x) + c 1 (x)E 0 (x) + c 0 (x)E app ] (15)
Conservation of matter in the interlayer reads in complex notation:

iωc 1 (x) = - dj dx (16)
Combining (15) and ( 16) one obtains a differential equation satisfied by E 1 (c 1 is its derivative, modulo a proportionality constant z/ 0 s , see equation ( 4)). One then replaces c 0 and E 0 by their expressions ( 5) and ( 7):

1 α 2 d 2 E 1 dx 2 - 2 α tan(αx) dE 1 dx - iω Dα 2 + 2 cos 2 (αx) E 1 (x) = 2 cos 2 (αx) E app (17)
The characteristic time that appears in (17),

τ D = 1/Dα
τ D = τ diff (αL) 2 (18)
The differential equation is solved using the parity of E 1 , together with the boundary condition satisfied by the evolution of the adsorption:

iωδΓ + = k 1 c 1 (L/2) -k 2 δΓ + (19)
or introducing the dimensionless parameter µ = τ χ /τ diff :

1 K d 1 + iωτ diff µ (αL) 2 1 + K d 2 sin αL αL δΓ + = sin αL 2αL × L c 1 (L/2) ( 20
)
where one has also used (1) and (9). µ characterizes the relative rates of chemical relaxation and diffusion: when µ 1 chemistry is "fast", and the density relaxation is dominated by diffusion (over the shortest distance between the screening length 1/α and the interlayer spacing L), while µ 1 corresponds to a relaxation of the density which is dominated by the chemical reaction at the surfaces.

After some algebra, detailed in Appendix A, one finds for the excess internal electric field:

                               E 1 (x) = - Dα 2 iω E app 1 cos 2 (αx) -ξ(ω)g(αx) g(u) = (γ -1) cos[(γ + 1)u] + (γ + 1) cos[(γ -1)u] cos 2 u γ = 1 - iω Dα 2 = √ 1 -iωτ D = 1 - iωτ diff (αL) 2 ξ(ω) = χ(ω) cos 2 (αL/2) [χ(ω)g(αL/2) + iωτ D sin(γαL/2) sin(αL/2)] χ(ω) = sin 2 (αL/2) + 1 K d 1 + iωτ χ 1 (αL) 2 1 + K d 2 sin αL αL (21)
The excess concentration c 1 can then be obtained by equation ( 4), but it is not necessary for computing the polarization. The frequency dependence of γ is reminiscent of that, unusual, of the Debye-Falkenhagen relaxation effect in bulk electrolytes [START_REF] Debye | [END_REF]. Inserting [START_REF] Andra | Tome 1: Contexte et objet[END_REF] in ( 14), we get the polarization:

P (ω) = 0 s E app 1 iωτ D tan(αL/2) αL/2 -ξ(ω) 1 αL 1 -1 g(u)du (22)
The last integral can be performed analytically, which gives the final result for the dielectric constant r :

r (ω) = s 1 iωτ D tan(αL/2) αL/2 1 -2 ξ(ω) sin γαL/2 sin αL/2 (23)
Equations ( 21) and ( 23) constitute the major findings of the present article, and will be the starting point of the discussion of the influence of the various parameters. We would like to underline here the fact that the (relative) simplicity of the solution, and the existence itself of a completely analytical solution were not obvious at all, considering the non-linearity of the problem.

Discussion

The frequency dependence of the dielectric constant is not trivial: it results from an interplay between electrical relaxation and chemical relaxation at the charged surfaces.

Chemistry is present in (23) only through the parameter ξ(ω) given by [START_REF] Andra | Tome 1: Contexte et objet[END_REF], which also depends on the electrostatics, via the dimensionless parameter αL. The frequency dependence of ξ(ω) is then incorporated in expression (23) which depends apart from ξ only on electrostatics. It is therefore quite difficult to have a clear feeling for what is the exact behavior of the complex permittivity, and one will have to compute it numerically as a function of frequency in order to gain insight in the influence of the thermodynamics (K d ) and on the dynamics (chemical vs. electrical relaxation rate, i.e. the value of µ) on the complex permittivity.

Interestingly, some important limits can be studied analytically. Among them, the static permittivity is of particular interest. It is defined as the low frequency limit of ( 23), and it can be shown to be real and finite, which is not obvious from equation ( 23). The result is simply:

r (ω → 0) = s 2 1 - tan(αL/2) αL/2 × cos 2 (αL/2) 1 + K d sin 2 (αL/2) ( 24
)
It is easily checked that the right hand side of equation ( 24) is positive, whatever the values of K d and αL. Two remarks can be made at that point. First, it does not depend on the chemical relaxation time τ χ , but only on the distribution constant K d , nor does it depend on the electrical relaxation time τ D = 1/Dα 2 , but only on the dimensionless structural parameter αL. Indeed, when the applied field varies with an infinitely low frequency, all relaxation phenomena (chemical and electrical) occur comparatively very fast, and all information about their dynamics becomes irrelevant. Second, equation (24) implies that the contribution of ionic motion perpendicular to the charged surface is systematically less important than (half of) that of the solvent. This very important point will need to be discussed when turning to the comparison with experiments. We would like to draw the reader's attention to the following point: since equation ( 24) gives only the ionic contribution to the static permittivity, the solvent contribution S should be added to give the full static permittivity of the interlayer, which satisfies from the previous discussion: s ≤ interlayer ≤ 1.5 s .

The second important limit is that of "slow chemistry". This limit corresponds to an infinitely long chemical relaxation time, and for a fixed K d (fixed effective charge), it implies that both forward (adsorption) and backward (desorption) reactions are extremely slow. The case of no chemical exchange at all is therefore encompassed in this limit (µ → ∞). The parameter ξ(ω) from equation ( 21) simply becomes:

1 ξ(ω) = cos 2 αL 2 g( αL 2 ) = 2γ cos γαL 2 cos αL 2 + 2 sin γαL 2 sin αL 2 (25)
The dielectric constant then reduces to:

r (ω) = 2 αL 1 iωτ D × γ tan(αL/2) -tan(γαL/2) γ + tan(αL/2) tan(γαL/2) ( 26
)
where γ is still given by ( 21). The polarizability of the interlayer, as described by its complex permittivity ( r = r -i r ) depends greatly on αL, or equivalently on L/λ, as can be seen from figure 3, where r / S has been plotted as a function of frequency (ωτ diff ) and of L/λ (both in logarithmic scales). Within this representation, the usual Bode diagrams ( r vs. log ω) are found along horizontal lines. There is no additional information contained in r (ω), because it can be directly deduced from r (ω) via the Kramers-Kronig relations. Each relaxation appears on the Bode diagram for the imaginary part as a peak at the corresponding frequency.

For small αL, the equilibrium profile is flat, and is not very sensitive to an applied field. Indeed, the permittivity is then:

r (ω) αL→0 ∼ (αL) 2 iωτ diff 1 - tan( √ -iωτ diff /2) √ -iωτ diff /2 (27)
and vanishes as (αL) 2 . One can show that the absorption peak (maximum of r (ω)) is located in this limit for ωτ diff ≈ π 2 (red curve of figure 3), independently of αL, because for low charge densities, the relevant distance is not the screening length but rather the interlayer spacing; hence the relevant timescale is not τ D but τ diff .

At the other extreme (of charge saturation), the ionic distribution is sharply peaked in the vicinity of the charged surfaces. The profile is deformed more easily by an applied field. The permittivity is in this limit:

r (ω) αL→π ∼ 2 iωτ diff √ π 2 -iωτ diff tan( √ π 2 -iωτ diff /2) (28)
and is finite. The absorption peak is located at ωτ diff ≈ π 3 , and is seen on figure 3 (green curve) to coincide at large αL with ωτ D ≈ π (blue curve). The relevant timescale is then the diffusion time over the screening length.

Therefore for the case of no exchange, the dynamics of the charge density is entirely determined by its structure (L/λ), and corresponds to the diffusion of the ions over the relevant length scale, which can be either the interlayer distance or the screening length, depending on the strength of the electrostatic interactions.

The last noteworthy limit is that of fast chemical exchange (µ → 0). In equation ( 21), ξ(ω) simplifies to:

1 ξ(ω) = cos 2 αL 2 g( αL 2 ) + iωτ D sin γαL 2 sin αL 2 × K d cos 2 αL/2 1 + K d sin 2 αL/2 (29)
In this limit too the chemical relaxation time has disappeared, and the only relevant timescale corresponds to the transport of mobile ions. Further simplifications are obtained in the weak and strong electrostatics limits. In the former, r (ω) also vanishes as (αL) 2 :

r (ω) αL→0 ∼ (αL) 2 iωτ diff     1 - tan( √ -iωτ diff /2) √ -iωτ diff /2 × 1 1 + K d √ -iωτ diff 4 tan( √ -iωτ diff /2)     (30)
The position of the absorption peak now depends on K d , but for weak adsorption (K d 1), one recovers the same result as in the slow exchange limit (equation ( 27)). This is because in the weak adsorption case, what happens at the surfaces is of secondary importance with respect to the transport of the mobile ions; in particular, the rate of the exchange does not significantly affects the overall dynamics. In the fast exchange, strong electrostatics limit (αL → π), one finds the same result as in the slow exchange, strong electrostatics limit (equation ( 28)).

We now turn to the more complex case where the chemical reaction takes place at a finite rate. The relaxation dynamics of the charge density now involves an interplay between chemical reaction at the surfaces and transport in the interlayer. In order to study simultaneously the influence of the chemical rate (or µ = τ χ /τ diff ) and of the relative abundance of bound/free ions (K d ), we have computed numerically the absorption spectrum (Bode diagram for r ) for a fixed surface charge density and fixed interlayer distance (σ 0 = 0.131 Cm -2 , L = 6 Å) for univalent ions in a solvent of dielectric constant S = 78.

These values were already used to illustrate the equilibrium profiles in section 2 (see figure 2b). As the affinity of the ions for the surface (K d ) increases, the effective charge of the surfaces decreases, resulting in weaker electrostatic interactions (αL decreases).

The results are reported in figure 4 for K d = 0.1 (a), 1.0 (b) and 10 (c), where r / S is represented as a function of log(ωτ diff ) and log µ. For each exchange rate, the usual Bode diagram is found along a horizontal line. Fast exchange is found at the bottom of each plot, while slow exchange is located at the top. It is obvious from figure 4 that whatever the adsorption, two relaxations are found when the exchange is slow (though with varying relative weights), and only one is left when the chemical exchange rate is very high. This feature is ubiquitous in spectroscopy: although the physical relaxations are different from the present case, it has long been used in Nuclear Magnetic Resonance to evaluate kinetic constants [33,34,35].

For slow exchange (µ 1), the dominant mechanism relaxation is diffusion for weak adsorption (4a), and exchange for strong adsorption (4c); for intermediate values of the distribution constant, both mechanisms give a significant contribution to the permittivity (4b). Indeed, the maximum of r is found not to depend on µ in this limit for the weak adsorption case. The secondary peak, which does depend on µ and corresponds to the chemical exchange, contributes to less than half of the former. The opposite trend is observed for strong adsorption (4c), the primary peak occurring at a frequency proportional to 1/τ χ .

For fast exchange (µ 1), the chemical reaction is diffusion limited, and the relaxation is governed by transport only. The position of the maximum in this limit is seen to depend on K d , as was explained before, and to be located at higher frequencies for smaller K d 's. For K d = 0.1 (4a), the relaxation frequency does not depend on µ, and therefore it coincides with that in the slow chemistry limit. This is because for this choice of the numerical parameters, we are in the range of "weak electrostatics" even for small K d .

Hence, in the fast exchange limit, the permittivity is given by equation (30), which was said to reduce to equation (27) for weak adsorption (small K d ). These numerical results therefore illustrate the announced limiting expressions.

It appears from the preceding discussion that both the thermodynamics (K d ) and the dynamics (τ χ ) of the adsorption / desorption processes have a significant influence on the dielectric permittivity of the interlayer. Therefore, one could in principle extract from experimental data the values of these parameters. Indeed, structural informations (bare charge, interlayer spacing, valency of the ions) are generally known from other experimental techniques (e.g. titration, X-ray diffraction), and one could deduce from dielectric spectroscopy measurements the remaining unknown parameters (K d , τ χ , τ diff ) by fitting the data with equation ( 23). This promising method suffers however from several limitations, which we now wish to discuss. Those limitations are of two kinds: theoretical on the one hand, and related to the experiments on the other hand.

Let us first focus on limitations of the theory exposed in the present article. The present treatment of the interlayer solution, and particularly of the solvent as a continuum, is valid for not too small interlayer distances. It is expected to break down for distances comparable to the water molecular size, but has been shown to give reasonable results for the density profiles for distances as short as 15 Å [19]. At shorter distances, the molecular character of the solvent and of the ions must be included to account for the oscillatory density profiles in the interlayer [36]; a complete description of the dynamical response in this case in order to calculate the dielectric permittivity would be a great challenge. In this respect, the numerical values chosen here to illustrate our discussion might be at the limit of validity of the present theory. However, the results obtained with another system with a larger interlayer distance and a smaller charge density (leaving their product L/λ unchanged), would be the same, and fully in the range of applicability of the model. This could easily be obtained e.g. with a less charged clay at higher water content. Therefore the model might not be applicable to the experimental situation that guided us in the choice of the numerical values of the parameters, but that does not rule out the validity of the theory. In the same spirit, the choice of the value of the dielectric constant S is not straightforward, since for small (a few molecular) layers, its value (if its very existence is admitted in that case) might be different from its bulk one. Under confinement down to a few Å, it has been shown that the dielectric constant of water has generally a lower value than in the bulk [37]. Thus, the results are expected to compare favorably with experiments only for sufficiently large interlayer distances, for which these problems are ruled out.

In order to compare with experimental data, one should be able to take into account the contribution of the lamellae, which were assumed to be charged planes in the present description. This simplification allowed each interlayer spacing to be treated independently. In real systems however, the finite extent of the lamellae, and their possible dielectric properties might give a contribution to the dielectric response of the whole system. In particular, a mismatch between the dielectric constant of the solvent and that of the lamellae would necessitate to take all the interactions with image charges into account [38]. These limitations are probably not too serious. First, assuming that each cell is independent of the others is valid as soon as the surfaces are charged enough, which results in a screening of the electrical interactions across the layers. Second, their dielectric properties can be considered -at least for inorganic lamellae (such as clays)-as frequency independent, in this frequency region (Hz to GHz), because electronic polarization occurs at much higher frequencies. Thus the contribution of the layers reduces to a single dielectric constant L , which can then be introduced in the calculation of the resulting dielectric permittivity of the lamellar system.

Finally, two difficulties arise not from the model, but from its predictions. As was mentioned earlier, the ionic contribution to the complex permittivity is less than half of that of the solvent (see equation (24)). This could make it difficult to separate from that of the solvent, although the latter occurs at higher frequencies. The other point is that in some cases (and only in these case), the predicted permittivity is not very sensitive to the value of the parameters. For example, close to the charge saturation (αL → π), the value of αL depends only little on the exact value of L/λ. Hence, extracting the latter from the former cannot be achieved with a great precision. Similarly, when the chemical exchange is either very fast or very slow, the resulting permittivity does not depend on the precise value of µ = τ χ /τ diff . Therefore, one can hope to extract precise microscopic information only if the system is such that the electrostatics is not too strong (L/λ not too large) and the chemical exchange is both not too fast and not too slow (µ not too large and not too small). These remarks do not minimize the range of the present model, but only the hope that we can have to deduce microscopic information in some particular cases.

Conclusion

A fully analytical expression of the complex permittivity of lamellar systems in absence of added salt has been derived. The chemical exchange governing the partition between free and bound ions was properly taken into account, and the influence of both thermodynamic (distribution constant) and kinetic characteristics of this equilibrium was discussed. The main conclusions are twofolds.

First, when chemical exchange is infinitely slow (no exchange, just a static partition between free and bound ions), the only relaxation process is due to transport. The study of the complex permittivity as a function of frequency reveals that depending on the "strength of electrostatics" (the ratio between the interlayer distance L and the Gouy length λ), the relevant length scale vary between the whole interlayer distance (weak interactions case) and the screening length only (strong interactions case). Accordingly, the characteristic relaxation time ranges from the diffusion time (over the interlayer distance)

to the Debye relaxation time (diffusion over the screening length α -1 ). Furthermore the magnitude of the considered contribution to complex dielectric permittivity also depends greatly on L/λ. When electrostatic interactions are weak, the ionic profile in the inter-layer is flat and not very deformable. Hence, the contribution of normal ionic fluxes to the complex permittivity is small, and other phenomenon must be accounted for. In particular, for small interlayer spacings, the ionic motion along the charged surfaces, which is usually not considered, must be taken into account. This motion is the subject of a forthcoming paper [39].

Second, when the rate of chemical exchange is finite, two mechanisms can give a significant contribution to the complex permittivity: to the diffusion relaxation mode now adds the chemical relaxation mode. The relative weight of both mechanisms depends on the distribution constant K d between bound and free ions, and on the rate of exchange.

For slow exchange, relaxation is dominated by transport for weak adsorption and exchange for strong adsorption; for intermediate adsorption both are important. For fast exchange, the reaction at the surface becomes limited by the transport of reactants (free ions) in the interlayer; the characteristic time then depends on L/λ (diffusion over L or the screening length α -1 ).

To our knowledge, this is the first time that the complex permittivity of charged lamellar systems without added salt was determined (and this could be done analytically), taking properly into account the chemical exchange at the surfaces. The conclusions of the present article will be useful in the analysis of experimental data on such systems, for which the predictions of the hitherto available models (with added salt) do not hold.

We have particularly in mind the study of clay gels (at sufficiently high water content).

Further refinements of the model are envisaged. First, the calculations of the permittivity of the whole system, including the contribution of the dielectric layers, in the spirit of Sørensen [12], would allow for a direct comparison with experiments even if the contribution of the latter is not negligible. Second, the generalization to other geometries would open the way for the study of other no added salt cases: the convex spherical case is that of colloidal suspensions, while the concave spherical case would provide a more realistic model of inclusions in rocks. Finally, it would be of particular interest to try to improve the model in order to investigate small interlayer separation; it will be necessary to account for the molecular character of the solvent and of the ions, which increases dramatically the complexity of the problem, but simulations have shown that even in that case, the relaxation of the ionic cloud does not differ significantly from the Debye relaxation time calculated with the bulk parameters [40], at least in a certain range of the latter.

with h(v) = f (u). Two independent solutions for -1 ≤ v ≤ 1 are: where P and Q are associated Legendre functions (see e.g. reference [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]). From these two solutions an even solution can be found by forming the appropriate linear combination, namely:

f 1 (u) = i √ 2π P
g(u) = cos(πγ/2) f 2 (u) -sin(πγ/2) f 1 (u) (36)
After some basic trigonometric manipulations this simplifies to the expression of g given by ( 21), which is obviously even.

Therefore, the full solution for E 1 (x) is of the form given by ( 21) where the integration constant ξ(ω) remains to be determined from the boundary conditions. Equation ( 20) can be rewritten in terms of E 1 using E 1 (L/2) = -δΓ + z/ 0 s and E 1 (L/2) = c 1 (L/2)z/ 0 s (see equation ( 4)). Using the analytic form ( 21) for E 1 (x), it is then straightforward to derive that ξ(ω) is given by [START_REF] Andra | Tome 1: Contexte et objet[END_REF].

  Figure 2a shows the evolution of αL with L/λ. For small charge densities (L/λ → 0) one has αL ∼ (L/λ) 1/2 , while αL → π for large values of L/λ. The concentration profiles corresponding to three values of the charge density have been reported in figure 2b. The values are chosen to reproduce typical charge density and interlayer spacing for a smectite clay at low water content (σ 0 = 0.131 Cm -2 , L = 6 Å, S = 78), and three values of the distribution constant K d (0.1, 1.0 and 10), at room temperature (k B T = 4.10 -21 J). The corresponding values of L/λ and αL are summarized in table

  1) cos[(γ + 1)(π/2 -u)] + (γ + 1) cos[(γ -1)(π/2 -u)] 1) sin[(γ + 1)(π/2 -u)] -(γ + 1) sin[(γ -1)(π/2 -u)] cos 2 u (35)

Table 1 :

 1 Ratio between the interlayer distance L and the Gouy length λ (resp. the screening length α -1 ) corresponding to the numerical values used for the discussion (see text and figure4).

Figure 1 :

 1 Figure 1: Schematic view of the interlayer. Ions can be either diffusing in the interlayer, or bound to the charged surfaces. Both forms are in chemical equilibrium, exchanging with rates k 1 and k 2 (see text). Adsorption of ions results in an effective charge of the surfaces.

Figure 2 :

 2 Figure 2: (a) Reduced screening parameter αL vs. the ratio L/λ between the interlayer distance and the Gouy length, which characterizes the strength of electrostatics. L/λ and hence αL decrease with increasing adsorption (K d ). At large fields, the screening length saturates (limit value αL = π). (b) Concentration profiles corresponding to αL = 0.63 ( ), 1.38 (•) and 1.74 (+).

Figure 3 :

 3 Figure3: Imaginary part of the reduced complex permittivity ( / S ) as a function of frequency (normalized by the diffusion time τ diff = L 2 /D), and of the ratio L/λ between the interlayer distance and the Gouy length, which controls the strength of electrostatics in the interlayer, in the case of infinitely slow exchange (no exchange). For L/λ 1, the concentration profile is very flat, and not very sensitive to an applied electric field, while for L/λ 1, electric interactions result in an ionic distribution sharply peaked at the charged surfaces, and deformable by an electric field. The characteristic time for the relaxation process corresponds in the former case to the diffusion over the interlayer distance (τ τ diff /π 2 -red curve) while in the latter case the relevant distance is the screening length α -1 (τ τ D /π -blue curve) which reaches its saturation value L/π when L/λ 1 (τ τ diff /π 3 -green curve).

Figure 4 :

 4 Figure 4: Imaginary part of the reduced complex permittivity ( / S ) for K d = 0.1 (a), 1.0 (b) and 10 (c), as a function of frequency (normalized by the diffusion time τ diff = L 2 /D), and of µ = τ χ /τ diff . For slow exchange (µ 1), the dominant mechanism is diffusion for weak adsorption (a), and exchange for strong adsorption (c); for intermediate values of the distribution constant, both mechanisms give a significant contribution (b). For fast exchange (µ 1), the relaxation is due to transport only, and the position of the maximum depends on K d .

Figure 1 :

 1 Figure 1: B. Rotenberg et. al

  2 , is analogous to the Debye relaxation time in bulk electrolytes, and depends on the distribution constant K d (because α does). It is conveniently expressed as a function of the diffusion time τ diff = L 2 /D which an ion takes to diffuse over the interlayer spacing L. One simply has:

Figure 2: B. Rotenberg et. al

Acknowledgments B.R. acknowledges financial support of the Ecole Normale Supérieure de Paris.

Appendix A

In this appendix we derive the solution [START_REF] Andra | Tome 1: Contexte et objet[END_REF] of the differential equation (17). We need to find an even solution satisfying the appropriate boundary condition (20). We first rewrite (17) using the dimensionless variables u = αx and γ = 1 -iω/Dα 2 as:

with f (u) = E 1 (x). It is easily checked that a particular solution of ( 31) is:

This solution is an even function. The general solutions can be found by introducing the change of variable v = sin u (0 ≤ |v| ≤ 1) and dropping the right hand side of equation (31), leading to:


Benjamin Rotenberg
email: rotenber@ccr.jussieu.fr

Virginie Marry

Jean-François Dufrêche

Eric Giffaut

Pierre Turq

A multiscale

Keywords: numbers clay, ion transport, diffusion, two-state, multiscale, coarse-graining

INTRODUCTION

Electrokinetic phenomena are related to the transport properties of mobile charged species (ions) near charged surfaces (e.g. colloids or macroscopic mineral surfaces). The presence of solid surfaces influences dramatically the ionic mobility, especially in porous media such as compacted clays, where geometrical confinement is reinforced by the interactions between the ionic species and the charged surfaces. This retention property has made clayey materials good candidates for the long term storage of chemical and nuclear waste, which often contain charged species.

The water/mineral interface plays a central role in the overall dynamics of ions, which is usually described using the sorption concept: An ion can be trapped temporarily, immobilized for some time, thereby reducing its mobility. The extensive literature on the related sorption data (empirical partitioning coefficients K d ) obtained by geologists and engineers provide a convenient tool for the modeling and prediction of ion mobility in these materials [START_REF] Andra | Matériaux à base d'argiles gonflantes[END_REF][START_REF] Daqing | [END_REF][3][4]. On the microscopic scale however, there is no clear-cut distinction between trapped and mobile ions. In compacted clays, the water content is very low, and there has been some evidence that the usual macroscopic equations (e.g. Navier-Stokes for solvent flow, Poisson-Boltzmann for ionic distributions) are at the limit of their validity [5]. A detailed description of ionic mobility is still at hand, however, using Molecular Dynamics (MD) simulations [6][7][8][9][10]. They provide an atomically resolved description of ion (and water) motion near the mineral surfaces, and allow for the computation of the short-time diffusion coefficient of ions.

The purpose of the present paper is to provide a link between the molecular simulations and the simplified diffusion-reaction picture, by an appropriate coarse-graining strategy. In section 2, we highlight some MD results on ionic diffusion in clay interlayers. These results are then used in section 4 to build a robust mesoscopic (Fokker-Planck) description, in which the interactions of the ions with the surfaces are introduced as an effective external potential, and the interactions of the ions with the solvent are described by a unique friction parameter.

In turn, the mesoscopic model is used in section 5 to build a diffusion-reaction scheme, by first determining the mobility of the ions in the interlayer as a function of the strength of their interaction with the surface. These results are eventually cast into a diffusion-reaction scheme, introducing in particular the fraction of mobile ions, or equivalently the distribution coefficient K d .

MOLECULAR MODELING OF COMPACTED CLAYS

The most detailed information on water and ion dynamics in clay interlayers can be obtained by molecular modeling, which provides an atomically resolved description of the system. Molecular simulations are done for a montmorillonite clay model whose unit cell formula is M 0.75 Si 8 (Al 3.25 Mg 0.75 ) O 20 (OH) 4 , where M is the counterion. Clay particles are formed by stacks of aluminosilicate layers. Each layer of montmorillonite consists of three sheets : An octahedral sheet of Al oxide between two tetrahedral sheets of Si oxide (see figure 1). Substitutions of Al III by Mg II give rise to a negative charge on the mineral compensated by cations (M). The above clay is an idealized montmorillonite which is close to the natural clay MX80 studied for its potential applications for the geological storage of nuclear waste [START_REF] Andra | Évaluation de la faisabilité du stockage géologique en formation argileuse[END_REF][START_REF] Bradbury | Near Field Sorption Data Bases for Compacted MX-80 Bentonite for Performance Assesment of a High-Level Radioactive Waste Repository in Opalinus Clay Host Rock[END_REF]. The simulations are done with Li + , Na + , K + and Cs + counterions.

Na + is one of the major natural cations in MX80 (with Ca 2+ ) and Cs + is a potential radionuclide. Li + and K + are also studied to complete the comparison of diffusion properties of the alkali cations.

The parameters used in the interaction potentials are taken from literature and detailed in previous articles [8,9]. For the counterions, we use the Lennard-Jones parameters of Koneshan [START_REF] Koneshan | [END_REF]. The simulation box contains two clay layers of eight unit cells each. The positions of the clay atoms are taken from X-ray data [START_REF] Brindley | Crystal Structures of Clay Minerals and their X-ray Identification[END_REF]. The lateral dimensions of the layer are then 20.72 × 17.94 Å2 , its thickness 6.54 Å. At low relative humidities, water can enter in the interlayer in order to solvate the counterions and form one, two or three layers of water in the interlayers, depending on the compaction, humidity and the type of counterion. Only the monohydrated state is simulated, here with N = 36 water molecules in each interlayer space, i.e. 6 water molecules per counterion. In this case, the space is too constrained for the cation to be fully hydrated by water. Hence cations are always in direct contact with the surface clay atoms. A snapshot of the simulation box is given in figure 1.

Prior to dynamic studies, the average interlayer spacings for each type of cation is evaluated by Monte Carlo simulations in the (N , σ zz =1 bar, T =298 K) ensemble. We found 12.25 Å, 12.3 Å, 12.4 Å and 12.7 Å for Li + , Na + , K + and Cs + respectively. These results agree with RX diffraction experiments [START_REF] Bérend | [END_REF][16][17][18] and other simulations on montmorillonites [19][20][21][22]. After equilibration, the interlayer structural properties are calculated by fixing interlayer spacings at their average value. Moreover, the hexagonal cavities formed by silicon and oxygen atoms of opposite surfaces are maintained face to face during the simulations.

This allows an easier interpretation of the preferential paths of the cations along the surface of the clay layer. All MD simulations were performed using the DLPOLY software. A timestep of 1 fs was used, and a set of 720 ps simulations was carried out for each cation (2 to 4 runs per cation). Simulations were done in the N V T ensemble, with T =298 K, using a Nosé-Hoover thermostat.

Water and cation distributions in the direction perpendicular to the clay surface are shown on figure 2 for one interlayer space. The peak of water oxygen atoms widens with increasing size of the cation and even starts to split in the case of K + and Cs + . The opposite is observed for the cation distributions. An obvious explanation is that the distance between surfaces is controled by the size of the largest species. In the case of small cations, the size of the water molecule is determining and the interlayer spacings remain around 12.3 Å : The cations are free to go from one surface to another according to their interactions with the surface atoms and their hydration state. In the case of large cations, the size of the cation becomes more important and the water moves perpendicular to the clay surface more easily.

It is already known that alkali cations have different behaviours towards water and clay surfaces. Experimentally, Na-and Li-montmorillonites can swell until the tri-hydrated state (three water layers in the interlayer space) under increasing humidity, although Cs-and K-montmorillonites can only reach the monohydrated state [START_REF] Bérend | [END_REF]. Simulations show that in the bihydrated case, Na + remains in the middle of the interlayer space, surrounded by six water molecules, like in its first hydration shell in bulk water [9]. On the contrary, Cs + remains close to the clay surface, even in a fictitious bihydrated state. The maps of their preferential positions along the clay surface in the monohydrated state are given on figure 3 with those of Li + and K + . They illustrate the different behaviour of the alkali cations. Large cations (Cs + and K + ) exhibit a site-to-site diffusion above surface silicon atoms. They complete their solvation shell with the three oxygen atoms linked to the silicon. Small cations (Na + and Li + ) prefer oxygen surface atoms to silicon and hexagonal cavity centers. A previous study showed that when the clay layers are allowed to move horizontally (in the two directions parallel to the surfaces) the hexagonal cavities silicon atoms are preferentially facing each other in the case of Cs-montmorillonite. It is not the case with Na-montmorillonite. However, even with another relative shift of the hexagonal cavities, a preference of Na + for surface oxygen atoms was already observed [8].

COARSE-GRAINING STRATEGY

Molecular simulations provide a picture of water and ion dynamics in clays on short timescales (ns) and small lengthscales (nm). On the macroscopic scale of observation, and for engineering purposes, these molecular details are not necessary, and a coarser representation at the macroscopic scale (>mm) is generally adopted. In order to link this microscopic information to macroscopic observables, and especially to the apparent diffusion coefficient of each species, one usually separates the geometrical confinement effect from the physical interactions with the mineral surfaces. For the geometrical part, a set of empirical parameters is introduced : Porosity , constrictivity δ and tortuosity τ . Porosity is the fraction of the volume accessible to the diffusive species. Constrictivity is a geometrical factor accounting for the closing/opening of diffusion pathways: δ > 1 if the pore width increases along the diffusion path, δ < 1 in the opposite case [START_REF] Kato | Scientific basis for nuclear waste man-12 agement[END_REF][START_REF] Ochs | [END_REF][25]. Tortuosity quantifies the length of the diffusion pathways in the medium with respect to a straight trajectory. Only and the ratio δ/τ 2 are available experimentally. In addition to geometrical effects, the interactions of the diffusing particles with the mineral surfaces are accounted for in the framework of a diffusion-reaction scheme. An ion can remain temporarily trapped at the mineral surface, thereby reducing its apparent diffusion coefficient, according to the "chemical" reaction :

mobile adsorbed ( 1
)
This mechanism is quantified by a partitioning coefficient K d , which is the ratio between the concentrations of bound and mobile forms of each chemical species. At the macroscopic scale, the apparent diffusion coefficient is deduced from the "free" diffusion coefficient D 0 using :

D a = D 0 × δ τ 2 × ( + ρK d ) (2)
with ρ the dry density of the medium.

As mentioned earlier, there is no clear-cut distinction on the molecular scale between mobile and bound ions. Therefore, the definition of a diffusion-reaction model from the microscopic dynamics (as obtained by MD simulations) is not without its ambiguities. We now suggest a tentative method to connect the atomic level of modeling developed in section 2 to the more approximate one just described. The general coarse-graining strategy involves two steps. The first one defines a mesoscopic description, calibrated on the molecular dynamics results, i.e. taking advantage of the microscopic insights gained from molecular modeling.

The solvent is described as a continuum, and the ion/surface interaction is included as an effective potential. The second step derives, from the mesoscopic model, the parameters defining the diffusion-reaction scheme, namely a dimensionless analogue of K d , and the kinetic rates of exchange. Section 4 describes the first step, while the second is postponed to section 5.

MESOSCOPIC DESCRIPTION

This section describes how to calibrate the continuous solvent description of the interlayer on MD simulations. Such an approach uses the actual results on the complex system (described in section 2) instead of an infinite dilution reference.

When the solute (ion) moves slower than the solvent (water) one can average over the solvent degrees of freedom [START_REF] Résibois | Electrolyte Theory[END_REF], and describe the solute/solvent interaction by a single friction parameter γ (in s -1 ) : The friction force on an ion is proportional to its velocity (F solv = -mγv), and the resulting diffusion coefficient is given by Einstein's relation D 0 = v 2 T /γ, with v T = k B T /m the thermal velocity of the ion. The value of the friction can be extracted from MD simulations using a fit of the velocity auto-correlation function (VACF) [START_REF] Hansen | Theory of Simple Liquids[END_REF]:

Z(t) = 1 3 v(t) • v(0) = k B T m exp(-γt) (3)
The right-hand side is the analytical form for Brownian motion, and the VACF is obtained by analysis of the MD trajectories. A typical value for Cs + is γ ∼ 5 ps -1 . However, the final results of our approach on the definition of K d , do not depend on the precise value of γ. This point will be discussed in section 5.

The interaction of the ion with the mineral surface is introduced as an external (free energy) potential V , which can be computed by an inversion procedure from the ionic density profiles obtained by molecular dynamics simulations. The solute dynamics is then characterized by the time evolution of the probability density function f (x, v, t) of finding an ion at a given position with a given velocity. The moments of f in velocity space are the ionic density and flux :

ρ(x, t) = f (x, v, t) dv (4) j(x, t) = f (x, v, t)v dv (5)
This probability density evolves, within this simplified mesoscopic description, according to the Fokker-Planck (FP) equation [START_REF] Risken | The Fokker-Planck Equation[END_REF] :

∂ t f + v • ∇f = ∇ v • γv + ∇V m + γk B T m ∇ v f (6)
The potential V encompasses all interactions involving the ion (except from the ion/water interaction, which is the friction): In addition to the ion/surface interaction, it can also include an external field (e.g. an electrical field), which is necessary to determine the mobility of the ions. To solve the FP equation numerically, we used a lattice method inspired by the Lattice-Boltzmann methodology: The lattice Fokker-Planck method [START_REF] Melchionna | [END_REF][30][31][32].

In order to investigate efficiently the effect of surface-ion interactions on ionic mobility (and later on the parameters of the two-state model), we did not use a potential directly derived from MD simulations, and used instead an analytical ansatz, which respects the characteristic features of the true potential (topology, barrier height, ...). The magnitude of the ion/surface interactions is quantified by the ratio V s /k B T = βV s , where V s is the barrier connecting two adjacent free energy minima. The precise form of the potential is given by:

V (x, y) = V s × 3 + 2 cos 2π a x + y √ 3 + 2 cos 2π a x - y √ 3 + 2 cos 2π a 2y √ 3 (7)
with a the distance between two maxima (a ∼ 5.2 Å for montmorillonite). is clearly room for improvement in the Na + case (a), although it already reflects its more diffuse behaviour. One should keep in mind that in the Li + and Na + cases the hydration tendency is rather different from the K + and Cs + cases, and that the results of section 2 were obtained for the monohydrated state only, with the constraint of face-to-face clay layers. In this situation, Li + (and to a smaller extent Na + ) tend to avoid the center of the interlayer (see the two peaks on figure 2), so that a strictly two dimensional modeling as the one undertaken here is not fully appropriate. The study of the bihydrated state, in which Li + and Na + are clearly located in the midplane of the interlayer [20,33] is under progress. This will open the way to the study of more hydrated states.

IONIC MOBILITY AND THE DIFFUSION-REACTION SCHEME

The second coarse-graining step consists in the definition of the two-state model from the mesoscopic (FP) description. The underlying idea is that in such a model, the reduction in mobility is traced back to the fixation of some ions on a surface: only a fraction f m is considered as mobile. The link between microscopic mobility and two-state model is thus provided by the equality f m = µ app /µ 0 , where µ 0 = 1/γ is the "free" mobility of the ions, that is when no interactions with the surfaces are present.

Thus we first need to determine the mobility of the ions interacting with the mineral surfaces. This is achieved by applying an external force F ext on the ions, in addition to the one arising from the clay layer itself. In practice, equation 6 is now solved with (∇V -F ext )/m on the right-hand side. At steady-state, the average velocity is proportional to the applied field, the proportionality coefficient being related to the apparent mobility :

v st = µ app /m × F ext .
The latter is of course lower than the "free" mobility µ 0 . The mobility µ app is reported as a function of βV s in figure 5. The value is rescaled by µ 0 . The two components of the mobility tensor, obtained by applying the force along each direction, are equal. Isotropy of the mobility is a known consequence of the symmetry of the honeycomb lattice [34,35]. As expected, the mobility decreases with increasing affinity for the surface (increasing βV s ). Two regimes are observed, with a crossover at approximately βV s ∼ 1.

Their significance will be discussed below, when interpreted in terms of K d .

Within the two-state diffusion-reaction scheme, the ratio µ app /µ 0 is precisely the fraction of mobile ions f m , that is the ions that are not bound to the surfaces. Translating this fraction of mobile ions in terms of the partioning coefficient K d , one obtains :

K d = 1 -f m f m = µ 0 -µ app µ app (8)
Note that the definition of K d is slightly different from the one given above, because no geometrical effects are included in the present description. In particular, this definition makes K d dimensionless. Our purpose is to give a microscopic basis to a simple two-state model, which is generally the starting point of even coarser models, such as the one presented in section 3. The simplicity of the latter allows to include other features of the real material (distribution of porosity sizes, presence of minerals other than clay, ...) which are out of the scope of the present study. Equation 8 clearly shows that K d only depends on the ratio µ app /µ 0 . In our simple model, which neglects ionic correlations, and in the high friction limit (the one relevant here), this ratio is in fact independent of the value of the friction γ, as was mentioned earlier.

Results for K d as a function of the interaction parameter βV s are represented on a logarithmic scale in figure 6. The two regimes appear now more obviously : For strong interactions with the walls or at low temperature, the diffusion follows an Arrhenius law, with an activation energy of V s , the barrier that an ion needs to overcome to move from one energy minimum to the next; for weak interactions or at high temperature, no such activated diffusion is seen. These results are consistent with the trajectories observed in MD simulations for sodium (weak interaction case) and cesium (strong interaction case).

From the knowledge of the amount of bound ions, it is also possible to determine: (a)

the localization of the bound and mobile states, and (b) the adsorption/desorption rates, thereby completely specifying the diffusion-reaction model [35]. The results, not shown, indicate that the border between fixed and mobile ions is independent of the precise value of βV s , provided that the latter is large enough (that is, for strong affinity for the surface).

This justifies the notion of trapping site, because it makes reference to the mineral surface only, and not to the ions moving in its vicinity. We refer the interested reader to a previous publication [35].

CONCLUSION

We have demonstrated how a two-state diffusion-reaction model can be developed from the microscopic dynamics of ions in compacted clays, following a two step coarse-graining strategy. Molecular dynamics simulations first provide the information which is necessary to calibrate a mesoscopic description of the clay/ion/water system. The latter allows for a determination of ionic mobility near the clay surface, which in turn can be translated in terms of the diffusion-reaction scheme via the fraction of mobile ions f m , or equivalently the distribution coefficient K d .

This coarse-graining strategy extends the continuous solvent description used in the theory of bulk electrolyte solution since Debye, Hückel and Onsager [START_REF] Barthel | Physical Chemistry of Electrolyte Solutions[END_REF] to non infinitely diluted reference states. It is of course not restricted to clays, and applies to any charged surfaces.

A generalization of the method is needed if one wants to consider semi-infinite systems (that is a surface facing a bulk fluid), as opposed to the confined, quasi two-dimensional system of a clay interlayer. However we expect no particular difficulty, except for the computational cost of solving the FP equation in three dimesions. The method used in this paper, lattice Fokker-Planck, is a method of choice to overcome this technical difficulty.

  The density profiles corresponding to this model potential are the steady state solutions of equation 6 with V given by equation 7. The solution statisfies the Boltzmann law: ρ(x, y) ∝ exp(-βV (x, y)). Results are shown in figure4. These profiles are to be compared to those obtained with the fully atomic description of figure3: The K + and Cs + behaviours are relatively well reproduced for the intermediate (b) and high (c) values of βV s . There

FIG. 1 :FIG. 2 :FIG. 3 :FIG. 4 :FIG. 5 : 14 FIG. 6 :

 12345146 FIG.1: Snapshot of the simulation box for Cs-montmorillonite for the case of a water monolayer.Oxygen atoms are in red, H in white, Si in yellow, Al in green and Cs in blue. The interlayer distance corresponds to a water content of 6 water molecules per cation, a temperature of 298 K and a pressure of 1 bar (normal to clay layers).

  Figure 2

Acknowledgments BR acknowlegdes financial support from ANDRA (Agence Nationale pour la Gestion des Déchets Radioactifs). Financial support of GdR PARIS (Physico-chimie des Actinides et autres Radioéléments aux Interfaces et en Solution) is gratefully acknowledged. LI2C is part of the SFR Sciences Chimiques pour l'Analyse et la Mesure de Paris-Centre.