Ionic Transport in Nano-Porous Clays with Consideration of Electrostatic Effects

The mineralogical and chemical properties of clays have been the subject of longstanding study in the research community—in fact, entire journals are devoted to the topic. In the field of hydrology where transport behavior is more routinely considered, clays and clay-rich rock were largely relegated to a minor role because of their low hydraulic conductivity. However, this very property explains in part the renewed interest in the behavior of clays and clay rocks in several important subsurface energy-related applications, including the long-term disposal of nuclear wastes in geological repositories and the storage of CO2 in subsurface geological formations. In these applications and environments, the low permeability of the clay-rich formations or engineered barriers provides at least part of the safety functions for radionuclide contaminants confinement and subsurface CO2 sequestration. From a geochemical and mineralogical point of view, the high adsorption capacity of clay minerals adds to the effect of low hydraulic conductivities by greatly increasing the retardation of radionuclides and other contaminants, making clays ideal where isolation from the biosphere is desired. The low permeability of clay-rich shales also explains why hydrocarbon resources are not easily exploited from these formations, thus requiring in many cases special procedures like hydraulic fracturing in order to extract them. Clay properties remain also topic of intensive research in the oilfield industry in connection with their swelling behavior, which has an adverse impact on drilling operations (Anderson et al. 2010; Wilson and Wilson 2014; De Carvalho Balaban et al. 2015).

While the low permeability and high adsorption capacity of clay minerals are widely acknowledged, it is clear nonetheless that there is a need for an improved understanding of how the chemical and mineralogical properties of clay rocks impacts transport through them. It is at the pore-scale that the chemical properties …


INTRODUCTION
The mineralogical and chemical properties of clays have been the subject of longstanding study in the research community-in fact, entire journals are devoted to the topic.In the fi eld of hydrology where transport behavior is more routinely considered, clays and clay-rich rock were largely relegated to a minor role because of their low hydraulic conductivity.However, this very property explains in part the renewed interest in the behavior of clays and clay rocks in several important subsurface energy-related applications, including the long-term disposal of nuclear wastes in geological repositories and the storage of CO 2 in subsurface geological formations.In these applications and environments, the low permeability of the clay-rich formations or engineered barriers provides at least part of the safety functions for radionuclide contaminants confi nement and subsurface CO 2 sequestration.From a geochemical and mineralogical point of view, the high adsorption capacity of clay minerals adds to the effect of low hydraulic conductivities by greatly increasing the retardation of radionuclides and other contaminants, making clays ideal where isolation from the biosphere is desired.The low permeability of clay-rich shales also explains why hydrocarbon resources are not easily exploited from these formations, thus requiring in many cases special procedures like hydraulic fracturing in order to extract them.Clay properties remain also topic of intensive research in the oilfi eld industry in connection with their swelling behavior, which has an adverse impact on drilling operations (Anderson et al. 2010;Wilson and Wilson 2014;De Carvalho Balaban et al. 2015).
While the low permeability and high adsorption capacity of clay minerals are widely acknowledged, it is clear nonetheless that there is a need for an improved understanding of how the chemical and mineralogical properties of clay rocks impacts transport through them.It is at the pore-scale that the chemical properties of clay minerals become important since their electrostatic properties can play a large role.For all above examples, it is necessary to predict gas, water and oil transport properties in clay materials porosity as a function of a range of physical and chemical contexts.For example, the prediction of diffusion and retention properties of clay-rocks and engineered bentonite barriers are of paramount importance for waste storage applications in order to assess the long-term safety of the storage systems under consideration.In these systems, advective fl ux is negligible as compared to diffusive fl ux since the clays and clayrich rock have permeability values as low as 10 -21 -10 -19 m 2 (Delay et al. 2006(Delay et al. , 2014)).Examples include the clay-rock from the Callovian-Oxfordian layers in Bure (France) (Homand et al. 2006) and the clay-rock from the Opalinus Clay formation in Mont-Terri (Pearson et al. 2003).Safety assessment calculations must be carried out over periods ranging from a few hundreds to hundreds of thousand years (Altmann 2008).This is despite the fact that the interpretation of laboratory results are typically obtained i) on time periods ranging from days to years and ii) on samples with centimeter dimension.This suggests the need for a good knowledge of the fundamental properties of clays at the pore-scale is required in order to extend these analyses to the longer time and space scales.Clay transport properties are however not simple to model, as they deviate in many cases from predictions made with models developed previously for "conventional" porous media such as permeable aquifers (e.g., sandstone).For example, water fl ow in clay cannot be described consistently with Darcy's law (Neuzil 1986).In addition, model predictions must also take into account spatial heterogeneities, which can be time-dependent due to physical and chemical perturbations.In this respect, a signifi cant advantage of modern reactive transport models is their ability to handle complex geometries and chemistry, heterogeneities and transient conditions (Steefel et al. 2014).Indeed, numerical calculations have become one of the principal means by which the gaps between current process knowledge and defensible predictions in the environmental sciences can be bridged (Miller et al. 2010).In this regard, the present article approaches the topic of pore-scale transport through clays and clay-rich rocks by adopting a reactive transport approach (Steefel and Maher 2009).
This article begins with a discussion of the classic theory of diffusion in porous media, the limitations of which for clay media are highlighted with selected diffusion experiment results from the literature.The second section introduces specifi cs of clay minerals and clay materials that explain some of the inconsistencies found in using the classical diffusion theory.In a third section, constitutive equations for diffusion in clay porous media are proposed and the link between the predictions made with these equations and the experimental results is achieved through simple examples.In a fourth section, reactive transport (RT) model applications to the selected examples provide the basis for a discussion on diffusion process understanding and the current limitations of the proposed approach.In a fi nal section, a summary is provided, together with perspectives on RT models and further code development that is needed.

CLASSICAL FICKIAN DIFFUSION THEORY Diffusion basics
Diffusion processes are most often treated in terms of Fick's laws.Fick's fi rst law states that the diffusive fl ux of a species i in solution (J i ) is proportional to its concentration (c i ) gradient (here in 1-D along x) (Steefel et al. 2014 where D e,i is the effective diffusion coeffi cient that is specifi c to the chemical species i.The diffusion coeffi cient includes a correction for the tortuosity () and the porosity () of the porous media: , , 0 , , where 0,i D is the diffusion coeffi cient of species i in water (or self-diffusion coeffi cient), and , ).The tortuosity is defi ned as the square of the ratio of the path length the solute would follow in water alone, L, relative to the tortuous Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 289 path length, it would follow in porous media, Le (Bear 1972): Note that the terminology of the diffusion coeffi cient terms is very diverse (Shackelford and Moore 2013).The terminology presented here is the most commonly used in geosciences.In particular the effective diffusion coeffi cient is defi ned here to include the porosity.
Fick's second law is derived from the mass conservation law that includes the divergence of the fl ux: tot, , where C tot,i is the concentration of species i in the porous media (i.e., the amount of species i in the solution and in the solid normalized to the solution and solid volumes).If the species i is in the solution only then: If the species i is also adsorbed on or incorporated into the solid phase, then it is possible to defi ne a rock capacity factor  that relates the concentration in the porous media to the concentration in solution: tot, .
The quantifi cation of adsorption processes is commonly translated into distribution ratio values, Rd (L•kg -1 ): surf , , where c surf, i is the concentration on the surface of the element of interest (mol•kg -1 ).If the concentration of species i on the solid is only due to adsorption processes, then Equation ( 6) can be combined with Equation (7), yielding: , where  d is bulk dry density of the material.In that case, Equation (4) transforms into: , .
When interpreting diffusion data, the distribution ratio is commonly assumed to be constant (the adsorption is linearly dependent on the concentration) and representative of an instantaneous and reversible adsorption process.Under these conditions, the Rd value is designated as the distribution coeffi cient, K D .If it is further assumed that the media is homogeneous, Equation (9) reduces to:

Anion, cation and water diffusion in clay materials
Diffusion parameters for cations, anions and water in clay materials have been extensively studied in the literature.Various experimental setups have been used to determine porosity, diffusion coeffi cients, tortuosity and -K D values.In the following, results obtained by Tachi and Yotsuji (2014), who performed through-diffusion experiments in order to study the diffusion of HTO, I -, Na + and Cs + in montmorillonite samples compacted at a bulk dry density of d   0.8 kg•dm -3 , are discussed.In the context of the geometry of their experimental setup (Fig. 1), all parameters of Equation ( 10) can be determined simultaneously by fi tting experimental points to the analytical solution: where: and where  m fulfi lls: C(0,0) is the initial concentration in the inlet reservoir with volume V in (mol m -3 ); V out is the volume of the outlet reservoir (mol m -3 ), A is the cross-sectional area of the sample (m 2 ) and L is the thickness of the sample (m).
Diffusion parameters derived from fi tting of the data presented in Figure 2 with Equation ( 11) are given in Table 1.The rock capacity factor for water (~0.77) is consistent with the total porosity value that can be obtained by considering the clay mineral "grain" density ( g ~2.7-2.8 kg dm -3 ) and the bulk dry density of the material according to: 1 .
This result is in agreement with the case where the entire porosity is assumed (or treated) as fully connected for diffusion of water as a tracer (tritiated water, or HTO).The higher  value for cations than for water can be related to their adsorption to the solid surfaces (K D > 0).The lower  values for anions than for water indicate that anions do not have access to all of the porosity.This result is a fi rst direct evidence of the limitation of the classic Fickian diffusion theory when applied to clay porous media: it is not possible to model the diffusion of water and anions with the same single porosity model.The observation of a lower  value for anions than for water led to the development of the important concept of anion accessible porosity (sometimes also improperly named 'geochemical' porosity) to be compared to the water saturated or total porosity (Pearson 1999).The ratio of the anion accessible porosity to the total porosity depends not only on the nature of the clay minerals in the material, but also on the chemical conditions, particularly the ionic strength (Glaus et al. 2010;Tournassat and Appelo 2011).Thus, changes in chemical conditions can lead to signifi cant modifi cations of anion diffusion properties (Fig. 3).
The effective diffusion coeffi cients normalized to the self-diffusion coeffi cient values depend on the nature of the aqueous species, with effective diffusion coeffi cient values in the order (D e /D 0 ) Cs + > (D e /D 0 ) Na + > (D e /D 0 ) HTO > (D e /D 0 ) I -when measured under similar experimental conditions (Table 1).Higher effective diffusion coeffi cient values for cations than for neutral species and higher values for neutral species than for anions have been reported repeatedly in the literature for clay materials (Nakashima 2002;Van Loon et al. 2003, 2004a,b, 2005;García-Gutiérrez et al. 2004;Appelo and Wersin 2007;Glaus et al. 2007Glaus et al. , 2010;;Descostes et al. 2008;Wersin et al. 2008;Birgersson and Karnland 2009;Melkior et al. 2009;Appelo et al. 2010;Gimmi and Kosakowski 2011;Wittebroodt et al. 2012).A second fundamental problem related to the classic Fickian diffusion theory arises when the pore diffusion coeffi cient, D p , i , Equation (2), is estimated for cations based on their effective diffusion coeffi cient values and a porosity value.Since  cannot be greater than the total porosity, it follows that , , / p i e i D D  .For Cs + diffusion in an experiment conducted at 0.01 mol•L -1 NaCl, D p,Cs + values were higher than 3.1 × 10 -9 m 2 •s -1 , a value that is higher than the self-diffusion coeffi cient of Cs + in pure water (2.07 × 10 -9 m 2 •s -1 , Li and Gregory, 1974).This result, also reported repeatedly in the literature (Van Loon et al. 2004b;Glaus et al. 2007;Wersin et al. 2008;Appelo et al. 2010;Gimmi and Kosakowski 2011), is not physically correct and points out the inconsistency of the classic Fickian diffusion theory for modeling diffusion processes in clay media.Again, the large changes of Cs + diffusion parameters as a function of chemical conditions (D e,Cs + decreases when the ionic strength increases, Table 1, Figure 2) highlight the need to couple the chemical reactivity of clay materials to their transport properties in order to build reliable and predictive diffusion models.

Diffusion under a salinity gradient
Most of reported diffusion experiments have been performed under spatially constant ionic strength conditions.Recently, Glaus et al. (2013) reported experimental results of 22 Na + diffusion under a gradient of NaCl concentration.The experimental setup was similar to the one depicted in Figure 1: the inlet and outlet reservoirs contained a 0.5 mol•L -1 and 0.1 mol•L -1 NaCl solution respectively.At time t = 0, both solutions were spiked with the same concentration of 22 Na + and the concentrations in both reservoir were monitored.From Fick's diffusion equation, it would have been expected that 22 Na + diffuses from both reservoirs at an equal rate into the clay material, eventually producing a zero concentration gradient inbetween the reservoirs (dashed lines on Figure 4).However, the experimental observations were completely different: 22 Na + accumulated in the high NaCl concentration reservoir as it was depleted in the low NaCl concentration reservoir, evidencing non-Fickian diffusion processes.

K D values obtained from static and diffusion experiments
The adsorption properties of a material can be evaluated using batch (static) experiments.Batch K D values can be evaluated independently from diffusion experiments and then compared with  parameters derived from diffusion results.Unfortunately, this comparison often leads to K D values that differ from the diffusion experiment-based values, calling into question the usefulness of batch K D measurements to predict transport parameters of adsorbing Tabl e 1. Diffusion parameters for HTO, I -, Na + and Cs + in a montmorillonite sample compacted at 0.8 kg•dm -3 (sample thickness L = 10 mm, and diameter d = 20 mm).Parameters were obtained from Tachi and Yotsuji (2014) and correspond to the fi tting lines shown in Figure 2 and obtained with Equation (11).D 0,i values were taken from Li and Gregory (1974).

Tracer i
NaCl conc.
V in V out  D e,i D e,i / D 0,i species.For some reasons these K D discrepancies have often been attributed to the fact that the conditions of batch adsorption tests "have long been known to be unrepresentative of those existing in compacted clays", i.e., the surface properties and/or the adsorption site accessibility depend on the compaction (Shackelford and Moore 2013).This statement fi nds its origin in the batch and diffusion K D discrepancies observed primarily for Cs + in a range of published studies, with batch K D values typically higher than diffusion K D values (Miyahara et al. 1991;Tsai et al. 2001;Jakob et al. 2009;Aldaba et al. 2010).Studies of Cs + that compared adsorption in loose and compacted clay material, without relying on diffusion experiments, led to apparently contradictory results.Oscarson et al. (1994) found that K D values for Cs + on bentonite decreased with increasing compaction.In contrast, Montavon et al. (2006)   Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 295

From classic diffusion theory to process understanding
The limitations of the classic Fickian diffusion theory must fi nd their origin in the fundamental properties of the clay minerals.In the next section, these fundamental properties are linked qualitatively to some of the observations described above.

