Relevance of diffuse-layer, Stern-layer and interlayers for diffusion in clays: A new model and its application to Na, Sr, and Cs data in bentonite

The diffusive flux of cations is enhanced and that of anions is decreased compared to that of water tracers in bentonite or other clays, as a result of electrostatic interactions between ions and the charged clay surfaces. Clays are often used or foreseen as barriers in waste disposal in order to protect the environment from hazardous materials. A consistent description of diffusion of various ions and uncharged species is important in this context, especially if long-term interactions between clays and other materials shall be predicted by reactive transport simulations. Here, diffusion of a suite of tracers (HTO, Cl, Na, Sr and Cs) in bentonite was investigated. Experimental through-diffusion data at different bentonite dry densities were described with models of increasing complexity. First, ‘standard ’ empirical single ion transport models (uncoupled, simple Fickian diffusion) were applied for each density. These models served as reference cases for comparisons. For sorbing tracers (Na, Sr, Cs), surface diffusion models were used in a next step, where average surface mobilities for the cations were determined based on comparisons with transport parameters from HTO diffusion. Finally, a more complex model was developed in order to describe anion and cation diffusion in a coupled way. This model accounts for locally parallel diffusion in different environments, namely in ‘free ’ water unaffected by surface charges, in diffuse (Donnan) layer water, within the Stern layer, and within interlayer water (D-S-I model). This coupled model requires additional parameters related to the bentonite microstructure as well as to cation mobilities in the Stern layer and in the interlayer. The latter were taken from literature. Microstructural parameters were constrained in a manner that overall anion exclusion matches anion accessible porosities found by the simple Fickian diffusion model. This was possible with a reasonable choice of microstructure parameters that are consistent with literature values. A good agreement between the experimental data and the simulated cation diffusion co-efficients of the D-S-I model was found, using constraints by the HTO and Cl diffusion data. The interlayer pathway was found to be most important for diffusion of Na, Sr and Cs through bentonite. Stern layer diffusion was significant and more important than diffusion through the diffuse layer for Cs and Sr, while for Na the two pathways were of equal importance.


