Diffusion calculations on reconstructed bentonite microstructures with anion exclusion effects

Abstract

investigated at the nanopore scale are upscaled using the Homogenization of Porous Media approach.Homogenization computations rely on a hierarchical description of bentonite that acknowledges the existence of pores networks at different scales.At the scale of montmorillonite layers, digitized TEM images have been employed to simulate diffusion of ionic solutes by considering electrostatic interactions in the vicinity of the negatively-charged clay platelets' surface.Finite Element microstructures are created after extraction of the contours of the layers using dedicated image processing algorithms.Local electric potential distribution, anion exclusion and cation inclusion are displayed by ion distribution maps.The effective diffusion tensor and the transport equation obtained through volume averaging are then used to simulate diffusion at the scale of a Wyoming bentonite sample composed of clay gels of variable density, mineral grains and micropores.Qualitative comparisons were made with existing diffusion data, and a particular attention is given to the anisotropy of the diffusion tensors at both the

Introduction
Due to their prevalence in the lithosphere and their high capability of sorbing pollutants, smectite clays play a foreground role in environmental pollution studies, waste management and soil science (Choi and Oscarson, 1996;Gonzales Sanchez et al., 2009;Suuronen et al., 2014;Whittaker et al., 2020).The use of compacted bentonite as a buffer material in radioactive or domestic waste storage is indeed motivated by its exceptional swelling and retention properties.At the nanoscopic scale, their microstructure is characterized by a regular stacking of montmorillonite layers of very high aspect ratio, with a permanent negative surface charge compensated by interlayer cations.The resulting arrangement will condition the microstructural properties (cationic exchange capacity, specific surface, swelling potential, retention of water and pollutants), and the resulting hydro-mechanical behavior of hydrated smectites at engineering scale (Hetzel et al., 1994;Pusch, 2001;Smith et al., 2004;Whittaker et al., 2020).
Diffusion measurements at the nanoscopic scale (Gonzales Sanchez et al., 2009) have shown the preponderance of surface effects, such as reduced water viscosity and electrostatic interactions between the diffusing species (water molecules, ions) and the clay surfaces and compensating cations, those effects being expressed by an electrostatic constraint factor.At the micrometer scale, the morphology and the connectivity of the pore network play a foreground role in the resulting diffusion through micropores, this contribution being expressed by an overall geometric factor or tortuosity.However, as pointed out in (Gonzales Sanchez et al., 2009), those factors cannot be related to the porosity or other easily derived properties of the porous medium, and must be regarded as empirical parameters.Monte Carlo and Molecular Dynamics Simulations have been employed in order to identify the interlayer molecular structure and local diffusion coefficients of tracers in the vicinity of charged montmorillonite mineral surfaces (Chang et al., 1995;Marry and Turq, 2003).The simulated values are then introduced into a macroscopic model of diffusion using the known porosity value and empirical parameters such as the constrictivity factor and the tortuosity.Although the resulting macroscopic diffusion coefficients are close to experimental data, Monte Carlo and Molecular Dynamics Simulations cannot take into account the geometry of compacted clays beyond the interlayer region (Marry and Turq, 2003).
Besides, as pointed out in (Pusch, 2001), most of the microstructural models of clays refer only to the detailed particle-to-particle interaction by assuming a parallel arrangement for the montmorillonite layers, and disregard the real pore distribution and the scale-dependent variation in density of the clay gels.
In complementarity with recent works, we try to model ionic diffusion using a realistic description of smectite microstructure.We aim to investigate the specific contributions of the surface charge and the interlayer porosity at the nanoscopic scale, as well as the part played by micropores and clay gels of varying density at the micrometric scale.Three levels of description are considered during the model derivation through homogenization.The microscopic level designates the scale of stacked montmorillonite lamellae and interlayer water whose physical properties are affected by the pore walls (Marry and Turq, 2003).The Homogenization of Periodic Media approach is employed to obtain the effective properties governing solute transport by diffusion at the scale of an assembly of montmorillonite layers and interlayer water.The mesoscopic level refers to the scale of clay gels, mineral grains and inter-aggregate pores (Keller et al., 2014;Tomioka et al., 2010).The effective properties computed by upscaling from the microscopic to the mesoscopic scales are affected to the clay gels, while diffusion within the inter-aggregate pores is assumed to take place without surface effects.The macroscopic level corresponds to samples subjected to diffusion tests in the laboratory.Upscaling is then performed from the mesoscopic scale to the scale of a laboratory sample of bentonite to compute the effective diffusion tensor and the effective coupling tensor of bentonite.

