INTRODUCTION

Nanoporous media consist of homogeneous or heterogeneous porous material in which a significant part of the pore size distribution lies in the nanometer range. Clayey rocks, sediments or soils are natural nanoporous media, and cementitious materials, the most widely used industrial materials in the world, are also nanoporous materials. Nanoporous materials also include compounds present at the interface between non-porous solids and a aqueous solution. For example, hydrous silica gel coatings are interfacial nanoporous media that control the weathering rates of silicate glasses and minerals (Bourg and Steefel 2012). The understanding of the reactive and transport properties of these interfacial phases is essential to our global understanding of the long-term evolution of natural systems, such as soil formation (Navarre-Sitchler et al. 2011), nutrient cycling in the oceans (Loucaides et al. 2010), and engineered applications, such as the prediction of radionuclides release in high-level radioactive waste disposals (Grambow 2006; Collin et al. 2018a,b; Frugier et al. 2018). In colloidal suspensions, the aggregation of nanoparticles can also lead to the formation of nanoporous aggregates in which the bulk properties of nanoparticles are strongly influenced by the surrounding nanopores. For example, the dynamics of contaminant retention in ferrihydrite aggregates can be slowed by diffusional processes in the into/out of aggregates (Beinum et al. 2005).

The bulk fluid transport properties and the in situ chemical reactivity properties of nanoporous media are notoriously difficult to characterize. Because of their large specific surface area, most of the fluid volume in nanoporous media is influenced by the close proximity of mineral1 surfaces, which explains the very low transmissivity of these materials. As a consequence, the experimental characterization of their permeability requires special techniques (Neuzil and Person 2017). Also, the large specific surface area of nanoporous material provides them with very high adsorption capacity. The strong adsorption and resulting retardation of many contaminants by nanoporous material make them ideal for use in natural or engineered barrier systems or in filtration technologies. A good understanding of their chemical reactivity coupled to their transport properties is necessary to predict the long-term evolution of these properties of interest as a function of a range of physical and chemical conditions and processes. In this regard, reactive transport modeling can help bridging the gap between current process knowledge and predictions of the long term evolution of natural and engineered nanoporous materials in geological and industrial settings. However, nanoporous media exhibit a remarkable array of macro-scale properties with marked departures from those observed in “conventional” porous media such as permeable aquifers, for the study of which reactive transport models and codes have been historically developed. These properties arise from the interactions of charged mineral surfaces with water and solutes present in the nanopores, which leads to coupling between flux terms. These couplings manifest themselves in macroscopic observations that have intrigued geologists for more than one century, such as geologic ultrafiltration, i.e., the accumulation of solutes on the inflow side of clay-rich lithologies (Lynde 1912; Neuzil and Person 2017).

The vast majority of published reactive transport studies dealing with clay and cement materials are related to the evaluation of the long-term stability of surface and underground radioactive waste storage systems (Claret et al. 2018 and references therein; Bildstein et al. 2019, this volume, and references therein). In these types of simulation, the modeling effort has been focused primarily on the reactivity of the nanoporous materials rather than on their transport properties. Traditionally, Fickian diffusion has been typically considered, i.e., without taking into account advection, and without taking into account the anomalous transport properties of nanoporous media. In the last decade, special capabilities have been developed in a limited number of reactive transport codes that make them able to model part of these unconventional properties, with consideration of coupled processes that go beyond the traditional coupling between advective flow, dispersion, diffusion and reactions. The present chapter reviews these recent developments and explores the need to develop additional code capabilities to encompass relevant properties of nanoporous media in a holistic model from the micro-continuum to macro-continuum scales. The micro-continuum scale is restricted here to the definition given by Steefel et al. (2015b), i.e., a scale with resolution intermediate between true pore scale models and macro-continuum models, and in which parameters and properties such as permeability or reactive surface area need to be averaged or upscaled in some fashion. This corresponds to the matrix domains in the hybrid micro-continuum scale description defined by Soulaine et al. (2016, 2018), in which flow in true pore scale domains are described with the Navier–Stokes equation, whereas flow in matrix is described with Darcy's Law. Clayey materials are, by far, the most studied nanoporous media using reactive transport modeling. Consequently, most of the examples described in this chapter deal with clays. Also, this review is limited to coupled processes in water-saturated nanoporous media.

This chapter begins with a description of nanoporous media and the semi-permeable behavior that arises from the combination of their microstructure and their surface properties. Then, non-coupled and coupled transport processes in porous media are defined. In the next section, a quantitative description of the ion concentration distribution in the porosity of nanoporous materials is given, together with the approximations used in reactive transport codes to calculate it effectively. The reactive transport treatment of coupled transport processes is then described in the two last sections, starting from a description of transport in the diffusive regime, and ending with the challenges associated with the consideration of advective flow in nanoporous media.

NANOPOROUS MEDIA: SMALL PORES, LARGE SURFACES AND MEMBRANE PROPERTIES

Surface and microstructural properties of nanoporous materials

Most nanoporous media found in geological and engineered settings are comprised of an assembly of layered minerals. Bentonite is used in geosynthetic clay liners for isolation of landfills, and is planned as a barrier material in radioactive waste repository sites (Gates et al. 2009; Guyonnet et al. 2009). Bentonite is composed primarily of clay minerals of the smectite group, particularly montmorillonite. Clayey rocks (claystone, mudstone, argillite, shales) are being investigated as cap-rocks for underground gas sequestration technologies, geological barriers for radioactive waste disposal and source rocks for unconventional hydrocarbons (Neuzil 2013; Bourg 2015). These rocks usually contain large amounts of illite and illite/smectite mixed-layer minerals (Gaucher et al. 2004). The fundamental structural unit of smectite and illite consists of layers made of a sheet of edge-sharing octahedra fused to two sheets (2:1 or TOT layers) of corner-sharing tetrahedra with a thickness of about 9.5 Å (Fig. 1). The cations in the octahedral sheet of illite and montmorillonite consist predominantly of Al3+. In montmorillonite, isomorphic substitutions take place mostly in the octahedral sheets where Al3+ is replaced by Fe3+ or a cation of lower charge (Mg2+, Fe2+), whereas in illite isomorphic substitutions take place in the octahedral sheet as in montmorillonite, and in the tetrahedral sheet where Si4+ is partially replaced by Al3+. Isomorphic substitutions by cations of lower charge results in a permanent negative layer charge (Brigatti et al. 2013). The permanent negative layer charge of clay minerals is compensated for the most part by the accumulation of compensating cations between the layers. These interlayer cations can be solvated, as in montmorillonite interlayer space, or not, as in the case of K+ in illitic interlayer space (Fig. 1). In the former case, interlayer cations can be readily exchanged with other cations present in the surrounding pore water, reflecting thus the composition of the pore water. Illite and montmorillonite particles are formed by a stacking of a small number of TOT layers. Illite particles typically consist of 5 to 20 stacked TOT layers, while for smectite the number of layers per particle increases with decreasing water chemical potential and with the valence of the charge-compensating cation (Banin and Lahav 1968; Shainberg and Otoh 1968; Schramm and Kwak 1982; Saiyouri et al. 2000). The nanometric dimension of a small number of stacked layers must be compared with the basal plane dimensions that range 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). Clay mineral particles are thus highly anisotropic, and their spatial organization in the matrix of a clayey materials can yield a wide distribution of pore sizes, the dimension of which range from less than one nanometer to several micrometers (Fig. 2). Inter-particle and inter-aggregate pores are usually elongated in the bedding direction or perpendicular to the compaction direction in remolded samples (Pusch 2001; Robinet et al. 2012; Gu et al. 2015; Gaboreau et al. 2016). These pores are usually represented by slit-shaped or cylindrical-shaped pores for modeling purposes. The smallest pores, i.e., the interlayer pores in swelling clay minerals, are logically usually described as slit-shaped pores (Birgersson and Karnland 2009; Appelo et al. 2010; Tournassat and Appelo 2011; Tournassat et al. 2016b; Appelo 2017; Wigger and Van Loon 2017).

Figure 1.

Examples of layered phases in clayey and cementitious materials. [Clay mineral structures reprinted from Tournassat et al. (2015a) Chapter 1—Surface properties of clay minerals. In: Tournassat C, Steefel CI, Bourg IC, Bergaya F (Eds.), Natural and Engineered Clay Barriers, Developments in Clay Science. Elsevier, p 5–31, Copyright 2015, with permission from Elsevier. C–S–H layer structure according to Grangeon et al (2016). AFm structure according to Marty et al. (2018).]

Figure 1.

Examples of layered phases in clayey and cementitious materials. [Clay mineral structures reprinted from Tournassat et al. (2015a) Chapter 1—Surface properties of clay minerals. In: Tournassat C, Steefel CI, Bourg IC, Bergaya F (Eds.), Natural and Engineered Clay Barriers, Developments in Clay Science. Elsevier, p 5–31, Copyright 2015, with permission from Elsevier. C–S–H layer structure according to Grangeon et al (2016). AFm structure according to Marty et al. (2018).]

Figure 2.

A: 3-D visualization of the microstructure of an illite sample compacted at 1.7 kg⋅dm−3 (backscattered electrons focused ion beam nanotomography stack images 350 × 350 × 180 voxels with voxel resolution size of 5 × 5 × 5 nm3); B: Segmented image of a TEM micrograph of the same illite sample. The blue color represents pores and light and dark gray are pore/solid mixtures and solids, respectively. C: Pore-width distribution obtained from a combination of bulk and microscopic characterization techniques. [Figure modified from Tournassat et al. (2016b).]

Figure 2.

A: 3-D visualization of the microstructure of an illite sample compacted at 1.7 kg⋅dm−3 (backscattered electrons focused ion beam nanotomography stack images 350 × 350 × 180 voxels with voxel resolution size of 5 × 5 × 5 nm3); B: Segmented image of a TEM micrograph of the same illite sample. The blue color represents pores and light and dark gray are pore/solid mixtures and solids, respectively. C: Pore-width distribution obtained from a combination of bulk and microscopic characterization techniques. [Figure modified from Tournassat et al. (2016b).]

Cementitious materials are complex assemblages of a range of solid phases. Among them, nanocrystalline calcium silicate hydrate (C–S–H) is often the main hydration product (Richardson 2008), and constitutes the “glue material” of the cement (Pellenq and Van Damme 2004; Masoero et al. 2012). C–S–H have a complex chemistry, with calcium to silicon (Ca/Si) atomic ratio that vary as a function of pore water chemical conditions, and especially pH. C–S–H have layered structures that can be described as nanocrystalline and defective tobermorite, in which the Si chains depolymerize as the Ca/Si atomic ratio increases (Grangeon et al. 2013, 2016, 2017; Richardson 2014). The depolymerization reaction creates silanol groups at the surface having amphoteric (protonation/deprotonation) properties and surface complexation capabilities (Fig. 1). The resulting surface charge of C–S–H can take positive or negative values as a function of pore water chemistry following layer structure changes and surface site reactivity (Haas and Nonat 2015). A C–S–H particle is composed of a stack of a small number of layers, which are separated by a hydrated interlayer space that may also contain cations (Grangeon et al. 2016; Roosz et al. 2016). The spatial distribution of C–S–H particles in the cement matrix creates a heterogeneous distribution of inter-particulate pore sizes (Ioannidou et al. 2016). Other major cement phases include the hydrated calcium aluminate phase (AFm), which is a member of the layered double hydroxides (LDHs) group. Its structure consists of stacked positively charged layers of Al and Ca polyhedra separated from each other by hydrated interlayer spaces, which contain charge compensating hydrated anions (Fig. 1). Similarly to interlayer cations in swelling clay minerals, interlayer anions in AFm are exchangeable with the surrounding pore water, making them an important contributor to the retention of anions that enter into contact with cement-based materials (Ma et al. 2017, 2018; Marty et al. 2018).

The layered nature of the major components of clayey and cementitious materials confer to these materials very high specific surface area, and, consequently, a large fraction of their porosity lies within the nanometer range. C–S–H, AFm, illite and smectite layers exhibit surface charges that are compensated in the adjacent porosity. Between the layers, these charges are mostly compensated by counter-ions (ions with charge sign opposite to that of the surface) in the interlayer porosity. The layer charge present at the outer surfaces of the particles (hereafter named basal surfaces), which border the interparticle pores, must also be compensated. This is also true for non-permanent charges present on the border surfaces of other solid phase (e.g., silica) that arise from the cleavage of the grains. The charge compensation mechanisms and their consequences for the water and ion transport in nanoporous media are explored in the next section.

Ion distribution in the vicinity of charged surfaces

The charge of the surfaces bordering the pores is responsible for the presence of a double layer or electrical double layer (EDL), i.e., the layers of interfacial water and electrolyte ions that screen the surface charge. The EDL is often conceptually subdivided into two regions, in agreement with spectroscopic results (Lee et al. 2010, 2012b), direct force measurements (Zhao et al. 2008; Siretanu et al. 2014) and molecular dynamics (MD) calculations (Marry et al. 2008; Tournassat et al. 2009a; Rotenberg et al. 2010; Bourg and Sposito 2011): the Stern layer located within the first water monolayers of the interface in which ions adsorb as inner-and outer-sphere surface complexes (ISSC, OSSC) and a diffuse layer (DL) located beyond the Stern layer in which a diffuse cloud or swarm of ions screens the remaining uncompensated surface charge (Fig. 3). The concentrations of ions in this region depend on the distance from the surface considered. At infinite distance from the charged surface, the solution is neutral and is commonly described as bulk or free solution (or water). The exact 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 (Delville 2000). It is important to note that the ion distribution in the EDL cannot be probed directly by measurements because EDL features are inherently disturbed by direct measurement techniques (Bourg et al. 2017). A range of experimental techniques have been used to characterize indirectly the EDL properties, but the interpretation of their results is always sensitive to an a priori choice of an EDL model. Fortunately, in the last two decades, atomistic level models, such as molecular dynamics, have provided a growing set of information on the properties of the EDL. It is now commonly acknowledged that the diffuse layer part of the EDL is well understood, while the characterization of the Stern layer remains a challenging issue (Bourg et al. 2017). In the following, we shall concentrate our attention on the diffuse layer, which has a large influence on transport properties in nanoporous media.

Figure 3.

A: Water and ion distribution at the vicinity of a montmorillonite surface according to MD simulations. Average density (mol⋅dm−3) of water (black), Na (dark blue), Ca (light blue) and Cl (orange) are plotted as a function of distance from the mid-plane of the montmorillonite layer at a NaCl/CaCl2 brine concentration of 0.34 mol⋅dm−3. B: Conceptual model of EDL structure. [Reprinted from Bourg and Sposito (2011) Molecular dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl–CaCl2) solutions. Journal of Colloid and Interface Science 360:701–715 Copyright 2011, with the permission of Elsevier.]

Figure 3.

A: Water and ion distribution at the vicinity of a montmorillonite surface according to MD simulations. Average density (mol⋅dm−3) of water (black), Na (dark blue), Ca (light blue) and Cl (orange) are plotted as a function of distance from the mid-plane of the montmorillonite layer at a NaCl/CaCl2 brine concentration of 0.34 mol⋅dm−3. B: Conceptual model of EDL structure. [Reprinted from Bourg and Sposito (2011) Molecular dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl–CaCl2) solutions. Journal of Colloid and Interface Science 360:701–715 Copyright 2011, with the permission of Elsevier.]

In the diffuse layer, water is often considered to have properties similar to bulk-liquid water and the distribution of ions predicted by molecular dynamic simulations is often in close agreement with predictions made with the Poisson–Boltzmann (PB) equation (Tinnacher et al. 2016). The Poisson equation relates the local charge imbalance at a position y in the direction perpendicular to the charged surface to the second derivative of the electrostatic potential (x in V) at the same position:

 
d2ψDL(y)dy2=1εΣiziFci,DL
(1)

where ε is the water dielectric constant (78.3 × 8.85419·10−12 F·m−1 at 298.15 K), F is the Faraday constant (96 485 C·mol−1), and zi is the charge of ion i. The concentration of species i in the DL, ci,DL, is given by the Boltzmann equation:

 
ci,DL=ci,0exp(ziFψDL(y)RT)
(2)