Introduction
Montmorillonite is the main component of bentonite.Pore water in montmorillonite and bentonite is to a large extent located in interlayers, but it occupies also pores on external surfaces of montmorillonite or other particles.Experimental porosity investigations with magnetic resonance (NMR) and small-angle X-ray scattering spectroscopy (SAXS) in MX80 bentonite and Ca-montmorillonite (Matusewicz et al., 2013;Muurinen et al., 2013;Matusewicz and Olin, 2019) suggest that the total pore space can be divided into an external or interparticle pore space and an interlayer pore space.Interlayer pore sizes of Na-and Camontmorillonite were estimated by X-ray diffraction (XRD) and NMR relaxometry (Ohkubo et al., 2016;Ohkubo et al., 2021) and XRD (Holmboe et al., 2012) for different degrees of compaction.Under high compaction, 1-3 hydrated (water) layers were found, while for less compacted Na-montmorillonite even a 4-layer hydration state has been reported.The swelling of a dry montmorillonite (d-value of ~1 nm) up to a 4-layer hydrate (d-value of ~2-2.2 nm) is generally considered to be of crystalline nature, that is, to be caused mainly by hydration of the cations and interaction of water molecules with the clay surfaces (e.g., Madsen and Müller-Vonmoos, 1989).Clay swelling above this hydration state is generally considered to be of osmotic origin (including repulsion by diffuse layers); it exist especially for Na clays.The interlayer pore size distribution was generally found to be dependent on clay dry density, prevailing cation, and ionic strength of the pore solution.
Adsorption of alkali and alkaline earth metals in bentonite and montmorillonite are usually described by cation exchange models (e.g., Bradbury and Baeyens, 2002;Missana and García-Gutiérrez, 2007;Van Loon and Glaus, 2008;Missana et al., 2014;Siroux et al., 2017).Missana and García-Gutiérrez (2007) found a 1-site cation exchange model sufficient for describing Ca and Sr sorption in FEBEX bentonite, while a 2site model was necessary to explain concentration dependent sorption of Cs (Missana et al., 2014).The second site was attributed to sorption on clay edges, which are dominant at low Cs concentrations.Siroux et al. (2017) developed a three-site cation exchange model in order to describe adsorption of Sr and Cs in Na-MX80 bentonite, where the third site accounts for increased sorption above a pH of about 8. Van Loon and Glaus (2008) found sorption of Cs in compacted bentonite to increase with increasing dry density.They explained this observation with changes in the thermodynamics of ion exchange at different degrees of mechanical compaction, that is, with a preference of cations with low hydration tendency over those with high hydration tendency in the smaller interlayer space at high bulk densities.Even though these models are able to quantitatively describe cation sorption in batch experiments, they do not give a detailed insight into cation distributions close to the surface and possible consequences on cation transport.
Theoretical models consider electrostatic interaction between the negatively charged clay surfaces and ions in pore solution.Diffuse double layer (DL) models (Tournassat et al., 2009) describe the excess of cations and the depletion of anions in the diffuse swarm using the Poisson-Boltzmann equation.Triple layer models (Davis et al., 1978;Leroy and Revil, 2004;Leroy et al., 2007) include an additional Stern layer (SL), where cations are more specifically sorbed as inner-sphere surface complexes (ISSC) or outer-sphere surface complexes (OSSC).The triple-layer structure has been confirmed in molecular dynamics (MD) simulation studies (Marry et al., 2008;Tournassat et al., 2009;Bourg and Sposito, 2011;Tinnacher et al., 2016;Lammers et al., 2017), where ion distributions showed distinct density peaks with distance from clay surfaces, which were attributed to ISSC/OSSC (Stern layer) and diffuse layer, respectively.In addition, MD simulations revealed that cations maintain a significant mobility in the Stern layer (Tournassat et al., 2009;Bourg and Sposito, 2011;Tinnacher et al., 2016) as well as in interlayers of smectites (Kosakowski et al., 2008;Bourg and Sposito, 2010;Holmboe and Bourg, 2014).
Diffusion of water, cations and anions in clays and clay rocks have been widely studied (Kozaki et al., 1998;Kozaki et al., 1999;Molera and Eriksen, 2002;Glaus et al., 2007;Van Loon et al., 2007;González et al., 2008;Melkior et al., 2009;Appelo et al., 2010;Glaus et al., 2010;Gimmi and Kosakowski, 2011;Bourg and Tournassat, 2015;Tinnacher et al., 2016;Wigger and Van Loon, 2018;Fukatsu et al., 2021).Results of these studies revealed that diffusion of cations is enhanced (i.e., shows a larger flux) while diffusion of anions is reduced in comparison to that of a water tracer.Therefore, when applying a simple Fickian diffusion model using concentration gradients of bulk aqueous phase species, comparably high (cation) and low (anion) effective diffusion coefficients were obtained.Diffusion experiments carried out under salt gradient conditions demonstrated that mobile surface-associated species may also act as a driving force in Fickian diffusion (Glaus et al., 2013).Diffusion coefficients of water, anions and cation tracers in bentonite all exhibit a decrease with increasing bentonite dry density, but to different degrees (Bourg and Tournassat, 2015).To explain the different behavior of cations, anions and water tracers, more sophisticated models were developed and applied.Based on pore structure and interlayer diffusivity Bourg et al. (2007) and Bourg and Sposito (2010) described diffusion of Na, Sr and Cs in bentonite clays.Anion exclusion on a pore scale was modeled for bentonite with different dry densities and at different ionic strengths by Tournassat and Appelo (2011) using microstructure information, considering DL effects and assuming interlayer being devoid of anions.Tachi and Yotsuji (2014) used an integrated sorption and diffusion model (Ochs et al., 2001) inlcuding DL diffusion and accounting for viscoelectric effects to model I, Na and Cs diffusion in Na-montmorillonite.Fukatsu et al. (2021) extended this model for Ca-montmorillonite by accounting for a heterogeneous pore structure.Soler et al. (2019) and Appelo et al. (2010) modeled anion and cation diffusion in Opalinus Clay with a dual continuum approach dividing the porosity in bulk and diffuse layer water.Interlayer diffusion is rarely considered explicitly in reactive transport models, but some studies contain interlayer diffusion modeling (e.g., Tournassat and Steefel (2015) using the reactive transport code PHREEQC and diffusion data from Tertre et al. (2015)).However, none of the above mentioned models explicitly accounted for diffusion in the Stern layer.Moreover, it is unclear how the dry density of the bentonite affects the relative importance of the different diffusion pathways, including the distribution of the different interlayer hydration states.
Here, through-diffusion experiments in bentonite at different bentonite dry densities (1.3, 1.6 and 1.9 kg L − 1 ) with HTO, Cl, Na, Sr and Cs tracers were investigated.Experimental results were evaluated with classical Fickian diffusion models for HTO and Cl and with a surface diffusion model (SD) for the cations Na, Sr and Cs.The parameters derived in this way varied with the dry density of the bentonite.In order to describe anion and cation diffusion in a more consistent way for all densities, a model including transport in the diffuse layer, the Stern layer, and the interlayer (D-S-I) besides transport in 'free' pore water was implemented in the continuum scale reactive transport code Flotran (Lichtner, 2007).The model is based on bentonite microstructure parameters combined with anion exclusion data and with Stern layer mobilities and interlayer mobilities taken from the literature.Predicted effective diffusion coefficients of this model for Na, Sr and Cs were compared to values determined from experimental data by the simple SD model.

Surface diffusion model
A fairly simple way to describe the experimental cation diffusion data in bentonite is the application of a surface diffusion model (Gimmi and Kosakowski, 2011;Krejci et al., 2021).This model incorporates, in addition to pore diffusion, the mass fluxes that result from enrichment of cations close to the clay surfaces.Their contributions are represented with one model parameter, the surface mobility μ s .At the local scale, pore and surface fluxes follow parallel pathways, but at the sample scale they follow also serial pathways, when considering the small particle sizes compared to a typical sample size.Equilibrium between both regions can be presumed because of their closeness.Then, the total mass flux for constant background concentrations can be written as: where j p and j s are the fluxes for the pore and surface regions, C is the concentration in solution, S is the sorbed concentration per mass of dry solid, D 0 is the cation diffusion coefficient in bulk water, ε is the total (water accessible) porosity, ρ bd is the bulk dry density, G is the so-called geometrical factor ('tortuosity') and accounts for the tortuous and constricted pathway in the pore water, ∂S/∂C is the derivative of the sorption isotherm, and ∂C/∂x is the concentration gradient in the pore water.The product of ρ bd and S equals the total sorbed amount (e.g., in moles) per bulk sample volume.
One problem when splitting the flux into separate fluxes of two or multiple transport domains as in Eq. ( 1) is the assignment of tortuosities (geometrical factors G) to each domain.As noted by Gimmi and P. Krejci et al.Kosakowski (2011), these tortuosities are only defined at the global scale of a sample, where we have a parallel and serial combination of the pathways, but not at the local scale for a single domain.We assume here that identical tortuosities for HTO and cations apply to the different pathways in Eq. (1), as in Gimmi and Kosakowski (2011), because the domains are considered to be in close contact and solutes to be in local equilibrium.This assumption may be an oversimplification and introduces uncertainty.A different tortuosity applies very likely for anions that are excluded from some pore regions.
The balance equation then is: Neglecting the second term on the right hand side in Eq. (1) and Eq.(2) one obtains the equations of the classical Fickian diffusion model with the effective diffusion coefficient D e being D e = εD p = εD 0 /G, and the pore diffusion coefficient The SD model was implemented in the reactive transport code Flotran using an implicit approach where pore and surface diffusion are combined in a single diffusion coefficient: The SD model, as written here, is based on the assumption that cations can diffuse in the same pore space as water plus in a surface layer, and therefore the porosity accessible to cations in the first term of Eq. ( 1) is the same as that for water ε = ε w .As mentioned, the geometrical factor G is a 'global' parameter at the sample scale and is typically derived from experimental diffusion data.

