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 the original work is properly cited.


The aqueous chemistry of water films confined between clay mineral surfaces remains an important unknown in predictions of radioelement migration from radioactive waste repositories. This issue is particularly important in the case of long-lived anionic radioisotopes (129I,99TcO4,36Cl) which interact with clay minerals primarily by anion exclusion. For example, models of ion migration in clayey media do not agree as to whether anions are completely or partially excluded from clay interlayer nanopores. In the present study, this key issue was addressed for Cl using MD simulations for a range of nanopore widths (6 to 15 Å) overlapping the range of average pore widths that exists in engineered clay barriers. The MD simulation results were compared with the predictions of a thermodynamic model (Donnan Equilibrium model) and two pore-scale models based on the Poisson-Boltzmann equation under the assumption that interlayer water behaves as bulk liquid water. The simulations confirmed that anion exclusion from clay interlayers is greater than predicted by the pore-scale models, particularly at the smallest pore size examined. This greater anion exclusion stems from Cl being more weakly solvated in nano-confined water than it is in bulk liquid water. Anion exclusion predictions based on the Poisson-Boltzmann equation were consistent with the MD simulation results, however, if the predictions included an ion closest approach distance to the clay mineral surface on the order of 2.0 ± 0.8 Å. These findings suggest that clay interlayers approach a state of complete anion exclusion (hence, ideal semi-permeable membrane properties) at a pore width of 4.2 ± 1.5 Å.


How does water confinement in nanopores impact anion exclusion? The answer to this question has important implications for the migration of salts and groundwater through clay-rich rocks (Fritz, 1986; Leroy et al., 2006), the design of engineered nanofluidic systems (Schoch et al., 2008; Bocquet and Charlaix, 2010), the performance of desalination membranes (Richards et al., 2012), and the transport of ions through channels in cell membranes (Nikaido and Rosenberg, 1981; Nikaido, 2003). It also is of importance to high-level radioactive waste (HLRW) management, where the performance of geologic repositories is determined, in many cases, by the migration of anionic, long-lived fission products (129I,99TcO4,36Cl) through nanoporous clay barriers (Tournassat et al., 2007; Altmann, 2008; Claret et al., 2010; Glaus et al., 2011; Altmann et al., 2012).

At present, extant models of anion adsorption and diffusion in nanoporous clayey media postulate either that anions are completely excluded from clay interlayer nanopores (Sposito, 1992; Wersin et al., 2004; Bourg et al., 2008; Churakov and Gimmi, 2011; Tournassat and Appelo, 2011) or that anions are partially excluded from these nanopores in a manner that can be predicted accurately by assuming nanopore water is the same as bulk liquid water, albeit with a lower electrostatic potential imposed by negatively charged clay surfaces (Leroy et al., 2006; Birgersson and Karnland, 2009; Glaus et al., 2013; Tachi et al., 2014). Efforts to test these two competing hypotheses using anion exclusion measurements have been inconclusive because the results of such measurements depend not only on the mechanistic details of anion exclusion in individual nanopores, but also on the pore-size distribution in saturated, compacted clayey media (Tournassat and Appelo, 2011; Tournassat et al., 2016; Tinnacher et al., 2016), which is not yet fully understood (Holmboe et al., 2012).

Molecular dynamics (MD) simulations offer a way out of this conundrum by directly examining the details of anion exclusion in individual clay interlayer nanopores. However, three recent studies that quantified anion exclusion in clay interlayers by MD simulation have yielded contradictory results: one concluded that interlayer water can be treated as having the same aqueous chemistry as bulk liquid water (Hedström and Karnland, 2012), whereas two others reached the opposite conclusion, that anion exclusion in clay interlayers is much greater than would be expected based on the aqueous chemistry of bulk liquid water (Rotenberg et al., 2007b; Hsiao and Hedström, 2015). Consequently, the correct way in which to account for anion exclusion in macroscopic models of ion migration in clay formations and engineered clay barriers remains uncertain (Muurinen et al., 2007; Birgersson and Karnland, 2009; Appelo et al., 2010; Tournassat and Appelo, 2011). The present study aimed to resolve the differences among previous MD simulations. The MD analysis was carried out using montmorillonite as a prototypical smectite clay mineral. This clay mineral contributes predominantly to the specific surface area and surface charge characteristics of most engineered clay barriers and many soil and sedimentary formations. It is also the most widely examined clay mineral in studies of anion diffusion through engineered clay barriers (Bourg and Tournassat, 2015). While the present study focuses on clay interlayers, the results should be relevant to ion and salt migration in other charged nanopores with pore widths below 2 nm.