where R is the gas constant (8.3145 J·mol−1·K−1), T is the temperature (K), and ci,0 is the concentration of species i in the bulk water at equilibrium with the diffuse layer. At infinite distance from the surface, the electrostatic potential vanishes and the concentration in the diffuse layer is equal to that in the bulk water, where electro-neutrality prevails. However, if pores are not large enough, there will be no regions within the pore that meet this electro-neutrality condition because of the overlap of the diffuse layers. The size of the diffuse layer depends primarily on the concentration of ions, and on the valence of the counter-ions (Carnie and Torrie 1984; Sposito 2004). As the ionic strength (I, here in mol⋅L−1) increases, the diffuse layer thickness decreases, and, as the valence of the counter-ions increases, the diffuse layer thickness decreases as well. So too does the overlap of the diffuse layers at two neighboring surfaces (Fig. 4, top and middle).

Figure 4.

Top and middle figures: Ion enrichment/depletion according to the Poisson–Boltzmann equation in a diffuse layer in contact with negatively charged surfaces (σ0 = 0.11 C⋅m−2) bordering a slit-shaped pore with a distance of 2 nm between the surfaces. A and D: NaCl electrolyte; B and E: CaCl2 electrolyte; C and F: diluted seawater. A, B and C correspond to an ionic strength of 0.72 mol L−1, while D, E and F correspond to a ionic strength of 0.072 mol L−1. The seawater composition was taken from Holland (1978) and was simplified into Na+: 481 mmol⋅kg−1, Mg2+: 55 mmol⋅kg−1; Ca2+: 11 mmol⋅kg−1; K+: 11 mmol⋅kg−1; Cl 564 mmol⋅kg−1; SO42−: 29 mmol⋅kg−1; HCO3: 2 mmol⋅kg−1. No aqueous complexation in bulk solution was considered for this calculation. In figures C and F, the enrichment/depletion are indicated for monovalent cations (Cat+), divalent cations (Cat+2), monovalent anions (An) and divalent anions (An−2). Calculations were done at 298 K. Bottom figures (G, H and I): Poisson–Boltzmann prediction of the electrostatic potential (in V) at the middle of a slit-shaped pore as a function of the inverse of the square root of ionic strength (x-axis) and half distance between the surfaces (y-axis) (σ0 = −0.11 C⋅m−2). The blue lines indicate 1, 2 and 3 times the Debye length (k−1, Eqn. 21), which is commonly assumed to be proportional to the thickness of the diffuse layer.

Figure 4.

Top and middle figures: Ion enrichment/depletion according to the Poisson–Boltzmann equation in a diffuse layer in contact with negatively charged surfaces (σ0 = 0.11 C⋅m−2) bordering a slit-shaped pore with a distance of 2 nm between the surfaces. A and D: NaCl electrolyte; B and E: CaCl2 electrolyte; C and F: diluted seawater. A, B and C correspond to an ionic strength of 0.72 mol L−1, while D, E and F correspond to a ionic strength of 0.072 mol L−1. The seawater composition was taken from Holland (1978) and was simplified into Na+: 481 mmol⋅kg−1, Mg2+: 55 mmol⋅kg−1; Ca2+: 11 mmol⋅kg−1; K+: 11 mmol⋅kg−1; Cl 564 mmol⋅kg−1; SO42−: 29 mmol⋅kg−1; HCO3: 2 mmol⋅kg−1. No aqueous complexation in bulk solution was considered for this calculation. In figures C and F, the enrichment/depletion are indicated for monovalent cations (Cat+), divalent cations (Cat+2), monovalent anions (An) and divalent anions (An−2). Calculations were done at 298 K. Bottom figures (G, H and I): Poisson–Boltzmann prediction of the electrostatic potential (in V) at the middle of a slit-shaped pore as a function of the inverse of the square root of ionic strength (x-axis) and half distance between the surfaces (y-axis) (σ0 = −0.11 C⋅m−2). The blue lines indicate 1, 2 and 3 times the Debye length (k−1, Eqn. 21), which is commonly assumed to be proportional to the thickness of the diffuse layer.

Semi-permeable properties

The presence of overlapping diffuse layers in charged nanoporous media is responsible for a partial or total repulsion of co-ions from the porosity. In the presence of a gradient of bulk electrolyte concentration, co-ion migration through the pores is hindered, as well as the migration of their counter-ion counterparts because of the electro-neutrality constraint. This explains the salt-exclusionary properties of these materials. These properties confer these media with a semi-permeable membrane behavior: neutral aqueous species and water are freely admitted through the membrane while ions are not, giving rise to coupled transport processes. Semi-permeable membrane properties of nanoporous materials have been extensively studied for decades. An overwhelming majority of these studies have dealt with clay materials, in the form of compacted powder samples of purified clay, or natural materials, such as soils or clayey rocks, with a strong emphasis on chemo-osmotic membrane properties (Lynde 1912; McKelvey and Milne 1960; Kemper 1960, 1961a,b; Olsen 1962; Kemper and Evans 1963; Young and Low 1965; Kemper and Rollins 1966; Bresler 1973; Hanshaw and Coplen 1973; Neuzil 2000; Cey et al. 2001; Garavito et al. 2007; Horseman et al. 2007; Rousseau-Gueutin et al. 2009; Gonçalvès et al. 2015; Sun et al. 2016).

NON-COUPLED AND COUPLED TRANSPORT PROCESSES

Transport equations in traditional reactive transport modeling

Single phase numerical reactive transport models aim at describing the fluxes Ji (mol⋅mmedium−2⋅s−1) of aqueous chemical species i in space and time, as well as the mineralogy and sorbed species evolution. The constitutive equations describing these fluxes at the continuum scale are thermodynamic and kinetic equations for chemical reactions coupled to advection (adv), mechanical dispersion (disp) and diffusion equations (diff) (Steefel and Maher 2009; Steefel et al. 2015a):

 
Ji=Jadv,i+Jdisp,i+Jdiff,i
(3)

Advection, dispersion and diffusion are non-steady, irreversible processes expressed in the form of partial differential field equations (Bear 1972). The advective flux is related to the average linear velocity, v (m⋅s−1), in the media:

 
Jadv,i=ϕCiv
(4)

where ϕ is the porosity (–), Ci (mol⋅mwater−3) is the aqueous concentration of species i. The average linear velocity in the porous media is commonly evaluated with Darcy's Law, which relates the volumetric flux of water (mfluid3⋅mmedium2⋅s−1), q, to the gradient in the hydraulic head, h (m), and the hydraulic conductivity of the medium, K (m⋅s−1):

 
q=ϕv=Kh
(5)

Molecular diffusion in porous media is usually described in terms of Fick's First Law:

 
Jdiff,i=De,iCi
(6)

where De,i is the effective coefficient of species i. The effective diffusion coefficient differs from the diffusion coefficient in water, because diffusion takes place only in the aqueous solution present in the medium. Consequently, the diffusion coefficient must be corrected using the porosity value and from a geometrical factor, τi (–), that accounts for the difference of paths lengths that the aqueous species would follow in water alone compared to the tortuous paths lengths it would follow in the porous medium (tortuosity), and for the constrictivity factor of the porous medium (Bear 1972):

 
De,i=ϕτiD0,i
(7)

where D0,i is the diffusion of species i in the solution alone. In the above formulation, the diffusive fluxes depend on the nature of the chemical species. This dependency would result in charge imbalance in conventional numerical reactive transport models if only a Fickian description of diffusion is used, so in practice, Equation (7) is usually simplified into:

 
JDiff,i=ϕτD0Ci
(8)

where τ and D0 are mean representative values.

The coefficient of hydrodynamic dispersion is defined as the sum of molecular diffusion and mechanical dispersion. Mechanical dispersion is a scale-dependent process, with larger dispersivities observed for larger observation scales (Steefel 2008). It is usually computed in reactive transport codes as the product of the fluid velocity and dispersivity with longitudinal and transverse components (Steefel et al. 2015a; Rolle and Le Borgne 2019, this volume and references therein). The low permeability of nanoporous media limits mass transport by advection, with the result that hydrodynamic dispersion is dominated in this case by molecular diffusion (Patriarche et al. 2004a,b; Mazurek et al. 2011). Consequently, the contribution of mechanical dispersion to the total flux is usually neglected for the modeling of nanoporous media (Gonçalvès et al. 2015; Neuzil and Person 2017).

In traditional numerical reactive transport models applied to nanoporous media, transport of solutes is thus described with the above two constitutive equations, Darcy's and Fick's laws. While Fick's and Darcy's Laws constitute a modeling framework that is adequate to deal with a wide range of reactive transport problems for various geological settings and processes such as those described in Maher and Navarre-Stichler (2019, this volume), Cama et al. (2019, this volume), and Lagneau et al. (2019, this volume) these laws are phenomenological equations rooted in the results of physical experiments carried out on systems involving salts in water and packed sand columns respectively (Fick 1855; Darcy 1856). These systems are very different from the nanoporous media we are addressing in the present contribution. Because of the charge imbalance in the solution present in the diffuse layer, the rigorous modeling of transport processes in nanoporous media is not possible by considering the advective flux equation based on the classical Darcy equation and the diffusive flux based on Fick's equations.

Coupled transport processes and thermodynamics of irreversible processes

The thermodynamics of irreversible processes (or non-equilibrium thermodynamics), founded by Onsager (1931a,b), provides the theoretical background for the existence of coupled flows (Kjelstrup and Bedeaux 2008).

Many observations of fluid and solute transport in geological media cannot be related to the Darcian and Fickian modeling framework because of the existence of coupling of forces of one type and flows of another type (Bear 1972). For example, chemical osmosis is a flow of water driven by a gradient of water chemical potential, which is usually created by a gradient of electrolyte concentration across a semi-permeable membrane. Water flow can also result from the presence of gradients in the electric potential (electro-osmosis) or temperature (thermo-osmosis). Osmotic processes have been extensively characterized on natural materials ranging from soils to biological membranes (Medved and Cern 2013).

Non-equilibrium thermodynamics is based on the second law of thermodynamics and the principle of microreversibility. In any system, the rate of entropy production per unit volume, σ, can be written:

 
σ=1VdSdt=ΣiJiXi
(9)

where S is the entropy of the system, V its volume, t the time, and Ji is a flux of type i associated with a driving force Xi for that particular flux (Lasaga 1998). The second law of thermodynamics predicts that:

 
ΣiJiXi0
(10)

With the assumption that all driving forces Xi needed to describe fully the evolution system are taken into account in Equation (10), each flux Ji is a function fi of the forces Xi, and takes a value of zero if all Xi are zero (reference point). Accordingly, the functions fi can be expanded in Taylor series of the form (Lasaga 1998):

 
Ji=fi(X1,X2,,Xn)=0+Σjfi(X1,X2,,Xn)Xj,X=0Xj+
(11)

where X = 0 indicates that the derivative of the function is evaluated at the reference point. Keeping only the linear terms of the Taylor series expansion, i.e., making the assumption that higher order terms are negligible, Equation (11) yields:

 
Ji=Σjfi(X1,X2,,Xn)Xj,X=0Xj=ΣjLi,jXj
(12)

where Li,j are termed phenomenological coefficients. Equation (12) indicates that a driving force Xj not only influences the conjugate flux Jj, but also the other non-conjugate fluxes Ji≠j as in the case of the various osmosis processes. The use of irreversible thermodynamics theory is thus dependent on the identification of a set of extensive independent variable Xj that defines completely the system of interest (de Groot and Mazur 1984; Kjelstrup and Bedeaux 2008). Under these conditions, Onsager (1931b) demonstrated that the following reciprocal relations apply:

 
Li,j=Lj,i
(13)

Transport problems investigated in Earth Sciences are focused mostly on the fluxes of fluids, chemical components, and heat. We shall see in the following section that in the case of aqueous solutions present in porous media in general, and in nanoporous media in particular, charge (or electrical) fluxes must also be considered. Ideally, all types of conjugate and coupled flux listed in Table 1 should be considered. However, Equations (12) and (13) do not give any information about the way fluxes and conjugate forces must be chosen, and the set of equations is not unique for a given system. Consequently, the problem one wants to describe dictates the form that is most convenient (Kjelstrup and Bedeaux 2008). Numerical methods for modeling macroscopic transport properties of nanoporous media, and especially clayey materials, with the explicit consideration of the presence of a diffuse ion swarm have received growing interest in diverse research communities in the past years (Revil 1999; Revil and Leroy 2004; Leroy et al. 2006; Gonçalvès et al. 2007, 2015; Jougnot et al. 2009; Revil et al. 2011). If we assign the value X1, X2, and X3 to the fluid pressure, temperature and electrical potential gradients respectively, and the value X4 to Xn+3 to each of the driving forces associated to the presence of n chemical species, then a complete description of coupled flux processes is given by:

 
J1=L1,1X1+L1,2X2+L1,3X3+L1,4X4++L1,n+3Xn+3J2=L2,1X1+L2,2X2+L2,3X3+L2,4X4++L2,n+3Xn+3J3=L3,1X1+L3,2X2+L3,3X3+L3,4X4++L3,n+3Xn+3J4=L4,1X1+L4,2X2+L4,3X3+L4,4X4++L4,n+3Xn+3=Jn+3=Ln+3,1X1+Ln+3,2X2+Ln+3,3X3+Ln+3,4X4++Ln+3,n+3Xn+3
(14)
Table 1.

Terminology of coupled fluxes (after de Marsily 1986)

FluxGradients
Fluid pressure ∇pTemperature ∇T/TElectrical potential ∇ψChemical potential∇μi
FluidDarcyThermo-osmosisElectro-osmosisChemical osmosis
HeatHeat filtration or mechano-caloric effectFourierPeltier effectDufour effect
ElectricalStreaming potential (Rouss effect)Seebeck effectOhmDiffusion current
Aqueous species i (relative to the fluid)UltrafiltrationThermal diffusion (Soret effect)ElectrophoresisFick
FluxGradients
Fluid pressure ∇pTemperature ∇T/TElectrical potential ∇ψChemical potential∇μi
FluidDarcyThermo-osmosisElectro-osmosisChemical osmosis
HeatHeat filtration or mechano-caloric effectFourierPeltier effectDufour effect
ElectricalStreaming potential (Rouss effect)Seebeck effectOhmDiffusion current
Aqueous species i (relative to the fluid)UltrafiltrationThermal diffusion (Soret effect)ElectrophoresisFick

It is clear that the evaluation of the phenomenological coefficients becomes increasingly more cumbersome as the numbers of fluxes and driving forces increase (Malusis et al. 2012), explaining why most of the effort has been put into predictive models restricted to simple electrolytes conditions (e.g., NaCl). If the electrolyte is not simple, then it is frequently assumed that the mobility of the different ions does not differ significantly from each other (Revil and Linde 2006). This approximation simplifies greatly the set of equations, but is a significant departure from reality, especially in systems at low or high pH, because the diffusion coefficient of OH and H3O+ are five to ten times greater than the diffusion coefficients of other ions (Li and Gregory 1974). The models and the formulation of the phenomenological coefficients are thus very dependent on the chemical conditions of the systems studied. Reactive transport modeling approaches, with their unique ability to deal with transient stages and chemically complex systems, may be able to overcome some of these limitations if it is possible to define equations linking the values of the phenomenological coefficients to the basic properties of the media of interest and to the chemical composition of the pore water. In the following, we shall explore in detail the sets of general equations that have been or could be implemented in reactive transport codes to deal with coupled transport processes in the presence of a diffuse layer. We have also highlighted that, in certain cases, the hindered mobility of co-ions controls the transport properties of all ions because of the requirement of charge balance at the macro-scale. A quantitatively correct description of ion concentration distribution in the diffuse layer volume is thus required to quantify ionic fluxes in nanoporous media. Consequently, the following section focuses on the methods used to compute the diffuse layer composition.

POISSON–BOLTZMANN EQUATION AND ION CONCENTRATION DISTRIBUTION CALCULATIONS IN REACTIVE TRANSPORT CODES

The mean electrostatic potential model

The Poisson–Boltzmann equation can be solved analytically, without approximations, in a very limited number of cases (Andrietti et al. 1976; Borkovec and Westall 1983; Chen and Singh 2002; Sposito 2004), and its solution for systems having complex electrolyte composition and/or geometries is possible only with numerical approaches (Leroy and Maineult 2018). In addition, nanoporous materials such as clayey and cementitious materials typically exhibit a range of pore size and shapes. For example, clay mineral particles are often segregated into aggregates delimited by inter-aggregate spaces, the size of which is usually larger than the inter-particle spaces inside the aggregates. In clayey 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; Robinet et al. 2012; Philipp et al. 2017) and these co-exist with pores having a width as narrow as the clay interlayer spacing, i.e., one nanometer.