Construction of the numerical microstructures 2.1 At the microscopic scale
The porous medium investigated at the nanoscale is composed of two phases: the solid phase of montmorillonite layers Ω s and the interstitial space saturated with the electrolyte solution Ω f .The microstructure illustrated in Figure 1 has been obtained by exploiting a high-resolution Transmission Electron Micrograph from (Hetzel et al., 1994), and represents one quarter of the complete Representative Volume Element (see Section 3.4).Using two-dimensional Diffusion calculations on reconstructed bentonite microstructures microstructures for describing solute transport around montmorillonite layers can be justified by typical lateral extensions observed for montmorillonite surfaces, which are more than a hundred times higher than their average thickness (Malikova et al., 2008;Pusch, 2001), as well as by experimental evidence (Gonzales Sanchez et al., 2009;Malikova et al., 2008;Marry and Turq, 2003).
The image analysis procedure developed using MATLAB is the result of Fig. 1 Finite Element Microstructure obtained after processing of a Transmission Electron Micrograph from (Hetzel et al, 1994).Image width:0.4 µm.Besides, it has been observed that a single threshold value leads generally to better final results than an adaptive thresholding.Contrast enhancement is then performed by bottom-hat and top-hat filtering.Next, a series of morphological opening and closing operations are performed with elongated elements of increasing size, oriented along the average orientation computed for the assembly of layers.Objects displaying an eccentricity smaller than the average eccentricity are then removed, as well as objects which area represents less than one percent of the maximal object area.Finally, the contours of the remaining objects are extracted using Moore-Neighboor's algorithm with Jacob's stopping criterion (Gonzalez et al., 2004).A dedicated procedure has been implemented to import the coordinates of the contour points and create the microstructure in the Finite Element-based code Castem (Le Fichoux, 2022).

At the mesoscopic scale
Numerical microstructures have been obtained at the mesoscopic scale by processing micrographs taken by Transmission Electron Microscopy (TEM) on resin-impregnated samples of Wyoming bentonite (Pusch, 2001).Four distinct phases are considered in the original TEM micrographs: non-smectite grains, dense clay gels, soft clay gels, and micropores.
The thresholding procedure is based on the algorithm of maximization of the gray level histogram's entropy developed by (Kapur et al., 1985).The threshold values obtained are validated by comparison of the resulting phase proportions with the cumulative distribution curve and the initial known phase proportions (Bouchelaghem and Pusch, 2017;Tomioka et al., 2010).The contours are extracted separately for each phase using Moore-Neighboor's algorithm with Jacob's stopping criterion, before being imported and used to create the Finite Element mesh in Castem.Such a microstructure is illustrated in Figure 6.The resolution is about 25 nanometers per pixel, which allows a good representation of the internal microstructure.
3 Homogenization at the microscopic scale

