This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.


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 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 flux is negligible as compared to diffusive flux since the clays and clay-rich rock have permeability values as low as 10−21–10−19 m2 (Delay et al. 2006, 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 flow 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 significant 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 specifics 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 final section, a summary is provided, together with perspectives on RT models and further code development that is needed.


Diffusion basics

Diffusion processes are most often treated in terms of Fick’s laws. Fick’s first law states that the diffusive flux of a species i in solution (Ji) is proportional to its concentration (ci) gradient (here in 1-D along x) (Steefel et al. 2014):



where De,i is the effective diffusion coefficient that is specific to the chemical species i. The diffusion coefficient includes a correction for the tortuosity (τ) and the porosity (φ) of the porous media:



where D0,i is the diffusion coefficient of species i in water (or self-diffusion coefficient), and Dp,i is the pore diffusion coefficient (Dp,i = τ ·D0,i ). The tortuosity is defined as the square of the ratio of the path length the solute would follow in water alone, L, relative to the tortuous path length, it would follow in porous media, Le (Bear 1972):



Note that the terminology of the diffusion coefficient terms is very diverse (Shackelford and Moore 2013). The terminology presented here is the most commonly used in geosciences. In particular the effective diffusion coefficient is defined here to include the porosity.

Fick’s second law is derived from the mass conservation law that includes the divergence of the flux:



where Ctot,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 define a rock capacity factor α that relates the concentration in the porous media to the concentration in solution:



The quantification of adsorption processes is commonly translated into distribution ratio values, Rd (L·kg−1):



where csurf,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 coefficient, KD. 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 coefficients, tortuosity and α–KD 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 fitting experimental points to the analytical solution:








and where θm fulfills:



C(0,0) is the initial concentration in the inlet reservoir with volume Vin (mol m−3); Vout is the volume of the outlet reservoir (mol m−3), A is the cross-sectional area of the sample (m2) and L is the thickness of the sample (m).

Diffusion parameters derived from fitting 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:



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 (KD > 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 first 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 significant modifications of anion diffusion properties (Fig. 3).

The effective diffusion coefficients normalized to the self-diffusion coefficient values depend on the nature of the aqueous species, with effective diffusion coefficient values in the order (De/D0)Cs+ > (De/D0)Na+ > (De/D0)HTO > (De/D0)I when measured under similar experimental conditions (Table 1). Higher effective diffusion coefficient 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. 2007, 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 coefficient, Dp,i, Equation (2), is estimated for cations based on their effective diffusion coefficient values and a porosity value. Since φ cannot be greater than the total porosity, it follows that Dp,iDe,i /φ. For Cs+ diffusion in an experiment conducted at 0.01 mol·L−1 NaCl, Dp,Cs+ values were higher than 3.1 × 10−9 m2·s−1, a value that is higher than the self-diffusion coefficient of Cs+ in pure water (2.07 × 10−9 m2·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 (De,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 22Na+ 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 22Na+ and the concentrations in both reservoir were monitored. From Fick’s diffusion equation, it would have been expected that 22Na+ diffuses from both reservoirs at an equal rate into the clay material, eventually producing a zero concentration gradient in-between the reservoirs (dashed lines on Figure 4). However, the experimental observations were completely different: 22Na+ accumulated in the high NaCl concentration reservoir as it was depleted in the low NaCl concentration reservoir, evidencing non-Fickian diffusion processes.

KD values obtained from static and diffusion experiments

The adsorption properties of a material can be evaluated using batch (static) experiments. Batch KD values can be evaluated independently from diffusion experiments and then compared with α parameters derived from diffusion results. Unfortunately, this comparison often leads to KD values that differ from the diffusion experiment-based values, calling into question the usefulness of batch KD measurements to predict transport parameters of adsorbing species. For some reasons these KD 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 finds its origin in the batch and diffusion KD discrepancies observed primarily for Cs+ in a range of published studies, with batch KD values typically higher than diffusion KD 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 KD values for Cs+ on bentonite decreased with increasing compaction. In contrast, Montavon et al. (2006) did not observe any significant differences in KD 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 KD values. Altogether, all recent Cs+ adsorption experiments showed no effect of compaction on the KD values, thus it is necessary to find a different reason for the discrepancy between KD values derived from batch and diffusion experiments.

From classic diffusion theory to process understanding

The limitations of the classic Fickian diffusion theory must find 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.


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 edge-sharing 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 first 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 Si4+ and Al3+ respectively. Their ideal structural formulae can be written Si2Al2O5(OH)4 for kaolinite and Si4Al2O10(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 Al3+ is replaced by Fe3+ or a cation of lower charge (Mg2+, Fe2+); in illite, a significant amount of the substitutions occur also in the tetrahedral sheet with Al3+ or Fe3+ replacing Si4+. 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.45 molc·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 NH4+). 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 specific area of basal surfaces (SSAbs) 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 m2·g−1 (Tournassat and Appelo 2011). From the specific surface area and the layer charge, it is possible to calculate a specific surface charge. For a montmorillonite with a typical layer charge of 0.32 molc·mol−1, the specific surface charge is −0.11 C·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 specific outer basal surface area (SSAbs_outer) depends on the average number of stacked layer (nst) in a single illite particle: SSAbs_outer = SSAbs / nst. 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 m2·g−1 considered as representative of SSAbs_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 specific surface area equal to SSAbs_inter = SSAbs ×(1-1/nst) and the inter-particle porosity in contact with a specific surface area equal to SSAbs_outer. It is not possible to give generic and representative values of SSAbs_inter and SSAbs_outer, since the value nst 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 field 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 infinite 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, κ−1:



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.3145 J·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, 1988; Sondi et al. 1996) are always model-dependent. The EDL can be conceptually subdivided into a Stern layer containing inner- and outer-sphere surface complexes, in agreement with spectroscopic results (Lee et al. 2010, 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 mineral-water 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 modified Gouy-Chapman (MGC) model prediction, where this model is applicable. In the MGC model, ion concentrations (DLci ) 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 (zi) and the electrostatic potential (ψDL) calculated from the Poisson equation:



where ci0 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 profile (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 influences the structure of the pore network and the pore size distribution. Pores as large as few micrometers are frequently observed (Keller et al. 2011, 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 difficult to calculate the anion accessible porosity from bulk sample data (e.g., specific surface area, total porosity and pore water ionic strength) (Tournassat and Appelo 2011).

Adsorption processes in clays

Adsorption processes on basal surfaces

The high specific 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 specifically adsorbed species are located on the same plane that corresponds also to the start of the diffuse layer. The quantification of the adsorbed cationic species Mei on basal surfaces sites, >B, is calculated 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 definition 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 (Ψ0) experienced by the species i. We do not consider any non-electrostatic surface activity coefficient 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 c>BTOT 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 modified Gouy–Chapman model for an infinite and flat surface in contact with an infinite 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 (Mei) 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-defined thermodynamic parameter. There is no unifying theory to calculate the activity coefficients 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 (Ei) and a surface activity coefficient (γi,exch), following the Gaines and Thomas convention (Gaines and Thomas 1953):



where exchci 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 jzjcexchj remains constant and is equal to the CEC (cCEC 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 coefficient term γi,exch is often, if not always, dropped owing to the difficulty to quantify it (Chu and Sposito 1981). Consequently, the chemical potential of exchanged species becomes:


μexchi=μexchi0+RTln Ei.

The equivalent fraction of species Mei on the exchanger sites are calculated according to the selectivity coefficients, exchKij, 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 specific edge surface area, SSAedge, depends on the lateral dimension of the TOT layers. The edge specific 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 SSAedge value ranges from 5 m2·g−1 to 30 m2·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 specific adsorption sites at frayed-edge sites on illite.

For metallic cations (e.g., Ni2+), 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 amphoteric sites may be modeled with surface complexation models, for which a generic reaction stoichiometry is written as:


>SOH+Meizi>B Meizi-1+H+,

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 field 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 coefficients 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 αci is responsible, in a large part, for the observed discrepancy between KD 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 affinities 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 significant fraction of the adsorption sites. The following simple example illustrates this problem. Na+ / Ca2+ cation exchange is considered to be the only reaction taking place at the clay mineral surface:


2NaX+Ca2+CaX2+2Na+log KNa/Ca=0.5.

The concentration of Ca2+ on the exchanger is obtained from:



where CCEC is the cation exchange capacity of the clay mineral in mol·L−1pore water. Equation (31) can be combined with the KD equation and yields:



Under the experimental conditions considered here, the Na+ concentration is constant and thus γCa2+ /[Na+]2 γNa+ is constant. If the concentration of Ca2+ on the exchanger is negligible as compared to the concentration of Na+, then [NaX] = CCEC and Equation (32) transforms into:



where qCEC is the CEC expressed in mol·kg−1clay. All the terms in Equation (33) are constant and Ca2+ adsorption is linear. As soon as [NaX] deviates from the CCEC value, however, Ca2+ 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 and Appelo 1999, 2013) in combination with i) a KD model or ii) a cation exchange model. The system modeled is similar to the one depicted in Figure 1. The Ca2+ 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 Ca2+ pore diffusion coefficient was set at 2 × 10−10 m2·s−1. According to Equation (33), and considering a qCEC value of 1 mol·kg−1clay and a NaCl background concentration of 0.1 mol·L−1, the KD value for Ca2+ at trace concentration is KD = 100 L·kg−1. While the KD model and the cation exchange models yielded identical results for the case with a Ca2+ concentration of 10−7 mol·L−1 in the inlet reservoir, the diffusion breakthrough curves were very different for the case with a Ca2+ 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 KD value to KD = 81 L·kg−1. Note that the KD value of Ca2+ at a concentration of 5 × 10−3 mol·L−1 is 42 L·kg−1 according to Equation (32) and thus a constant-KD 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 findings 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 define 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 defined 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 flux 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 flux of ions i in the bulk water bJi can be written as:



where bui is the mobility in bulk water (m2·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:



Combining Equation (36) with Equation (35), we obtain the Nernst–Planck equation:



with bDi = buiRT / |zi|F, the diffusion coefficient of species i in the bulk water (Steefel et al. 2014).

In the absence of an external electric field, there is no electrical current and so:



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:



Consequently, it is possible to express the Nernst-Planck equation with known parameters only, i.e., concentrations and activity coefficients (Boudreau et al. 2004).



Diffusive flux in the diffuse layer

Flux equation

According to Equation (17), it is possible to write:



The profiles of concentrations obtained from the Poisson–Boltzmann equation are representative of equilibrium between the solution in the diffuse layer and the solution at an infinite distance from the diffuse layer, i.e., the bulk water:



The combination of Equations (41) and (42) gives:



And so, by following the MGC model we implicitly assume:





The flux 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 field and thus electrical current in the diffuse layer along the axis x, we obtain:



with DDLj=uDLjRT|zj|F, the diffusion coefficient of species i in the diffuse layer water along the x-axis. The fluxes of species i in the diffuse layer along the x-axis become:



It can be noted that Equations (40) and (49) have the same form (Appelo and Wersin 2007; Appelo et al. 2010):



with the ‘DL enrichment factor’ Ai=e-ziFψDL(x,y)RT for the flux in the diffuse layer (DL) and Ai = 1 for the flux in the bulk water (b).

According to Equation (50), the diffusive flux can be split in three contributions for both diffuse layer and bulk waters:

  • the concentration gradient: -Db/DLiAibcix.

  • the activity coefficient gradient: -Db/DLiAicbilnγbix.

  • the diffusion potential: ziDb/DLicbiAijzjDb/DLjAj[cbjx+cbjlnγbjx]jzjD2b/DLjcbjAj.

Calculation 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 simplified model can be considered where only average ion concentrations in the diffuse layer, cDLi¯, are calculated:



where LDL 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 cDLi¯.

The average ion concentrations in the EDL are related to the surface charge, σ (in C·m−2), according to:



Note that hDL, the considered width of the diffuse layer, can be different from LDL, 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 flux 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:



and the diffusive flux 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, cexchi=EicCECzi, and so:



The diffusion potential term in the interlayer, ψexchdiffx, is obtained by considering the zero charge flux condition:





Approximations for Nernst-Planck equation for bulk and EDL water

Ion activity coefficient gradients

Ions activity coefficients are often calculated according to an extended Debye-Huckel model:



with I the ionic strength:



It follows:



In principle, this expression can be substituted into Equation (50), but in practice the consideration of the gradient of ionic strength can be difficult 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) reduces to:



Alternatively, by considering:



it is possible to further simplify Equation (50) into:



Equation (63) is for example the flux equation embedded in PHREEQC. It must however be remembered that Equation (62) is a simplification which neglects the cross-coupling terms since lnγbix=jlnγbicbjcbjx. The simplification made in Equation (62) leads to calculation results strictly equivalent to Equation (50) as long as it is evaluated numerically and as long as the concentration gradient cbix does not go to zero. If cbix=0, it is not possible to evaluate numerically lnγbilncbi, the denominator for which is zero when discretized between two numerical cells and one must assume that [1+lnγbilncbi]cbix=0. In this very special case Equation (63) is equivalent to Equation (61) and the activity gradient is neglected, even though the term cbilnγbix 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 coefficient 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. Their activity coefficients are the same and thus if follows that ctracerln γx=cbackground1000ln γx. However applying Equation (64), it follows that ctracerln γtracerx=cbackground106ln γbackgroundx because ctracerctracerx=cbackground106cbackgroundx. 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 coefficients

In PHREEQC, the diffusion coefficients in the diffuse layer are set proportional to the diffusion coefficients in the bulk water. This simplifying assumption is however not strictly necessary, and diffusion coefficients in the diffuse layer could be set to values which are not correlated with the diffusion coefficients values in the bulk water, as done in CrunchFlowMC. The latter case would imply that the tortuosity, and thus the effective diffusion coefficient, need not be the same in the bulk and EDL porosity. One can however justify the PHREEQC simplification by noting that the decrease of the diffusion coefficient 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 flux 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 flux 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 22Na and 36Cl 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 coefficients are set to values ten times lower in the porous medium than in pure water, i.e., Dp,Na+ = 1.33 × 10−10 m2·s−1 and Dp,Cl = 2.03 × 10−10 m2·s−1 in both bulk water and diffuse layer porosity. For simplicity, the Davies activity coefficient model is used (ADH = 0.5, Bai0 = 1, bi = −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 molc·kg−1, that is to say, values representative of an illite-rich clay material. The specific surface area is set to SSA = 100 m2·g−1, a value in good agreement with the properties of illite particles. The chosen thickness of the diffuse layer, hDL, is assumed = 3 × 10−9 m. The total porosity is φtot=1-ρdryρillite=1-12.8=0.64, the diffuse layer porosity is φDL = hDL ×SSA×ρdry = 0.3, and the volumetric charge compensated in the diffuse layer is q=-CEC×ρdryθDL=-0.33mol·L-1. Based on these parameters, it is possible to calculate the mean potential in the diffuse layer using Equation (52) as ψM = −0.033 V.

Example 1: Constant ionic strength

In the first 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 coefficient 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 jzjD2b/DLjcbjAj is high as compared to jzjDb/DLjAj[cbjx+cbjlnγbjx]. While the 36Cl flux is higher than the 22Na+ flux in the bulk water because of the higher of diffusion coefficient 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 22Na+ and the depletion of 36Cl 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 fluxes 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 flux in the bulk water is increased as a result of this term, while the Cl diffusive flux 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 coefficient term is minor as compared to the two other terms, but it represents a negative contribution of up to 12% of the total 36Cl flux in the bulk water (up to 8% for 22Na+, Fig. 11). Moreover the contribution of this term is not constant along the x-axis and its negative contribution to the total fluxes 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 22Na+ and 36Cl diffusive fluxes remain equal in the diffuse layer, the 22Na+ diffusive flux 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 coefficient term. For 36Cl, the concentration gradient and the diffusion potential terms have almost the same positive contribution to the diffusive flux. At x = 1 cm, Na+ and Cl diffusive fluxes 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 36Cl 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 flux intensity. Similar to the bulk water, the calculated 22Na+ and 36Cl fluxes 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 flux Equation (50). If this term is dropped as in Equation (61), then the fluxes increase by up to 4% for the conditions considered here (Fig. 11).

Example 3: Gradient in ionic strength and no tracer gradient

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 22Na+ and 36Cl 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 coefficient gradient and diffusion potential terms are the only remaining contributors to the total diffusive flux. In those conditions, the diffusive flux is positive for 22Na+ in the bulk water and negative for 36Cl while the reverse is true in the diffuse layer (Fig. 12). In these conditions, the diffusion potential term dominates over the activity coefficient gradient term (Fig. 13).

Links to experimental diffusion results

In Example 1 above, the fluxes 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 coefficients are the same in the EDL and bulk water. Under the same conditions, the diffusive flux 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 fluxes at steady state are in the order JCations > JHTO > JAnions (Glaus et al. 2010; Tachi and Yotsuji 2014).

Example 2 points out the extent of the approximation made in the calculation of the diffusive flux if the activity coefficient gradient is dropped in the solution of the Nernst–Planck equation. The relative error made in the calculation of the diffusive flux is less in the EDL water than in the bulk water, thereby justifying this simplification if the contribution of the EDL to the total flux dominates the flux 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 flux is two orders of magnitudes higher in the EDL than in the bulk water (Fig. 12). For anions, the direction of total diffusive flux 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 flux for the tracers in this case (Fig. 13). Again, this observation helps to justify the simplifying assumption of ignoring the activity coefficient gradient term in the Nernst–Planck equation. Up-gradient diffusion of 22Na has been recently reported in an experimental system similar to that described in Example 3 (Glaus et al. 2013).


The analysis presented above of the different contributions of the terms in the flux equation highlights the feasibility of explaining qualitatively experimental observations of water and ion fluxes with the Nernst Planck model developed from Equations (34) to (64). Quantitative simulations of those fluxes including transient conditions can be achieved only by implementing the flux 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 Ci 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 decay), then Ci = φci, where φ is the total porosity and ci is the total concentration in the porosity:



In the presence of a diffuse layer and of an interlayer space, the term φci can be split into three contributions (bulk – b, diffuse layer – DL, and interlayer, or exchangeable, – exch):



The same approach (and nomenclature) can be used for the flux term Ji. 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 flux terms includes porosity and tortuosity contributions within their effective diffusion coefficient terms:



If the porous media is not homogeneous ( φb/DL/exchx0), then it is difficult to define 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 flux between the bulk water and the diffuse layer water should be considered. This problem is exemplified in the following section by considering simple illustrative examples.

Summation of bulk and diffuse layer diffusive fluxes over an interface

Normalized flux term

We define here the terms, i.e., fluxes normalized to a porosity value of one:



Simple case

We will first 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 bJi* and



Complex case 1: Same total porosity, different diffuse layer porosity

In this first complex case (Fig. 15), the calculation of the diffuse layer and bulk water contributions to the overall diffusion flux is not straightforward. A simple method can be the use of an arithmetic mean, which is also consistent with Equation (71) if bφ1 = bφ2 and DLφ1 = DLφ2:



Complex case 2: Same total porosity, one zone without diffuse layer porosity

In the second complex case exemplified in Figure 16, the calculation of the diffuse layer and bulk water contributions to the overall diffusion flux 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 flux summation becomes:



Consequently, Equation (72) lacks a condition of continuity. This problem can be solved by considering the more complex equation:


Ji*=φb2φJbi*+φDL1φJDLi*+φ-φDL1-φb2φJb to DLi*,

where b to DLJi* is a flux 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) fulfills 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 fluxes in the diffuse layer and in the bulk water can be obtained theoretically from the consideration of the Nernst–Planck equation coupled to the modified Gouy–Chapman model (or its simplified form, the mean potential model). In contrast, the summation of the flux at complex interfaces between two numerical cells does not follow rules that are dictated by any firmly 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 filter (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 flux at interface between two numerical grid cells

The numerical differentiation of the flux poses another challenge. For a grid cell n, at position xn, the numerical equivalent of Equation (65) is given by:



where superscripts new and old refer to two consecutive time steps, and where xn+1/2 and xn-1/2 are the positions at the interface between cells n and n + 1, and n and n − 1, respectively.

Ji,xn+1/2 and Ji,xn–1/2 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 coefficient). For the bulk and diffuse layer contributions of, Ji,xn+1/2 the equations are represented numerically by finite differences as:



While the term cbi,xn+1-cbi,xnxn+1-xn is a well-defined quantity, the terms bcj,xn+1/2 (the concentration in the bulk water at the interface between cells n and n+1) and Aj,xn+1/2 (the 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 Aj,xn+1/2 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 exemplified in Figure 18. In that case an arithmetic or harmonic mean may be inaccurate for estimating the Ai value at the interface. The same kind of problem may occur for a gradient of surface charge between two grid cells.


Code limitations

While a range of reactive transport codes handle the Nernst–Planck equations for diffusive fluxes 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. 2008, 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 specific surface area of 750 m2·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 KNa = 0.7 (in agreement with the MD results from Tournassat et al. 2009) and log KCs = 2.5:


>Surf-+Na+>SurfNa KNa=c>SurfNaaNa+·c>Surf-exp(FΨ0RT),


>Surf-+Cs+>SurfCs KNa=c>SurfCsaCs+·c>Surf-exp(FΨ0RT).

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 fitted for each species and are reported in Table 2. Results are plotted in Figure 19. The tortuosity values follow the order τI < τHTO < τ22Na+ < τ137Cs+, 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 flux calculation make it possible to derive a physically feasible value (< 1) for τ137Cs+, 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 KD value of 430 L·kg−1 at a ionic strength of 0.1, i.e., a value in agreement with KD 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 flux 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 steady-state conditions and/or the estimation of apparent diffusion coefficients 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. 2008, 2010), but still remains the exception rather than the rule for the interpretation of diffusion data.

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 filters was explicitly taken into account by assuming a tortuosity of 0.25 and a porosity of 0.32 for the filters. This corresponds to an effective diffusion coefficient of ~10−10 m2·s−1 for Na+ in the filter, 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 fluxes 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 coefficient is about (1.5–1.7)·10−11 m2·s−1 in similar conditions (Glaus et al. 2010, 2013). The modeling results are shown on Figure 20 and agree almost perfectly with the experimental data.

Glaus et al. (2013) have also modeled their data using two types of model. The first 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 fit 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 fits 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 flat surfaces). The exchange sites in the vermiculite were saturated with Ca2+. The interlayer width in the vermiculite corresponded to a bi-layer hydrate. Following contact with a NaCl solution, a release of Ca2+ was observed in solution. This experiment showed unambiguously that interlayer diffusion exists, and it made it possible also to quantify the diffusion coefficient for Ca2+ 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 Ca2+ interlayer diffusion coefficient 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.

The data from Tertre et al. (2015) were reinterpreted using the interlayer diffusion option of PHREEQC (Fig. 21b). It was possible to reproduce the data with the same level of quality as the Brownian dynamics calculations by considering an interlayer Ca2+ diffusion coefficient value of 0.8 × 10−11 m2·s−1, together with an interlayer Na+ diffusion coefficient of 4 × 10−11 m2·s−1, 1 × 10−11 m2·s−1, 0.1 × 10−11 m2·s−1, and 0.1 × 10−11 m2·s−1, for the experiments at NaCl concentration of 1, 0.1, 0.05 and 0.003 mol·L−1 respectively. The decrease of the Na+ interlayer diffusion coefficient with NaCl concentration is in agreement with MD results from Tertre et al. (2015), which show a decrease of its value with an increase of the Ca2+/Na+ occupancy ratio in the interlayer (Fig. 21c). At 1 mol·L−1 NaCl, the PHREEQC and MD results fully agree. At 0.05 and 0.003 mol·L−1 NaCl, the interlayer diffusion coefficient fitted with PHREEQC correspond to the lowest value obtained with MD (within the range of the error bands). Consequently, the apparent decrease in Ca2+ interlayer diffusion coefficient can also be interpreted as a result of the decrease of the Na+ interlayer diffusion coefficient that arises from the coupling between the diffusion of these two species through the diffusion potential term in Equation (57).


Diffusion processes through clay materials is the result of a complex interplay of transport and non-linear adsorption processes under the influence of electrostatic fields. In this context, the classical Fickian diffusion model applied in the framework of a single pore diffusion model with linear adsorption processes (KD model) appears to be inappropriate for describing diffusion data without making less than satisfying modeling assumptions, such as i) different definitions 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 coefficients (or tortuosity values).

It is clear that RT codes have been capable of solving the problem of the transfer of KD values from batch experiments to compact porous media systems for some time. Much of the difficulty 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 fits 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,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 finer details (e.g., adsorption processes, activity coefficients)

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.


This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and by the European Atomic Energy Community Seventh Framework Programme [FP7 – Fission – 2009] under grant agreement n∘624 249624 Collaborative Project Catclay. C. Tournassat acknowledges funding from L’Institut Carnot BRGM for his visit to the Lawrence Berkeley National Laboratory.