In practice, the information about ion concentrations in the diffuse layer must be upscaled so that calculations can be carried out at the continuum scale. A common upscaling approach relies on the use of a mean electrostatic potential model (Revil and Leroy 2004; Appelo and Wersin 2007; Appelo et al. 2010; Revil et al. 2011; Tournassat and Steefel 2015; Tournassat et al. 2016b). This model averages ion concentrations in the diffuse layer Ci,DL¯ by scaling them to a mean electrostatic potential (MEP, ψM) that applies to the diffuse layer volume (VDL):

 
ci,DL¯=1VDLDLzici,0exp(ziFψDL(x,y,z)RT)dxdydzci,0exp(ziFψMRT)
(15)

The mean electrostatic potential value can be deduced from the charge balance in the diffuse layer:

 
ΣiziFci,DL¯=ΣiziFci,0exp(ziFψMRT)=QDL
(16)

where QDL(in C⋅m−3) is the volumetric charge that must be balanced in the diffuse layer. The accuracy of this method compared to a full resolution of the Poisson–Boltzmann equation can be appreciated by calculating the difference of average concentrations in the diffuse layer volume computed by the two methods in a simple geometry such as a slit-shaped pore (Fig. 5). When applied to the entire pore volume, the mean electrostatic potential model underestimates the average concentration of co-ions, and the accuracy of its predictions degrades as the size of the pore increases and as the ionic strength of the solution decreases (Fig. 5 A, B, C, D). The prediction of the average concentration of counter-ions is good for the simple NaCl and CaCl2 electrolytes. It must be noted, however, that the average concentration value is governed primarily by the requirement of charge balance with the surface in the diffuse layer. Because the surface charge is mostly balanced by the increase of concentration of counter-ions in the diffuse layer, it is not surprising to find a good agreement between the two models. In the presence of a more complicated electrolyte composition (seawater dilution), the mean electrostatic potential model overestimates the contribution of the monovalent counter-ions in the diffuse layer, and underestimates the contribution of the divalent counter-ions. The mean electrostatic potential model is often referred as a Donnan model in the literature (Revil and Leroy 2004; Appelo and Wersin 2007; Birgersson and Karnland 2009; Appelo et al. 2010; Tournassat and Appelo 2011; Alt-Epping et al. 2015, 2018; Gimmi and Alt-Epping 2018). The Donnan equilibrium model, however, 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 Donnan model gives strictly equivalent results to the mean electrostatic potential model only if two important conditions are met: first, the mean electrostatic potential model must be applied to the entire pore volume, and second, in the Donnan equilibrium model, the activity coefficients of chemical species present in the Donnan volume must be set at the same value as those in bulk water (Tournassat et al. 2016a). The consideration of a Donnan equilibrium model in clay media faces two challenges. First, the activity of chemical species in the Donnan volume cannot be measured, and thus the choice to set the activity coefficients of species in the Donnan volume to the same values as in bulk water must be considered as arbitrary. This choice can be relaxed by choosing other relationships for the calculation of the activity coefficients, but experimental data are lacking to calibrate such relationships (Birgersson 2017). Second and more importantly, a strict application of Henry's Law (and, hence, of the Donnan equilibrium model) 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). This second hypothesis is not in agreement with the observation of cations condensation at clay mineral surfaces, using diffraction techniques or molecular dynamics simulations (Marry et al. 2002; Schlegel et al. 2006; Tournassat et al. 2009a; Lee et al. 2010, 2012b; Holmboe and Bourg 2014; Bourg et al. 2017). The Poisson–Boltzmann model does not rely on the validity of Henry's Law. Consequently, the mean electrostatic model predictions remain valid within the limits of the underlying molecular-scale hypotheses and approximations of the model, while the Donnan model does not (Tournassat et al. 2016a).

Figure 5.

Difference of average diffuse layer concentration values computed with the Poisson–Boltzmann (cPB) equation and the mean electrostatic potential (cMEP) model (Eqns. 15 and 16). The difference, given by the color scale, is computed with the equation (cPBcMEP) / cPB) as a function of the inverse of the square root of ionic strength (I−1/2 on x-axis) and as a function of the half pore width (y-axis). The equivalence between the two models corresponds to (cPBcMEP) / cPB)=0. The geometry consists in a diffuse layer in contact with negatively charged surfaces (σ0 =0.11 C⋅m−2) bordering a slit-shaped pore. In figure C, D, G and H, the seawater composition was taken from Holland (1978) and was simplified into Na+: 481 mmol⋅kg−1, Mg2+: 55 mmol⋅kg−1; Ca2+: 11 mmol⋅kg−1; K+: 11 mmol kg−1; Cl 564 mmol⋅kg−1; SO42−: 29 mmol⋅kg−1; HCO3: 2 mmol⋅kg−1. This composition corresponds to the lowest value if I−1/2 on the figures. Other values of I−1/2 correspond to the dilution of this seawater composition with pure water. No aqueous complexation in bulk solution was considered for this calculation. The calculation results are shown for monovalent anions (An) in C, divalent anions (An−2) in D, monovalent cations (Cat+) in G and divalent cations (Cat+2) in H. Calculations were done at 298 K.

Figure 5.

Difference of average diffuse layer concentration values computed with the Poisson–Boltzmann (cPB) equation and the mean electrostatic potential (cMEP) model (Eqns. 15 and 16). The difference, given by the color scale, is computed with the equation (cPBcMEP) / cPB) as a function of the inverse of the square root of ionic strength (I−1/2 on x-axis) and as a function of the half pore width (y-axis). The equivalence between the two models corresponds to (cPBcMEP) / cPB)=0. The geometry consists in a diffuse layer in contact with negatively charged surfaces (σ0 =0.11 C⋅m−2) bordering a slit-shaped pore. In figure C, D, G and H, the seawater composition was taken from Holland (1978) and was simplified into Na+: 481 mmol⋅kg−1, Mg2+: 55 mmol⋅kg−1; Ca2+: 11 mmol⋅kg−1; K+: 11 mmol kg−1; Cl 564 mmol⋅kg−1; SO42−: 29 mmol⋅kg−1; HCO3: 2 mmol⋅kg−1. This composition corresponds to the lowest value if I−1/2 on the figures. Other values of I−1/2 correspond to the dilution of this seawater composition with pure water. No aqueous complexation in bulk solution was considered for this calculation. The calculation results are shown for monovalent anions (An) in C, divalent anions (An−2) in D, monovalent cations (Cat+) in G and divalent cations (Cat+2) in H. Calculations were done at 298 K.

Dual continuum representation of the pore space

We have highlighted in a previous paragraph that the semi-permeable membrane properties of nanoporous media originate from the exclusion of co-ions in the diffuse layer. It is thus desirable to build a model that reproduces as well as possible the mean co-ion concentration in the diffuse layer. The mean electrostatic potential model can be adapted to do so, following the basic principle of a subdivision of the pore space into two compartments, one being electroneutral, and the other being influenced by a non-zero mean electrostatic potential value (Fig. 6). This type of representation will be referred to as “Dual Continuum” in the following.

Figure 6.

Comparison of the Poisson–Boltzmann (left figure) and MEP (right figure) models in a slit-shaped pore. The specific surface charge is −0.11 C⋅m−2 and the bulk concentration of NaCl electrolyte is 0.1 mol⋅L−1. The half pore width (hpore) is 20 Å. The width of the diffuse layer in the dual continuum model is set to hDL = fDL × hpore, and the value of fDL is adjusted to obtain the same average pore concentrations of Cl and Na+ with the two models. Dashed curves = electrostatic potential, ψ; plain lines = Na+ concentration; dotted line = Cl concentration.

Figure 6.

Comparison of the Poisson–Boltzmann (left figure) and MEP (right figure) models in a slit-shaped pore. The specific surface charge is −0.11 C⋅m−2 and the bulk concentration of NaCl electrolyte is 0.1 mol⋅L−1. The half pore width (hpore) is 20 Å. The width of the diffuse layer in the dual continuum model is set to hDL = fDL × hpore, and the value of fDL is adjusted to obtain the same average pore concentrations of Cl and Na+ with the two models. Dashed curves = electrostatic potential, ψ; plain lines = Na+ concentration; dotted line = Cl concentration.

Using fDL as the volumetric fraction of the porosity to which the mean electrostatic potential is applied, Equations (15) and (16) transform into:

 
ci,DL¯=1VDLDLzici,0exp(ziFψDL(x,y,z)RT)dxdydzci,0exp(ziFψMRT)
(17)
 
ΣiziFci,DL¯=ΣiziFci,0exp(ziFψMRT)=QDL
(18)

where cCl,pore¯ is the average concentration of species i in the pore. Note that (1fDL)FΣizici,0 is equal to zero because bulk water is electroneutral. The value of fDL can be adjusted as a function of the conditions (pore volume, surface charge, electrolyte composition) to achieve consistency with the theoretical PB predictions or with experimental constraints. It should be also noted here that this model refinement does not imply necessarily that an electroneutral bulk water is present at the center of the pore in reality. This can be appreciated in Figure 6, which shows that the Poisson–Boltzmann predicts an overlap of the diffuse layers bordering the two neighboring surfaces, while the dual continuum model divides the same system into a bulk and a diffuse layer water volume in order to obtain an average concentration in the pore that is consistent with the Poisson–Boltzmann model prediction. Consequently, the pore space subdivision into free and DL water must be seen as a convenient representation that makes it possible to calculate accurately the average concentrations of ions, but it must not be taken as evidence of the effective presence of bulk water in a nanoporous medium. This modified MEP model was introduced in PHREEQC to model diffusion properties of Opalinus Clay at the Mont Terri underground research laboratory (URL) and of Callovian-Oxfordian argillite at the Laboratoire Souterrain de Meuse Haute Marne (LSMHM) URL (Appelo and Wersin 2007; Appelo et al. 2008). The method to calculate fDL on the basis of anion (Cl) exclusion experimental data in clay samples was later developed in Appelo et al. (2010). Average Cl concentrations in pore water (cCl,pore) can be obtained from Cl concentration measured in leaching experiments (cCl,leach), in which a mass msample of clay sample is dispersed in a volume Vleach of milli-Q water. As cCl,leach is proportional to the ratio msample/Vleach in clay samples, it follows that (Tournassat et al. 2015b):

 
ccl,pore¯=ccl,leachVleachmsample×ωw
(19)

where ωW is the water content of the sample (in kg⋅kg−1). The value of cCl,pore¯ is further compared to the Cl concentration value obtained from core squeezing experiments or seepage water composition from in situ equipped borehole, which is assumed to be representative of the bulk concentration cCl,0 (Vinsot et al. 2008a,b; Pearson et al. 2011; Tournassat et al. 2011; Fernández et al. 2014; Mazurek et al. 2015). The value of cCl,pore¯ is always lower than that of cCl,0, demonstrating the depletion of Cl in the diffuse layer. From simple mass balance consideration, it follows that:

 
ccl,pore¯ccl,0=(1fDL)+fDLexp(zclFψMRT)
(20)

in which the unknown parameters are fDL and ψM. If the bulk concentrations of the other species present in solution are also known, e.g., from analysis of squeezing or seepage waters or from (bulk) pore water composition modeling (Gaucher et al. 2009; Pearson et al. 2011; Beaucaire et al. 2012; Gailhanou et al. 2017), Equation (18) provides the necessary constraints to calculate ΨM as a function of the volumetric charge that is compensated in the diffuse layer, QDL. The value of QDL can be itself estimated on the basis of macroscopic measurements (e.g., cation exchange capacities of clays) or mineralogical considerations. However, not all of the surface charge is usually compensated in the diffuse layer, and it may be necessary to consider the contribution of the Stern layer to the charge compensation (Appelo et al. 2010). Geochemical calculation codes capabilities make the calculation possible with the use of surface complexation models (Steefel et al. 2015a), although then this surface complexation model must be coupled (loosely or tightly) to the calculation of the diffuse layer charge. In absence of additional constraints, the experimental value of cCl,pore/cCl,0 can thus be fitted by adjusting the values of fDL, QDL, or both (Leroy and Revil 2004; Leroy et al. 2007; Appelo et al. 2010; Appelo 2017). The contribution of the Stern layer to the total charge compensation can be estimated with molecular dynamic simulations, which provides a growing set of independent information that can be used to calibrate models of ion distribution and mobility in nanopores (Rotenberg et al. 2007, 2014; Jardat et al. 2009; Tournassat et al. 2009b; Bourg and Sposito 2011; Obliger et al. 2014; Tinnacher et al. 2016; Bourg et al. 2017). The dual continuum model was first developed in PHREEQC, but is also now available in CrunchClay (Tournassat and Steefel 2015, 2019; Soler et al. 2019).

Influence of ionic strength on the dual continuum model predictions

Changes in ionic strength influence the ion concentration distribution in the diffuse layers (Fig. 4), and thus fDL changes also as a function of ionic strength. In the absence of overlap, the effective thickness of the diffuse layer is proportional to the Debye length (k−1 in m) (Fig. 4, bottom) (Sposito 2004; Hunter 2013):

 
κ1=εRT2F21000I
(21)

In PHREEQC and CrunchClay, the volume of the diffuse layer (VDL in m3), and hence the fDL value, can be defined as a multiple of the Debye length in order to capture this effect of ionic strength on fDL:

 
VDL=αDLκ1SfDL=VDLVpore
(22)

where Vpore is the total volume of the pore (m3), S is the surface area that borders the pore (m2), and αDL is an empirical multiplying factor. This relationship has, however, severe limitations. The value of fDL is linearly related to the Debye length for a limited range of conditions only. For example, if αDL=2, a strict equality between the dual continuum model and the Poisson–Boltzmann equation predictions is observed in slit-shaped pores that are filled with NaCl electrolyte, the width of which is large compared to 4κ−1 (Sposito 2004; Tournassat and Appelo 2011; Tournassat et al. 2015a; Appelo 2017). As the width of these pores decreases, the value of fDL approaches the limiting value of one, as shown by the good agreement of the original MEP model with the PB model under these conditions (Fig. 5). Also, it is obvious that fDL cannot exceed 1. Equation (22) must then be seen as an approximation, the validity of which may be limited to small variations of ionic strength compared to the conditions at which fDL is determined experimentally. This can be appreciated by looking at the results obtained with a simple model where:

 
αDL=2 if 4κ1VporeS and,
(23)
 
fDl=1 otherwise.
(24)

In a homogeneous nanoporous medium represented by a single slit-shaped pore, these model parameters make it possible to improve estimations of ion concentration distribution compared to the MEP model: concentrations of counter-ions are calculated very accurately (less than 2 % error, except for monovalent cations in the “seawater dilution” case on Fig. 7) while the prediction of co-ion concentrations is notably improved compared to the MEP model applied to the entire pore volume (Fig. 7). Still, some conditions cannot be modeled satisfactorily, especially at low ionic strength (large values of I−1/2 on Fig. 7).

Figure 7.

Difference of average diffuse layer concentration values computed with the Poisson–Boltzmann (cPB) equation and the Dual continuum model (cDual) model described in Figure 6 and Equations (23) and (24). The difference, given by the color scale, is computed with the equation (cPB - cDual) / cPB) as a function of the inverse of the square root of ionic strength (I−1/2 on the x-axis) and as a function of the half pore width (y-axis). The equivalence between the two models corresponds to (cPB - cDual) / cPB) = 0. Same caption as in Fig. 5 otherwise.

Figure 7.

Difference of average diffuse layer concentration values computed with the Poisson–Boltzmann (cPB) equation and the Dual continuum model (cDual) model described in Figure 6 and Equations (23) and (24). The difference, given by the color scale, is computed with the equation (cPB - cDual) / cPB) as a function of the inverse of the square root of ionic strength (I−1/2 on the x-axis) and as a function of the half pore width (y-axis). The equivalence between the two models corresponds to (cPB - cDual) / cPB) = 0. Same caption as in Fig. 5 otherwise.

It is frequently argued that the Poisson–Boltzmann equation, and hence the mean electrostatic potential model or the dual continuum model, are not always accurate to model the diffuse layer composition, particularly in the presence of complex electrolytes or high surface charges (Torrie and Valleau 1982; Carnie and Torrie 1984; Luo et al. 2006; Wang and Chen 2007; Lima et al. 2008; Lee et al. 2012a). In PHREEQC, corrections can be applied by assigning an enrichment factor in the diffuse layer to each of the chemical species (Appelo et al. 2010). This correction makes it possible to fit any target ion concentration distribution in the diffuse layer, but it must be in turn calibrated with additional experimental or theoretical inputs.