CLAY MINERAL SURFACES AND RELATED PROPERTIES
Electrostatic properties, high surface area, and anion exclusion Crystallographic origin of clay mineral electrostatic properties.The fundamental structural unit of phyllosilicate clay minerals consists of layers made of a sheet of edgesharing octahedra fused to one sheet (1:1 or TO layers) or two sheets (2:1 or TOT layers) of corner-sharing tetrahedra (Fig. 5).The metals in the octahedral sheet of clay minerals consist predominantly either of divalent or trivalent cations.In the fi rst case, all octahedral sites are occupied (trioctahedral clay minerals) whereas in the second case, only two-thirds of the octahedral sites are occupied (dioctahedral clay minerals).The clay minerals smectite and illite may constitute ~30 % of the material making up sedimentary rocks (Garrels and Mackenzie 1971).Those minerals have 2:1 layer structures and they are frequently dioctahedral.Kaolinite is also a very common clay mineral, the structure of which is dioctahedral and made up of 1:1 layers.For these three minerals, tetrahedral and octahedral cations are primarily Si 4+ and Al 3+ respectively.Their ideal structural formulae can be written Si 2 Al 2 O 5 (OH) 4 for kaolinite and Si 4 Al 2 O 10 (OH) 2 for dioctahedral illite and smectite.In the following, most examples consider illite, the principal constituent of most clay-rocks, and montmorillonite, a smectite that is the main constituent of bentonite, which is the most studied material in diffusion experiments.Illite and montmorillonite layers differ by the nature and amount of the isomorphic substitutions taking place in their octahedral and tetrahedral sheets: in montmorillonite, most of the substitutions occur in the octahedral sheet where Al 3+ is replaced by Fe 3+ or a cation of lower charge (Mg 2+ , Fe 2+ ); in illite, a signifi cant amount of the substitutions occur also in the tetrahedral sheet with Al 3+ or Fe 3+ replacing Si 4+ .Isomorphic substitutions by cations of lower charge results in negative layer charge.Smectites have negative layer charges ranging between 0.2 and 0.6 molc•mol -1 (the layer charge of montmorillonite is commonly in the range 0.2-0.45molc•mol -1 ), while illite has negative layer charges values between 0.6 and 0.9 molc•mol -1 (Brigatti et al. 2013).
Morphology of illite and montmorillonite particles.Clay mineral particles are made of layer stacks and the space between two adjacent layers is named the interlayer space (Fig. 5).Illite particles typically consist of 5 to 20 stacked TOT layers (Sayed Hassan et al. 2006).The interlayer spaces of illite particles are occupied by non-solvated cations (K + and sometimes NH 4 + ).In contrast, the interlayer spaces of montmorillonite particles are occupied by cations and variable amounts of water.The number of layers per montmorillonite particle depends on the water chemical potential and on the nature and external concentration of the layer charge compensating cation (Banin and Lahav 1968;Shainberg and Otoh 1968;Schramm and Kwak 1982a;Saiyouri et al. 2000).The variable amount of interlayer water in montmorillonite leads to variations in interlayer distance, i.e., swelling, with discrete basal spacing (crystalline swelling) from 11.8-12.6Å (one-layer hydrate), 14.5-15.6Å (two-layer hydrate), and up to 19-21.6 Å (four-layer hydrate) (Holmboe et al. 2012;Lagaly and Dékány 2013).In the case of Na + -and Li + -smectites, swelling can result in even larger basal spacing values in a continuous manner (osmotic swelling) (Méring 1946;Norrish 1954).
The TOT layer thickness from the center-to-center of oxygen atoms is approximately 6.5 Å.As a rough estimation, the total thickness of the TOT layer can be obtained by summing this value to twice the ionic radius of external oxygen atoms: hTOT = 6.5 + 2 × 1.5 = 9.5 Å.This dimension can be compared with the lateral dimension of the TOT layers: from 50 to 100 nm for illite (Poinssot et al. 1999;Sayed Hassan et al. 2006) and from 50 to 1000 nm for montmorillonite (Zachara et al. 1993;Tournassat et al. 2003;Yokoyama et al. 2005;Le Forestier et al. 2010;Marty et al. 2011).Consequently, illite and montmorillonite particles have high aspect ratios (from 2.5 to 1000) and their surface area is dominated by the contribution of the basal surfaces corresponding to the plane of TOT layer external oxygen atoms.The contribution of the surface area corresponding to the layer terminations (the edge surface) to the total surface area is minor (in the case of illite) if not negligible (in the case of montmorillonite).
The specifi c area of basal surfaces (SSA bs ) does not depend on the TOT layer lateral dimension, and it can be calculated from the structural formula and the crystallographic parameters, amounting to approximately 750 m 2 •g -1 (Tournassat and Appelo 2011).From the specifi c surface area and the layer charge, it is possible to calculate a specifi c surface charge.For a montmorillonite with a typical layer charge of 0.32 mol c •mol -1 , the specifi c surface charge is -0.11C•m -2 .With illite particles, only part of this surface, corresponding to the outer basal surfaces, is in contact with the porosity and interacts with water, while the other part that corresponds to the interlayer basal surfaces has no contact with water because of the collapse of the interlayer space.The specifi c outer basal surface area (SSA bs_outer ) depends on the average number of stacked layer (n st ) in a single illite particle: SSA bs_outer = SSA bs / n st .The outer basal surface area of illite particles can be measured by atomic force microscopy (AFM) and by gas adsorption experiments using the derivative isotherm summation (DIS) method, with a typical value of 100 m 2 •g -1 considered as representative of SSA bs_outer (Sayed Hassan et al. 2006).With montmorillonite particles, all the basal surfaces are in contact with water because interlayer spaces are hydrated.However, the number of stacked TOT layers in montmorillonite particles dictates the distribution of water in two distinct types of porosity: the interlayer porosity in contact with a specifi c surface area equal to SSA bs_inter = SSA bs ×(1-1/n st ) and the inter-particle Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 297 porosity in contact with a specifi c surface area equal to SSA bs_outer .It is not possible to give generic and representative values of SSA bs_inter and SSA bs_outer , since the value n st changes as a function of conditions, such as the degree of compaction, salt background nature and concentration (Banin and Lahav 1968;Schramm and Kwak 1982a,b;Bergaya 1995;Benna et al. 2002;Melkior et al. 2009).
From particle structure and morphology to anion exclusion.The negative charge of the clay layers is responsible for the presence of a negative electrostatic potential fi eld at the clay mineral basal surface-water interface.The concentrations of ions in the vicinity of basal planar surfaces of clay minerals depend on the distance from the surface considered.In a region known as the electrical double layer (EDL), concentrations of cations increase with proximity to the surface, while concentrations of anions decrease.This leads to a progressive screening of the surface charge by a solution having an opposite charge.At infi nite distance from the surface, the solution is neutral and is commonly described as bulk or free solution (or water).This spatial distribution of anions and cations gives rise to the anion exclusion process that is observed in diffusion experiments.Measurements of anion exclusion and electrophoretic mobility in aqueous dispersions of clay mineral particles indicate that the EDL has a thickness on the order of several nanometers with a strong dependence on ionic strength (Sposito 1992).As the ionic strength increases, the EDL thickness decreases, with the result that the anion accessible porosity increases as well.The EDL thickness, where more than 90 % of the surface charge is screened, is commensurable with 2-3 Debye lengths, where I is the non-dimensional ionic strength (Solomon 2001),  is the water dielectric constant (78.3 × 8.85419 × 10 -12 F•m -1 at 298 K), F is the Faraday constant (96 485 C•mol -1 ), R is the gas constant (8.3145J•mol -1 •K -1 ) and T is the temperature (K).
The ionic charge distribution in the EDL is related to the potentials of mean force for the various ions and those potentials are for the most part related to the local magnitude of the electrostatic potential.Unfortunately, there is no experimental method to measure directly the electrostatic potential: the values derived from experiments such as electrophoretic measurements (Delgado et al. 1986(Delgado et al. , 1988;;Sondi et al. 1996) are always model-dependent.The EDL can be conceptually subdivided into a Stern layer containing inner-and outersphere surface complexes, in agreement with spectroscopic results (Lee et al. 2010(Lee et al. , 2012)), and a diffuse layer (DL) containing ions that interact with the surface through long-range electrostatics (Leroy et al. 2006;Gonçalvès et al. 2007), in agreement with direct force measurements (Zhao et al. 2008;Siretanu et al. 2014).Molecular dynamics (MD) calculations can also provide information on the Stern layer and diffuse layer structure at the clay mineralwater interface (Marry et al. 2008;Tournassat et al. 2009a;Rotenberg et al. 2010;Bourg and Sposito 2011).The ion distribution in the diffuse layer obtained from MD simulation is in close agreement with the prediction of the simple modifi ed Gouy-Chapman (MGC) model prediction, where this model is applicable.In the MGC model, ion concentrations ( DL i c ) at a position y from the starting position of the diffuse layer follow a Boltzmann distribution: that is related to the charge of the ions (z i ) and the electrostatic potential ( DL  ) calculated from the Poisson equation: where 0 i c is the concentration of species i in the bulk water.
An equivalent anion accessible porosity can be estimated from the integration of the anion concentration profi le (Fig. 6) from the surface to the bulk water (Sposito 2004).In compacted clay material, the pore sizes may be small as compared to the EDL size.In that case, it is necessary to take into account the EDLs overlap between two neighboring surfaces.If the pores have all the same size, this calculation is straightforward.Clay mineral particles are, however, often segregated into aggregates delimiting inter-aggregate spaces whose size is usually larger than inter-particle spaces inside the aggregates.In clay-rocks, the presence of non-clay minerals (e.g., quartz, carbonates, pyrite) also infl uences the structure of the pore network and the pore size distribution.Pores as large as few micrometers are frequently observed (Keller et al. 2011(Keller et al. , 2013) ) and these co-exist with pores having a width as narrow as the clay mineral interlayer spacing, i.e one nanometer.In practice, this complex distribution of pore sizes makes it diffi cult to calculate the anion accessible porosity from bulk sample data (e.g., specifi c surface area, total porosity and pore water ionic strength) (Tournassat and Appelo 2011).

Adsorption processes in clays
Adsorption processes on basal surfaces.The high specifi c basal surface area and their electrostatic properties give rise to adsorption processes in the diffuse layer, but also in the Stern layer.The composition of the Stern layer can be calculated according to various models such as the double layer model (DDL), the Triple layer-or plane-model (TLM or TPM), and the charge distributed model (CD) etc., depending on the required level of details.In the following, the double layer model is selected.In the DDL model, all specifi cally adsorbed species are located on the same plane that corresponds also to the start of the diffuse layer.The quantifi cation of the adsorbed cationic species Me i on basal surfaces sites, >B -, is calculated Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 299 according to the surface complexation reactions (Dzombak and Hudson 1995;Appelo and Wersin 2007;Appelo et al. 2010;Tournassat et al. 2011): and their related equilibrium constants: where a is an activity term.The defi nition of surface species activity has always been problematic (Kulik 2009).Here, we consider that the activity of surface species is related to their concentrations on the surface and to the surface potential (  ) experienced by the species i.
We do not consider any non-electrostatic surface activity coeffi cient term, in accordance with most of geochemical codes conventions.The surface concentration term is calculated based on a coverage mole fraction convention in order to avoid thermodynamic inconsistencies with heterovalent reactions (Parkhurst and Appelo 1999;Kulik 2009;Tournassat et al. 2013;Wang and Giammar 2013): where c is a concentration term (mol•L -1 ) and TOT B c  is the total concentration of adsorption sites >B on the basal surface.The surface potential 0  is calculated according to: Equation ( 23) corresponds to the exact analytical solution of the modifi ed Gouy-Chapman model for an infi nite and fl at surface in contact with an infi nite reservoir of 1:1 electrolyte.As such, it must be remembered that Equation ( 23) is only an approximation for systems where the solution contains multi-valent species, for which the surface potential is lower, all other parameters (surface charge, ionic strength) remaining equal (Fig. 6).

Interlayer basal surface-solution interaction.
The interlayer space can be seen as an extreme case where the diffuse layer vanishes leaving only the Stern layer of the adjacent basal surfaces.For this reason, the interlayer space is often considered to be completely free of anions (Tournassat and Appelo 2011), although this hypothesis is still controversial (Rotenberg et al. 2007c;Birgersson and Karnland 2009).The reactions between the species in the interlayer space can be accounted for in the framework of the ion exchange theory (Vanselow 1932;Gapon 1933;Gaines and Thomas 1953;Sposito 1984).In this theory, negatively charged sites (X -) are fully compensated by counter cations (Me i ) in the vicinity of the sites according to the reaction and species on the exchanger sites can be exchanged with other species in solution: The sum of the sites X -is referred as the cation exchange capacity (CEC).As for the adsorbed surface species, the activity of exchangeable cations is an ill-defi ned thermodynamic parameter.
There is no unifying theory to calculate the activity coeffi cients of surface or exchange species.In the following, we will rely on the equivalent ratio form of activity in cation exchange theory, where activity is related to the site molar ratio that is occupied by the species i (E i ) and a surface activity coeffi cient ( i,exch ), following the Gaines and Thomas convention (Gaines and Thomas 1953): where exch c i is the concentration of species i one the exchanger in mol•L -1 of interlayer water.Note that there is no electrostatic potential term in Equation ( 25); the surface charge is considered to be fully compensated by the exchangeable cations leaving no surface potential.Usually, the equivalent fraction convention of Gaines and Thomas is preferred over the mole fraction convention of Vanselow (1932) in geochemical and reactive transport codes because the term  remains constant and is equal to the CEC (c CEC in mol•L -1 of interlayer water) whatever the composition of the exchanger.This choice is not dictated by any theoretical reasons and it must be seen as being arbitrary, the Gaines and Thomas convention being easier to implement numerically.In addition the activity coeffi cient term ,exch i  is often, if not always, dropped owing to the diffi culty to quantify it (Chu and Sposito 1981).Consequently, the chemical potential of exchanged species becomes: The equivalent fraction of species Me i on the exchanger sites are calculated according to the selectivity coeffi cients, exch i j K  , of binary reactions (Eqn.24): It can be emphasized that Equation ( 27) is equivalent to the combination of two Equations ( 22) for species i and j with the consideration of an electrostatic potential value of 0. This last condition corresponds indeed to the condition of full compensation of the surface charge by adsorbed cations.In this respect, the exchange model can be seen as a limiting case of the surface complexation model given above.
Adsorption processes at edge surfaces.Unlike basal surface area, the specifi c edge surface area, SSA edge , depends on the lateral dimension of the TOT layers.The edge specifi c surface area can also be measured by AFM and DIS methods (Bickmore et al. 2002;Tournassat et al. 2003;Le Forestier et al. 2010;Marty et al. 2011;Reinholdt et al. 2013).Average reported values of SSA edge value ranges from 5 m 2 •g -1 to 30 m 2 •g -1 .Although the edge surface area is quantitatively less important than the basal surface area, edge surfaces dominate in determining the surface complexation properties of clay minerals (Tournassat et al. 2013) and they cannot be ignored for the modeling of strongly interacting species such as divalent and trivalent metallic cations, lanthanides or actinides that interact with amphoteric sites at the clay mineral layer edges, or Cs + that interacts with size specifi c adsorption sites at frayed-edge sites on illite.
For metallic cations (e.g., Ni 2+ ), at least two adsorption sites-a high energy and a low energy amphoteric adsorption sites-are necessary for modeling adsorption isotherms as a function of pH and as function of concentration.Adsorption reactions on layer edge Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 301 amphoteric sites may be modeled with surface complexation models, for which a generic reaction stoichiometry is written as: where >SOH is an amphoteric surface site.A range of surface complexation models are available for adsorption processes on clay mineral layer edges (Bradbury and Baeyens 2005;Tertre et al. 2009;Gu et al. 2010).The most successful models in terms of range of conditions of applicability are the non-electrostatic models despite the abundant evidence of the presence of an electrostatic fi eld at clay mineral layer edges (Tournassat et al. 2013).
For Cs + , frayed edge sites reactivity is commonly modeled with a cation exchange model with high values of selectivity coeffi cients for Cs + and K + as compared to Na + (Brouwer et al. 1983;Poinssot et al. 1999;Bradbury and Baeyens 2000;Zachara et al. 2002;Steefel et al. 2003;Gaboreau et al. 2012;Chen et al. 2014).

Effect of nonlinearity of adsorption processes on experimental diffusion parameters.
The superposition of various adsorption processes can result in highly non-linear adsorption isotherms.In such cases, adsorption needs to be described with multi-site models that include multiple cation exchange sites (as described above) and/or another set of more selective adsorption sites.However, interpretations of diffusion data for adsorbing species almost invariably rely on Equation ( 10), which assumes that adsorption is linear and rapid relative to the time-scale of diffusion.The consideration in Equation ( 6) of a non-linearity in the adsorption process yields the following alternative to Equation (9): , .
According to Appelo et al. (2010), the term i c   is responsible, in a large part, for the observed discrepancy between K D values obtained from batch and diffusion experiments, and for which a difference in Cs + adsorption capacities from dispersed to compact system is often invoked (Miyahara et al. 1991;Tsai et al. 2001;Van Loon et al. 2004b;Maes et al. 2008;Wersin et al. 2008).It is also possible to interpret diffusion data with numerical reactive transport software, in which case there is no special restriction to linear and/or equilibrium treatments of adsorption.Non-linear adsorption is not restricted to the case where a range of adsorption sites are present with different affi nities for the tracer of interest.Even in the case of a single adsorption site, a non-linear adsorption isotherm can occur if the tracer concentration is high enough that the ions occupy a signifi cant fraction of the adsorption sites.The following simple example illustrates this problem.Na + / Ca 2+ cation exchange is considered to be the only reaction taking place at the clay mineral surface: 2 NaX + Ca 2+ CaX 2 + 2 Na + log K Na/Ca = 0.5.

(30)
The concentration of Ca 2+ on the exchanger is obtained from: where C CEC is the cation exchange capacity of the clay mineral in mol•L -1 pore water .Equation ( 31) can be combined with the K D equation and yields: Under the experimental conditions considered here, the Na + concentration is constant and thus  Ca 2+ /[Na + ] 2  Na + is constant.If the concentration of Ca 2+ on the exchanger is negligible as compared to the concentration of Na + , then   CEC NaX C  and Equation (32) transforms into: where q CEC is the CEC expressed in mol•kg -1 clay .All the terms in Equation ( 33) are constant and Ca 2+ adsorption is linear.As soon as [NaX] deviates from the C CEC value, however, Ca 2+ adsorption becomes non-linear.
In the following, the effect of the non-linearity of adsorption is explored by comparing the simulation of diffusion breakthrough curves using reactive transport (RT) modeling with PHREEQC (Parkhurst andAppelo 1999, 2013) in combination with i) a K D model or ii) a cation exchange model.The system modeled is similar to the one depicted in Figure 1.The Ca 2+ concentration was held constant in the inlet reservoir, either at 10 -7 mol•L -1 , at 10 -3 mol•L -1 , or at 5 × 10 -3 mol•L -1 , and the outlet concentration was simulated as a function of time.The inlet and outlet reservoir had a volume of 1 L and the diameter and length of the sample were set at 2 cm and 1 cm, respectively.The dry bulk density was set at 0.8 kg•dm -3, and the porosity value was 0.72.The reference Ca 2+ pore diffusion coeffi cient was set at 2 × 10 -10 m 2 •s -1 .According to Equation (33), and considering a q CEC value of 1 mol•kg -1 clay and a NaCl background concentration of 0.1 mol•L -1 , the K D value for Ca 2+ at trace concentration is K D = 100 L•kg -1 .While the K D model and the cation exchange models yielded identical results for the case with a Ca 2+ concentration of 10 -7 mol•L -1 in the inlet reservoir, the diffusion breakthrough curves were very different for the case with a Ca 2+ concentration of 10 -3 mol•L -1 and 5 × 10 -3 mol•L -1 (Fig. 7).At 5 × 10 -3 mol•L -1 , the two breakthrough curves could only be matched by decreasing the K D value to K D = 81 L•kg -1 .Note that the K D value of Ca 2+ at a concentration of 5 × 10 -3 mol•L -1 is 42 L•kg -1 according to Equation (32) and thus a constant-K D model calibrated on this value would have led to an incorrect position of the breakthrough curve.This simple example highlights the need to take into account the non-linearity of adsorption for the interpretation of experimental diffusion data.This result echoes the previous fi ndings of various authors (Melkior et al. 2005;Jakob et al. 2009;Appelo et al. 2010) who showed that Cs + diffusion breakthrough curves can only be adequately modeled when the non-linearity of its adsorption isotherm obtained from batch experiment is correctly taken into account.Non-linear adsorption isotherms are not easily introduced in analytical solutions of the diffusion equations, whereas RT codes can handle easily various and complex adsorption models such as cation exchange and surface complexation reactions (Steefel et al. 2014).