Modeling assumptions and local description of the problem
The modeling work starts with equations that have classically been employed to generalize Fick's model of diffusion, in order to account for chemical and electrostatic couplings when ionic species are transported within a charged porous medium (Cheng and Hendry, 2014;Moyne and Murad, 2006;Smith et al., 2004).The molar flux density is expressed by Nernst-Planck equation in the stationary mass balance equation written for each chemical species i of concentration c i and valence z i : In Eq. (1), D 0,i is the self-diffusion coefficient of solute i in the pore fluid, F is Faraday's constant, R the universal gas constant and T the constant absolute temperature.
∂ ∂xi e i is the nabla operator of differentiation with respect to the Cartesian spatial coordinates.
Poisson's equation with a source term representing the net charge density (Cheng and Hendry, 2014;Moyne and Murad, 2006) for N solutes is written to describe the distribution of the electrical potential ψ under steady-state conditions: where the permittivity ε is the product of the permittivity of vacuum ϵ 0 and the relative permittivity of water ϵ w : ε = ϵ 0 ϵ w .
Eqns. ( 1) and ( 2) are completed with boundary conditions in the vicinity of the montmorillonite platelets' surface Γ: (3) Diffusion calculations on reconstructed bentonite microstructures n, the normal unit vector on Γ, points towards Ω s .
The voltage gradient flux is imposed at the surface of the montmorillonite layers in order to express the relationship between ψ and the known surface charge density σ (Cheng and Hendry, 2014;Liu et al., 2019): In what follows, to simplify the model presentation, we consider a binary solu- Following (Cheng and Hendry, 2014;Moyne and Murad, 2006;Sasidhar and Ruckenstein, 2003), ψ is decomposed into two components: where ψ ex is associated with the exterior solution and the induced membrane potential, while ψ in is related with the interior diffuse double-layer (DDL) and the negative charge on the montmorillonite layers.ψ in is effective within the Debye length, so that within the DDL, the concentration of cations is dominant relative to the concentration of anions, but both tend to approach their respective equilibrium concentrations in the exterior solution.This implies a Boltzmann distribution of concentration for cations (c + ) and anions (c − ): with c ± ex the electrolyte concentration associated with the exterior solution.
As a result, Eqn.
(2) can be decomposed in two equations:

Scaling analysis
To normalize the equations, proper estimates of the physical variables are required.The subscript c is used to denote a characteristic value, such that each quantity f (f = c + , c − , ψ, ...) is written as: f = f c f ′ , with f ′ the corresponding dimensionless variable.
Following an approach detailed in (Sanchez-Palencia, 1980), two characteristic lengths scales are introduced: -a microscopic length l ′ associated with heterogeneities present within the microstructure.As in (Cheng and Hendry, 2014;Moyne and Murad, 2006), we take l ′ as the Debye length τ D : . τ D represents the scale of variation of ψ or the thickness of the double-layer.τ D , which depends on the electrolyte concentration c c = (c ± ex ) c (typical values being comprised between 1 and 100 mol • m −3 (Liu et al., 2019;Whittaker et al., 2020)), varies generally between one and ten nanometers.
-a mesoscopic length l associated with typical dimensions of the microstructure (i.e. the characteristic dimension of a montmorillonite particle, several hundreds of nanometers).
The perturbation parameter ε ′ is then defined by the following ratio: , such that ε ′ ≪ 1 guarantees the separation of scales required for homogenization (Sanchez-Palencia, 1980).
(2) the characteristic value for the electrolyte concentration is taken to be and Hendry, 2014;Liu et al., 2019;Moyne and Murad, 2006).
Introducing the dimensionless space variable x ′ i = xi l (i = 1, 2 in 2D), using the previous estimates of c c , ψ c and σ c and accounting for Eqn.(7), the system of Eqs.(1),(3),( 4) and ( 8) is rewritten in dimensionless form as follows: ) By assuming a local periodicity for the microstructure and all physical variables, c + ex , c − ex , ψ in and ψ ex can be written as asymptotic developments with respect to the mesoscopic coordinates x ′ = x l and the local coordinates ε (to simplify the notation, in the following we drop the ′ symbol in the expression of dimensionless quantities): ) where each term f (i) (x, y) in the developments above is Ω ′ -periodical.
As a result, differentiation with respect to x, ∇ x , is replaced by differentiation with respect to x and y, ∇ x + 1 ε ∇ y (Bouchelaghem, 2018;Sanchez-Palencia, 1980).Substituting asymptotic developments ( 17)-(20) into Eqns.( 9)-( 16) and accounting for the new rules of differentiation, local boundary-value problems are then derived at successive orders of approximation with respect to ε ′ .The local problems allow to identify the mesoscopic quantities (that depend only on x) and the local distributions of concentration and electric potential, as well as to obtain the effective properties by averaging over the RVE.