The dual continuum model introduced in the present section is currently the only model available to calculate the diffuse layer composition in PHREEQC and CrunchClay, i.e., the two reactive transport codes that consider explicitly the presence of a diffuse layer in aqueous species transport calculation. Recently, the MEP model has been added to an in-house version of Flotran by Gimmi and Alt-Epping (2018) with the consideration of fixed charges in solution, as opposed to surface charge in PHREEQC and CrunchClay.

In the following section, we focus on the transport equations in the presence of a diffuse layer. Additional work is certainly needed to improve the accuracy of the prediction of average concentration in the diffuse layer. However, the type of method used for this calculation does not influence the concepts and equations used in transport equations, and in the following, the concentration ci of the species i in bulk or DL water is defined as a function of its concentration in bulk water and of an accumulation (or depletion) factor Ai, which takes the value of 1 in bulk water:

 
ci=ci,0Ai
(25)

and with the following charge balance constraints, which is equivalent to Equation (18):

 
ΣiziFci,pore¯=F(1fDL)Σizici,0+FfDLΣizici.0Ai,DL=FfDLΣizici.0Ai,DL=QDL
(26)

This notation makes it possible to treat the transport equations without a priori choice of a diffuse layer model.

COUPLED TRANSPORT PROCESSES IN REACTIVE TRANSPORT CODES—THE ISOTHERMAL NO-FLOW CONDITION

Nernst–Planck equation in the isothermal, no-flow, and no external electric field condition

Tournassat and Steefel (2015) reviewed the basics of the application of the Nernst–Planck equation to the diffuse layer bordering charged surfaces in reactive transport codes in the isothermal, no-flow and no current condition. The basic hypothesis of this application is that, at each time step, the diffuse layer equilibrates infinitely rapidly with a bulk water composition in each of the grid cells that describe the system. As explained in the previous section, the bulk water present in the grid cell can be real or fictitious, depending on the microstructural meaning of the dual continuum representation (Fig. 6). The size of a grid cell, which is at least a few micrometers for the reactive transport calculations with the highest spatial resolution, is far larger than the size of a diffuse layer, which is usually less than a few nanometers. This difference of length scale justifies the assumption of rapid equilibrium between the bulk water and the diffuse layer. The equilibrium condition is written as:

 
μi,0=μi,DLμi,0=μio+RTlnγici,0+ziFψe
(27)

where μi,0 is the electro-chemical potential of species i in bulk water, μi,DL is the electro-chemical potential of species i in the diffuse layer, μio is the standard electro-chemical potential of species i (in J⋅mol−1), γi is its activity coefficient (here in L⋅mol−1), and ψe is the electrical potential. Note that ψe does not correspond to the potential in the diffuse layer, but to a diffusion potential or an external electric potential. According to the Nernst–Planck equation, in the absence of fluid flow, the flux of species i is (in one dimension):

 
Ji=uiciμix=uiciRTln(γici,0)xuiziFciψex
(28)

where ui is the mobility of species i. Equation (28) applies to the bulk water and to the diffuse layer, and introducing Ai, the accumulation factor defined in Equation (25), it follows:

 
Ji=uici,0AiRTln(γici,0)xuiziFci,0Aiψex
(29)

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

 
ΣiziJi=0=Σi(uizici,0AiRTln(γici,0)x+uizi2Fci,0Aiψex)
(30)

It follows:

 
ψex=Σiuizici,0AiRTln(γici,0)xΣiuizi2Fci,0Ai
(31)

The mobility of ion i in the porosity is:

 
ui=ϕτiD0,iRT
(32)

and Equation (29) becomes:

 
Ji=ϕτiD0,ici,0Ailn(γici,0)x+ϕτiD0,izici,0AiΣjτjD0,jzjcj,0Ajln(γjcj,0)xΣkτkD0,kzk2ck,0Ak
(33)

Equation (33) can be further rearranged to yield:

 
Ji=ϕτiD0,ici,0Ai(1τiD0,izici,0zi2ΣkτkD0,kzk2ck,0Ak)ln(γici,0)x+ΣjiϕτiτjD0,iD0,jzizjci,0cj,0AiAjln(γjcj,0)xΣkτkD0,kzk2ck,0Ak
(34)

With isothermal, no-flow, and no external electric field conditions assumed here, the only driving forces are chemical potential gradients, and the phenomenological coefficients can be obtained directly from Equation (34):

 
Li,i=ϕτiD0,ici,0Ai(1τiD0,izici,0zi2ΣkτkD0,kzk2ck,0Ak)Li,ji=ϕτiτjD0,iD0,jzizjci,0cj,0AiAjΣkτkD0,kzk2ck,0Ak
(35)

The Onsager relationships Li,j = Lj,i are verified for these phenomenological coefficients. These coefficients apply to the bulk water and to the diffuse layer water defined in the dual continuum model. These coefficients can also be used to characterize the entire pore volume by substituting an average accumulation factor Ai¯ to the bulk and diffuse layer accumulation factors Ai:

 
Ai¯=1fDL+fDLAi
(36)

This substitution requires however that average geometrical factors τi can also be defined in some fashion.

In classical isotopic tracer diffusion experiments (Molera and Eriksen 2002; Van Loon et al. 2003, 2005, 2007; Descostes et al. 2008; Wersin et al. 2008; Tachi and Yotsuji 2014; Glaus et al. 2015; Dagnelie et al. 2018), the concentration of the diffusing species i is very low compared to the constant electrolyte background concentration. Under these conditions, Equation (35) reduces to:

 
Li,iϕτiD0,ici,0AiLi,ji0
(37)

In addition, the activity coefficients are also constant, and thus Equation (33) reduces to Fick's law, augmented by the accumulation term Ai that accelerate the diffusive flux of counterions and decelerate that of co-ions. In the presence of an electrolyte background concentration gradient, i.e., a condition that is encountered in salt diffusion experiments (Shackelford 1991; Shackelford and Daniel 1991; Malusis et al. 2015), the coupling terms between the ionic species cannot be neglected, explaining thus the need to distinguish inter-diffusion (or counter-diffusion, or mutual-diffusion) coefficients from salt diffusion coefficients. The experimental observation that effective salt-diffusion coefficients in clay membrane are dependent on both the salt concentration gradient and the effective stress applied to the clay (Malusis et al. 2015) is directly traceable from Equation (35). The dependence on salt concentration gradient is related to the non-zero values of Li,j terms, while the dependence on effective stress is related to the variations of Ai values as a function of changes in pore size distribution following the compaction of the material. For diffusion processes, it is thus necessary to take into account the electrophoretic processes corresponding to the coupling between the charge imbalance in the diffuse layer, the differences in concentration gradients for each species having different charges and thus also different accumulation factors in the diffuse layer, together with the diffusion coefficient of each individual species (Eqn. 31).

Reactive transport codes such as PHREEQC and CrunchClay are able to solve Equation (33) (without the consideration of the activity coefficient gradient in CrunchClay) under transient and stationary conditions with a diffuse layer dual continuum model, making it thus possible to calculate diffusion membrane properties from species dependent self-diffusion coefficients and from the intrinsic properties of the material, i.e., diffusion pathway geometry and ion concentration distribution. In theory, these last two parameters could be calculated from the exact knowledge of microstructure and mineralogical composition. In practice, the geometrical factor values of diffuse layer and bulk water, as well as the concentration distribution between the bulk porosity and the diffuse layer in the dual continuum model, are fitted on experimental data. Reactive transport calculations have been essential to decipher and quantify the mechanisms explaining the contrasting diffusional behavior of cations, anions and neutral species in clay materials (Appelo and Wersin 2007; Appelo et al. 2010; Glaus et al. 2013, 2015; Tournassat and Steefel 2015; Tinnacher et al. 2016; Soler et al. 2019). The fact that the accumulation of cations in the diffuse layer enhances their effective diffusion coefficient according to De,i = ϕτiD0,iAi, when present at trace concentrations, is an important result for the prediction of the fate of radionuclides such as radioactive Cs+ or Sr2+ migrating from an underground radioactive waste disposal (Altmann et al. 2012). It is therefore surprising that reactive transport codes have not yet been used extensively to interpret experimental data obtained with geosynthetic clay liners used for the isolation of landfills.

Counter-intuitive results obtained with cationic tracer experiments in the presence of a gradient of salinity can also be explained and quantified using the dual-continuum model. In our opinion, the best example of such counter-intuitive results was published by Glaus et al. (2013). The authors carried out a 22Na+ diffusion experiment through a compacted Na-montmorillonite in the presence of a NaClO4 (stable isotopes) background electrolyte concentration gradient from 0.5 mol⋅L−1 to 0.1 mol⋅L−1. At time t0, after pre-equilibration of the clay with the two background electrolytes at each end of the clay sample, the two reservoirs were spiked with the same 22Na+ concentration. A simple application of Fick's Law would predict that 22Na+ diffuses into the clay and then 22Na+ concentration stabilizes rapidly at the same value in both reservoirs. The actual results evidenced a continuous increase of the 22Na+ concentration in the NaClO4 0.5 mol⋅L−1 reservoir. This result can be explained by the difference of accumulation factor A22Na+ in the porosity at both ends of the sample. The value of A22Na+ was anti-correlated with the concentration of stable Na+:

 
A22Na+=ANa+c22Na+cNa+
(38)

and the gradient in accumulation factor A22Na+ drove the seemingly “up-hill” diffusion of A22Na+ through the clay sample (Glaus et al. 2013) (Fig. 8).

Figure 8.

22Na+ “Up-hill” diffusion experiment of Glaus et al. (2013). A: experimental setup. The clay sample (1) is in contact with a low salinity reservoir (2) and high salinity reservoir (3) by liquid lines. Blue colors represent the concentration of the background electrolyte; the red bars indicate that the experiment started with equal initial concentrations of the 22Na+ tracer in the two reservoirs. B: Experimental (symbols) and modeling (lines) results obtained with CrunchClay for the experiment carried out with with a ∼ 5 mm thick clay plug and reservoir volumes of ~ 250 mL [Modified after Tournassat and Steefel (2015) Rev Mineral Geochem 80:287–330 CCBY and with permission from Glaus et al. (2013) Seeming steady-state uphill diffusion of 22Na+ in compacted montmorillonite Environmental Science & Technology 47:11522–11527, copyright 2013 American Chemical Society.]

Figure 8.

22Na+ “Up-hill” diffusion experiment of Glaus et al. (2013). A: experimental setup. The clay sample (1) is in contact with a low salinity reservoir (2) and high salinity reservoir (3) by liquid lines. Blue colors represent the concentration of the background electrolyte; the red bars indicate that the experiment started with equal initial concentrations of the 22Na+ tracer in the two reservoirs. B: Experimental (symbols) and modeling (lines) results obtained with CrunchClay for the experiment carried out with with a ∼ 5 mm thick clay plug and reservoir volumes of ~ 250 mL [Modified after Tournassat and Steefel (2015) Rev Mineral Geochem 80:287–330 CCBY and with permission from Glaus et al. (2013) Seeming steady-state uphill diffusion of 22Na+ in compacted montmorillonite Environmental Science & Technology 47:11522–11527, copyright 2013 American Chemical Society.]

Nernst–Planck equation in the isothermal and no-flow condition in the presence of an external electric field

Electro-migration experiments consist in applying an external electric potential with two electrodes at both ends of a diffusion cell in order to accelerate the migration of ionic species in very resistive porous media such as cementitious or clayey materials (Goto and Roy 1981; Maes et al. 1999; Truc et al. 2000a,b; Samson et al. 2003; Friedmann et al. 2004, 2008a; Beauwens et al. 2005; Narsilio et al. 2007; Shi et al. 2011). In cementitious materials, the diffusion coefficient of chloride is the most commonly studied parameter. The interpretation of the results obtained with electro-migration tests is not straightforward, and a vast literature has been dedicated to the modeling of these experimental tests (Krabbenhøft and Krabbenhøft 2008). When the diffusion parameters of a single species are sought, a common practice consists in applying the so-called single-species model, in which only the transport of the ion of interest is considered, and the electric field is approximated by a constant equal to the electric potential change over the sample divided by the sample length. More elaborate models rely on the solution of the Poisson–Nernst–Planck equation. According to Krabbenhøft and Krabbenhøft (2008), while the constant field approximation may be a very poor approximation, the effective diffusivities predicted by this model are often in agreement with natural diffusion tests, i.e., diffusion experiments carried out in the absence of an external electric field. The modeling of both the Cl breakthrough together and the current measured during an electro-migration experiment is however more challenging (Appelo 2017). Discrepancies between model predictions with the Poisson–Nernst–Planck equation and actual measurements were apparent, and the influence of the diffuse layer on the transport properties was put forward by various authors as a possible explanation of these discrepancies (Goto and Roy 1981; Friedmann et al. 2004, 2008a,b; Krabbenhøft and Krabbenhøft 2008). However, no modeling framework other than those based on phenomenological coefficients was available to decipher and quantify the possible mechanisms, until PHREEQC was adapted to consider the application of an external electrical field (Appelo 2017).

The “no external electric field” condition used to solve Equation (29) in reactive transport codes can be relaxed to model electro-migration experiments. In the presence of an external electrical field, an electrical flux Jc (A⋅m−2 or C⋅m−2⋅s−1) of ionic nature appears, which is related to the fluxes of aqueous species. The following method can be used to calculate the electrical flux in 1D (Appelo 2017):

 
Jc=ΣjzjFJj=Jc,d+Jc,e
(39)

Jc,d is the current contribution from the chemical activity gradient:

 
Jc,d=ΣjzjFϕτjD0,jcj,0Ajln(γjcj,0)x
(40)

and Jc,e is the current contribution from the electrical potential gradient:

 
Jc,e=Σjzj2F2ϕτjD0,jcjAjRTψex=ΣjϕτjΛm,jcjAjψex
(41)

where:

 
Λm,j=zj2F2D0,jRT
(42)

is the molar conductivity of species i (A⋅V−1⋅m2⋅mol−1). If the subscript 1 is attributed to the driving force due to the presence of an electrical potential gradient and to the electrical flux, and the other subscripts are attributed to the driving forces due to the chemical activity gradients and the conjugated ionic fluxes, it follows:

 
L1,1=ΣjϕτjΛm,jcjAj
(43)
 
Lj1,j1=ϕτjD0,jcj,0Aj
(44)
 
L1,j1=zjFϕτjD0,jcj,0Aj=Lj1,1
(45)
 
L1,1=ΣjϕτjΛm,jcjAj
(46)

As a consequence, the model proposed by Appelo (2017) follows Onsager reciprocal relationships for bulk water as well as for diffuse layer water. Equation (46) means that there is no effect of the activity gradient of chemical species i on the flux of a chemical species j other than the influence mediated by the electrical potential gradient. Appelo (2017) introduced in his model the influence of the ionic strength on the self-diffusion coefficient of ionic species to model the change of molar conductivities with concentration and composition of the solution. This effect can be included in the D0,j terms present in Equations (43) to (45).

The numerical resolution of Equation (29) in the framework of a reactive transport approach requires that the local gradient of electrical potential ψex is correctly evaluated. This evaluation should rely on the solution of the Nernst–Planck equation together with the Poisson equation (i.e., the Poisson–Nernst–Planck equation):

 
d2ψedx21εΣi(ziFci+q)=0
(47)

where q = QDL in the diffuse layer and q = 0 in bulk water. However, the order of magnitude of F/ε in Equation (47) is 1015 V⋅m⋅mol−1. Consequently, the second derivative of the potential (first term) in Equation (47) is negligible in comparison to the charge balance (second term) except for very small length scales or very large electric potentials (Truc et al. 2000b,c; Krabbenhøft and Krabbenhøft 2008). The length scale that is usually considered in reactive transport calculations, i.e., the size of a grid cell, is typically larger than a few micrometers. As a consequence, this approximation must be correct. In addition, Appelo (2017) pointed out that the precision of any numerical model for calculating chemical reactions is usually not better than 10−8 mol⋅m−3, and a calculation of the second derivative of the potential from charge imbalance could not be more precise than 1.4 × 106 V⋅m−2, a value that is usually not reached in electro-migration experiments (Appelo 2017). Therefore, the Poisson equation reduces to a local electro-neutrality condition, and there is no transient local charge build-up: Jc is thus constant throughout the sample.

Rearranging Equation (41) yields:

 
ψexJc,eΣiϕτjΛm,jcjAj=RcJc,e
(48)