Combined D-S-I model
In order to gain a more detailed insight into the diffusive processes of cations in bentonite, and also to be able to predict the dependence of the diffusion behavior for a variation of the bentonite dry densities or the solution ionic strength, a more sophisticated model is necessary.Here, a new model is developed that accounts for cation fluxes in free (bulk, uncharged) water, in the diffuse layer, in the Stern layer, and in the interlayer space of smectites.Conceptually, the model splits the total pore space into an external pore space, which includes bulk, diffuse layer and Stern layer space, and an interlayer pore space, which may be characterized by a distribution of widths or hydration states.Again, it is assumed that the bulk, diffuse and Stern layer fluxes are parallel and the pore environments are locally in equilibrium.Local equilibrium is also assumed between the interlayer pore space and the pore water of the external pores.The total mass flux of the D-S-I model is then written as: where ε free and ε DL are the porosities of the pore space not influenced by the surface and in the DL, respectively, and G DL,SL,IL is the geometrical factor of the D-S-I model.Stern layer and interlayer fluxes are described via their derivatives of the sorption isotherms ∂S SL,IL,i /∂C and their specific mobilities μ SL,IL,i .The sum for the interlayer contributions originates from considering a distribution of interlayer hydration states i (i.e., average number of water layers).Note that the product of ρ bd S SL,IL,i equals the sorbed amount (e.g., moles) in the Stern layer or an interlayer i per bulk sample volume, and that this product could equally be expressed as ε SL,IL,i C SL,IL,i , i.e., as porosity of the corresponding pore environment multiplied by the corresponding concentration.Here, no explicit Stern layer pore volume is considered, and it is assumed that the total porosity ε tot is the sum of ε free , ε DL and all ε IL,i (and possibly, see Section 4, an additional porosity of interlayer-equivalent pores).A combined diffusion coefficient for the D-S-I model can be derived from Eq. (4): The cation concentrations in the diffuse layer C DL are determined by the amount of charge compensated in the diffuse layer.C DL can be approximated by a Donnan approach where ion concentrations are averaged over the volume of the diffuse layer and expressed via a Boltzmann type distribution: where Ψ D is the Donnan potential, z is the charge of the ion, F is the Faraday constant, R is the universal gas constant and T is the temperature.The distribution (enrichment or depletion) factor is defined as ξ = C DL /C.The Donnan potential Ψ D can be obtained from a charge balance: where Q is the net surface charge compensated in the DL per Donnan volume: Here CEC DL is the fraction of the total cation exchange capacity that is compensated by ions in the diffuse layer.
Anion diffusion is modeled by the same equation (Eq.( 5)), but assuming zero anion concentrations in the Stern layer and in the interlayer pore space (consisting of a few water layers only at the considered dry densities).Also, a different geometrical factor than for water or cations has to be taken into account.This is a consequence of the fact that G or the tortuosity is a property at the sample scale that depends on the network of pathways, which is definitively different for anions compared to cations or water tracers.
The D-S-I model is implemented in the reactive transport code Flotran.It is basically a combination of a Donnan approach (Gimmi and Alt-Epping, 2018) for the DL and of a multi-site SD-model (Krejci et al., 2021) for SL and IL, both verified in the mentioned papers.The implementation here is essentially a combination of the Donnan approach and the SD-model implementations.In order to model diffusion in free and DL water in parallel, an additional ('virtual') dimension is added to the simulation set up.Free and DL cells are treated to be in equilibrium with each other.In the DL cells an immobile negative charge is added to mimic the surface charge compensated in the diffuse layer (Gimmi and Alt-Epping, 2018), and the flux of ion concentrations are calculated using the Nernst-Planck equation.This leads then to the desired Donnan behavior of this cell.Equilibrium concentrations in the Stern layer and in interlayers are modeled via cation exchange on different mineral sites using the multi-site SD-model.The Stern layer mineral site is placed within the diffuse layer cell, while the interlayer mineral sites are placed in the free pore water cells.Overall 'sorption', that is, the overall distribution of ions in DL, SL and IL environments of the D-S-I model should be equal to that of a cation exchange model for any comparisons.As DL sorption does not reflect the ion specific selectivities reported for ion exchange models (when assuming identical activity coefficients in DL and free pore water, as generally done), selectivity coefficients for the Stern layer were adapted such that the total of DL, SL and IL cation concentrations match the distribution of a representative cation exchange model (Appelo et al., 2010;Glaus et al., 2020).Site specific mobilities of cations in the Stern layer and in the interlayer water, which allow to consider a mobility of these 'sorbed' ions, are incorporated in the overall mass flux according to Eq. ( 4), i.e., using the same approach as in the SD model.The above introduced model set up basically describes a two-dimensional domain for a one-dimensional transport P. Krejci et al. problem.Two parallel fluxes, namely, the flux in the cells of free pore water plus IL contribution and the flux in the cells of DL plus SL contribution are calculated.As the fluxes are parallel and free and DL cells are in equilibrium, these two fluxes add up to the overall mass flux (Eq.4).Alternatively, in situations of constant background electrolyte concentrations, a combined effective diffusion coefficient of the D-S-I model can be directly calculated from (Eq. 5), when all parameters are known.The verification of the implementation of the D-S-I model (combination of the Donnan approach and the multi-site SD model in a two-dimensional set up) in the code is found in the Supplementary Material S1.