Local problems
At the lowest level, it is straightforward to show that c are mesoscopic quantities, which do not depend on the local variable y (see Appendix A).This allows to identify , i = 1, 2, as the mesoscopic force components of the imposed electrolyte concentration gradient, and as the macroscopic force of the induced membrane potential gradient (Cheng and Hendry, 2014).
If we consider Eq. ( 12) at order O(ε ′0 ) = O(1) and Eq. ( 16) at order O(1), in verifies the following dimensionless problem within the unit cell Ω ′ : At the next order of approximation, we obtain the problems verified by ψ The problems verified by β l (y), l = 1, 2, are detailed in Appendix A. After some developments (see Appendix A), we find that c The local problems verified by c + l , c − l , d + l and d − l (l = 1, 2) are detailed in Appendix A.
By writing Eqs. ( 9) and ( 10 In Eqs.( 26), the components of the mesoscopic diffusion tensors D ± and the tensors D ± ψ expressing the effect on ion diffusion of the membrane potential are defined as follows (i, l = 1, 2): For the sake of simplicity, as in (Cheng and Hendry, 2014;Smith et al., 2004), a single binary monovalent salt (NaCl), completely dissolved in water, is considered, such that z + = z − = 1 and c in and the osmotic pressure computed using the microstructure illustrated in Figure 1.The swelling or osmotic pressure that develops within the interlayer pores can be expressed as a function of the first order approximation of ion concentrations (Cheng and Hendry, 2014):

Electric potential and Ion concentration maps
ex Neumann boundary condition ( 22) applied on Γ ′ constitutes the main source  term for ψ (0) in , as attested by higher isovalues in the vicinity of Γ.However, through the source term in Eq. ( 21), the membrane potential depends also strongly on each ion concentration c ±(0) ex weighted by its valence z ± .As a result, the spatial distribution of cations and anions in interlayer pores, and the effective diffusion tensor will be affected by the external salt concentration, as has been experimentally observed in Glaus et al. (2007) for more compacted montmorillonites.

Effective diffusion tensor
According to Eq. ( 24), the cation concentration c Electron Micrograph from (Hetzel et al., 1994).A qualitative comparison can be made with experimental diffusivities obtained for various cations and anions with Na-Wyoming montmorillonite gels (Nakashima, 2003), which are comprised between 0.5D ± and D ± .In (Kemper and van Schaik, 1966;Kozaki et al., 2001) the experimental diffusivities reported for NaCl in montmorillonite are lower, and vary between 0.20D ± and 0.25D ± .The numerical diffusivities computed fall within the range of the experimental values obtained under the assumption of isotropy.More details would be required about the microstructural features of the smectite clays employed in the experiments to push the comparison further.The effective diffusion tensor obtained at the mesoscopic scale is then used to simulate diffusion at the scale of a Wyoming bentonite sample comprising clay gels (79.7 %), mineral grains (11.3 %) and micropores (9 %) (Pusch, 2001).Figure 6 illustrates the local fields ψ 1 and ψ 2 that have to be derived with respect to the mesoscopic spatial coordinates in order to obtain the macroscopic diffusion tensor (see Bouchelaghem (2018)).The diffusion coefficients computed at the scale of montmorillonite layers are assigned to the clay gels, and their contribution is predominant for macroscopic diffusion.The  (Sato and Suzuki, 2003;Smith et al., 2004) for smectite clays.As already observed in (Sato and Suzuki, 2003;Suzuki et al., 2004), there is no definite influence of the compaction level on the diffusion tensor anisotropy.This is mainly due to the presence of incompressible mineral grains and large inter-particle pores which contribute to orient the diffusion pathways in all directions, thereby smearing out the anisotropy present at the scale of montmorillonite layers (Suuronen et al., 2014).

Conclusion
A simple model has been proposed to describe ion transport by diffusion through montmorillonite, by relying on a realistic description of the interlayer space and taking into account the coupling between diffusion and electrostatic At the scale of clay gels, mineral grains and inter-particle pores, the microstructures need to be improved in order to be able to simulate diffusion in three-dimensional conditions.Other perspectives for the work include accounting for advection in order to compute the hydrostatic pressure gradient ∇ x p (0) ex and the membrane efficiency ∇xp (0) ex 2 RT ∇xc (0) ex (Cheng and Hendry, 2014).Finally, even at the level of the montmorillonite layers, it would be interesting to apply the model to 3D microstructures, as recent research indicates large ranges of variation for montmorillonite particles's sizes, that imply a three-dimensional network of interconnected pores (Whittaker et al., 2020).
Appendix A Development of the local problems Taking Eq. ( 11) at order O((ε ′ ) −2 ) and Eq. ( 15) at order O((ε ′ ) 0 ) = O(1), we have, using Einstein summation convention on repeated indices: Taking Eq. ( 11) at the order O((ε ′ ) −1 ) and Eq. ( 15) at the order O(ε ′ ), we have: ex (x i ), the previous system can be simplified as follows: ex depends linearly on ∂ψ (0) ex ∂xi .We can therefore introduce two functions that depend only on y i , β l (y) for l = 1, 2 in 2D, such that: ex (A1) Diffusion calculations on reconstructed bentonite microstructures and, by inserting Eq. (A1) in the system verified by ψ (1) ex , we have: By taking Eq. ( 9) at order O((ε ′ ) −2 ) and Eq. ( 13) at order O((ε ′ ) −1 ), we have: ex (x i ) does not depend on y i , we obtain simply: Similarly, Eq. ( 10) at order O((ε ′ ) −2 ) and Eq. ( 14) at order O((ε ′ ) −1 ) lead to Taking Eq. ( 9) at order O((ε ′ ) −1 ) and Eq. ( 13) at order O(1), and taking into account that ψ ex (x i ), we can write: According to Eq. (A1), ψ ex is a linear function of ex ∂xi , and we can assume: with c + l and d + l , l = 1, 2, Ω ′ -periodic characteristic functions.Similarly, from Eq. ( 10) at the order O((ε ′ ) −1 ) and Eq. ( 14) at the order O(1) and the same approach, we have: Eq. ( B11) is integrated in Ω ′ F .Using Green's theorem, the Ω ′ -periodicity of the physical variables and the Ω ′ -antiperiodicity of the outward normal vector n on Γ ′ , and boundary conditions (B12) and (B13), all the derivatives with respect to y i disappear.And we have: trial and error.Each grayscale TEM image is thresholded first to obtain a binary image, without prior filtering.The threshold value is found based on the histogram of the image, and the choice is confirmed by visual inspection.
varies within the microstructure under the action of the induced mesoscopic potential gradient ∇ x ψ (0) ex on the montmorillonite layers' surface (we use Einstein summation convention on repeated indices):

+
) at order O((ε ′ ) 0 ) = O(1) and averaging over the RVE, it is possible to obtain transport equations for c clay gels in the following forms (see Appendix B): 28)3.4 Numerical Resolution with the Finite ElementMethodOn the outer boundary of the RVE c i and ψ are assumed to be periodical.Using symmetry considerations, the periodicity boundary conditions are transformed into more classical Dirichlet and Neumann boundary conditions.We assume that the domain occupied by the montmorillonite layers is symmetrical with respect to the coordinate axes and the origin of the coordinate system is a center of symmetry.As a result, the computations are performed on a quarter of the Representative Volume Element (RVE) in 2D.The model is implemented in the Finite Element-based software Castem (http://www-cast3m.cea.fr).

Figure 3
Figure 3 illustrates the concentrations c +(0) of Na + ions and c −(0) of Cl − ions obtained at the first order of approximation when electrical couplings are taken into account.As expected from the assumed Boltzmann concentration distribution, cation inclusion and anion exclusion effects are clearly observed, particularly within the more densely stacked montmorillonite layers.
by parts in Ω ′ F while accounting for the associated boundary condition on Γ ′ , the Ω ′ -periodicity of ∂ψ (0) ex ∂yi and the Ω ′ -antiperiodicity of the normal vector n, we obtain that ∂ψ boundary condition on Γ ′ , the Ω ′ -periodicity of ∂c +(0) ex ∂yi and the Ω ′ -antiperiodicity of the normal vector n, we obtain that ∂c +