where Rc is the resistivity of the nanoporous media. After discretization on a grid cell, Equation (48) becomes:

 
ψe(x+Δx)ψe(x)=ΔxRc(x)Jc,e(x)
(49)

and:

 
ψe(xn+1)ψe(x0)=Σk=0nΔxkRc(xk)Jc,e(xk)
(50)

where x0 and xn+1 correspond to the positions where the potential ψe(xn + 1) − ψe(x0) is imposed. and ΔxkRc(xk) are known quantities. Combining Equations (50) and (39), with Jc being constant throughout the sample, yields:

 
ψe(xn+1)ψe(x0)=Σk=0nΔxkRc(xk)(JcJc,d(xk))
(51)

Jc,d (xk) is a known quantity from Equation (40) and the value of Jc can thus be computed.

The implementation of these equations in PHREEQC made it possible to unravel some of the diffusion properties of cementitious materials. In particular, the role of the diffuse layer and its average charge in the establishment of the measured current and ion breakthrough during electro-migration experiments was quantified, demonstrating that reactive transport can be used as an essential tool to understand and quantify highly coupled transport processes in nanoporous media (Appelo 2017).

BEYOND DIFFUSION, COUPLINGS WITH ADVECTIVE FLOW

Advective displacement method and reactive transport calculations

The “advective displacement” method consists in forcing flow through a confined sample by applying an hydraulic gradient in order to collect time series of small aliquots of displaced pore water (Fig. 9). It is usually assumed that early extracts closely resemble the in situ composition of the pore water. In addition, injection of tracers and the record of their breakthrough, or the elution of elements present initially in the sample, make it possible to derive a range of parameters for transport, adsorption/desorption, and dissolution/precipitation reactions using a reactive transport approach (Grambow et al. 2014; Montavon et al. 2014; Mäder and Waber 2017; Mäder 2018). Reactive transport modeling of advective displacement experiments have been carried out in a handful of studies with the consideration of the presence of a diffuse layer and the dual continuum model described in a previous section (Grambow et al. 2014; Alt-Epping et al. 2015, 2018). In all of these studies, a simplifying approximation was made: the porosity affected by the advective flow was set to the same exact value as the initial bulk porosity defined with the dual continuum model (Fig. 10A). While a significant portion of the diffuse layer porosity may remain indeed unaffected by the advective flow, it is unlikely that this portion is always equal to that of the diffuse layer porosity, which must be seen itself as a convenient simplified representation of the sample microstructure (see previous section and Fig. 6). Because the water density and viscosity in the diffuse layer is similar to that of bulk water, an increase of the diffuse layer porosity from ϕDL,1 to ϕDL,2, following e.g., a decrease of ionic strength or an increase of the ratio of monovalent counter-ions relative to multi-valent counter-ions in the bulk water, should lead to an advective displacement of a fraction of diffuse layer porosity equivalent to ϕDL,2 – ϕDL,1. This displacement should result in the appearance of coupled flux because of charge imbalance (Fig. 10B). Yet, modeling results obtained with a fixed value of ϕDL led to satisfying agreement between experimental data and modeling predictions (Fig. 11). Several explanations can be put forward to explain this good agreement including: the change of the input solution may have been too limited to change significantly the size of the diffuse layer; or the concentrations of the elements of interest were mostly controlled by reaction processes (dissolution/precipitation, cation exchange in the absence of formation of a diffuse layer); or the fitted transport parameters may also lump together several coupled processes. In this last case, the parameters must be seen as being apparent rather than intrinsic to the material studied.

Figure 9.

Advective displacement experimental setup. Reproduced from Mäder (2018) Advective displacement method for the characterisation of pore water chemistry and transport properties in claystone. Geofluids ID 8198762 CCBY.

Figure 9.

Advective displacement experimental setup. Reproduced from Mäder (2018) Advective displacement method for the characterisation of pore water chemistry and transport properties in claystone. Geofluids ID 8198762 CCBY.

Figure 10.

Advective flow in the presence of a diffuse layer and related problematics of charge balance. A: The advective flow is limited to the bulk water porosity. An electro-neutral solution displaces an electro-neutral solution and the charge balance is maintained throughout the sample. B: The advective flow displaces the bulk water porosity as well as a portion of the diffuse layer porosity. An electro-neutral solution displaces a charged solution: coupled processes must occur to maintain the local charge balance in the sample porosity and in the output solution.

Figure 10.

Advective flow in the presence of a diffuse layer and related problematics of charge balance. A: The advective flow is limited to the bulk water porosity. An electro-neutral solution displaces an electro-neutral solution and the charge balance is maintained throughout the sample. B: The advective flow displaces the bulk water porosity as well as a portion of the diffuse layer porosity. An electro-neutral solution displaces a charged solution: coupled processes must occur to maintain the local charge balance in the sample porosity and in the output solution.

Figure 11.

Example of experimental and modeling results obtained with an advective displacement experiment using a Callovian–Oxfordian argillite sample from the underground research laboratory in Bure (Meuse-Haute Marne, France). The type of model used in the study corresponds to the model shown in Fig. 10A. [Reproduced from Grambow et al. (2014) Nuclear waste disposal: I Laboratory simulation of repository properties Applied Geochemistry 49:237–246, Copyright 2014, with permission from Elsevier.]

Figure 11.

Example of experimental and modeling results obtained with an advective displacement experiment using a Callovian–Oxfordian argillite sample from the underground research laboratory in Bure (Meuse-Haute Marne, France). The type of model used in the study corresponds to the model shown in Fig. 10A. [Reproduced from Grambow et al. (2014) Nuclear waste disposal: I Laboratory simulation of repository properties Applied Geochemistry 49:237–246, Copyright 2014, with permission from Elsevier.]

Advective displacement studies can be seen as the adaptation to intact clayey rocks of earlier experiments carried out on compacted bentonite or montmorillonite samples (McKelvey and Milne 1960; Kemper and Maasland 1964; Milne et al. 1964; Hanshaw and Coplen 1973). Ultrafiltration processes were clearly evidenced with steady state output solutions having a lower electrolyte concentration than input solutions. It is thus not possible to model these experiments with the model shown in Fig. 10A.

Osmotic membrane efficiency

The influence of the diffuse layer on water and ion water transport processes in nanoporous media is most often taken into account implicitly using empirical phenomenological coefficients (Soler 2001; Malusis et al. 2012; Neuzil and Person 2017). As shown in a previous section, this modeling approach is restricted to systems in which a simple electrolyte such as NaCl is considered. In this approach, a key parameter is the coefficient of osmotic efficiency ω (–), which quantifies exclusion of ions from the total pore volume because of the presence of a diffuse layer. In the absence of thermal or electrical potential gradients, the ion flux of a single salt is modeled as the result of four contributions (Soler 2001; Sun et al. 2016; Neuzil and Person 2017):

 
Ji=Jadv,i+Jdiff,i+Jco,i+JHYP,i
(52)

where JCO,i and JHYP,i are the chemical osmosis and hyperfiltration fluxes. Chemical osmosis corresponds to the viscous flow of water down the water activity gradient (or up the osmotic pressure gradient), which drags solutes along. The osmotic pressure of an aqueous solution, π, may be expressed in the Raoult-Lewis equation (Lewis 1908):

 
παπ22=RTVolnaw
(53)

where α is the compressibility coefficient of water, Vo is its molecular volume, and aw is its activity. Neglecting the water compressibility, the following equation is obtained for the osmotic pressure gradient:

 
πx=RTΣjcjx
(54)

For aqueous solutions with molar concentration below 1 mol⋅L−1, the osmotic pressure can be approximated with van't Hoff's Law:

 
π=RTΣjcj
(55)

and the osmotic pressure gradient is:

 
πx=RTΣjcjx
(56)

The corresponding chemical osmosis flux is:

 
Jco,i=ciωKΠx
(57)

where Π = π/ρg is the osmotic pressure head (in m; ρ is the density of water and g the acceleration due to gravityd in m⋅s−2). The hyperfiltration flux corresponds to a correction to the advective and chemical osmosis flux due to the exclusion of solutes from the diffuse layer:

 
JHYP,i=ciωK(hxωΠx)
(58)

The diffusion coefficient is usually taken at the same value for anions and cations (salt diffusion coefficient) and corrected from the membrane efficiency coefficient. It follows:

 
Ji=Di(1ω)cix+(1ω)ci(Khx+ωKΠx)
(59)

A value of ω = 1 means that the membrane is ideal and there is no ionic flux through it. A value of ω = 0 means that the media has no membrane properties.

Equation (59) or similar equations have been widely used in the literature to characterize membrane properties of clay materials (Soler 2001; Malusis and Shackelford 2002, 2004; Malusis et al. 2003, 2012; Neuzil and Person 2017). The coefficient of membrane efficiency is however dependent on the electrolyte nature and concentration, and as noted by Malusis and Shackelford (2012), the development of concentration gradients within the membrane preclude the calculation of intrinsic w values from experimental data using Equation (59).

Nernst–Planck equation in the isothermal, no external electric field, no current condition, and in the presence of advective flow

In the presence of advective flow, the Nernst–Planck equation becomes:

 
Ji=ϕτiD0,ici,0Ailn(γici,0)xziFϕτiD0,ici,0AiRTψex+ciAiϕν
(60)

where v again is the average pore fluid linear velocity (m⋅s−1). The no current condition imposes:

 
ΣiziJi=0=ΣiziϕτiD0,ici,0Ailn(γici,0)xΣizi2FϕτiD0,ici,0AiRTψex+Σizici,0Aiϕν
(61)

By rearranging Equation (61), it is possible to express the value of ψex as a function of chemical potential gradients, average linear pore fluid velocity, accumulation factors and concentrations in a multi-component system:

 
ψex=QDLFνΣiziτiD0,ici,0Ailn(γici,0)xΣizi2FτiD0,ici,0AiRT
(62)

The first term of the numerator of the right part in Equation (62) corresponds to a streaming potential, while the second term of the right part corresponds to a diffusion potential. Equation (62) is then reinserted into Equation (60):

 
Ji=ϕτiD0,ici,0Ailn(γici,0)xziϕτiD0,ici,0Ai(QDLFνΣjzjτjD0,jcj,0Ajln(γjcj,0)xΣkzk2τkD0,kck,0Ak)+ci,0Aiϕν
(63)

Equation (63) can be further rearranged to yield:

 
Ji=ϕτiD0,ici,0Ai(1zi2τiD0,ici,0Σkzk2τkD0,kck,0Ak)ln(γici,0)x+ΣjiϕτiτjD0,iD0,jzizjci,0cj,0AiAjln(γjcj,0)xΣkzk2τkD0,kck,0Ak+ci,0ϕ(Ai+ziτiD0,iAiQDLFΣkzk2τkD0,kck,0Ak)ν
(64)

which is a generalization of the single salt equation given in Revil and Leroy (2004). The species dependent coefficient of filtration efficiency, ωi becomes:

 
ωi=1AiziτiD0,iAiQDLFΣkzk2τkD0,kck,0Ak
(65)

Equation (65) makes clear that, in multicomponent systems, the coefficient of osmotic efficiency is specific to each chemical species. As expected, Equation (65) indicates that the nanoporous medium is a perfect membrane with regards to the species ii = 1) if Ai = 0, i.e., if the species i is completely excluded from the porosity. Equation (65) also indicates that the coefficient of osmotic efficiency is equal to 0 if Ai = 1 for all species, because ΣiziFci,0Ai=QDL. In this last case, the porous media has no membrane properties. Equation (64) has not been implemented in a reactive transport code yet. To do this, it would be necessary to relate the linear fluid velocity to hydrostatic and osmotic pressure terms in a multi-component model framework.

From dual continuum to multi-continuum

If the entire porosity is affected by the advective fluid flow, the coefficients of osmotic efficiency can be calculated on the basis of the average pore concentrations computed in the previous section. Similarly, the osmotic pressure gradient can be calculated with Equation (56) combined with one of the models used to estimate the concentration in a charged pore. However, the concentration used for calculating the osmotic pressure may be different from the average concentration in the pore. In smectite osmotic swelling models, the reference concentration for this calculation is usually taken as the mid-pore concentration (Langmuir 1938; Warkentin and Schofield 1958; Blackmore and Miller 1961; Madsen and Müller-Vonmoos 1989; Liu 2013; Massat et al. 2016), i.e., a parameter the value of which cannot be easily computed from a mean electrostatic potential modeling approach. It would be surprising in any case if the advective flow affects uniformly the entire porosity, especially in clayey rocks and cementitious materials that exhibit a large range of pore sizes from nanometer size (interlayer porosity) to few micrometers (macropores) (Curtis et al. 2012; Kuila and Prasad 2013; Hemes et al. 2015, 2016; Gaboreau et al. 2016). It may be thus necessary to subdivide the porosity into different domains, in which specific surface charge and fluid flow are distributed heterogeneously. Multiple-porosity model for fluid flow in shale reservoirs are currently developed to capture the effect of the complexity of the shales pore structure on fluid flow (Yan et al. 2016). The consideration of ionic transport in modeling work adds one more level of complexity because of the presence of charged surfaces that make it necessary to adopt a multi-continuum approach. Even in the absence of fluid flow, a multi-continuum instead of a dual-continuum description may be necessary to model the diffusion behavior of some species such as Cs+, where interlayer diffusion may be important compared to diffusion in bulk and diffuse layer water (Appelo et al. 2010).

CONCLUSIONS

The current implementations in reactive transport codes of coupled transport processes in the presence of a diffuse layer are in agreement with the theory of irreversible thermodynamics processes. These implementations available in only two codes—CrunchClay and PHREEQC—are currently limited to the modeling of purely diffusive systems, i.e., without advective flow where there is a need to consider flow with that portion of the pore structure affected by the diffuse layer. It is also possible to model electro-migration in PHREEQC. The multi-component capability of these codes makes it possible to overcome many limitations encountered with simplified approaches that are limited to single salt transport. In particular, the transport of electrolyte background ions and tracer ions can be simulated simultaneously, demonstrating the existence of highly coupled processes during their transport through nanoporous media. In addition, reactive transport codes are not limited to the simulation of systems with fixed surface charge and fixed diffuse layer dimensions, and can incorporate microstructural information in dual- to multi-continuum descriptions of the porosity. Applications of these capabilities to the simulation of clayey and cementitious material properties highlighted the prominent role of the diffuse layer in the transport properties of ions. This review provides some insights to advance reactive transport modeling in the direction of a full implementation of multi-component advective flux in the presence of a diffuse layer. In particular, the osmotic efficiency parameter must be revisited so that ion and neutral molecule specific properties can be simulated simultaneously. Holistic models of flow and ion transport in nanoporous charged media still face significant but ultimately energizing scientific and computational challenges.

Acknowledgments

This work was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy (BES-DOE) under Contract No. DE-AC02-05CH11231, the TRANSREAC project from the Défi NEEDS – MIPOR, and by the TelluS Program of CNRS/INSU. Carl I. Steefel acknowledges funding from L'Institut Carnot for his visit to the BRGM. We gratefully acknowledge, Thomas Gimmi and Josep Soler for their constructive review of the manuscript.

1
In the following, we will use the term “mineral” as a generic term designing true minerals (e.g., clay minerals), but also other types of solid phases (e.g., cementitious phases).

REFERENCES

Alt-Epping
P
,
Tournassat
C
,
Rasouli
P
,
Steefel
CI
,
Mayer
KU
,
Jenni
A
,
Mäder
U
,
Sengor
SS
,
Fernandez
R
(
2015
)
Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects
.
Comput Geosci
 
19
:
535
550
Alt-Epping
P
,
Gimmi
T
,
Wersin
P
,
Jenni
A
(
2018
)
Incorporating electrical double layers into reactive-transport simulations of processes in clays by using the Nernst–Planck equation: A benchmark revisited
.
Appl Geochem
 
89
:
1
10
Altmann
S
,
Tournassat
C
,
Goutelard
F
,
Parneix
J-C
,
Gimmi
T
,
Maes
N
(
2012
)
Diffusion-driven transport in clayrock formations
.
Appl Geochem
 
27
:
463
478
Andrietti
F
,
Peres
A
,
Pezzotta
R
(
1976
)
Exact solution of the unidimensional Poisson–Boltzmann equation for a 1: 2 (2: 1) electrolyte
.
Biophys J
 
16
:
1121
1124
Appelo
CAJ
(
2017
)
Solute transport solved with the Nernst–Planck equation for concrete pores with “free” water and a double layer
.
Cem Concr Res
 