Diffusion experiments and model input parameters
The Fickian and SD models were applied to experimental data of HTO, Cl, Na and Sr and Cs diffusion in bentonite reported by Glaus et al. (2017).The mineral composition and the pore water compositions at different dry densities for the here investigated Volclay KWK (Südchemie, Germany) bentonite, a Wyoming-type clay, can be found in Table 1 and Tables S1-S3.Its exchangeable cations are Na, Ca, Mg and K and its cation exchange capacity (CEC) is reported to be in the range of 0.76-1.20 eq kg − 1 (Van Loon et al. (2007).Experiments were carried out using a through-diffusion set-up.The cylindrical bentonite specimen had a diameter of 25.6 mm and a height of 10.4 mm for HTO, Cl and Na experiments and 5.3 mm for Sr and Cs experiments.Initial tracer concentrations in the upstream reservoir are found in Table S4.The detailed experimental set-up and description can be found in Van Loon et al. (2003) and Glaus et al. (2017).Tracer concentrations in the upstream reservoir, the tracer flux at the downstream boundary as well as the total concentration profile for the cations Na, Sr and Cs at the end of the experiment were measured.
The experimental fluxes in the through-diffusion experiments were modeled by considering the upstream reservoir, the filters and the bentonite clay specimen.Downstream, a zero-concentration boundary condition was applied for HTO and Cl.For the cations Na, Sr and Cs the measured concentrations were used as downstream boundary conditions, because a zero-concentration boundary condition may lead to systematically underestimated effective diffusion coefficients (Glaus et al., 2015).Chemical reactions (aquatic complexes and ion exchange reactions) were calculated based on activities a = γC/C 0 , where C 0 is the reference concentration equal to unity and the activity coefficient γ being described by the model of Davies (Stumm and Morgan, 1996).Ion specific bulk diffusion coefficients used for modeling are listed in Table 2. Filter diffusion coefficients were adapted within the value range of fresh and used filters given by Glaus et al. (2008) (Table S5).
Sorption for all cations is described by a one-site cation exchange model.A general cation exchange reaction between cations A and B with their respective charges z A and z B onto a surface sorption site ≡ S can be written as: The cation exchange capacity was calculated from the smectite content of the bentonite (71%) and the CEC of montmorillonite being 0.85 eq kg − 1 (Bradbury and Baeyens, 1997).Selectivity coefficients for the cations present in the pore water are taken from literature data, except for Sr, where selectivity was fitted to the concentration profiles (Table 3).The derivative of the sorption isotherm ∂S/∂C was calculated from the one-site cation exchange model numerically.A detailed description of this calculation can be found in Krejci et al. (2021).The same parameters were used for all simulations, except for the different Cs selectivity coefficients that account for a reported bentonite density dependence.

Results of HTO and Cl diffusion
From upstream reservoir concentration and downstream flux data, porosities and geometrical factors of HTO and Cl were determined (Table 4 and Figs. 1, S2 and S4).Reasonable matches between experimental data and the simple Fickian diffusion model were achieved.Accessible porosities of HTO and Cl decreased while geometrical factors increased with increasing bentonite dry density.Steady-state diffusion data of HTO showed some fluctuations, which made a precise determination of porosity and geometrical factor difficult.However, for the two replicate HTO experiments per bentonite dry density the same porosity ε HTO was determined.The values are in good agreement with porosities calculated from mineralogy (see Section 4, Eq. ( 11), Table 8) except for the highest bentonite dry density (about 20% difference).For further analysis with the SD model and D-S-I model the porosity based on the mineralogy was used at the highest bentonite dry density and geometrical factors were recalculated to match the steady-state mass flux (Table 4 values in brackets).

Results of Na, Sr and Cs diffusion
The porosity ε HTO and the geometrical factor G HTO determined by  Flury and Gimmi (2018).** from Applin and Lasaga (1984).HTO diffusion are used as input parameter for the surface diffusion model.So, the only unknowns are the surface mobilities of Na, Sr and Cs, which were determined by fitting the SD model to the experimental data (upstream reservoir concentration, downstream flux and sorbed concentration profile) (Figs. 2, S3 and S5).The three different pore waters exhibited significant amounts of NaSO − 4 (10%, 9% and 6% of total Na) and SrSO 4 (56%, 52% and 42% of total Sr) for the dry densities 1.3, 1.6 and 1.9, respectively, which were considered in the modeling with different diffusion coefficients (Table 2).For all cations fitted mobilities decreased with increasing bentonite dry density (Table 5).Mostly good, partly reasonable, agreement between experimental data and the fitted SD model was achieved (Figs. 2,S3 and S5).The discontinuities in the flux curves are a result of the downstream concentration boundary condition.At all bentonite dry densities sorbed Cs concentrations showed a significant decrease close to the filters.This decrease is a result from a decreased bentonite dry density close to the filters (Van Loon et al., 2007).Such effects are not included in the SD model.Other than that, ∂S Cs /∂C Cs is in very good agreement with distribution coefficients reported by Van Loon and Glaus (2008).

Application of D-S-I model to experimental results and discussion
The anion accessible porosities derived from the experiments by applying Fick's law are in good agreement with values from other experiments (Van Loon et al., 2007), and the mobilities of Na, Sr and Cs derived by applying the SD model are in good agreement with average mobilities reported by Gimmi and Kosakowski (2011).The cation mobilities seem well constrained as data and model show good agreement.However, their physical meaning, and therefore their predictive capability, is limited; they provide only a phenomenological description of experimental data at each bulk dry density, but cannot explain the variations between different densities.Such variations could be linked to variable contributions of different diffusion pathways (e.g., diffuse layer, Stern layer, interlayer) depending on the bulk dry density.Accordingly, a more predictive model should be able to quantify the significance of different pathways on the overall solute flux.The D-S-I model explicitly accounts for these different contributions.
In order to compare the D-S-I model with the experimental data or, more precisely, with the effective diffusion coefficients derived from the experimental data by the SD model, the combined effective diffusion coefficient of the D-S-I model (Eq.( 5)) is directly calculated by assuming equilibrium between concentrations in all water phases.For the calculation we make use of the fact that the overall sorption, and therefore the overall derivative of the sorption isotherm ∂S tot /∂C, is given by the onesite cation exchange model.This overall derivative must equal the sum of the sorption contributions of DL, SL and IL.The derivatives of the sorption isotherm for the various interlayer environments ∂S IL,i /∂C were calculated by multiplying ∂S tot /∂C with the fraction of CEC compensated in each interlayer environment.Thereby, it was assumed that the DL and SL sites (exposed to the external pore space) and the different interlayer sites have the same selectivity for the cations Na and Sr.For Cs, the same selectivity as in the SD model was used for the DL and SL sites, but, as Cs exhibits a preference for smaller interlayers expressed by higher selectivity coefficient (Van Loon and Glaus, 2008), individual selectivities were attributed to the different interlayer environments (Table S6).The selectivities of Cs for the different interlayer environments were calculated in a way that the weighted (by individual interlayer abundance) sum of all interlayer selectivities equals the overall selectivity used in the SD model at each bentonite dry density (Table 3).For the pore water at the different bentonite dry densities the distribution factor ξ and the Stern layer fraction (∂S SL /∂C)/(∂S tot /∂C) were derived from Donnan equilibrium calculations including a cation exchange site representing the Stern layer.The selectivities of this cation exchange site were adapted (Table S7) so that the overall adsorption in DL and SL matched the distribution of cations predicted by the one-site cation exchange model.Here, it was assumed that the activity coefficient in the diffuse layer is the same as in bulk water.This assumption is, however, arbitrary (Tournassat and Steefel, 2019) as the activity coefficients in the diffuse layer cannot be measured.
The determination of the distribution factor ξ and the fractions of Stern layer and the interlayer environments requires detailed knowledge of the bentonite microstructure, namely the distribution of the porosity and the charge over the different pore environments.

Bentonite pore structure model
The two key microstructural parameters to describe the bentonite structure are the stack size n c and the distribution of the different interlayer hydration states.They determine the external and internal surface area and porosity, respectively.However, there exists no such data specific to the compacted bentonite investigated here.Therefore, literature data from montmorillonite, the main component of bentonite,  is taken.From these data combined with anion exclusion data for bentonite (Van Loon et al., 2007), a model of the pore space is derived.From mass balance considerations the anion accessible porosity can be calculated from the sum of the free bulk porosity ε free and the DL porosity ε DL multiplied by the distribution factor of Cl ξ Cl .
The pore space model was constrained so that Eq. ( 10) matches the Cl accessible porosity determined with the Fickian diffusion model (Table 4).Stack sizes and interlayer pore size distributions were taken from the range of literature values so that this condition was met.
Table 6 shows a summary of the results of the pore space model and of one alternative configuration, which are both derived in the next subsections.

Bentonite microstructure
Pusch (2001) reported stack sizes for Na-montmorillonite around 3-5. Higher stack sizes n c of 6-8 were found by Matusewicz et al. (2013) for Ca-montmorillonite, independent of the clay dry density.In the bentonite investigated here, pore water is dominated by Na ions, therefore the stack size used in the model was chosen to be 4.5.Under high compaction 1 WL (water layer), 2 WL and 3 WL hydration states are present in montmorillonite (Holmboe et al., 2012;Ohkubo et al., 2016;Ohkubo et al., 2021), and their distribution strongly depends on the clay Fig. 2. Flux at downstream boundary (blue, first column), concentration in upstream reservoir (red, first column) and total (sorbed + aqueous) concentration profile (black, second column) for Na (first row), Sr (second row), and Cs (third row) at 1.6 kg L − 1 bentonite dry density.Symbols with error bars: experimental data; lines: surface diffusion model.Results for the other dry densities are shown in the Supplementary Material in Figs.S3 and S5.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)P. Krejci et al. dry density, the prevailing cations in the pore water and to some extent on the ionic strength.Bourg et al. (2006) accounted for the dry density effects based on data of Kozaki et al. (1998) by a linear transition from 3 WL to 2 WL between dry densities of Na-montmorillonite of 1.3 kg L − 1 and 1.6 kg L − 1 .Tournassat and Appelo (2011) added an ionic strength dependency to this model estimated from data of Kozaki et al. (1998), Kozaki et al. (2001) and Kozaki et al. (2008).More recent studies (Ohkubo et al., 2016;Ohkubo et al., 2021), however, suggested that ionic strength effects may be less significant under high compaction.At very high compaction of 1.8 kg L − 1 , Holmboe et al. (2012) found about 33% 1 WL interlayers in Na-montmorillonite.An exact quantification of the interlayer water distribution is difficult considering the nonhomoionic pore water of the bentonite and also pore size distributions for Na-montmorillonite estimated by XRD measurements at a given dry density (Holmboe et al., 2012;Ohkubo et al., 2021) differ significantly.For the lowest bentonite dry density, a configuration of the interlayer water distribution was chosen (Table 7) that lies between a pure 3 WL configuration (100% 3 WL, 0% 2 WL) and a combination of 55% 3 WL and 45% 2 WL that would have been estimated with the model of Tournassat and Appelo (2011).For the medium bentonite dry density fractions of 85% 2 WL and 15% 3 WL were used, similar to the findings of Ohkubo et al. (2016).At the highest bentonite dry density, 40% 1 WL and 60% 2 WL were taken, which is between the findings of Holmboe et al. (2012) and estimated values of Van Loon et al. (2007).