Donnan equilibrium: the thermodynamic view

The most widely applied thermodynamic model of anion exclusion in clayey media, the Donnan equilibrium (DE) model, relies on two fundamental assumptions (Babcock, 1960): (1) that chemical equilibrium exists between two distinct aqueous phases, a hydrated clay phase and a bulk aqueous solution phase, and (2) that Henry’s Law applies at infinite dilution to all of the ions in both phases. The first assumption implies that anions must always be present in the hydrated clay phase at equilibrium, whereas the second one implies that ions in the hydrated clay phase are completely dissociated from the clay portion at infinite dilution. As shown in the Supplemental Materials section (deposited with the Editor in Chief and available at, the two assumptions of the DE model yield the following expression for the concentration ratio of the anions in the case of an aqueous NaCl solution in equilibrium with a hydrated clay (Babcock, 1960, 1963):  
mi and γi are the molality (mol kg−1) and activity coefficient of ion i,  
is the ratio of anion charge concentration in aqueous solution to that of the clay anionic charge (both in moles of charge per kg of water), and the superscripts clay and soln refer to the hydrated clay phase and the bulk aqueous solution phase, respectively. Thus, the anion concentration ratio of the two phases is determined solely by the anion-to-clay charge ratio y and the electrolyte activity coefficient ratio θ. Note that θ = 1 if the electrolyte in each phase has the same mean ionic activity coefficient. The left side of Equation 1 may be very small, indicating anion exclusion, or may be larger than 1, indicating positive anion adsorption in the hydrated clay phase, depending on the values of y and θ. (In general, θ will depend on y.) Consequently, the DE model predicts quite generally that, not only are anions always present in clay interlayers, anions can even swarm there under certain circumstances (Babcock, 1960). This prediction, of course, is a direct consequence of the two model assumptions and thus cannot be considered as a theoretical proof of the presence of anions in clay interlayer nanopores in equilibrium with bulk aqueous solution. Moreover, Equation 1 cannot be applied in a predictive manner, because no accurate model of γiclay exists. Studies of anion exclusion and semi-permeable membrane properties in clayey media have, therefore, routinely used Equation 1 along with the simplifying approximation that γiclay=1 or θ = 1 for all ions (Kharaka and Berry, 1973; Fritz, 1986; Keijzer et al., 1999; Sherwood and Craster, 2000; Revil et al., 2011), an approximation that forms the basis of the well known Teorell-Meyer-Sievers model of ion migration in charged semi-permeable membranes (Wang et al., 1995; Schoch et al., 2008). Brownian dynamics calculations, however, strongly suggest that θ < 1 for a three-layer hydrate of montmorillonite equilibrated with a NaCl solution, whereas θ ≈ 1 in larger interlayer pores (Jardat et al., 2009).

Gouy-Chapman model: the molecular scale view

An alternative to thermodynamic modeling of anion exclusion is to describe, at the molecular scale, the structure of the electrical double layer (EDL), i.e. the distribution of aqueous ions in the vicinity of a charged surface. The simplest approach derives from the Gouy-Chapman model, which is based on six assumptions: first, that the Poisson equation, which relates the electrostatic potential in a swarm of ions to the net charge density in the swarm, remains physically meaningful on molecular length scales; second, that the surface from which x is measured is a uniform, infinite plane of charge characterized by a surface charge density σ (in C m−2); third, that the charged species in the nearby aqueous solution are fully dissociated from the planar surface lying at x = 0; fourth, that these ions interact among themselves and with the surface exclusively through the coulomb force; fifth, that water can be represented as a uniform continuum characterized solely by the dielectric constant ε; and, finally, that the potential of mean force W(x) is proportional to the mean electrostatic potential, ψEDL(x), at a distance x from the charged surface (Sposito, 1984, 2004). Under these assumptions, the Poisson equation can be transformed into the Poisson-Boltzmann (PB) equation of the Gouy-Chapman model:  
where F is Faraday’s constant (96 485 C mol−1), R is the ideal gas constant (8.314 J mol−1 K−1), T is absolute temperature, ε0 is the permittivity of vacuum (8.845×10−12 C V−1m−1), ε is the dielectric constant of water, and ci is the molality of ion i (mol dm−3). The electrostatic potential ψEDL(x) has no thermodynamic significance. In particular, it cannot be substituted for ψclay in Equation 1 because it cannot be measured without relying on ad hoc non-thermodynamic concepts and models (Babcock, 1963; Sposito, 1984).
The Gouy-Chapman model is routinely used to describe anion exclusion and migration in clay-water systems and other charged nanoporous media (Smith et al., 2004; Schoch et al., 2008; Cheng and Hendry, 2014). Two specialized forms of the model are widely used. The first, known as the modified Gouy-Chapman (MGC) model, is based on the assumptions listed above along with the additional condition that ions can be modeled as hard spheres and, therefore, have a distance of closest approach a to the clay surface (Sposito, 1992; Kim and Netz, 2006). The MGC model is consistent with experimental data on anion exclusion in dilute clay suspensions if a ≈ 1.8 Å (Sposito, 1992). The second specialized form, known as the mean electrostatic potential (MEP) model, is based on the same assumptions as the Gouy-Chapman model along with the simplifying approximation that ψEDL(x) does not vary with distance from the charged surface within the EDL. This MEP approximation allows Equation 2 to be solved analytically for complex electrolytes or for systems where the diffuse layers on different charged surfaces overlap. As shown in the Supplemental Materials section, the MEP model yields the expression:  
where cClEDL¯ is the mean chloride concentration in the EDL,  
σ is again the surface charge density (C m−2), and hEDL is the thickness of the region occupied by the EDL (in m). A distance of closest approach of ions to the charged surface can be incorporated into the MEP model by substituting (hEDLa) for hEDL, thus defining an “ion-accessible EDL thickness” in Equation 3.


Equations 1 and 3 obviously have very similar forms. Because of this similarity, the terms “Donnan model,” or “Donnan volume,” are sometimes erroneously used to describe the MEP model despite its molecular-scale assumptions (Appelo and Wersin, 2007; Schoch et al., 2008; Birgersson and Karnland, 2009; Tertre et al., 2011; Tournassat et al., 2011; Glaus et al., 2013). The DE model, however, is based on chemical thermodynamics and, therefore, is fundamentally different from the MGC and MEP models, or any other model involving molecular-scale concepts.

Whether a Donnan equilibrium model or a model based on the Poisson-Boltzmann equation is best suited for predicting anion concentrations in clay mineral interlayers cannot be readily answered by comparing model predictions with macroscopic experimental results. Macroscopic data can be fitted to the DE model by adjusting the clay-phase activity coefficients in θ. Since neither the individual activity coefficient values nor the Donnan potential is accessible experimentally, θ cannot be constrained independently. Similarly, macroscopic anion-exclusion data can be fitted with the MGC or MEP model by adjusting the pore-size distribution or the distance of closest approach of ions to the clay surface (Muurinen et al., 2007; Tournassat and Appelo, 2011; Tournassat et al., 2016). Independent information must, therefore, be sought that can be used to test these modeling approaches. In the following, MD simulations are shown to be insightful in this respect.


Simulations of clay-water systems with three interlayer nanopores in contact with a larger mesopore were carried out with the program LAMMPS (Plimpton, 1995). Water molecules were assumed to be rigid according to the SHAKE algorithm (Ryckaert et al., 1977). The clay mineral structure was also assumed to be rigid, as in previous simulations of this type (Rotenberg et al., 2007b; Hedström and Karnland, 2012; Hsiao and Hedström, 2015), but clay H atoms were allowed to move. Recent MD simulations have shown that the equilibrium distributions of interlayer water and ions (but not their dynamics) are insensitive to the motions of clay atoms (Holmboe and Bourg, 2014). Interatomic interactions were described with the SPC/E model of liquid water (Berendsen et al., 1987), the CLAYFF model of mineral-water interactions (Cygan et al., 2004), and the ion-water interaction parameters of Smith and Dang (1994). Van der Waals interactions between unlike atomic species were derived using the Lorentz-Berthelot combining rules. The CLAYFF model includes Van der Waals and Coulomb interaction parameters for atoms in the bulk clay mineral structure and for certain functional groups on the clay edge surfaces (>Si–OH, >Al2–OH), but not for other clay mineral edge groups (>AlOH, >Al–O–Si<). For the present study, CLAYFF-like interatomic parameters were derived for edge >AlOH and >Al–O–Si< groups as described in the Supplemental Materials section.

The cis-vacant structure of montmorillonite was taken from Tsipursky and Drits (1984). Individual clay mineral layers (hereafter referred to as TOT layers) were truncated at the (110) and ( 11¯0) structural planes in order to create edge surfaces, which were then healed by adding −OH or −H groups to under-coordinated edge atoms in order to obtain edge surfaces with zero net proton surface charge (Churakov, 2006). The resulting edge surface structures, which correspond to the AC chains of White and Zelazny (1988), are thought to be the most stable or one of the most stable smectite edge structures (White and Zelazny, 1988; Bickmore et al., 2003; Churakov, 2006; Liu et al., 2014; Newton and Sposito, 2015; Newton et al., 2016). Each TOT layer had a thickness of ~10 Å (ambiguities associated with the definition of the TOT layer thickness are discussed below), a finite length of 45.7 Å in the direction normal to the edge surfaces, and a width of 63.4 Å (effectively an infinite width, because of the periodic boundary conditions of the simulation cell). Each simulation cell contained three TOT layers, three interlayer nanopores, and a larger mesopore (Figure 1). The cell had an initial length of 90 Å in the direction normal to the clay mineral edges (before pressure equilibration), a width of 63.4 Å, and a height of 46.8 Å, 56.7 Å, or 76.5 Å in the direction normal to the clay basal surfaces for simulations of the two-, three-, and five-layer hydrates (referred to hereafter as 2WL, 3WL, and 5WL systems, respectively). Isomorphic substitutions of Mg for Al were randomly distributed on cis octahedral sites. In order to minimize uncertainties associated with the parameterization of edge surface sites, isomorphic substitutions were not allowed to occur near the edge surfaces. The mean number of isomorphic substitutions (108 substitutions randomly distributed between the three TOT layers) was selected in order to obtain a mean basal surface charge density near −0.1 C m−2, a value typical of montmorillonite (Sposito, 1984).

Each simulation cell included 148 Na+ ions and 40 Cl ions as required to balance the negative surface charge of the TOT layers and establish a NaCl concentration in the mesopore of ~0.5 mol dm−3, ~0.4 mol dm−3, or ~0.3 mol dm−3 in the 2WL, 3WL, and 5WL systems, respectively. The NaCl concentration in the mesopore was chosen such that the MEP model calculations of Hedström and Karnland (2012) (based on Equation 2, with hEDL equal to one-half the width of a clay interlayer nanopore and with a closest approach distance of a = 0) would predict one Cl ion in each interlayer nanopore, on average, during the entire simulation. Previous calculations have indicated that, at the ionic strengths examined in the present study, the region where the electrostatic potential spills over between the nanopore and mesopore extends ~1 nm from the clay mineral edges into the interlayer nanopores and into the mesopore (Kraepiel et al., 1999; Jardat et al., 2009). These simulated systems should, therefore, be sufficiently large to contain interlayer and mesopore regions that are not influenced by proximity to the clay edge surfaces.

As in the simulations of Hedström and Karnland (2012), a single Cl ion (hereafter denoted Clint) was initially placed in the center of each interlayer nanopore. The three Clint ions were held fixed during equilibration, then allowed to move freely during the simulation. Production runs were carried out in the NVT ensemble (constant composition, volume, and temperature) at T = 350 K with a 1 fs time step. Simulation runs were preceded by 0.5 ns of equilibration in the NVT ensemble, 0.5 ns of equilibration in an NPyT ensemble (constant composition, stress in the direction normal to the edge surfaces, and temperature) at Py = 1 bar and T = 350 K, and another 2 ns of equilibration in the NVT ensemble. Electrostatic and dispersion interactions beyond 15 Å were computed by the Ewald sum method during part of the simulation (22 ns, 17 ns, and 15 ns for the 2WL, 3WL, and 5WL systems, respectively), then by the similarly accurate particle-particle particle-mesh (PPPM) method (Eastwood et al., 1980; Isele-Holder et al., 2012) during the remaining part of a simulation (up to 53 ns, 40 ns, and 38 ns for the 2WL, 3WL, and 5WL systems, respectively).

The contributions of various electrostatic interactions to the potential energy of all Cl and Na+ ions were calculated using LAMMPS. The contributions of electrostatic interactions with water molecules, clay mineral atoms, and other ions (Na+, Cl) to the total per-atom potential energy of Na+ and Cl ions were calculated by canceling the charges of all water atoms, then of all water and montmorillonite atoms during different reexaminations of a trajectory. For the last calculation, where only Na+ and Cl mutual electrostatic interactions were calculated, the system was not neutral (an excess of Na+ was present due to the TOT layer surface charge) and the computation of a relevant total electrostatic energy (including pairwise and k-space energy contributions) was not possible. Consequently, only pairwise energies were computed with a cutoff at 25 Å.


Cl access to the interlayer

All Clint atoms left the 2WL and 3WL interlayer space within three nanoseconds of simulation. During the same period of time, Cl ions from the mesopore entered the 3WL nanopores but did not enter the 2WL nanopores. In the 5WL system, not all of Clint atoms left the interlayer space within the same period of time. According to the concentration maps that were computed after 3 ns of simulation, Cl has access to 3WL and 5WL interlayers, but not to 2WL interlayers (Figure 2). Note that in the 2WL system, the reported Cl concentration in the interlayer space is truly zero, not simply a low value that is unresolved by the color concentration scale.

Na, Cl, and water distribution profiles in the interlayer nanopores

Concentration profiles from molecular dynamics

Total Cl and Na+ concentration profiles and water density profiles (Figure 3) were computed with a spatial resolution of 0.5 Å on the z-axis in the y [65–85 Å] interval, i.e. in the middle of the interlayer nanopores. The calculated water-density profiles correspond to the concentration of O atoms in water molecules multiplied by the mass of one mole of water. The Na+ ions were concentrated in a single peak at the center of the 2WL interlayer space, while two main peaks were present in the 3WL and 5WL interlayer spaces (Figure 3) in agreement with previous MD and Monte Carlo (MC) simulations of Na-montmorillonite with structural charge arising from isomorphic substitutions in the octahedral sheet (Marry et al., 2002; Tambach et al., 2004; Holmboe and Bourg, 2014). The Na density maxima were located at a distance of ~4.5 Å from the siloxane surface O atoms of the TOT layer and separated from the TOT layer by a hydration layer characterized by a water O density peak. Peak positions were consistent with previous studies of clay mineral interlayer nanopores and larger pores (Marry et al., 2002; Tournassat et al., 2009; Bourg and Sposito, 2011). The height of the Na density peaks increased with the number of isomorphic substitutions in the nearby TOT layer, as expected from charge balance considerations. The height of the main water density peak was correlated, but was not strictly proportional, to that of the Na peak, indicating that the structure of the interfacial water layer was significantly influenced by both the clay surface and the adsorbed cations. Few Na+ ions were engaged in inner sphere surface complexes. In the case of Cl, density profiles showed significant anion exclusion even at the mid-plane of the nanopores (Figure 3). The distance of closest approach of Cl to the plane of O atoms on the basal surface was ~ 3.5 Å. This value coincides with the distance of closest approach of Na+ ions that were not bound in innersphere surface complexes.

Na+ and Cl mean interlayer concentrations according to MD calculations

At the scale of an individual interlayer nanopore, average concentrations calculated from the MD concentration profiles depend on an arbitrary choice of the boundary between the pore space and the solid phase. The DE model calculations presented in a previous section are based on a molality convention (mol per kilogram of water) while the MGC and MEP model calculations are based on a molarity convention (mol dm−3) with the underlying assumption of a water density of 1 kg dm−3 (note that the effect of temperature on water density was neglected). Consequently, the location of the clay mineral-water boundary was identified based on the choice that yields an average interlayer water density of 1 kg dm−3. The resulting effective TOT layer thickness equals 9.78 ± 0.45 Å for all three systems investigated in the present study (Table 1). Accordingly, average interlayer ion densities were calculated based on interlayer nanopore widths of 5.7, 9.1, and 15.6 Å in the 2, 3, and 5WL systems, respectively (Table 1).

Concentration profiles predicted with the MGC and MEP models

For comparison with the mean interlayer concentrations predicted by MD simulation, interlayer concentrations of Na+ and Cl were calculated using the MGC and MEP models. The MGC model predictions were obtained by solving the PB equation (Equation 2) numerically using the values of σ and cisoln from the MD simulation as input parameters (Table 2). The distance of closest approach of ions to the solid surface was identified based on the MD concentration profiles, using either the location of the ‘foot’ of the Na+ density profiles or the location of the Na+ density peak (Figure 3 and Table 1). In the case of the 2WL system, the second option was not available because only one Na+ density peak was present. The MEP model calculations assumed that the EDL filled the entire interlayer pore space and involved the same assumptions as the MGC model regarding the distance of closest approach of ions to the clay mineral surface (Table 2). Predicted Na+ and Cl density profiles were converted to average concentrations per mass of interlayer water in order to ensure consistency between different results.

Maps of per-atom electrostatic energy

The per-atom electrostatic energy corresponds to the coulomb energy attributed to each individual Cl (Na+) ion in the system. For the 2WL system (Figure 4), the computation of the per-atom energy for Cl in the interlayer space was not possible because chloride did not enter the interlayer space during the simulation. On the other hand, the total per-atom electrostatic energy was the same in the mesopore as in the interlayer in the 3WL and 5WL systems (Figures 5 and 6). The most remarkable outcome of this calculation is that, in the 3WL and 5WL systems, the electrostatic energy arising from Cl solvation was much higher (less favorable interaction) in interlayer nanopores than in the mesopore. The opposite effect was not observed in the case of Na+; electrostatic interaction of Na+ ions with water molecules was approximately equally favorable in the interlayer nanopores as in the mesopores. A similar asymmetry in Na+ and Cl solvation energies has been observed in model polar nanopores and ascribed to hindered orientations of confined water molecules hydrating Cl in nanopores (Crozier et al., 2001). As expected, interactions with other ions in the interlayers were attractive in the case of Cl and repulsive in the case of Na+ (lower row in Figures 4, 5, and 6).


2WL hydrate

In the case of the 2WL hydrate, no Cl ion entered the interlayer space during the course of the simulation, in agreement with the modeling results of Rotenberg et al. (2007b), but in disagreement with those of Hsiao and Hedström (2015). The Cl concentrations expected from the energy barriers predicted by Rotenberg et al. (2007b) or Hsiao and Hedström (2015) for bulk NaCl concentrations similar to those obtained in the present study are given in Table 3, together with the average number of Cl ions expected in the interlayer space during the course of an entire simulation. The simulations were carried out at the same temperature (350 K) as the simulations of Hsiao and Hedström (2015) and with similar simulation times (50 ns vs. 100–200 ns) and volumes (27×104 Å3vs. 15×104 Å3), thus ensuring roughly equally reliable output statistics. The fact that Cl ions did not enter the interlayer space cannot, therefore, be attributed to a lack of convergence in the present simulation, as Hsiao and Hedström have postulated to explain the difference between their results and those of Rotenberg et al. (2007b).

The mapping of the per-atom electrostatic energy contributions evinced the major role played by solvation in inhibiting the entry of Cl into the interlayer space. This effect was particularly pronounced in the 2WL system, perhaps because a significant fraction of the water molecules associated with the water density peaks nearest to the clay surfaces (i.e. most of the interlayer water in the 2WL hydrate) were engaged in solvating Na+ outer-sphere surface complexes while also donating a hydrogen bond to the siloxane surface O atoms (Marry et al., 2008; Bourg and Sposito, 2011). In short, water located in the first monolayer on the clay surfaces had a very different structure from bulk liquid water (Low, 1979), and this difference became acute in the case of the 2WL hydrate. These findings are consistent with the hypothesis formulated by Hsiao and Hedstrom (2015) that differences in solvation energy play an important role in inhibiting the entry of Cl in the interlayer space.

Molecular dynamics simulations are well known to be sensitive to the choice of parameters used to describe interatomic interactions. Consequently, differences in the extent of anion exclusion predicted by Rotenberg et al. (2007b) and Hsiao and Hedström (2015) may reflect differences in the parameterization of interatomic interactions in these two studies. The present work used the same interatomic potential models as Hsiao and Hedstrom (2015), namely, the CLAYFF and SPC/E models. To the best of the present authors’ knowledge, the only major methodological difference between the present simulations and those of Hsiao and Hedstrom (2015) is that these latter authors treated the negative structural charge of the clay mineral particles as uniformly distributed over all octahedrally coordinated atoms in the clay mineral structure. In the present simulations, the negative charge arising from each Al → Mg isomorphic substitution was assigned to the site of substitution and its six surrounding O atoms in accordance with the original CLAYFF model. The treatment of clay mineral negative structural charge used in the present study was derived from quantum mechanical calculations of the electrostatic potential near clay mineral surfaces (Cygan et al., 2004), which predicts interlayer water and ion distributions and dynamics in excellent agreement with experimental results (Bourg and Sposito, 2010; Ferrage et al., 2011; Marry et al., 2011). The present results do not allow an accurate quantification of the height of the energy barrier opposing Cl entry into the interlayer space of the 2WL hydrate, but strongly suggest that this barrier is greater than that predicted by Hsiao and Hedstrom (2015).

3WL and 5WL hydrates

The numerical solution of the PB equation for the 3WL and 5WL hydrates was consistent with the average interlayer Na+ concentration and with the interlayer Na+ concentration profile predicted by the MD simulations (Figure 3). This close correspondence was particularly true if the distance of closest approach of ions to the clay mineral surface was identified with the foot of the Cl density profiles, which yields a = 2.0±0.8 Å. This predicted a value was in good agreement with the range of values regarded as a typical ion size in studies of ion adsorption (a =1.8 to 2 Å) (Sposito, 1992; Kim and Netz, 2006) and with the hydrated radii of Na+ and Cl (1.7 and 2.0 Å, respectively) (Jardat et al., 2009). The agreement between PB calculations and MD simulation predictions was somewhat worse in the case of the Cl concentration profiles than in the case of the Na+ profiles (Figure 3), perhaps reflecting the poorer statistics for interlayer Cl concentrations or the influence of short-range ion-ion interactions (and possibly ion-water interactions, as noted above) that are not accounted for in the PB equation. Nevertheless, reasonable quantitative agreement was found (Table 2). These calculations further showed that, for the systems examined, the average interlayer Na+ and Cl concentrations predicted with the MGC and MEP models were essentially identical, i.e. a full numerical resolution of the PB equation was not necessary to predict average ion concentrations in the nanopores. Instead, predicted interlayer anion concentrations were primarily sensitive to the distance of closest approach used in the MGC or MEP models, the best fit being obtained with a distance of closest approach intermediate between those obtained using the ‘closest approach’ and ‘Na peaks’ approximations. The similarity of the MGC and MEP models in the case of monovalent ions is consistent with recent predictions of the density profiles of Na+, Cl, and Br (but not Ca2+) in a 31.5 Å-wide clay interlayer nanopore (Tinnacher et al., 2016).

According to Hsiao and Hedstrom (2015), solving the Poisson-Boltzmann equation in order to model the behavior of ions in clayey media is not necessary, because the use of the Donnan equilibrium (DE) model is warranted for compacted clay systems. A strict application of Henry’s Law (and, hence, of the DE model) to clay systems would require that the counterions be dissociated fully from the mineral surface, such that mean ionic concentrations in the hydrated clay phase be unambiguously determined (Babcock, 1963). The present MD simulations were carried out at a relatively high salinity and, therefore, the simulations cannot determine whether Henry’s Law is valid for counter-ions in a clay system. However, MGC model calculations (Sposito, 1992), X-ray reflectivity experiments (Schlegel et al., 2006; Lee et al., 2010, 2012), and MD simulation studies (Marry et al., 2002; Rotenberg et al., 2007a; Holmboe and Bourg, 2014) have provided strong evidence that cations remain condensed at clay mineral–water and mica–water interfaces, even under bulk solution conditions approaching infinite dilution. The results presented here indicate MGC model calculations and MD simulations are mutually consistent, even in highly compacted clays, and, therefore, suggest that the DE model cannot be valid, even in highly compacted clay.

The PB model does not rely on the validity of Henry’s Law. Consequently, the MEP model predictions are similar to those of the DE model and it remains valid within the limits of the underlying molecular-scale hypotheses and approximations. The interlayer pore scale of the 3WL system can be considered as the lower limit for water content in this model. Below this scale, the mean force potential that applies to Cl ions entering the interlayer region cannot be approximated using mean field theories, because Cl is significantly more weakly solvated by confined water than by bulk liquid water. Within the context of the MEP (or MGC) model, this effect can be modeled in a straightforward manner by accounting for the distance of closest approach of ions to the clay surface, a ~ 2 Å.


The present MD simulations demonstrated that anion exclusion from clay interlayer nanopores is significantly greater than expected from mean field theories, such as the modified Gouy-Chapman (MGC) or the simplified mean electrostatic potential (MEP) models. The calculations showed that this effect can be viewed as arising from the weaker solvation of Cl by confined water than by bulk liquid water and it can be modeled by accounting for the closest approach distance of ions to the clay surface, a ~ 2 Å. Calculations reported here further suggest the existence of a threshold value in nanopore width below which Cl cannot readily access clay interlayer nanopores.

Comparison of the MD simulation results with the underlying hypothesis of the Donnan model established that the Donnan model cannot be applied to clay-water mixtures because these systems do not satisfy the two basic assumptions made about the behavior of ions in a Donnan system. This problem is related to the geometry of the charged surfaces (infinite plane), which leads to counterion condensation in the vicinity of the surface (Sposito, 1992), thereby preventing the full dissociation of Na+ ions from the surface at infinite dilution required by the Donnan model. Although a mean electrostatic potential model derived from the Poisson-Boltzmann equation may be equivalent numerically to the Donnan model, it is fundamentally different in terms of the underlying hypotheses.

The ability of Cl to enter clay interlayer nanopores decreases sharply as pore size decreases from 9.1 Å (3WL hydrate) to 5.7 Å (2WL hydrate). This finding should be taken into account when predicting the membrane properties of compacted clay materials. The present results suggest that clay materials may behave as nearly ideal semi-permeable membranes under conditions where only 2WL interlayer spaces (or inter-particle pores of similar width) connect the diffusion paths between larger pores, in agreement with recent experimental results by Chagneau et al. (2015).


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 DE-AC02-05CH11231. C. Tournassat acknowledges funding from L’Institut Carnot BRGM for his visit to the Lawrence Berkeley National Laboratory, and from the CNRS-NEEDS project TRANSREAC. G. Sposito acknowledges funding from the University of California at Berkeley under the auspices of a Chancellor’s Professorship. The research reported in this paper used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC-02-05CH11231.

This paper is published as part of a special issue on the subject of ‘Computational Molecular Modeling’. Some of the papers were presented during the 2015 Clay Minerals Society-Euroclay Conference held in Edinburgh, UK.