101
:
102
113
Appelo
CAJ
,
Wersin
P
(
2007
)
Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide, and sodium in Opalinus clay
.
Environ Sci Technol
 
41
:
5002
5007
Appelo
CAJ
,
Vinsot
A
,
Mettler
S
,
Wechner
S
(
2008
)
Obtaining the porewater composition of a clay rock by modeling the in- and out-diffusion of anions and cations from an in-situ experiment
.
J Contam Hydrol
 
101
:
67
76
Appelo
CAJ
,
Van Loon
LR
,
Wersin
P
(
2010
)
Multicomponent diffusion of a suite of tracers (HTO, Cl, Br, I, Na, Sr, Cs) in a single sample of Opalinus clay
.
Geochim Cosmochim Acta
 
74
:
1201
1219
Babcock
KL
(
1960
)
Some characteristics of a model Donnan system
.
Soil Sci
 
90
:
245
252
Babcock
KL
(
1963
)
Theory of the chemical properties of soil colloidal systems at equilibrium
.
Hilgardia
 
34
:
417
542
Banin
A
,
Lahav
N
(
1968
)
Particle size and optical properties of montmorillonite in suspension
.
Israel J Chem
 
6
:
235
250
Bear
J
(
1972
)
Dynamics of Fluids in Porous Media
.
Courier Dover Publications
Beaucaire
C
,
Tertre
E
,
Ferrage
E
,
Grenut
B
,
Pronier
S
,
Madé
B
(
2012
)
A thermodynamic model for the prediction of pore water composition of clayey rock at 25 and 80 °C—Comparison with results from hydrothermal alteration experiments
.
Chem Geol
 
334
:
62
76
Beauwens
T
,
De Canniere
P
,
Moors
H
,
Wang
L
,
Maes
N
(
2005
)
Studying the migration behaviour of selenate in Boom Clay by electromigration
.
Eng Geol
 
77
:
285
293
Beinum
W van
,
Hofmann
A
,
Meeussen
JC
,
Kretzschmar
R
(
2005
)
Sorption kinetics of strontium in porous hydrous ferric oxide aggregates: I The Donnan diffusion model
.
J Colloid Interfac Sci
 
283
:
18
28
Bildstein
O
,
Claret
F
,
Frugier
P
(
2019
)
RTM for waste repositories
.
Rev Mineral Geochem
 
85
:
419
457
Birgersson
M
(
2017
)
A general framework for ion equilibrium calculations in compacted bentonite
.
Geochim Cosmochim Acta
 
200
:
186
–(
200
)
Birgersson
M
,
Karnland
O
(
2009
)
Ion equilibrium between montmorillonite interlayer space and an external solution-Consequences for diffusional transport
.
Geochim Cosmochim Acta
 
73
:
1908
1923
Blackmore
AV
,
Miller
RD
(
1961
)
Tactoid size and osmotic swelling in calcium montmorillonite
.
Soil Sci Soc Proc
 
25
:
169
173
Borkovec
M
,
Westall
J
(
1983
)
Solution of the Poisson–Boltzmann equation for surface excesses of ions in the diffuse layer at the oxide electrolyte interface
.
J Electroanal Chem
 
150
:
325
337
Bourg
IC
(
2015
)
Sealing shales versus brittle shales: A sharp threshold in the material properties and energy, technology uses of fine-grained sedimentary rocks
.
Environ Sci Technol Lett
 
2
:
255
259
Bourg
IC
,
Sposito
G
(
2011
)
Molecular dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl–CaCl2) solutions
.
J Colloid Interfac Sci
 
360
:
701
715
Bourg
IC
,
Steefel
CI
(
2012
)
Molecular dynamics simulations of water structure and diffusion in silica nanopores
.
J Phys Chem C
 
116
:
11556
11564
Bourg
IC
,
Lee
SS
,
Fenter
P
,
Tournassat
C
(
2017
)
Stern layer structure and energetics at mica–water interfaces
.
J Phys Chem C
 
121
:
9402
9412
Bresler
E
(
1973
)
Anion exclusion and coupling effects in nonsteady transport through unsaturated soils: I Theory
.
Soil Sci Soc Am J
 
37
:
663
669
Brigatti
MF
,
Galán
E
,
Theng
BKG
(
2013
)
Chapter 2—Structure and mineralogy of clay minerals
. In:
Bergaya
F
,
Lagaly
G
(Eds.),
Handbook of Clay Science, Developments in Clay Science
 .
Elsevier
, p
21
81
Cama
J
,
Soler
JM
,
Ayora
C
(
2019
)
Acid water–rock–cement interaction and multicomponent reactive transport modeling
.
Rev Mineral Geochem
 
85
:
459
498
Carnie
SL
,
Torrie
GM
(
1984
)
Advance in chemical physics
. In:
Prigogine
I A RS
(Eds.),
Advances in Chemical Physics
 .
John Wiley & Sons, Inc.
, p
141
253
Cey
BD
,
Barbour
S
,
Hendry
MJ
(
2001
)
Osmotic flow through a Cretaceous clay in southern Saskatchewan, Canada
.
Can Geotech J
 
38
:
1025
1033
Chen
Z
,
Singh
RK
(
2002
)
General solution for Poisson–Boltzmann equation in semiinfinite planar symmetry
.
J Colloid Interfac Sci
 
245
:
301
306
Claret
F
,
Marty
N
,
Tournassat
C
(
2018
)
Modeling the long-term stability of multi-barrier systems for nuclear waste disposal in geological clay formations
. In:
Xiao
Y
,
Whitaker
F
,
Xu
T
,
Steefel
C
(Eds.),
Reactive Transport Modeling: Applications in Subsurface Energy, and Environmental Problems
 .
John Wiley & Sons, Ltd
Chichester, UK
, p
395
451
Collin
M
,
Fournier
M
,
Frugier
P
,
Charpentier
T
,
Moskura
M
,
Deng
L
,
Ren
M
,
Du
J
,
Gin
S
(
2018a
)
Structure of international simple glass and properties of passivating layer formed in circumneutral pH conditions
.
npj Materials Degradation
 
2
:
4
Collin
M
,
Gin
S
,
Dazas
B
,
Mahadevan
T
,
Du
J
,
Bourg
IC
(
2018b
)
Molecular dynamics simulations of water structure and diffusion in a 1 nm diameter silica nanopore as a function of surface charge and alkali metal counterion identity
.
J Phys Chem C
 
122
:
17764
17776
Curtis
ME
,
Sondergeld
CH
,
Ambrose
RJ
,
Rai
CS
(
2012
)
Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging
.
AAPG Bull
 
96
:
665
677
Dagnelie
R
,
Rasamimanana
S
,
Blin
V
,
Radwan
J
,
Thory
E
,
Robinet
J-C
,
Lefèvre
G
(
2018
)
Diffusion of organic anions in clay-rich media: Retardation and effect of anion exclusion
.
Chemosphere
 
213
:
472
480
Darcy
H
(
1856
)
Note II Observation relative à l'écoulement de l'eau dans l'aqueduc du Rosoir
. In:
Les fontaines publiques de la ville de Dijon
 .
Victor Dalmont
, p
638
639
de Groot
SR
,
Mazur
P
(
1984
)
Non-equilibrium Thermodynamics. Courier Corporation de Marsily G (1986) Quantitative Hydrogeology, Groundwater Hydrology for Engineers
.
Academic Press
,
New York
Delville
A
(
2000
)
Electrostatic interparticle forces: from swelling to setting
. In:
PRO 13: 2nd International RILEM Symposium on Hydration and Setting—Why does cement set? An interdisciplinary approach
.
RILEM Publications
, p
37
Descostes
M
,
Blin
V
,
Bazer-Bachi
F
,
Meier
P
,
Grenut
B
,
Radwan
J
,
Schlegel
ML
,
Buschaert
S
,
Coelho
D
,
Tevissen
E
(
2008
)
Diffusion of anionic species in Callovo-Oxfordian argillites and Oxfordian limestones (Meuse/Haute-Marne, France)
.
Appl Geochem
 
23
:
655
677
Fernández
AM
,
Sánchez-Ledesma
DM
,
Tournassat
C
,
Melón
A
,
Gaucher
EC
,
Astudillo
J
,
Vinsot
A
(
2014
)
Applying the squeezing technique to highly consolidated clayrocks for pore water characterisation: Lessons learned from experiments at the Mont Terri Rock Laboratory
.
Appl Geochem
 
49
:
2
21
Fick
A
(
1855
)
V On liquid diffusion
.
The London, Edinburgh Dublin Phil Mag J Sci
 
10
:
30
39
Friedmann
H
,
Amiri
O
,
Ait-Mokhtar
A
,
Dumargue
P
(
2004
)
A direct method for determining chloride diffusion coefficient by using migration test
.
Cem Concr Res
 
34
:
1967
1973
Friedmann
H
,
Amiri
O
,
Ait-Mokhtar
A
(
2008a
)
Shortcomings of geometrical approach in multi-species modelling of chloride migration in cement-based materials
.
Mag Concr Res
 
60
:
119
124
Friedmann
H
,
Amiri
O
,
Ait-Mokhtar
A
(
2008b
)
Physical modeling of the electrical double layer effects on multispecies ions transport in cement-based materials
.
Cem Concr Res
 
38
:
1394
1400
Frugier
P
,
Minet
Y
,
Rajmohan
N
,
Godon
N
,
Gin
S
(
2018
)
Modeling glass corrosion with GRAAL
.
npj Mater Degrad
 
2
:
35
Gaboreau
S
,
Robinet
J-C
,
Prêt
D
(
2016
)
Optimization of pore network characterization of compacted clay materials by TEM and FIB/SEM imaging
.
Micropor Mesopor Mater
 
224
:
116
128
Gailhanou
H
,
Lerouge
C
,
Debure
M
,
Gaboreau
S
,
Gaucher
EC
,
Grangeon
S
,
Grenèche
J-M
,
Kars
M
,
Madé
B
,
Marty
NCM
,
Warmont
F
,
Tournassat
C
(
2017
)
Effects of a thermal perturbation on mineralogy and pore water composition in a clay-rock: an experimental and modeling study
.
Geochim Cosmochim Acta
 
197
:
193
214
Garavito
AM
,
De Cannière
P
,
Kooi
H
(
2007
)
In situ chemical osmosis experiment in the Boom Clay at the Mol underground research laboratory
.
Phys Chem Earth, Parts A/B/C
 
32
:
421
433
Gates
WP
,
Bouazza
A
,
Churchman
GJ
(
2009
)
Bentonite clay keeps pollutants at bay
.
Elements
 
5
:
105
110
Gaucher
E
,
Robelin
C
,
Matray
J-M
,
Negrel
G
,
Gros
Y
,
Heitz
JF
,
Vinsot
A
,
Rebours
H
,
Cassabagnere
A
,
Bouchet
A
(
2004
)
ANDRA underground research laboratory: Interpretation of the mineralogical and geochemical data acquired in the Callovian–Oxfordian Formation by investigative drilling
.
Phys Chem Earth, Parts A/B/C
 
29
:
55
77
Gaucher
EC
,
Tournassat
C
,
Pearson
FJ
,
Blanc
P
,
Crouzet
C
,
Lerouge
C
,
Altmann
S
(
2009
)
A robust model for pore-water chemistry of clayrock
.
Geochim Cosmochim Acta
 
73
:
6470
6487
Gimmi
T
,
Alt-Epping
P
(
2018
)
Simulating Donnan equilibria based on the Nernst–Planck equation
.
Geochim Cosmochim Acta
 
232
:
1
13
Glaus
MA
,
Birgersson
M
,
Karnland
O
,
Van Loon
LR
(
2013
)
Seeming steady-state uphill diffusion of 22Na+ in compacted montmorillonite
.
Environ Sci Technol
 
47
:
11522
11527
Glaus
M
,
Aertsens
M
,
Appelo
C
,
Kupcik
T
,
Maes
N
,
Van Laer
L
,
Van Loon
L
(
2015
)
Cation diffusion in the electrical double layer enhances the mass transfer rates for Sr2+, Co2+ and Zn2+ in compacted illite
.
Geochim Cosmochim Acta
 
165
:
376
388
Gonçalvès
J
,
Rousseau-Gueutin
P
,
Revil
A
(
2007
)
Introducing interacting diffuse layers in TLM calculations: A reappraisal of the influence of the pore size on the swelling pressure and the osmotic efficiency of compacted bentonites
.
J Colloid Interfac Sci
 
316
:
92
99
Gonçalvès
J
,
Adler
PM
,
Cosenza
P
,
Pazdniakou
A
,
de Marsily
G
(
2015
)
Chapter 8—Semipermeable membrane properties and chemomechanical coupling in clay barriers
. In:
Tournassat
C
,
Steefel
CI
,
Bourg
IC
,
Bergaya
F
(Eds.),
Natural and Engineered Clay Barriers, Developments in Clay Science
 .
Elsevier
, p
270
327
Goto
S
,
Roy
DM
(
1981
)
Diffusion of ions through hardened cement pastes
.
Cem Concr Res
 
11
:
751
757
Grambow
B
(
2006
)
Nuclear waste glasses-How durable?
Elements
 
2
:
357
364
Grambow
B
,
Landesman
C
,
Ribet
S
(
2014
)
Nuclear waste disposal: I Laboratory simulation of repository properties
.
Appl Geochem
 
49
:
237
246
Grangeon
S
,
Claret
F
,
Lerouge
C
,
Warmont
F
,
Sato
T
,
Anraku
S
,
Numako
C
,
Linard
Y
,
Lanson
B
(
2013
)
On the nature of structural disorder in calcium silicate hydrates with a calcium/silicon ratio similar to tobermorite
.
Cem Concr Res
 
52
:
31
37
Grangeon
S
,
Claret
F
,
Roosz
C
,
Sato
T
,
Gaboreau
S
,
Linard
Y
(
2016
)
Structure of nanocrystalline calcium silicate hydrates: insights from X-ray diffraction, synchrotron X-ray absorption and nuclear magnetic resonance
.
J Appl Crystallogr
 
49
:
771
783
Grangeon
S
,
Fernandez-Martinez
A
,
Baronnet
A
,
Marty
N
,
Poulain
A
,
Elkaim
E
,
Roosz
C
,
Gaboreau
S
,
Henocq
P
,
Claret
F
(
2017
)
Quantitative X-ray pair distribution function analysis of nanocrystalline calcium silicate hydrates: a contribution to the understanding of cement chemistry
.
J Appl Crystallogr
 
50
:
14
21
Gu
X
,
Cole
DR
,
Rother
G
,
Mildner
DF
,
Brantley
SL
(
2015
)
Pores in Marcellus shale: a neutron scattering and FIB SEM study
.
Energy Fuels
 
29
:
1295
1308
Guyonnet
D
,
Touze-Foltz
N
,
Norotte
V
,
Pothier
C
,
Didier
G
,
Gailhanou
H
,
Blanc
P
,
Warmont
F
(
2009
)
Performance-based indicators for controlling geosynthetic clay liners in landfill applications
.
Geotextiles Geomembranes
 
27
:
321
331
Haas
J
,
Nonat
A
(
2015
)
From C–S–H to C–A–S–H: Experimental study and thermodynamic modelling
.
Cem Concr Res
 
68
:
124
138
Hanshaw
BB
,
Coplen
TB
(
1973
)
Ultrafiltration by a compacted clay membrane—II Sodium ion exclusion at various ionic strengths
.
Geochim Cosmochim Acta
 
37
:
2311
2327
Hemes
S
,
Desbois
G
,
Urai
JL
,
Schröppel
B
,
Schwarz
J-O
(
2015
)
Multi-scale characterization of porosity in Boom Clay (HADES-level, Mol, Belgium) using a combination of X-ray µ-CT, 2D BIB-SEM and FIB-SEM tomography
.
Micropor Mesopor Mater
 
208
:
1
20
Hemes
S
,
Desbois
G
,
Klaver
J
,
Urai
JL
(
2016
)
Microstructural characterisation of the Ypresian clays (Kallo−1) at nanometre resolution, using broad-ion beam milling and scanning electron microscopy
.
Netherlands J Geosci
 
95
:
293
313
Holland
HD
(
1978
)
The Chemistry of the Atmosphere and Oceans
.
Wiley-interscience
,
New York
Holmboe
M
,
Bourg
IC
(
2014
)
Molecular dynamics simulations of water and sodium diffusion in smectite interlayer nanopores as a function of pore size and temperature
.
J Phys Chem C
 