From real porosity distributions to reactive transport model representation
The pores in clay material and clay rocks are usually fully saturated as evidenced by the agreement between the porosity values derived from water loss measurements, density measurement (wet, dry, and grain densities), and water diffusion-accessible porosity measurements (Fernández et al. 2014).While these results imply that the pore network is fully connected, a full characterization of the connected pore geometry (particularly the pore throats) is still beyond the scope of what is possible, at least in the case of clay materials.Continuum reactive transport codes do not handle in full this complexity and average properties of the porosity must be considered (Steefel et al. 2014).Still it is possible to defi ne three porosity domains, or water domains, that can be handled separately: the bulk water, the diffuse layer water and the interlayer water, the properties for which can be each defi ned independently.One limitation of reactive transport models currently is that it is necessary to consider a non-zero volume for the bulk water and that the bulk water volumes are connected from one cell to the other.This representation of the system can be at variance with a real system in which the macropores (with bulk water) are inter-connected with only small pores with only EDL or interlayer water, i.e., a system where there is a discontinuity of the diffusion path in macropores.

Diffusive fl ux in bulk water
Fick's law as presented in Equation ( 1) is a strictly phenomenological relationship that is more rigorously treated with the Nernst-Planck equation.In the following, we will consider a pseudo 2-D Cartesian system in which diffusion takes place along the x axis only (Fig. 8).In absence of an external electric potential, the electrochemical potential in the bulk water can be expressed as (Ben-Yaakov 1981;Lasaga 1981) The gradient in chemical potential along the x axis is the driving force for diffusion and the fl ux of ions i in the bulk water b J i can be written as: , the diffusion coeffi cient of species i in the bulk water (Steefel et al. 2014).
In the absence of an external electric fi eld, there is no electrical current and so: b 0.
The combination of Equations ( 37) and ( 38) provides an expression for the gradient of the diffusion potential along the x-axis in the bulk water:

Diffusive fl ux in the diffuse layer
Flux equation.According to Equation ( 17), it is possible to write: The profi les of concentrations obtained from the Poisson-Boltzmann equation are representative of equilibrium between the solution in the diffuse layer and the solution at an infi nite distance from the diffuse layer, i.e., the bulk water: The combination of Equations ( 41) and ( 42) gives:

Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 305
And so, by following the MGC model we implicitly assume: The fl ux along x in the diffuse layer can be expressed in the same way as in the bulk water (Appelo and Wersin 2007): Since the chemical potential in the diffuse layer is the same as in the bulk water, it follows that: Again, in absence of an external electric fi eld and thus electrical current in the diffuse layer along the axis x, we obtain:

 
, the diffusion coeffi cient of species i in the diffuse layer water along the x-axis.The fl uxes of species i in the diffuse layer along the x-axis become: It can be noted that Equations ( 40) and ( 49 with the 'DL enrichment factor' for the fl ux in the diffuse layer (DL) and 1 i A  for the fl ux in the bulk water (b).
According to Equation (50), the diffusive fl ux can be split in three contributions for both diffuse layer and bulk waters: Calcu lation of  DL .A numerical resolution of the full Poisson-Boltzmann (PB) equation is possible, but it is often too demanding computationally for RT codes applications that aim to solve the Nernst-Planck equation together with a full geochemical reaction network.If the details of the ion distribution in the diffuse layer as a function of the position along y are not needed, a simplifi ed model can be considered where only average ion concentrations in the diffuse layer, DL i c , are calculated: where L DL is the length of the DL along the y axis and M  is the 'mean potential' in the diffuse layer.This model is often referred as a Donnan model (Muurinen et al. 2004;Appelo and Wersin 2007;Birgersson and Karnland 2009;Appelo et al. 2010;Tournassat and Appelo 2011) although it is based on an approximation of the PB equation and its underlying hypotheses.Note that Equation ( 51) can be solved for each individual ion, but that the resulting value of  M that is calculated can be different depending on the ions considered.Since only a single value of  M is desired, the right side of Equation ( 51) can only be an approximation of DL i c .
The average ion concentrations in the EDL are related to the surface charge,  (in C•m -2 ), according to: Note that h DL , the considered width of the diffuse layer, can be different from L DL , the actual width of the diffuse layer, in order to have a better agreement between the left and right terms of Equation ( 51).Re-substituting Equation ( 52) into Equation ( 49) makes it possible to express the species fl ux in the diffuse layer along the x-axis independently of the position within the diffuse layer along the y-axis, since an average electrostatic potential is used.

Interlayer diffusion
According to Equation ( 26), the gradient of chemical potential in the interlayer is: exch R , and the diffusive fl ux in the interlayer can be expressed as: The concentration of species i in the interlayer is related to the cation exchange capacity and to the equivalent fraction of species i on the exchanger, exch , and so: The diffusion potential term in the interlayer, exch diff x    , is obtained by considering the zero charge fl ux condition:

Approximations for Nernst-Planck equation for bulk and EDL water
Ion activity coeffi cient gradients.Ions activity coeffi cients are often calculated according to an extended Debye-Huckel model: with I the ionic strength: b 2 0 0.5 .
In principle, this expression can be substituted into Equation ( 50), but in practice the consideration of the gradient of ionic strength can be diffi cult as a function of the numerical scheme used for solving the transport and chemistry equations because the ionic strength is a function of all concentrations terms.Nonetheless, it is possible to do so within the standard RT numerical frameworks and there are some cases where it is necessary to do so (Molins et al. 2012;Steefel et al. 2014).
In most reactive transport codes, ionic strength and activity gradients are neglected and Equation ( 50 .
Equation ( 63) is for example the fl ux equation embedded in PHREEQC.It must however be remembered that Equation ( 62 . In this very special case Equation ( 63) is equivalent to Equation ( 61) and the activity gradient is neglected, is not equal to zero if there is a gradient of salinity.
If Equation ( 62) is solved analytically, one obtains: From Equation ( 64), it is clear that the approximation made in Equation ( 62) leads to activity coeffi cient gradient terms which are underestimated (in absolute value).As an example let's consider tracers whose charges are the same as the salt background ions, but whose concentrations are 1000 times lower than the background salt concentration.
. This result implies that for tracers, Equation (62) neglects the activity gradient terms and thus is very similar to Equation ( 63).

Diffuse layer and interlayer effective diffusion coeffi cients.
In PHREEQC, the diffusion coeffi cients in the diffuse layer are set proportional to the diffusion coeffi cients in the bulk water.This simplifying assumption is however not strictly necessary, and diffusion coeffi cients in the diffuse layer could be set to values which are not correlated with the diffusion coeffi cients values in the bulk water, as done in CrunchFlowMC.The latter case would imply that the tortuosity, and thus the effective diffusion coeffi cient, need not be the same in the bulk and EDL porosity.One can however justify the PHREEQC simplifi cation by noting that the decrease of the diffusion coeffi cient in the diffuse layer is similar for all species according to molecular dynamics calculations (Tournassat et al. 2009b;Bourg and Sposito 2011;Holmboe and Bourg 2014), although such calculations are incapable of taking into account the heterogeneous, potentially hierarchical structure of natural porous media.

Model system
An important consideration is that the relative contributions of the three terms of Equation ( 50) to the total diffusive fl ux are not the same in the diffuse layer as in the bulk water.In the following sections, we explore the coupling of these different terms using simple examples in which the fl ux equations are solved numerically so as to explain experimental observations which have been reported repeatedly for diffusion in clay materials.
We consider a one centimeter long system with a constant salt background concentration along the x-axis (0.1 mol•L -1 NaCl), but with a linear gradient of 22 Na and 36 Cl tracer concentrations along x from 10 -9 mol•L -1 at x = 0 to 0 mol•L -1 at x = 1 cm (Fig. 9).The pore diffusion coeffi cients are set to values ten times lower in the porous medium than in pure water, i.e., = 2.03 × 10 -10 m 2 •s -1 in both bulk water and diffuse layer porosity.For simplicity, the Davies activity coeffi cient model is used (A DH = 0.5, 0 i Ba = 1, b i = -0.3 for all ions).The surface charge value is chosen to be representative of a clay plug packed at a bulk density of  dry = 1 kg•dm -3 and with a cation exchange capacity (CEC) of 0.1 mol c •kg -1 , that is to say, values representative of an illite-rich clay material.The specifi c surface area is set to SSA = 100 m 2 •g -1 , a value in good agreement with the properties of illite particles.The chosen thickness of the diffuse layer, h DL, is assumed = 3 × 10 -9 m.The total porosity is Based on these parameters, it is possible to calculate the mean potential in the diffuse layer using Equation ( 52) as  M = -0.033V.

Example 1: Constant ionic strength
In the fi rst example, the ionic strength is taken as constant across the length of the system, which is equilibrated with a 0.1 mol•L -1 NaCl solution.The tracers concentrations are assumed to show a linear gradient along x from 10 -9 mol•L -1 at x = 0 to 10 -11 mol•L -1 at x = 1 cm.The contribution of the tracers to the ionic strength is negligible and thus, the activity coeffi cient gradient is negligible.Under these conditions, the diffusion potential term is also negligible in the bulk water and in the diffuse layer water because the term . While the 36 Cl -fl ux is higher than the 22 Na + fl ux in the bulk water because of the higher of diffusion coeffi cient values of the former, the reverse is true in the diffuse layer (Fig. 9).This effect is due to the term A in Equation ( 50) which is related to the accumulation of 22 Na + and the depletion of 36 Cl -in the diffuse layer water as compared to the bulk water.
Example 2: Gradient in ionic strength and tracer concentration In a second example, we add a linear gradient of NaCl concentrations along the x-axis from 0.1 mol•L -1 at x = 0 to 0.001 mol•L -1 at x = 1 cm, i.e., the tracers gradient and the NaCl gradient are now the same once scaled to their respective maximum concentrations.The diffusion fl uxes in the bulk water are dominated by the concentration gradient term (Fig. 10).In this case, the diffusion potential term is not negligible.The Na + diffusive fl ux in the bulk water is increased as a result of this term, while the Cl -diffusive fl ux is decreased because of the higher diffusivity of Cl -as compared to Na + coupled to the electroneutrality condition in the bulk water.In this example, the activity coeffi cient term is minor as compared to the two other terms, but it represents a negative contribution of up to 12% of the total 36 Cl -fl ux in the bulk water (up to 8% for 22 Na + , Fig. 11).Moreover the contribution of this term is not constant along the x-axis and its negative contribution to the total fl uxes decreases over the length of the system as the total salt background concentration decreases.Consequently, the initial linear gradient of NaCl concentration is not truly representative of steady state conditions.
In the diffuse layer water, the contributions of the different terms are radically different compared to the bulk water.While 22 Na + and 36 Cl -diffusive fl uxes remain equal in the diffuse layer, the 22 Na + diffusive fl ux is dominated by the concentration gradient (positive) term that is counteracted primarily by the diffusion potential (negative) term (recall that, in contrast, the diffusion potential term is positive in the bulk water for Na + ) and to a lesser extent by the activity coeffi cient term.For 36 Cl -, the concentration gradient and the diffusion potential terms have almost the same positive contribution to the diffusive fl ux.At x = 1 cm, Na + and Cl -diffusive fl uxes in the diffuse layer drop to a very low value because of the low ionic strength (I = 0.001), which in turn is responsible for the very low value of the enrichment factor A for Cl -and 36 Cl -in the EDL.Because of the surface charge compensation condition, the enrichment factor A remains high for Na + in the diffuse layer, but the charge conservation conditions and the related diffusion potential term cancel its effect on the total fl ux intensity.Similar to the bulk water, the calculated 22 Na + and 36 Cl -fl uxes in the diffuse layer are not constant as a function of the position along the x-axis and thus the system is not representative of steady state conditions because of the activity gradient term in the fl ux Equation (50).If this term is dropped as in Equation ( 61), then the fl uxes increase by up to 4% for the conditions considered here (Fig. 11).In this last example, we consider a system where the NaCl concentration follows the same linear concentration gradient as in Example 2. However, in contrast to Example 2, we consider 22 Na + and 36 Cl -concentrations gradients that are initially zero.This last condition cancels the concentration gradient term in the diffusion equation for the tracers in both bulk and diffuse layer waters, and the activity coeffi cient gradient and diffusion potential terms are the only remaining contributors to the total diffusive fl ux.In those conditions, the diffusive fl ux is positive for 22 Na + in the bulk water and negative for 36 Cl -while the reverse is true in the diffuse layer (Fig. 12).In these conditions, the diffusion potential term dominates over the activity coeffi cient gradient term (Fig. 13).

Links to experimental diffusion results
In Example 1 above, the fl uxes are independent of the position along the x-axis.Therefore, this system is representative of a through-diffusion experiment that has reached steady state conditions.The presence of a diffuse layer at a negatively charged surface increases the diffusivity of cationic species and decreases the diffusivity of anionic species as compared to bulk water, at least in the case where the diffusion coeffi cients are the same in the EDL and bulk water.Under the same conditions, the diffusive fl ux of a neutral species (such as HTO, a water tracer) is not impacted by the presence of a diffuse layer (the enrichment term A is equal to 1 and there is no diffusion potential term for neutral species).As such, this very simple example captures essential features of diffusion experiments performed on clay samples where the observed fl uxes at steady state are in the order J Cations > J HTO > J Anions (Glaus et al. 2010;Tachi and Yotsuji 2014).  2Na + and 36 Cl -concentrations and fl uxes in the bulk and in the diffuse layer waters.The salt background NaCl concentration follows a linear gradient from 0.1 to 0.001 mol•L -1 along the x-axis in this case where there is no tracer concentration gradient.
Example 2 points out the extent of the approximation made in the calculation of the diffusive fl ux if the activity coeffi cient gradient is dropped in the solution of the Nernst-Planck equation.The relative error made in the calculation of the diffusive fl ux is less in the EDL water than in the bulk water, thereby justifying this simplifi cation if the contribution of the EDL to the total fl ux dominates the fl ux in the bulk water.(Fig. 11) Example 3 shows that the presence of a salt background concentration gradient can lead to up-gradient as well as to down-gradient diffusion of charged tracers in the bulk and in the EDL water.If the contribution of the diffuse layer to the total volume of water is not negligible, cationic tracers experience up-gradient diffusion, because the diffusive fl ux is two orders of magnitudes higher in the EDL than in the bulk water (Fig. 12).For anions, the direction of total diffusive fl ux depends on the relative contributions of bulk and EDL waters.Diffusion does not follow the activity gradients, since the diffusion potential term governs the direction of the diffusive fl ux for the tracers in this case (Fig. 13).Again, this observation helps to justify the simplifying assumption of ignoring the activity coeffi cient gradient term in the Nernst-Planck equation.Up-gradient diffusion of 22 Na has been recently reported in an experimental system similar to that described in Example 3 (Glaus et al. 2013).

FROM DIFFUSIVE FLUX TO DIFFUSIVE TRANSPORT EQUATIONS
The analysis presented above of the different contributions of the terms in the fl ux equation highlights the feasibility of explaining qualitatively experimental observations of water and ion fl uxes with the Nernst Planck model developed from Equations ( 34) to (64).Quantitative simulations of those fl uxes including transient conditions can be achieved only by implementing the fl ux equation in transport equations.

Diffusive transport equation for porous medium with interlayer and EDL water
The 1D Cartesian diffusion equation in a saturated porous media has the form: where i C is the concentration of species i in the porosity (in mol•dm -3 ).If the species i is not adsorbed on a surface, and if it does not participate in dissolution-precipitation reactions or other processes which would reduce its concentration in solution (for example, radioactive . In the presence of a diffuse layer and of an interlayer space, the term i c  can be split into three contributions (bulk -b, diffuse layer -DL, and interlayer, or exchangeable, -exch): b b DL DL exch exch .
The same approach (and nomenclature) can be used for the fl ux term i J .Thus, the transport equation becomes: .
The solution of Equation ( 68) can pose several problems.Firstly, the porosity terms on the left hand side of the equation are not constant if the salinity changes as a function of time, and in this case the derivative of the porosity terms as a function of time must be solved (for example, an increase of salinity contracts the diffuse layer, and thus the bulk porosity value increases proportionally to the decrease of the diffuse layer porosity).Secondly, each of the fl ux terms includes porosity and tortuosity contributions within their effective diffusion coeffi cient terms: If the porous media is not homogeneous ), then it is diffi cult to defi ne unambiguously summation rules at the interface for the different terms.In particular, if one part of the modeled system contains only bulk water and the other part contains bulk and diffuse layer (and a fortiori interlayer water), it is not clear whether or not a fl ux between the bulk water and the diffuse layer water should be considered.This problem is exemplifi ed in the following section by considering simple illustrative examples.

Summation of bulk and diffuse layer diffusive fl uxes over an interface
Normalized fl ux term.We defi ne here the terms, i.e., fl uxes normalized to a porosity value of one: Simple case.We will fi rst consider an interface between two adjacent numerical domains (porous media 1 and 2) whose bulk water porosities are the same, and whose diffuse layer porosities are also the same (Fig. 14).In this particular case, we can safely assume that the surface available for bulk and DL water diffusion are proportional to the bulk and diffuse layer water porosities.Consequently, can be expressed as the weighed sum of b * i J and

Tournassat & Steefel
Complex case 1: Same total porosity, different diffuse layer porosity.In this fi rst complex case (Fig. 15), the calculation of the diffuse layer and bulk water contributions to the overall diffusion fl ux is not straightforward.A simple method can be the use of an arithmetic mean, which is also consistent with Equation ( 71 Complex case 2: Same total porosity, one zone without diffuse layer porosity.In the second complex case exemplifi ed in Figure 16, the calculation of the diffuse layer and bulk water contributions to the overall diffusion fl ux is again not straightforward.One can see this case as a limiting case of Complex case 1, and so one should write: However, if we consider that the bulk water of the porous medium 2 is in contact with the diffuse layer water of porous medium 1, the scheme of Figure 16 is equivalent to the scheme in Figure 14, but with the consideration that the diffuse layer 2 volume has the same properties as bulk 2 water (i.e., the electrostatic potential is zero: this condition can be easily considered for the solution of Equation ( 63)).In that case the fl ux summation becomes: Consequently, Equation (72) lacks a condition of continuity.This problem can be solved by considering the more complex equation: , where b to DL * i J is a fl ux term that is calculated with Equation ( 63) with the consideration of an electrostatic potential of zero on one site of the interface and a non-zero electrostatic potential on the other side of the interface.As such, Equation (75) fulfi lls the continuity conditions.

Complex case 3. Variable total porosity, differing diffuse layer and bulk porosities.
A still more complex case where the total porosity is not constant along the x-axis (Complex case 3, Fig. 17) cannot be unambiguously treated with Equation (75).In this case, the priority that should be given to the connectivity between diffuse layer or bulk water volumes (or interlayer porosity) is not clear and the convention that is chosen must then be seen as arbitrary.

Summary of cases.
The calculation of the fl uxes in the diffuse layer and in the bulk water can be obtained theoretically from the consideration of the Nernst-Planck equation coupled to the modifi ed Gouy-Chapman model (or its simplifi ed form, the mean potential model).In contrast, the summation of the fl ux at complex interfaces between two numerical cells does not follow rules that are dictated by any fi rmly grounded theory.The choices of summation rules which are given above, therefore, must be considered as intuitive choices.One should note that the complex case 2 illustrates conditions similar to a boundary condition between a fi lter (without diffuse layer) and a clay plug (with a diffuse layer), so it is far from an academic scenario.The Complex case 1 can be seen as representative of a medium with homogeneous surface charge properties, but one with a gradient in salinity.

Differentiation of the fl ux at interface between two numerical grid cells
The numerical differentiation of the fl ux poses another challenge.For a grid cell n, at position x n , the numerical equivalent of Equation ( 65) is given by: where superscripts new and old refer to two consecutive time steps, and where x n+1/2 and x n-1/2 are the positions at the interface between cells n and n + 1, and n and n -1, respectively.1/2 , n i x J  and 1/2 , n i x J  can be calculated according to Equations (57) (for the interlayer contribution) and (61) (for the bulk and diffuse layer contributions and neglecting the gradient of activity coeffi cient).For the bulk and diffuse layer contributions of the equations are represented numerically by fi nite differences as:  diffuse layer enrichment factor water concentration related to the electrostatic potential of the diffuse layer at the interface between cells n and n+1) are not known quantities, and their value must be approximated.A simple arithmetic averaging is frequently used for the bulk concentration that is strictly valid for an equally spaced grid: Alternatively, a harmonic mean can be used if the numerical grid is heterogeneous.The averaging of the diffuse layer enrichment factor 1/2 , n j x A  is more problematic if there is a gradient of electrostatic potential from one cell to the next, as in the case of a gradient of ionic strength.A linear gradient of ionic strength results in a highly non-linear gradient of electrostatic potential as exemplifi ed in Figure 18.In that case an arithmetic or harmonic mean may be inaccurate for estimating the A i value at the interface.The same kind of problem may occur for a gradient of surface charge between two grid cells.

APPLICATIONS Code limitations
While a range of reactive transport codes handle the Nernst-Planck equations for diffusive fl uxes in the bulk water, very few of them can handle transport processes in the diffuse layer (Steefel et al. 2014).Available publications with reactive transport simulations considering diffusion in the diffuse layer are limited to PHREEQC and CrunchFlowMC application studies (Appelo and Wersin 2007;Appelo et al. 2008Appelo et al. , 2010;;Alt-Epping et al. 2014).To our knowledge, interlayer diffusion processes are handled by PHREEQC only, and the interlayer diffusion option has been applied in only two published studies (Appelo et al. 2010;Glaus et al. 2013).In the following PHREEQC and CrunchFlowMC are used to illustrate the importance of considering coupled diffusion/surface reaction effects in order to understand and to predict migration processes and associated parameters in charged porous media, especially clays.

Simultaneous diffusion calculations of anions, cations, and neutral species
The data from Tachi and Yotsuji (2014), which were presented previously (Fig. 2), were modeled with CrunchFlowMC (Steefel et al 2014).The total porosity of the montmorillonite plug was set at 0.71, in agreement with the bulk dry density of the material (0.8 kg•dm -3 ).A typical total specifi c surface area of 750 m 2 •g -1 was assumed for montmorillonite.The total surface charge was set at 1 mol•kg -1 .Constants of the surface adsorption reactions were set at log K Na = 0.7 (in agreement with the MD results from Tournassat et al. 2009) and log K Cs = 2.5: Half of the total porosity was attributed to the diffuse layer, so that the mean Cl -(or I -) accessible porosity was ~0.41, in close agreement with the value given in Table 1.The tortuosity values were fi tted for each species and are reported in Table 2. Results are plotted in Figure 19.The tortuosity values follow the order  I -<  HTO <  22 Na + <  137 Cs + , indicating that the tortuous diffusion pathways are not the same for all of these species.Or alternatively, that the adsorbed species in the Stern layer, considered in the calculations as immobile, are in fact mobile.Nevertheless, the contribution of the diffuse layer to the diffusion fl ux calculation make it possible to derive a physically feasible value (< 1) for  137 Cs + , in contrast with the results obtained with a single porosity diffusion model.Note also that the assumed adsorption constant for Cs + is responsible for a K D value of 430 L•kg -1 at a ionic strength of 0.1, i.e., a value in agreement with K D values obtained from batch experiments from Tachi and Yotsuji (2014).
The bulk + diffuse layer water diffusion model presented here makes it possible to calculate the diffusive fl ux of neutral species, anions and cations with the same conceptual model and with physically feasible parameters.In this respect, the advantage of the RT codes is their ability to test these kinds of models under transient conditions and in the context of complex geometries in order to derive porosity and tortuosity values as well as adsorption parameters.Model benchmarking is thus not restricted to the comparison of data under steadystate conditions and/or the estimation of apparent diffusion coeffi cients where the adsorption and the diffusion parameters are lumped together.The effectiveness of this type of integrated approach has previously been put forward in the geochemical literature (Appelo and Wersin 2007;Appelo et al. 2008Appelo et al. , 2010)), but still remains the exception rather than the rule for the interpretation of diffusion data.
Figure 18.Averaging methods for the bulk concentration (left) and the diffuse layer enrichment factor, A i (right) for a monovalent cationic tracer at the interface between two grid cells whose centers are located at x n = 0 and x n+1 = 0.01 m.The variation of A i is due to a linear gradient of ionic strength from 0.1 (left of the system) to 0.001 (right of the system).

Diffusion under a salinity gradient
The experiments conducted and described by Glaus et al. (2013) were also modeled using CrunchFlowMC.The resistance to the diffusion by the fi lters was explicitly taken into account by assuming a tortuosity of 0.25 and a porosity of 0.32 for the fi lters.This corresponds to an effective diffusion coeffi cient of ~10 -10 m 2 •s -1 for Na + in the fi lter, a value in agreement with published values (Glaus et al. 2008).The total porosity in the montmorillonite plug was set to a value of 0.3, in agreement with the 1.9 kg•dm -3 dry bulk density of the material (Glaus et al. 2010).At this high degree of compaction, the actual presence of bulk water is not certain.However, it is necessary to have a non-zero bulk water volume to run CrunchFlowMC (or PHREEQC).This bulk water volume was set to a very low porosity value (0.02), so that the overall fl uxes were not impacted by the diffusion in this volume.The same surface adsorption model as was used for the modeling of Tachi and Yotsuji's data (see previous section) was used in this case.The tortuosity parameter was set to a value of 0.014 in the montmorillonite, in close agreement with the value measured for water, whose effective diffusion coeffi cient is about (1.5-1.7)•10 -1 m 2 •s -1 in similar conditions (Glaus et al. 2010(Glaus et al. , 2013)).The modeling results are shown on Figure 20 and agree almost perfectly with the experimental data.  2Na + , Cs + , HTO and I -through a montmorillonite plug equilibrated with a 0.1 mol•L -1 NaClO 4 solution in the experimental conditions from Tachi and Yotsuji (2014).See Figure 2 for the reference data and Table 1 for the experimental conditions.Glaus et al. (2013) have also modeled their data using two types of model.The fi rst model they considered was an interlayer diffusion model implemented in PHREEQC.The second model was based on that of Birgersson and Karnland (2009) (BK-model), in which the entire porosity was represented by a diffuse layer with a homogeneous mean electrostatic potential.The authors were able to obtain an equally good fi t of the data with both models.In Figure 20, a third alternate and equally good model (from the point of view of how well it fi ts the data) is given.It should be noted that the three models are totally different from a conceptual point of view: i) the interlayer diffusion model makes the assumption of a complete screening of the surface charge by the exchanged cations and the macroscopic driving force for cations mobility is the gradient of activity of the exchanged cations; ii) the BK-model makes the assumption of no screening of surface charge other than in the diffuse layer, and the macroscopic driving force for the diffusion is the gradient of concentration in the diffuse layer; and iii) the present model makes the assumption of partial screening of the surface charge by cations adsorbed in the Stern layer, where the mobility of cations is assumed to be negligible, and the macroscopic driving force of the diffusion is the diffusion potential in the diffuse layer, as shown in Figure 12.

Interlayer diffusion
Recently, Tertre et al. (2015) published data from their diffusion experiments carried out with centimeter-size mono-crystal of vermiculite.Vermiculite is a swelling clay mineral that has a high CEC (~1.8 mol•kg -1 ) originating primarily from isomorphic substitutions in the tetrahedral sheet.The size of the sample and the experimental setup (Fig. 21a) are ideal for probing self-diffusion processes in the interlayer porosity, as the setup makes it possible to eliminate the tortuosity parameter in the diffusion equation (the interlayer spaces are sandwiched between two fl at surfaces).The exchange sites in the vermiculite were saturated with Ca 2+ .The interlayer width in the vermiculite corresponded to a bi-layer hydrate.Following contact with a NaCl solution, a release of Ca 2+ was observed in solution.This experiment showed unambiguously that interlayer diffusion exists, and it made it possible also to quantify the diffusion coeffi cient for Ca 2+ in the interlayer.By immersing the vermiculite in a NaCl solution with a high salinity (0.1 mol•L -1 ), the authors were able demonstrate that the Ca 2+ interlayer diffusion coeffi cient was in good agreement with the value that they obtained from MD simulations.For experiments performed at lower ionic strength, they observed a discrepancy between the two values that they attributed to the diffusion of the solute species at the interface between the NaCl reservoir and the vermiculite interlayers.Their conclusion was based on the results of Brownian dynamics calculations.  2Na + under a salinity gradient for the experimental conditions considered by Glaus et al. (2013) using CrunchFlowMC.The blue plain curves correspond to the ratio C/C 0 of the classical Fickian diffusion model applied in the framework of a single pore diffusion model with linear adsorption processes (K D model) appears to be inappropriate for describing diffusion data without making less than satisfying modeling assumptions, such as i) different defi nitions of the porosity as a function of the nature of the tracer of interest (see anion accessible porosity) and as a function of the conditions (anion accessible porosity change with time; ii) change of the adsorption parameters from batch to diffusion experiment; and iii) physically unrealistic pore diffusion coeffi cients (or tortuosity values).
It is clear that RT codes have been capable of solving the problem of the transfer of K D values from batch experiments to compact porous media systems for some time.Much of the diffi culty stems from the non-linearity of the adsorption process for strongly adsorbed species on clay mineral surfaces, which can be handled readily by numerical reactive transport codes (Steefel et al 2014).Perhaps surprisingly, this has been done only recently for Cs + diffusion data (Appelo et al. 2010), although it was done some time in the past for column percolation experiments (Steefel et al. 2003).
Recent developments of selected RT codes that can handle diffusion processes in diffuse layer and interlayer porosities made it possible to model the diffusion data of neutral, anionic and cationic species within the same conceptual framework.Also, the contribution of the diffuse layer and/or the interlayer to the overall diffusion of ions makes it possible to explain the origin of the apparent acceleration of cation diffusion as compared to water, which otherwise would require unrealistic tortuosity values for cations in classical Fickian diffusion models.RT modeling has also helped to improve our understanding of anomalous diffusion behavior such as that of up-hill diffusion.
Despite the successes of these new RT modeling approaches, it must be stressed that the model and its parameters derived from diffusion experiments are not always unique.Two examples given above highlight the fact that several different conceptual models can provide equally good fi ts of the data.As such, the modeling effort is typically under constrained, a fact that explains the multitude of conceptual and numerical models available in the literature that describe the ionic transport properties of clay media (Leroy et al. 2006;Appelo and Wersin 2007;Gonçalvès et al. 2007;Birgersson and Karnland 2009;Gimmi and Kosakowski 2011;Tachi et al. 2014).
The transport properties of clay nanopores have been the matter of intensive research using advanced computational methods such as molecular dynamics, lattice-Boltzmann etc. (Bourg and Sposito 2010;Rotenberg et al. 2010;Obliger et al. 2013).The input of these techniques is clearly needed to constrain the macroscopic diffusion models (Fig. 21).Upscaling strategies have been developed to derive macroscopic diffusion parameters from microscopic information (Rotenberg et al. 2007a(Rotenberg et al. ,b, 2014;;Jardat et al. 2009;Bourg and Sposito 2010;Churakov and Gimmi 2011;Churakov et al. 2014) and these can be complemented by RT modeling that takes into account the complexity of the chemical reactivity of the material in fi ner details (e.g., adsorption processes, activity coeffi cients) However, it is noteworthy that the interpretation of diffusion data is complicated by complex microstructures that are dependent on physical and chemical conditions, and that these microstructures have not been yet characterized down to the scale of the smallest pores, i.e., the interlayer pores.The connectivity of the pore network connectivity across the full range of pore sizes has not been successfully determined for any of the investigated systems as well.These properties cannot be probed easily.For example, microscopic and tomography techniques still fail at imaging the porous network at the correct resolution (Hemes et al. 2013).The development of new observation techniques is thus necessary to make a step forward in the understanding of transport processes in clay materials.

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Example of a through diffusion cell setup: (a) inlet reservoir, (b) peristaltic pump, (c) throughdiffusion cell, and (d) outlet reservoir.Arrow heads indicates the circulation of water from the reservoir to the fi lter in order to homogenize the inlet and outlet solutions compositions.[Figure from Tachi and Yotsuji (2014) Diffusion and sorption of Cs + , Na + , I -and HTO in compacted sodium montmorillonite as a function of porewater salinity: Integrated sorption and diffusion model.Geochimica et Cosmochimica Acta Vol.132, p 75-93.Reproduced with the permission from Elsevier.] did not observe any signifi cant differences in K D values for differing degrees of compaction under otherwise similar experimental conditions.Van Loon et al. (2009) reached the same conclusion when comparing Cs + adsorption on crushed clay-rock from the Opalinus Clay formation (Mont-Terri, Switzerland) dispersed in water with Cs + adsorption on intact samples.Chen et al. (2014) concluded also that there was no effect of compaction on Cs + adsorption on i) the clay mineral fraction of a natural clay-rock (Callovian-Oxfordian clay-rock from Bure, France) and ii) on samples of the clay-rock itself.Whether the samples were i) powdered clay-rock samples dispersed in water, or ii) re-compacted powdered clay-rock samples, or iii) intact clay-rock samples had no impact on the measured K D values.Altogether, all recent Cs + adsorption experiments showed no effect of compaction on the K D values, thus it is necessary to fi nd a different reason for the discrepancy between K D values derived from batch and diffusion experiments.

