Reactive Transport Modeling of Coupled Processes in Nanoporous Media

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Reactive Transport Modeling of Coupled Processes in Nanoporous Media Christophe Tournassat, Carl Steefel


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., 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 microcontinuum 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, noncoupled 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 if 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 cornersharing tetrahedra with a thickness of about 9.5 Å (Fig. 1).The cations in the octahedral sheet of illite and montmorillonite consist predominantly of Al 3+ .In montmorillonite, isomorphic substitutions take place mostly in the octahedral sheets where Al 3+ is replaced by Fe 3+ or a cation of lower charge (Mg 2+ , Fe 2+ ), whereas in illite isomorphic substitutions take place in the octahedral sheet as in montmorillonite, and in the tetrahedral sheet where Si 4+ is partially replaced by Al 3+ .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 interlayers 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 smallest dimension of which range from less than one nanometer to several micrometers (Fig. 2).Inter-particle and interaggregate 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 variy 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, where the Si chains depolymerize as the Ca/Si atomic ratio increases (Grangeon et al., 2013(Grangeon et al., , 2016(Grangeon et al., , 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(Ma et al., , 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 nonpermanent 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(Lee et al., , 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 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 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 ( DL  in V) at the same position: where  is the water dielectric constant (78.3 × 8.85419•10 -12 F•m -1 at 298 K), F is the Faraday constant (96 485 C•mol -1 ), and zi is the charge of ion i.The concentration of species i in the DL, ci,DL, is given by the Boltzmann equation: where R is the gas constant (8.3145J•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).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, while D, E and F correspond to a ionic strength of 0.072.The seawater composition was taken from Holland (1978) and was simplified into Na + : 481 mmol⋅kg -1 , Mg 2+ : 55 mmol⋅kg -1 ; Ca 2+ : 11 mmol⋅kg -1 ; K + : 11 mmol⋅kg -1 ; Cl -564 mmol⋅kg -1 ; SO4 2-: 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 (xaxis) and half distance between the surfaces (y-axis) (σ0=-0.11C⋅m -2 ).The blue lines indicate 1, 2 and 3 times the inverse of the Debye length (Eq.( 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 saltexclusionary properties of these materials.These properties confer these media with a semipermeable 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, 1960Kemper, , 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).

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): 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: 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 (mfluid 3 ⋅mmedium 2 ⋅s -1 ), q, to the gradient in the hydraulic head, h (m), and the hydraulic conductivity of the medium, K (m⋅s -1 ): Molecular diffusion in porous media is usually described in terms of Fick's First Law: 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): 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, Eq. ( 7) is usually simplified into: 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 Leborgne, 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 chapters 12, 15, and 16 in the present 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: 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: With the assumption that all driving forces Xi needed to describe fully the evolution system are taken into account in Eq. ( 10), each flux Ji is a function fi of the forces Xi, and takes a value of zero if all Xi are zero.Accordingly, the functions fi can be expanded in Taylor series of the form (Lasaga, 1998): ,0 , ,..., , ,..., 0 ...
Keeping only the linear terms of the Taylor series expansion, i.e. making the assumption that higher order terms are negligible, Eq. ( 11) yields: where Li,j are termed phenomenological coefficients.Eq. ( 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: 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, Eqs. ( 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., 2007Gonçalvès et al., , 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: 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(Keller et al., , 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 ( , i DL c ) by scaling them to a mean electrostatic potential (MEP, M  ) that applies to the diffuse layer volume (VDL): ( ) The mean electrostatic potential value can be deduced from the charge balance in the diffuse layer: 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, in the mean electrostatic model, 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., 2010Lee et al., , 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).15) and ( 16)).The difference, given by the color scale, is computed with the equation (cPB-cMEP)/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 (cPB-cMEP)/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 , Mg 2+ : 55 mmol⋅kg -1 ; Ca 2+ : 11 mmol⋅kg -1 ; K + : 11 mmol⋅kg -1 ; Cl -564 mmol⋅kg -1 ; SO4 2-: 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.Using fDL as the volumetric fraction of the porosity to which the mean electrostatic potential is applied, Eqs. ( 15) and ( 16) transform into: ( ) where , i pore c is the average concentration of species i in the pore.Note that ( ) 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 Fig. 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) where w  is the water content of the sample (in kg⋅ kg -1 ).The value of , Cl pore c 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 , Cl pore c 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: ( ) 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), Eq. ( 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 , ,0 / Cl pore Cl cc 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(Rotenberg et al., , 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;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 (κ -1 in m) (Fig. 4, bottom) (Sposito, 2004;Hunter, 2013): In PHREEQC and CrunchClay, the volume of the diffuse layer (VDL in m 3 ), 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: where Vpore is the total volume of the pore (m 3 ), S is the surface area that borders the pore (m 2 ), 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. Eq. ( 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: 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).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 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: and with the following charge balance constraints, which is equivalent to Eq. ( 18): ( ) 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: where ,0 i  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, o i  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): ( ) where ui is the mobility tensor of species i. Eq. ( 28) applies to the bulk water and to the diffuse layer, and introducing Ai, the accumulation factor defined in Eq. ( 25), it follows: ( ) In the absence of an external electric field, there is no electrical current and so: ( ) It follows: The mobility of ion i in the porosity is: and Eq. ( 29) becomes: Eq. ( 33) can be further rearranged to yield: With isothermal, no-flow, and no external electric field conditions assumed here, the only driving forces are chemical potential gradients, and the phenomenological coefficient can be obtained directly from Eq. ( 34): 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 i A to the bulk and diffuse layer accumulation factors Ai: 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, Eq. ( 35) reduces to: 0, , ,0 , 0 In addition, the activity coefficients are also constant, and thus Eq. ( 33) reduces to Fick's law, augmented by the accumulation term Ai that accelerate the diffusive flux of counter-ions 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 Eq. ( 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 (Eq.( 31)).
Reactive transport codes such as PHREEQC and CrunchClay are able to solve Eq. ( 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., 2013Glaus et al., , 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 , when present at trace concentrations, is an important result for the prediction of the fate of radionuclides such as radioactive Cs + or Sr 2+ 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 22 Na + 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 22 Na + concentration.A simple application of Fick's Law would predict that 22 Na + diffuses into the clay and then 22 Na + concentration stabilizes rapidly at the same value in both reservoirs.The actual results evidenced a continuous increase of the 22 Na + concentration in the NaClO4 0.5 mol⋅L -1 reservoir and a decrease in the NaClO4 0.1 mol⋅L -1 reservoir.This result can be explained by the difference of accumulation factor A 22 Na + in the porosity at both ends of the sample.The value of A 22 Na + was anti-correlated with the concentration of stable Na + : and the gradient in accumulation factor A 22 Na + drove the seemingly "up-hill" diffusion of A 22 Na + through the clay sample (Glaus et al., 2013) (Fig. 8).obtained with CrunchClay for the experiment carried out with with a ∼5 mm thick clay plug and reservoir volumes of ∼250 mL (modified after Glaus et al., 2013 andTournassat andSteefel, 2015).

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., 2004Friedmann et al., , 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 socalled 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., 2004Friedmann et al., , 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 Eq. ( 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): is the molar conductivity of species i (A⋅V -1 ⋅m 2 ⋅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: As a consequence, the model proposed by Appelo (2017) follows Onsager reciprocal relationships for bulk water as well as for diffuse layer water.Eq. ( 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 Eqs. ( 43) to (45).
The numerical resolution of Eq. ( 29) in the framework of a reactive transport approach requires that the local gradient of electrical potential ( ) where q=QDL in the diffuse layer and q=0 in bulk water.However, the order of magnitude of F ε in Eq. ( 47) is 10 15 V⋅m⋅mol -1 .Consequently, the second derivative of the potential (first term) in Eq. ( 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 × 10 6 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 Eq. (41) yields: where Rc is the resistivity of the nanoporous media.After discretization on a grid cell, Eq. ( 48) becomes:  50) and ( 39), with Jc being constant throughout the sample, yields: x is a known quantity from Eq. ( 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 electromigration 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).

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 counterions 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).In the most recent, 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  (Soler, 2001;Sun et al., 2016;Neuzil and Person, 2017) 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): 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: For aqueous solutions with molar concentration below 1 mol⋅L -1 , the osmotic pressure can be approximated with van't Hoff's Law: The corresponding chemical osmosis flux is: where g  = is the osmotic pressure head (in m; ρ is the density of water and g the acceleration due to gravity 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: 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: 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.
Eq. ( 59) or similar equations have been widely used in the literature to characterize membrane properties of clay materials (Soler, 2001;Malusis andShackelford, 2002, 2004;Malusis et al., 2003Malusis et al., , 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 ω values from experimental data using Eq. ( 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: ( ) where v again is the average pore fluid linear velocity (m⋅s -1 ).The no current condition imposes: By rearranging Eq. ( 61), it is possible to express the value of The first term of the numerator of the right part in Eq. ( 62) corresponds to a streaming potential, while the second term of the right part corresponds to a diffusion potential.Eq. ( 62) is then reinserted into Eq.( 60 which is a generalization of the single salt equation given in Revil and Leroy (2004).The species dependent coefficient of filtration efficiency, ωi becomes: Eq. ( 65) makes clear that, in multicomponent systems, the coefficient of osmotic efficiency is specific to each chemical species.As expected, Eq. ( 65) indicates that the nanoporous medium is a perfect membrane with regards to the species i (ωi = 1) if Ai = 0, i.e. if the species i is completely excluded from the porosity.Eq. ( 65) also indicates that the coefficient of osmotic efficiency is equal to 0 if Ai = 1 for all species, because ,0 F 0 . In this last case, the porous media has no membrane properties.Eq. ( 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 Eq. ( 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., 2015Hemes et al., , 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 PHREEQCare 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 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.

Figure 2 .
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 nm 3 ); 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 Gaboreau et al. (2016) and Tournassat et al. (2016b).

Figure 3 .
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.Modified from Bourg and Sposito (2011).

Figure 4 .
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

Figure 6 .
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.11C⋅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 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 Fig. 6 and Eqs.(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 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.5otherwise.

Figure 8 .
Figure8.22 Na + "Up-hill" diffusion experiment ofGlaus 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 22 Na + 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 afterGlaus et al., 2013 and Tournassat and Steefel, 2015).
is the current contribution from the chemical activity gradient: Jc,e is the current contribution from the electrical potential gradient: 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): xn+1 correspond to the positions where the potential is imposed.

Figure 10 .
Figure10.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 .
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.Modified from Grambow et al. (2014), with the permission of the Publisher.
average linear pore fluid velocity, accumulation factors and concentrations in a multi-component system:

Table 1 .
Terminology of coupled fluxes (after deMarsily, 1986) . Average Cl -concentrations in pore water ( , (Tournassat et al., 2015bl - 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):