118
:
1001
1013
Horseman
S
,
Harrington
J
,
Noy
D
(
2007
)
Swelling and osmotic flow in a potential host rock
.
Phys Chem Earth, Parts A/B/C
 
32
:
408
420
Hunter
RJ
(
2013
)
Zeta Potential in Colloid Science: Principles and Applications
.
Academic Press
Ioannidou
K
,
Krakowiak
KJ
,
Bauchy
M
,
Hoover
CG
,
Masoero
E
,
Yip
S
,
Ulm
F-J
,
Levitz
P
,
Pellenq
RJ-M
,
Del Gado
E
(
2016
)
Mesoscale texture of cement hydrates
.
PNAS
 
113
:
2029
–(
2034
)
Jardat
M
,
Dufreche
JF
,
Marry
V
,
Rotenberg
B
,
Turq
P
(
2009
)
Salt exclusion in charged porous media: a coarse-graining strategy in the case of montmorillonite clays
.
Phys Chem Chem Phys
 
11
:
2023
–(
2033
)
Jougnot
D
,
Revil
A
,
Leroy
P
(
2009
)
Diffusion of ionic tracers in the Callovo–Oxfordian clay-rock using the Donnan equilibrium model and the formation factor
.
Geochim Cosmochim Acta
 
73
:
2712
2726
Kemper
WD
(
1960
)
Water and ion movement in thin films as influenced by the electrostatic charge and diffuse layer of cations associated with clay mineral surfaces
.
Soil Sci Soc Proc
 
10
16
Kemper
W
(
1961a
)
Movement of water as effected by free energy and pressure gradients: II Experimental analysis of porous systems in which free energy and pressure gradients act in opposite directions
.
Soil Sci Soc Am J
 
25
:
260
265
Kemper
W
(
1961b
)
Movement of water as effected by free energy and pressure gradients: I Application of classic equations for viscous and diffusive movements to the liquid phase in finely porous media
.
Soil Sci Soc Am J
 
25
:
255
260
Kemper
W
,
Evans
N
(
1963
)
Movement of water as effected by free energy and pressure gradients III Restriction of solutes by membranes
.
Soil Sci Soc Am J
 
27
:
485
490
Kemper
WD
,
Maasland
DEL
(
1964
)
Reduction in salt content of solution on passing through thin films adjacent to charged surfaces 1
.
Soil Sci Soc Am J
 
28
:
318
323
Kemper
WD
,
Rollins
JB
(
1966
)
Osmotic efficiency coefficients across compacted clays
.
Soil Sci Soc Am J
 
30
:
529
534
Keller
LM
,
Holzer
L
,
Wepf
R
,
Gasser
P
(
2011
)
3D geometry and topology of pore pathways in Opalinus clay: Implications for mass transport
.
Appl Clay Sci
 
52
:
85
95
Keller
LM
,
Holzer
L
,
Schuetz
P
,
Gasser
P
(
2013
)
Pore space relevant for gas permeability in Opalinus clay: Statistical analysis of homogeneity, percolation, and representative volume element
.
J Geophys Res
 
118
:
2799
2812
Kjelstrup
S
,
Bedeaux
D
(
2008
)
Non-equilibrium Thermodynamics of Heterogeneous Systems
.
World Scientific Krabbenhøft
K
,
Krabbenhøft
J
(
2008
)
Application of the Poisson–Nernst–Planck equations to the migration test
.
Cem Concr Res
 
38
:
77
88
Kuila
U
,
Prasad
M
(
2013
)
Specific surface area and pore-size distribution in clays and shales
.
Geophys Prospect
 
61
:
341
362
Lagneau
V
,
Regnault
O
,
Descostes
M
(
2019
)
Industrial deployment of reactive transport simulation: An application to uranium in situ recovery
.
Rev Mineral Geochem
 
85
:
499
528
Langmuir
I
(
1938
)
The role of attractive and repulsive forces in the formation of tactoids, thixotropic gels, protein crystals and coacervates
.
J Chem Phys
 
6
:
873
896
Lasaga
AC
(
1998
)
Kinetic Theory in the Earth Sciences
.
Princeton University Press
Le Forestier
L
,
Muller
F
,
Villiéras
F
,
Pelletier
M
(
2010
)
Textural and hydration properties of a synthetic montmorillonite compared with a natural Na-exchanged clay analogue
.
Appl Clay Sci
 
48
:
18
25
Lee
SS
,
Fenter
P
,
Park
C
,
Sturchio
NC
,
Nagy
KL
(
2010
)
Hydrated cation speciation at the muscovite (001)–water interface
.
Langmuir
 
26
:
16647
16651
Lee
JW
,
Nilson
RH
,
Templeton
JA
,
Griffiths
SK
,
Kung
A
,
Wong
BM
(
2012a
)
Comparison of molecular dynamics with classical density functional and Poisson–Boltzmann Theories of the electric double layer in nanochannels
.
J Chem Theory Comput
 
8
:
2012
2022
Lee
SS
,
Fenter
P
,
Nagy
KL
,
Sturchio
NC
(
2012b
)
Monovalent ion adsorption at the muscovite (001)–solution interface: Relationships among ion coverage and speciation, interfacial water structure, and substrate relaxation
.
Langmuir
 
28
:
8637
8650
Leroy
P
,
Maineult
A
(
2018
)
Exploring the electrical potential inside cylinders beyond the Debye–Hückel approximation: a computer code to solve the Poisson–Boltzmann equation for multivalent electrolytes
.
Geophys J Int
 
214
:
58
69
Leroy
P
,
Revil
A
(
2004
)
A triple-layer model of the surface electrochemical properties of clay minerals
.
J Colloid Interfac Sci
 
270
:
371
380
Leroy
P
,
Revil
A
,
Coelho
D
(
2006
)
Diffusion of ionic species in bentonite
.
J Colloid Interfac Sci
 
296
:
248
255
Leroy
P
,
Revil
A
,
Altmann
S
,
Tournassat
C
(
2007
)
Modeling the composition of the pore water in a clay-rock geological formation (Callovo-Oxfordian, France)
.
Geochim Cosmochim Acta
 
71
:
1087
1097
Lewis
GN
(
1908
)
The osmotic pressure of concentrated solutions, and the laws of the perfect solution
.
J Am Chem Soc
 
30
:
668
683
Li
Y-H
,
Gregory
S
(
1974
)
Diffusion of ions in sea water and in deep-sea sediments
.
Geochim Cosmochim Acta
 
38
:
703
714
Lima
E
,
Horinek
D
,
Netz
R
,
Biscaia
E
,
Tavares
F
,
Kunz
W
,
Boström
M
(
2008
)
Specific ion adsorption and surface forces in colloid science
.
J Phys Chem B
 
112
:
1580
1585
Liu
L
(
2013
)
Prediction of swelling pressures of different types of bentonite in dilute solutions
.
Colloids Surf A
 
434
:
303
318
Loucaides
S
,
Behrends
T
,
Van Cappellen
P
(
2010
)
Reactivity of biogenic silica: Surface versus bulk charge density
.
Geochim Cosmochim Acta
 
74
:
517
530
Luo
G
,
Malkova
S
,
Yoon
J
,
Schultz
DG
,
Lin
B
,
Meron
M
,
Benjamin
I
,
Vansek
P
,
Schlossman
ML
(
2006
)
Ion distributions near a liquid-liquid interface
.
Science
 
311
:
216
218
Lynde
C
(
1912
)
Osmosis in soils: Soils act as semi-permeable membranes
.
Agron J
 
4
:
102
108
Ma
B
,
Fernandez-Martinez
A
,
Grangeon
S
,
Tournassat
C
,
Findling
N
,
Carrero
S
,
Tisserand
D
,
Bureau
S
,
Elkaïm
E
,
Marini
C
,
Aquilanti
G
,
Koishi
A
,
Marty
NCM
,
Charlet
L
(
2018
)
Selenite uptake by Ca–Al LDH: a description of intercalated anion coordination geometries
.
Environ Sci Technol
 
52
:
1624
1632
Ma
B
,
Fernandez-Martinez
A
,
Grangeon
S
,
Tournassat
C
,
Findling
N
,
Claret
F
,
Koishi
A
,
Marty
NC
,
Tisserand
D
,
Bureau
S
,
Salas-Colera
E
(
2017
)
Evidence of multiple sorption modes in layered double hydroxides using Mo as structural probe
.
Environ Sci Technol
 
51
:
5531
5540
Mäder
U
(
2018
)
Advective displacement method for the characterisation of pore water chemistry and transport properties in claystone
.
Geofluids
 :
8198762
Mäder
UK
,
Waber
HN
(
2017
)
Characterization of pore water, ion transport and water-rock interaction in claystone by advective displacement experiments
.
Procedia earth and planetary science
 
17
:
917
920
Madsen
FT
,
Müller-Vonmoos
M
(
1989
)
The swelling behaviour of clays
.
Appl Clay Sci
 
4
:
143
156
Maes
N
,
Moors
H
,
Dierckx
A
,
De Cannière
P
,
Put
M
(
1999
)
The assessment of electromigration as a new technique to study diffusion of radionuclides in clayey soils
.
J Contam Hydrol
 
36
:
231
247
Maher
K
,
Navarre-Stichler
A
(
2019
)
Reactive transport processes that drive chemical weathering: From making space for water to dismantling continents
.
Rev Mineral Geochem
 
85
:
349
380
Malusis
MA
,
Shackelford
CD
(
2002
)
Theory for reactive solute transport through clay membrane barriers
.
J Contam Hydrol
 
59
:
291
316
Malusis
MA
,
Shackelford
CD
(
2004
)
Predicting solute flux through a clay membrane barrier
.
J Geotech Geoenviron Eng
 
130
:
477
487
Malusis
MA
,
Shackelford
CD
,
Olsen
HW
(
2003
)
Flow and transport through clay membrane barriers
.
Eng Geol
 
70
:
235
248
Malusis
MA
,
Shackelford
CD
,
Maneval
JE
(
2012
)
Critical review of coupled flux formulations for clay membranes based on nonequilibrium thermodynamics
.
J Contam Hydrol
 
138
:
40
59
Malusis
MA
,
Kang
JB
,
Shackelford
CD
(
2015
)
Restricted salt diffusion in a geosynthetic clay liner
.
Environmental Geotechnics
 
2
:
68
Marry
V
,
Turq
P
,
Cartailler
T
,
Levesque
D
(
2002
)
Microscopic simulation for structure and dynamics of water and counterions in a monohydrated montmorillonite
.
J Chem Phys
 
117
:
3454
3463
Marry
V
,
Rotenberg
B
,
Turq
P
(
2008
)
Structure and dynamics of water at a clay surface from molecular dynamics simulation
.
Phys Chem Chem Phys
 
10
:
4802
4813
Marty
NCM
,
Cama
J
,
Sato
T
,
Chino
D
,
Villiéras
F
,
Razafitianamaharavo
A
,
Brendlé
J
,
Giffaut
E
,
Soler
JM
,
Gaucher
EC
,
Tournassat
C
(
2011
)
Dissolution kinetics of synthetic Na-smectite. An integrated experimental approach
.
Geochim Cosmochim Acta
 
75
:
5849
5864
Marty
NC
,
Grangeon
S
,
Elkaim
E
,
Tournassat
C
,
Fauchet
C
,
Claret
F
(
2018
)
Thermodynamic and crystallographic model for anion uptake by hydrated calcium aluminate (AFm): an example of molybdenum
.
Sci Rep
 
8
:
7943
Masoero
E
,
Del Gado
E
,
Pellenq
R-M
,
Ulm
F-J
,
Yip
S
(
2012
)
Nanostructure and nanomechanics of cement: polydisperse colloidal packing
.
Phys Rev Lett
 
109
:
155503
Massat
L
,
Cuisinier
O
,
Bihannic
I
,
Claret
F
,
Pelletier
M
,
Masrouri
F
,
Gaboreau
S
(
2016
)
Swelling pressure development and inter-aggregate porosity evolution upon hydration of a compacted swelling clay
.
Appl Clay Sci
 
124
:
197
210
Mazurek
M
,
Alt-Epping
P
,
Bath
A
,
Gimmi
T
,
Niklaus Waber
H
,
Buschaert
S
,
Cannière
PD
,
De Craen
M
,
Gautschi
A
,
Savoye
S
,
Vinsot
A
,
Wemaere
I
,
Wouters
L
(
2011
)
Natural tracer profiles across argillaceous formations
.
Appl Geochem
 
26
:
1035
1064
Mazurek
M
,
Oyama
T
,
Wersin
P
,
Alt-Epping
P
(
2015
)
Pore-water squeezing from indurated shales
.
Chem Geol
 
400
:
106
121
McKelvey
JG
,
Milne
IH
(
1960
)
The flow of salt solutions through compacted clay
.
Clay Clay Mineral
 
9
:
248
259
Medved
I
,
Cern
R
(
2013
)
Osmosis in porous media: A review of recent studies
.
Micropor Mesopor Mater
 
170
:
299
317
Milne
I
,
McKelvey
J
,
Trump
R
(
1964
)
Semi-permeability of bentonite membranes to brines
.
AAPG Bull
 
48
:
103
105
Molera
M
,
Eriksen
T
(
2002
)
Diffusion of 22Na+, 85Sr2+, 134Cs+ and 57Co2+ in bentonite clay compacted to different densities: experiments and modeling
.
Radiochim Acta
 
90
:
75
760
Montavon
G
,
Sabatié-Gogova
A
,
Ribet
S
,
Bailly
C
,
Bessaguet
N
,
Durce
D
,
Giffaut
E
,
Landesman
C
,
Grambow
B
(
2014
)
Retention of iodide by the Callovo-Oxfordian formation: An experimental study
.
Appl Clay Sci
 
87
:
142
149
Narsilio
G
,
Li
R
,
Pivonka
P
,
Smith
D
(
2007
)
Comparative study of methods used to estimate ionic diffusion coefficients using migration tests
.
Cem Concr Res
 
37
:
1152
1163
Navarre-Sitchler
A
,
Steefel
CI
,
Sak
PB
,
Brantley
SL
(
2011
)
A reactive-transport model for weathering rind formation on basalt
.
Geochim Cosmochim Acta
 
75
:
7644
7667
Neuzil
C
(
2000
)
Osmotic generation of “anomalous” fluid pressures in geological environments
.
Nature
 
403
:
182
Neuzil
CE
(
2013
)
Can shale safely host US nuclear waste?
Eos, Trans Am Geophys Union
 
94
:
261
262
Neuzil
CE
,
Person
M
(
2017
)
Reexamining ultrafiltration and solute transport in groundwater
.
Water Resour Res
 
53
:
4922
4941
Obliger
A
,
Jardat
M
,
Coelho
D
,
Bekri
S
,
Rotenberg
B
(
2014
)
Pore network model of electrokinetic transport through charged porous media
.
Phys Rev E
 
89
:
043013
Olsen
HW
(
1962
)
Hydraulic flow through saturated clays
.
Clays Clay Minerals
 
131
161
Onsager
L
(
1931a
)
Reciprocal relations in irreversible processes
.
I Phys Rev
 
37
:
405
Onsager
L
(
1931b
)
Reciprocal relations in irreversible processes
.
II Phys Rev
 
38
:
2265
Patriarche
D
,
Ledoux
E
,
Michelot
J-L
,
Simon-Coinçon
R
,
Savoye
S
(
2004a
)
Diffusion as the main process for mass transport in very low water content argillites: 2. Fluid flow and mass transport modeling
.
Water Resour Res
 
40
:
W01517
Patriarche
D
,
Michelot
J-L
,
Ledoux
E
,
Savoye
S
(
2004b
)
Diffusion as the main process for mass transport in very low water content argillites: 1. Chloride as a natural tracer for mass transport—Diffusion coefficient and concentration measurements in interstitial water
.
Water Resour Res
 
40
:
W01516
Pearson
FJ
,
Tournassat
C
,
Gaucher
EC
(
2011
)
Biogeochemical processes in a clay formation in situ experiment: Part E - Equilibrium controls on chemistry of pore water from the Opalinus Clay, Mont Terri Underground Research Laboratory, Switzerland
.
Appl Geochem
 
26
:
990
1008
Pellenq
RJ-M
,
Van Damme
H
(
2004
)
Why does concrete set?: The nature of cohesion forces in hardened cement-based materials
.
MRS Bull
 