Porosity distribution
The total porosity (Table 8) of the bentonite at the different dry densities is calculated by: where ρ g is the bentonite grain density, estimated (from experimental data in Glaus et al. (2017)) to be 2.84 kg L − 1 .The total specific surface area of the montmorillonite fraction in the bentonite was taken to be 770 m 2 g − 1 , in agreement with values calculated from the montmorillonite crystal structure (Tournassat and Appelo, 2011;Holmboe et al., 2012;Leroy et al., 2015).The interlayer surface area can be calculated according to Tournassat and Appelo (2011): where a × b is the surface area of the montmorillonite unit cell with a value of 46.7 Å 2 (Holmboe et al., 2012), N A is the Avogadro constant (6.022⋅10 23 molecules mol − 1 ) and M W is the molecular weight of montmorillonite taken to be 733 g mol − 1 .This results in an interlayer surface area of about 590 m 2 g − 1 .The external basal surface area is the difference between the total and the interlayer surface area (180 m 2 g − 1 ).As adsorption of Na, Sr and Cs was well explained without considering edge sorption (c.f.3.3), the edge surface area was neglected.
The total porosity can be described as the sum of external ε ext and interlayer ε IL porosity: Here, the external porosity is further subdivided into the free bulk porosity ε free , the diffuse layer porosity ε DL and an external, but anion inaccessible (interlayer-equivalent) porosity ε ext,inacc .The derivation of these porosities is shown in the following paragraphs.
The total interlayer porosity (Table 8) is calculated as the sum over the different interlayer hydration states according to their fractions p i (Table 7): where the interlayer distances t IL are 0.32 nm, 0.64 nm and 0.92 nm for 1 WL, 2 WL and 3 WL interlayers, respectively, and w mont is the montmorillonite weight fraction in the bentonite.In the pore model presented here, it is assumed that anions do not enter the pore space between the layers.This can be justified by MD simulations of Tournassat et al. (2016), which showed that anions do not enter 2 WL interlayers and are depleted 20-fold in 3 WL interlayers, which is rather low, such that, it hardly leads to a significant effect on anion mass fluxes.
Van Loon et al. ( 2007) measured anion exclusion for variable NaCl solutions (ranging from 0.01 M to 1 M) using different techniques for the same bentonite densities as used here.From highest NaCl concentration (1 M) they derived a maximum anion accessible porosity ε max,Cl (Table 8).Diffuse double layer effects at such high ionic strength are negligible, as most charge is compensated in the Stern layer and diffuse layer thickness is very small.Therefore, one can expect that this value of ε max,Cl represents the whole external pore space.
In addition, a free (bulk) porosity ε free (not influenced by surface charge) (Table 8) is approximated from results at lowest NaCl concentrations of 0.01 M, where the diffuse double layer effect is high.
The diffuse layer porosity ε DL is then taken to be the difference between the maximum accessible porosity and the free bulk porosity:

Table 6
Diffuse layer parameters used to fulfill the mass balance in Eq. ( 10) (parameters of the alternative pore space model in brackets).However, the experimentally determined external porosity (assumed to be equal to ε max,Cl ) and the interlayer porosity ε IL (calculated from microstructure data) do not sum up to the total porosity ε tot (Table 8).
This discrepancy was also described by several other authors (Tournassat and Appelo, 2011;Matusewicz et al., 2013;Wigger and Van Loon, 2017;Wigger and Van Loon, 2018).Matusewicz et al. (2013) found stack sizes determined by XRD and SAXS being independent of dry density in Ca-montmorillonite, while stack sizes (and therefore interlayer porosity) determined by anion exclusion significantly increased at high dry density.Their explanation was that at high densities the distance between disordered platelets becomes small enough to exclude anions and this effect cannot be captured by XRD measurements.Tournassat and Appelo (2011) hypothesized that under high compaction rearrangement of clay platelets occurs leading to distances between adjacent surfaces in the order of interlayer spacings and accounted for that by an increased stack size.Wigger and Van Loon (2017) and Wigger and Van Loon (2018) had to introduce a so-called interlayer-equivalent porosity or bottleneck porosity, respectively in order to match experimental data on anion exclusion in Opalinus Clay.Here, in our model the stack size was not increased as it is rather unclear how disordered platelets in the external pore space rearrange with increasing density, e. g., parallel orientation of clay stacks with an interlayer-like distance or creation of bottleneck situations blocking interconnectivity of bigger pores.To account for the discrepancy an external, but anion-inaccessible porosity ε ext,inacc = ε tot − ε IL − ε max,Cl was introduced and calculated from the difference between the total porosity and ε IL and ε max,Cl .The external porosity is then ε ext = ε max,Cl + ε ext,inacc .A surface area proportional to its porosity fraction is attributed to ε ext,inacc .Without any knowledge about its detailed nature, its diffusion behavior for HTO, Na, Sr and Cs is assumed to be the same as that for a 2 WL interlayer pore.Alternatively, ε ext,inacc may be interpreted as a Stern layer porosity devoid of anions.In this case the whole external surface area is attributed to the diffuse and the Stern layer, as these compartments are considered to occur next to each other and not in series.In the following both a pore structure model with ε ext,inacc and with a Stern layer porosity (alternative pore model) are investigated.