Figure 4 .
Figure 4. 22 Na + diffusion under a gradient of salinity.The dashed lines indicate the expected concentration profi les as a function of time in inlet and outlet reservoirs.Symbols indicate the measured concentrations.[Reprinted with permission from Glaus MA, Birgersson M, Karnland O, Van Loon LR (2013) Seeming steady-state uphill diffusion of 22 Na + in compacted montmorillonite.Environmental Science & Technology, Vol.47, p 11522-11527, Fig. 1.Copyright 2013.]

Figure 5 .
Figure 5. Structure of a (dioctahedral) TOT layer and scheme of TOT layer and compensating cations in a clay mineral particle.[Figure adapted from Tournassat C, Bizi M, Braibant G, Crouzet C (2011) Influence of montmorillonite tactoid size on Na-Ca cation exchange reactions.Journal of Colloid and Interface Science, Vol.364, p 443-454, Fig 1. with permission from Elsevier.]

Figure 6 .
Figure 6.Electrostatic potential and chloride concentration in the diffuse layer calculated according to the MGC model, as a function of the distance from a clay mineral surface (x-axis), for two ionic strengths (left: I = 0.015; and right: I = 0.15) and with two different electrolytes (NaCl: blue plain lines; CaCl 2 : red dashed line).The specifi c surface charge is -0.1 C•m -2 .