29
:
319
323
Philipp
T
,
Amann-Hildenbrand
A
,
Laurich
B
,
Desbois
G
,
Littke
R
,
Urai
J
(
2017
)
The effect of microstructural heterogeneity on pore size distribution and permeability in Opalinus Clay (Mont Terri, Switzerland): insights from an integrated study of laboratory fluid flow and pore morphology from BIB-SEM images
.
Geol Soc
 ,
London
Spec Publ
454
:
85
106
Poinssot
C
,
Baeyens
B
,
Bradbury
MH
(
1999
)
Experimental and modelling studies of caesium sorption on illite
.
Geochim Cosmochim Acta
 
63
:
3217
3227
Pusch
R
(
2001
)
The microstructure of MX−80 clay with respect to its bulk physical properties under different environmental conditions
.
SKB, TR-01−08
Revil
A
(
1999
)
Ionic diffusivity, electrical conductivity, membrane and thermoelectric potentials in colloids and granular porous media: a unified model
.
J Colloid Interfac Sci
 
212
:
503
522
Revil
A
,
Leroy
P
(
2004
)
Constitutive equations for ionic transport in porous shales
.
J Geophys Res-Solid Earth
 
109
Revil
A
,
Linde
N
(
2006
)
Chemico-electromechanical coupling in microporous media
.
J Colloid Interfac Sci
 
302
:
682
694
Revil
A
,
Woodruff
W
,
Lu
N
(
2011
)
Constitutive equations for coupled flows in clay materials
.
Water Resour Res
 
47
:
1
21
Richardson
I
(
2008
)
The calcium silicate hydrates
.
Cem Concr Res
 
38
:
137
158
Richardson
IG
(
2014
)
Model structures for C-(A)-SH (I)
.
Acta Crystallogr Section B
 
70
:
903
923
Robinet
J-C
,
Sardini
P
,
Coelho
D
,
Parneix
J-C
,
Prêt
D
,
Sammartino
S
,
Boller
E
,
Altmann
S
(
2012
)
Effects of mineral distribution at mesoscopic scale on solute diffusion in a clay-rich rock: Example of the Callovo-Oxfordian mudstone (Bure, France)
.
Water Resour Res
 
48
:
W05554
Rolle
M
,
Le Borgne
T
(
2019
)
Mixing and reactive fronts in the subsurface
.
Rev Mineral Geochem
 
85
:
111
142
Roosz
C
,
Gaboreau
S
,
Grangeon
S
,
Prêt
D
,
Montouillout
V
,
Maubec
N
,
Ory
S
,
Blanc
P
,
Vieillard
P
,
Henocq
P
(
2016
)
Distribution of water in synthetic calcium silicate hydrates
.
Langmuir
 
32
:
6794
6805
Rotenberg
B
,
Marry
V
,
Dufrêche
J-F
,
Giffaut
E
,
Turq
P
(
2007
)
A multiscale approach to ion diffusion in clays: Building a two-state diffusion-reaction scheme from microscopic dynamics
.
J Colloid Interfac Sci
 
309
:
289
295
Rotenberg
B
,
Marry
V
,
Malikova
N
,
Turq
P
(
2010
)
Molecular simulation of aqueous solutions at clay surfaces
.
J Phys Condens Matter
 
22
:
284114
Rotenberg
B
,
Marry
V
,
Salanne
M
,
Jardat
M
,
Turq
P
(
2014
)
Multiscale modelling of transport in clays from the molecular to the sample scale
.
C R Geosci
 
346
:
298
306
Rousseau-Gueutin
P
,
De Greef
V
,
Gonçalvès
J
,
Violette
S
,
Chanchole
S
(
2009
)
Experimental device for chemical osmosis measurement on natural clay-rock samples maintained at in situ conditions: Implications for formation pressure interpretations
.
J Colloid Interfac Sci
 
337
:
106
116
Saiyouri
N
,
Hicher
PY
,
Tessier
D
(
2000
)
Microstructural approach and transfer water modelling in highly compacted unsaturated swelling clays
.
Mech Cohesive-frictional Mater
 
5
:
41
60
Samson
E
,
Marchand
J
,
Snyder
KA
(
2003
)
Calculation of ionic diffusion coefficients on the basis of migration test results
.
Mater Struct
 
36
:
156
165
Sayed Hassan
M
,
Villieras
F
,
Gaboriaud
F
,
Razafitianamaharavo
A
(
2006
)
AFM and low-pressure argon adsorption analysis of geometrical properties of phyllosilicates
.
J Colloid Interfac Sci
 
296
:
614
623
Schlegel
ML
,
Nagy
KL
,
Fenter
P
,
Cheng
L
,
Sturchio
NC
,
Jacobsen
SD
(
2006
)
Cation sorption on the muscovite (0 0 1) surface in chloride solutions using high-resolution X-ray reflectivity
.
Geochim Cosmochim Acta
 
70
:
3549
3565
Schramm
LL
,
Kwak
JCT
(
1982
)
Influence of exchangeable cation composition on the size and shape of montmorillonite particles in dilute suspension
.
Clay Clay Mineral
 
30
:
40
48
Shackelford
CD
(
1991
)
Laboratory diffusion testing for waste disposal – A review
.
J Contam Hydrol
 
7
:
177
217
Shackelford
CD
,
Daniel
DE
(
1991
)
Diffusion in saturated soil. I: Background
.
J Geotech Eng
 
117
:
467
484
Shainberg
I
,
Otoh
H
(
1968
)
Size and shape of montmorillonite particles saturated with Na/Ca ions (inferred from viscosity and optical measurements)
.
Israel J Chem
 
6
:
251
259
Shi
X
,
Yang
Z
,
Liu
Y
,
Cross
D
(
2011
)
Strength and corrosion properties of Portland cement mortar and concrete with mineral admixtures
.
Construct Building Mater
 
25
:
3245
3256
Siretanu
I
,
Ebeling
D
,
Andersson
MP
,
Stipp
SLS
,
Philipse
A
,
Stuart
MC
,
Ende
D
van den
,
Mugele F
(
2014
)
Direct observation of ionic structure at solid–liquid interfaces: a deep look into the Stern Layer
.
Sci Rep
 
4
:
4956
4956
Soler
JM
(
2001
)
The effect of coupled transport phenomena in the Opalinus Clay and implications for radionuclide transport
.
J Contam Hydrol
 
53
:
63
84
Soler
JM
,
Steefel
CI
,
Gimmi
T
,
Leupin
OX
,
Cloet
V
(
2019
)
Modeling the ionic strength effect on diffusion in clay. The DR-A experiment at Mont Terri
.
ACS Earth Space Chem
 
3
:
442
451
Soulaine
C
,
Tchelepi
HA
(
2016
)
Micro-continuum approach for pore-scale simulation of subsurface processes
.
Transp Porous Media
 
113
:
431
456
Soulaine
C
,
Roman
S
,
Kovscek
A
,
Tchelepi
HA
(
2018
)
Pore-scale modelling of multiphase reactive flow: application to mineral dissolution with production of CO2
.
J Fluid Mech
 
855
:
616
645
Sposito
G
(
2004
)
The Surface Chemistry of Natural Particles
.
Oxford University Press
,
New York
Steefel
CI
(
2008
)
Geochemical kinetics and transport
. In:
Kinetics of Water–Rock Interaction
 .
Springer
, p
545
589
Steefel
CI
,
Maher
K
(
2009
)
Fluid-rock interaction: A reactive transport approach
.
Rev Mineral Geochem
 
70
:
485
532
Steefel
CI
,
Appelo
CAJ
,
Arora
B
,
Jacques
D
,
Kalbacher
T
,
Kolditz
O
,
Lagneau
V
,
Lichtner
PC
,
Mayer
KU
,
Meeussen
JCL
,
Molins
S
,
Moulton
D
,
Shao
H
,
Šimunek
J
,
Spycher
N
,
Yabusaki
SB
,
Yeh
GT
(
2015a
)
Reactive transport codes for subsurface environmental simulation
.
Comput Geosci
 
19
:
445
478
Steefel
CI
,
Beckingham
LE
,
Landrot
G
(
2015b
)
Micro-continuum approaches for modeling pore-scale geochemical processes
.
Rev Mineral Geochem
 
80
:
217
246
Sun
X
,
Wu
J
,
Shi
X
,
Wu
J
(
2016
)
Experimental and numerical modeling of chemical osmosis in the clay samples of the aquitard in the North China Plain
.
Environ Earth Sci
 
75
:
59
Tachi
Y
,
Yotsuji
K
(
2014
)
Diffusion and sorption of Cs+, Na+, I and HTO in compacted sodium montmorillonite as a function of porewater salinity: Integrated sorption and diffusion model
.
Geochim Cosmochim Acta
 
132
:
75
93
Tinnacher
RM
,
Holmboe
M
,
Tournassat
C
,
Bourg
IC
,
Davis
JA
(
2016
)
Ion adsorption and diffusion in smectite: molecular, pore, and continuum scale views
.
Geochim Cosmochim Acta
 
177
:
130
149
Torrie
GM
,
Valleau
JP
(
1982
)
Electrical double layers. 4. Limitations of the Gouy–Chapman theory
.
J Phys Chem
 
86
:
3251
3257
Tournassat
C
,
Appelo
CAJ
(
2011
)
Modelling approaches for anion-exclusion in compacted Na-bentonite
.
Geochim Cosmochim Acta
 
75
:
3698
3710
Tournassat
C
,
Steefel
CI
(
2015
)
Ionic transport in nano-porous clays with consideration of electrostatic effects
.
Rev Mineral Geochem
 
80
:
287
330
Tournassat
C
,
Steefel
CI
(
2019
)
Modeling diffusion processes in the presence of a diffuse layer at charged mineral surfaces. A benchmark exercise
.
Comput Geosci
, accepted
Tournassat
C
,
Neaman
A
,
Villiéras
F
,
Bosbach
D
,
Charlet
L
(
2003
)
Nanomorphology of montmorillonite particles: Estimation of the clay edge sorption site density by low-pressure gas adsorption and AFM observations
.
Am Mineral
 
88
:
1989
1995
Tournassat
C
,
Chapron
Y
,
Leroy
P
,
Boulahya
F
(
2009a
)
Comparison of molecular dynamics simulations with Triple Layer and modified Gouy–Chapman models in a 0.1 M NaCl–montmorillonite system
.
J Colloid Interfac Sci
 
339
:
533
541
Tournassat
C
,
Gailhanou
H
,
Crouzet
C
,
Braibant
G
,
Gautier
A
,
Gaucher
EC
(
2009b
)
Cation exchange selectivity coefficient values on smectite and mixed-layer illite/smectite minerals
.
Soil Sci Soc Am J
 
73
:
928
942
Tournassat
C
,
Alt-Epping
P
,
Gaucher
EC
,
Gimmi
T
,
Leupin
OX
,
Wersin
P
(
2011
)
Biogeochemical processes in a clay formation in situ experiment: Part F - Reactive transport modelling
.
Appl Geochem
 
26
:
1009
1022
Tournassat
C
,
Bourg
IC
,
Steefel
CI
,
Bergaya
F
(
2015a
)
Chapter 1 - Surface properties of clay minerals
. In:
Tournassat
C
,
Steefel
CI
,
Bourg
IC
,
Bergaya
F
(Eds.),
Natural and Engineered Clay Barriers, Developments in Clay Science
 .
Elsevier
, p
5
31
Tournassat
C
,
Vinsot
A
,
Gaucher
EC
,
Altmann
S
(
2015b
)
Chapter 3 - Chemical conditions in clay-rocks
. In:
Tournassat
C
,
Steefel
CI
,
Bourg
IC
,
Bergaya
F
(Eds.),
Natural and Engineered Clay Barriers, Developments in Clay Science
 .
Elsevier
, p
71
100
Tournassat
C
,
Bourg
IC
,
Holmboe
M
,
Sposito
G
,
Steefel
CI
(
2016a
)
Molecular dynamics simulations of anion exclusion in clay interlayer nanopores
.
Clay Clay Mineral
 
64
:
374
388
Tournassat
C
,
Gaboreau
S
,
Robinet
J-C
,
Bourg
IC
,
Steefel
CI
(
2016b
)
Impact of microstructure on anion exclusion in compacted clay media
.
CMS Workshop lecture series
 
21
:
137
149
Truc
O
,
Ollivier
J
,
Carcassès
M
(
2000a
)
A new way for determining the chloride diffusion coefficient in concrete from steady state migration test
.
Cem Concr Res
 
30
:
217
226
Truc
O
,
Ollivier
J-P
,
Nilsson
L-O
(
2000b
)
Numerical simulation of multi-species transport through saturated concrete during a migration test—MsDiff code
.
Cem Concr Res
 
30
:
1581
1592
Truc
O
,
Ollivier
J-P
,
Nilsson
L-O
(
2000c
).
Numerical simulation of multi-species diffusion
.
Mater Struct
 
33
:
566
573
Van Loon
LR
,
Baeyens
B
,
Bradbury
MH
(
2005
)
Diffusion and retention of sodium and strontium in Opalinus clay: Comparison of sorption data from diffusion and batch sorption measurements, and geochemical calculations
.
Appl Geochem
 
20
:
2351
2363
Van Loon
LR
,
Soler
JM
,
Bradbury
MH
(
2003
)
Diffusion of HTO
,
36Cl and 125I in Opalinus Clay samples from Mont Terri: Effect of confining pressure
.
J Contam Hydrol
 
61
:
73
83
Van Loon
LR
,
Glaus
MA
,
Müller
W
(
2007
)
Anion exclusion effects in compacted bentonites: Towards a better understanding of anion diffusion
.
Appl Geochem
 
22
:
2536
2552
Vinsot
A
,
Appelo
CAJ
,
Cailteau
C
,
Wechner
S
,
Pironon
J
,
De Donato
P
,
De Cannière
P
,
Mettler
S
,
Wersin
P
,
Gäbler
HE
(
2008a
)
CO2 data on gas and pore water sampled in situ in the Opalinus Clay at the Mont Terri rock laboratory
.
Phys Chem Earth, Parts A/B/C
 
33
,
S54
S60
Vinsot
A
,
Mettler
S
,
Wechner
S
(
2008b
)
In situ characterization of the Callovo-Oxfordian pore water composition
.
Phys Chem Earth, Parts A/B/C
 
33
,
S75
S86
Wang
M
,
Chen
S
(
2007
)
Electroosmosis in homogeneously charged micro- and nanoscale random porous media
.
J Colloid Interfac Sci
 
314
:
264
273
Warkentin
BP
,
Schofield
RK
(
1958
)
Swelling pressures of dilute Na-montmorillonite pastes
.
Clay Clay Mineral
 
7
:
343
349
Wersin
P
,
Soler
JM
,
Van Loon
L
,
Eikenberg
J
,
Baeyens
B
,
Grolimund
D
,
Gimmi
T
,
Dewonck
S
(
2008
)
Diffusion of HTO, Br, I, Cs+, 85Sr2+ and 60Co2+ in a clay formation: Results and modelling from an in situ experiment in Opalinus Clay
.
Appl Geochem
 
23
:
678
691
Wigger
C
,
Van Loon
LR
(
2017
)
Importance of interlayer equivalent pores for anion diffusion in clay-rich sedimentary rocks
.
Environ Sci Technol
 
51
:
1998
2006
Yan
B
,
Wang
Y
,
Killough
JE
(
2016
)
Beyond dual-porosity modeling for the simulation of complex flow mechanisms in shale reservoirs
.
Comput Geosci
 
20
:
69
91
Yokoyama
S
,
Kuroda
M
,
Sato
T
(
2005
)
Atomic force microscopy study of montmorillonite dissolution under highly alkaline conditions
.
Clay Clay Mineral
 
53
:
147
154
Young
A
,
Low
PF
(
1965
)
Osmosis in argillaceous rocks
.
AAPG Bull
 
49
:
1004
1007
Zachara
JM
,
Smith
SC
,
McKinley
JP
,
Resch
CT
(
1993
)
Cadmium sorption on specimen and soil smectites in sodium and calcium electrolytes
.
Soil Sci Soc Am J
 
57
:
1491
1501
Zhao
H
,
Bhattacharjee
S
,
Chow
R
,
Wallace
D
,
Masliyah
JH
,
Xu
Z
(
2008
)
Probing surface charge potentials of clay basal planes and edges by direct force measurements
.
Langmuir
 
24
:
12899
12910
Open access: Article available to all readers online.

Supplementary data