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).

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.

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).

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)

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).

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.

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).

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).

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.

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).
Open access: Article available to all readers online.