Figure 7 .
Figure 7. Effect of the non-linearity of adsorption on the K D parameter derived from diffusion breakthrough curves.

u
is the mobility in bulk water (m 2 •s -1 •V) and where b diff  is the diffusion potential that arises because of the diffusion of charged species at different rates.The gradient in chemical potential along x is:

Figure 8 .
Figure 8. Pseudo-2-D Cartesian system with diffusion along the x-axis and electrostatic potential developing along the y-axis due to the negative charge at clay mineral surfaces.

.
) is a simplifi cation which neglects the cross-The simplifi cation made in Equation (62) leads to calculation results strictly equivalent to Equation (50) as long as it is evaluated numerically Ionic Transport in Nano-Porous Clays Considering Electrostatic Eff ects 309 and as long as the concentration gradient

Figure 9 .
Figure 9. Ex ample 1. Calculation of 22 Na + and 36 Cl -concentrations and fl uxes in the bulk and in the diffuse layer waters according to Equation (50).Ionic strength is constant along the x-axis (I = 0.1).The dashed line and the open green circles in the upper left Figure correspond to the gradient of the activity coeffi cient calculated according to Equation (60) or Equation (62), respectively.

Figure 10 .
Figure 10.E xample 2. Calculation of 22 Na + and 36 Cl -concentrations and fl uxes in the bulk and in the diffuse layer waters.The salt background NaCl concentration follows a linear gradient from 0.1 to 0.001 mol•L -1 along the x-axis.The dashed line and the open green circles in the upper left Figure correspond to the gradient of activity coeffi cient calculated according to Equation (60) or Equation (62) respectively.The blue open circles on the bottom-right Figure corresponds to the J 22 Na + values in the diffuse layer for comparison with the J 36 Cl -values.