Charge distribution
Charge distribution between diffuse and Stern layer can be approximated by a surface complexation reaction (Wersin et al., 2008;Appelo et al., 2010;Tinnacher et al., 2016): with logK SuNa being the respective surface complexation constant.Results of MD simulations of a montmorillonite nanopore in Na dominated pore water at ionic strength of I = 0.1 M (Tinnacher et al., 2016) showed about 40%-50% of surface charge being compensated in the Stern layer.Therefore, logK SuNa was chosen to be 0.9, which results in 45% charge compensation in the Stern layerat I = 0.1 M.
From the Na concentration in the pore water and the surface complexation constant the fraction of charge compensated in the diffuse layer (Table 9) can be calculated: The total cation exchange capacity is distributed equally over the total specific surface area as external CEC and interlayer CEC: The external CEC is further subdivided between the DL, SL and anion-inaccessible pore space according to the charge fraction f DL and the respective volumes, which leads to: In the case where the anion-inaccessible pore space is interpreted as Stern layer porosity (alternative pore model), ε ext,inacc is set to zero and the ratio ε max,Cl /ε ext is set to one in Eq. ( 19).The interlayer CEC is distributed over the different interlayer porosities according to their fractions p i (Table 7): 4.1.4.Site-specific mobilities for cations and water tracers Specific mobilities for the Stern layer and different interlayers were taken from literature when available (Table 10 (bracket values)).Mobility values of 1 WL, 2 WL and 3 WL interlayer for HTO and Na were taken from Holmboe and Bourg (2014), for Sr and Cs from Bourg and Sposito (2010).The latter literature data has however a large uncertainty, and therefore, mobilities of Sr and Cs were allowed to vary in the given range of uncertainty in order to achieve a better match with the experimentally determined mobilities.The Stern layer mobility of Na was taken from Tinnacher et al. (2016).For Sr a similar mobility as for Ca (Bourg and Sposito, 2011) was used, as both cations have similar diffusion coefficient, similar water chemistry, and similar adsorption behavior.No literature data is available for the Stern layer mobility of Cs; it was therefore treated to have the same value as in a 3 WL interlayer, because for Na and Sr Stern layer and 3 WL interlayer mobilities show very similar values.The mobility μ ext,inacc is set to that of a 2 WL interlayer in the basic pore model (cf.4.1.2).In the case where ε ext,inacc is interpreted as Stern layer porosity (alternative pore model) the mobility μ ext,inacc in Eq. ( 21) is set to the Stern layer mobility of water of 0.66 (Tinnacher et al., 2016).

Geometrical factor
Because water has a lower mobility in bentonite interlayer (Table 10) than in bulk water, the geometrical factor G HTO derived from HTO diffusion and used in the SD model is overestimated, as it does not account for the decreased water mobility in the interlayer environment.The D-S-I model uses a corrected geometrical factor G DL,SL,IL , which is written as: where μ are the mobilities of water in the different pore environments and f are the respective porosity fractions.The mobility of water in the external (anion accessible) pore space is assumed to be equal to that of bulk water (μ ext = 1, as for cations).

Comparison of experimental data with predictions based on the D-S-I model
The comparison between experimental data, represented by the combined effective diffusion coefficients of Na, Sr and Cs fitted from the diffusion experiments (Eq.(3), Table 5), and the ones predicted by the D-S-I via Eq.( 5) model shows good agreement (Fig. 3).Predicted D DL,SL,IL values are within a 15% error of the values fitted from the diffusion experiments with the SD or the Fickian model, except for Cs at highest bentonite dry density, where the predicted value from the D-S-I model is

Table 9
Fraction of charge compensated in the diffuse layer at the different bentonite dry densities for the different pore chemistries (Tables S1-S3).
ρ bd [kg L − 1 ] 1.3 1.6 1.9 fDL 0.46 0.42 0.36 P. Krejci et al. higher by a factor of about two and a factor of three for the alternative pore model.For the alternative pore model Sr diffusion at highest bentonite density was also overestimated by 50%.The reason for this big difference of Cs diffusion coefficients is not clear.The lower clay bulk density towards the filter boundaries is not considered in the evaluation of the experimental data, and may affect the diffusive behavior.It also has to be noted that the available mobility data for Cs (Table 10) has large uncertainty.Another reason that may play a role is that Cs sorption on edge sites should also be taken into account (Missana et al., 2014), even though the experimental data here (Figs. 2, S3, S5) could be well described with a 1-site cation exchange model.The results for Na are remarkable, because its site-specific mobilities were well constrained by MD simulation results.Diffusion coefficients were highest for Cs and lowest for Na for all bentonite densities.4, and Tables S8-S10 for both pore models.As expected from the large smectite content of the bentonite, for all cations at all densities interlayer diffusion mostly turned out to be dominant in the model with ε ext,inacc and no Stern layer porosity, except for Cs at highest dry density where Stern layer diffusion was dominant.For Cs, Stern layer diffusion became increasingly important with increasing bentonite dry density.Stern layer diffusion is also important for Sr.It makes up for about 20% of the total diffusivity at all dry densities, because Sr is comparably concentrated within the Stern layer (Table 11) and also has a relatively high mobility there.Diffusion in the DL was least important for Cs, as most Cs is located in the Stern layer.Stern layer adsorption of Na (Table 11) increases with increasing bentonite dry density, because the fraction of charge compensated in the DL decreases with dry density (Table 9).As Sr and Cs are much more selective, they are hardly influenced by this decrease of the fraction of charge compensated in the diffuse layer.
For the alternative pore model, Stern layer diffusion and diffuse layer diffusion become more important for all cations, because the surface charge of the total external surface area (Eq.( 19)) is attributed to diffuse and Stern layer.Consequently, Stern layer diffusion was most important for Cs at intermediate and highest bentonite density, and diffuse layer and Stern layer diffusion was most important for Sr at highest bentonite density.
Comparison of effective diffusion coefficients of the two alternative pore models shows very good agreement at low and intermediate bentonite density.However, the D-S-I model with the alternative pore model overestimates Sr and Cs diffusion coefficients at highest bentonite dry density.The two different pore models may be viewed as variation of pore size distribution of external pore space.The first pore model including ε ext,inacc accounts for a part of the external porosity that has pore sizes similar to interlayer pores and a part with somewhat bigger pores, while the alternative pore model does not include interlayer-like pores, but instead attributes this volume to the Stern layer.The comparison of the effective diffusion coefficients of Sr and Cs calculated for both pore models at highest bentonite density suggests that the pore size distribution of the external pore space is indeed not uniform and significant parts of the external pore space has pore sizes equivalent to interlayer pore sizes.Attributing this volume to the Stern layer leads then to an overestimation of the diffusion coefficients because of the larger cation mobilities within the Stern layer as compared to the otherwise assumed properties of a 2WL interlayer (Table 10).However, exact knowledge of the distribution of the external pore space in bentonite clays and its properties would be required for a more detailed analysis.Finally, we would like to stress that here, in the experiments with trace ions at constant background concentrations, a constant overall diffusion coefficients applies.This would be different in cases with varying background concentrations.Then, diffusion coefficients of the D-S-I model would become space and time dependent, and this would prevent the application of a simple Fickian model.Applying the D-S-I model to such cases would of course be the next step to test its applicability.

Summary and conclusions
Experimental through-diffusion data of HTO, Cl, Na, Sr and Cs in bentonite at different dry densities were modeled.A Fickian diffusion model was applied to determine water and anion accessible porosities as well as geometrical factors (tortuosities).A surface diffusion model was used to derive an average surface mobility of Na, Sr and Cs from experimental data, as well as the overall Fickian diffusion coefficient.The experimental data could be well described with a Fickian or a surface diffusion model.However, these models have no predictive capability, as their input parameters are not linked to structural and chemical parameters that may vary in different systems.A new model was developed and implemented in the reactive transport code Flotran.It accounts for diffusion in the diffuse layer, in the Stern layer, and in the interlayer space of clays.This model is mostly based on microstructural information (porosities and charge distribution in these pore environments) and values of cation mobilities in the Stern layer and interlayer environments that were derived from the available literature.Deriving the cation exchange selectivities for the different pore environments from literature data is more difficult, but they can be reasonably constrained.The pore structure model was further constrained by the maximum anion accessible porosity obtained form an anion diffusion experiment.One parameter that is difficult to determine independently is the tortuosity.Here, we proceeded by using the tortuosity of a water tracer (HTO) as input for the cations.With this set of input data, the D-S-I model was able to predict effective diffusion coefficients of Na, Sr and Cs that match the experimentally derived diffusion coefficients at three different bentonite dry densities.A successful (but challenging) link between diffusion of water, anions and cations was made in this way.In bentonite at the given densities, overall diffusion of the three cations is mostly dominated by the interlayers.Diffusion in the Stern layer is significant and more important than that in the diffuse layer for Cs and Sr, while these two pathways are of about equal importance for Na.Overall, the D-S-I model can be a powerful predictive tool for anion and cation diffusion when detailed information about the clay microstructure and cation mobilities are available.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Flux at downstream boundary (blue) and concentration in upstream reservoir (red) for HTO (first two columns) and Cl (third column) at 1.6 kg L − 1 bentonite dry density.Symbols with error bars: experimental data; lines: Fickian diffusion model.Results for the other dry densities are shown in the Supplementary Material in Figs.S2 and S4.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) While the fitted parameters D e of a Fickian diffusion model (or equally the fitted μ s of the SD-model) vary for the cations with background electrolyte concentrations and bentonite density (Table 5), the D-S-I model based on a single pore structure model and a set of mobility parameters leads to a consistent description of the experimental data at nearly all conditions.Moreover, as Na and Cs are both singly charged, the large difference between Na and Cs diffusion coefficients cannot be explained by a model including a diffuse layer only, unless largely different activities within this diffuse layer are assigned.The results clearly show the advantages of the D-S-I model over classical Fickian, SD or DL only models: it exhibits significantly larger predictive capability compared to these models.The derived effective diffusion coefficient of the D-S-I model can be seen as much more fundamental as it is solely based on independent diffusion data of HTO and Cl, microstructural data and mobility data of MD simulations.With the combined D-S-I model, we can furthermore compare relative contributions of the different pore environments to the overall effective diffusion coefficient of each cation.These comparisons are shown in Fig.

Fig. 3 .
Fig. 3. Comparison of effective diffusion coefficients for Na, Sr and Cs determined from diffusion experiments with a simple Fickian (or an SD) diffusion model (black dots) and predicted by the D-S-I model (blue diamonds and magenta squares for the alternative pore model).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Relative contributions (in %) to the effective diffusion coefficient (Eq.(5)) of the D-S-I model for Na, Sr and Cs.The values of the relative contributions (also for the alternative pore model) are found in the Supplementary Material TablesS8-S10.

Table 2
Species specific bulk diffusion coefficients used in modeling.

Table 3
Selectivity coefficients used in the cation exchange model with respect to Na.

Table 4
HTO and Cl accessible porosities and geometrical factors.Values in brackets represent porosity based on mineralogy and corresponding tortuosity.

Table 5
Surface mobilities of Na, Sr and Cs determined with the SD model and derivatives of the sorption isotherms based on the one site cation exchange model.The last lines show the effective diffusion coefficients that result from the fitted parameters of the SD model, or equally from a direct fit of the data with the simple Fickian diffusion model.

Table 8
Porosities derived for the different pore environments.In the alternative model, ε ext,inacc is considered to represent the Stern layer porosity.

Table 10
Site specific mobilities used in the D-S-I model, literature values in brackets.The asterisk indicates the Stern layer mobility of HTO used in the alternative pore model.

Table 11
Percentage of Na, Sr and Cs distributed in the Stern layer on external surfaces (alternative pore model in brackets).
P.Krejci et al.