Figure 11 .
Figure 11.Example 2. Calculation of 36 Cl -fl uxes in the bulk and in the diffuse layer waters according to Equation (50) (full equation, lines) and Equation (61) (no activity gradient term, open circles).Results for

Figure 12 .
Figure12.E xample 3. Calculation of22 Na + and 36 Cl -concentrations and fl uxes in the bulk and in the diffuse layer waters.The salt background NaCl concentration follows a linear gradient from 0.1 to 0.001 mol•L -1 along the x-axis in this case where there is no tracer concentration gradient.

Figure 13 .
Figure 13.E xample 3. Calculation of 36 Cl -fl uxes in the bulk and in the EDL porosity according to Equation (50) (full equation, lines) and Equation (61) (no activity gradient term, open circles).

Figure 14 .Figure 15 .
Figure 14.Scheme of an interface between two numerical domains.Simple case with no gradient in porosity properties.Porous medium 1Porous medium 2 co ncentration in the bulk water at the interface between cells n and n+1) and 1

Figure 16 .
Figure 16.Scheme of an interface between two numerical domains.Complex case 2.Porous medium 1Porous medium 2

Figure 17 .
Figure 17.Scheme of an interface between two numerical domains.Complex case 3.

Figure 19 .
Figure19.Modeling with CrunchFlowMC and a bulk + diffuse layer water diffusion model of the diffusion of22 Na + , Cs + , HTO and I -through a montmorillonite plug equilibrated with a 0.1 mol•L -1 NaClO 4 solution in the experimental conditions from Tachi and Yotsuji (2014).See Figure2for the reference data and Table1for the experimental conditions.

Figure 20 .
Figure20.Modeling of the diffusion of22 Na + under a salinity gradient for the experimental conditions considered byGlaus et al. (2013) using CrunchFlowMC.The blue plain curves correspond to the ratio C/C 0 of the diffuse layer porosity is DL

Table 2 .
Tortuosity values calculated with CrunchFlow for each species in the diffusion experiment from Tachi and Yotsuji (2014), and according to a bulk + diffuse layer water diffusion model.