A model based on the code CrunchClay is presented for a fracture-clay matrix system that takes electrostatic effects on transport into account. The electrostatic effects on transport include those associated with the development of a diffusion potential as captured by the Nernst-Planck equation, and the formation of a diffuse layer bordering negatively charged clay particles within which partial anion exclusion occurs. The model is based on a dual continuum formulation that accounts for diffuse layer and bulk water pore space, providing a more flexible framework than is found in the classical mean electrostatic potential models. The diffuse layer model is obtained by volume averaging ion concentrations in the Poisson-Boltzmann equation, but also includes the treatment of longitudinal transport within this continuum. The calculation of transport within the bulk and diffuse layer porosity is based on a new formulation for the Nernst-Planck equation that considers averaging of diffusion coefficients and accumulation factors at grid cell interfaces. Equations for function residuals and the associated Jacobian matrix are presented such that the system of nonlinear differential-algebraic equations can be solved with Newton’s method. As an example, we consider a 2D system with a single discrete fracture within which flow and advective transport occurs that is coupled to diffusion in the clay-rich matrix. The simulation results demonstrate the lack of retardation for anions (e.g., 36Cl−) of the contaminant plume within the fracture flow system because they are largely excluded from the charged clay rock, while the migration of cations (e.g., 90Sr++) is more strongly attenuated. The diffusive loss of divalent cations in particular from the fracture is accentuated by their accumulation in the diffuse layer within the clay-rich matrix.

Molecular diffusion of dissolved species is a fundamental mass transport process affecting many environmental and technical processes. Whereas diffusive transport of single tracers can be described by Fick's law, a multicomponent approach based on the Nernst-Planck equation is required for charge-coupled transport of ions. The numerical solution of the Nernst-Planck equation requires special attention with regard to properties that are required at interfaces of numerical cells when using a finite difference or finite volume method. Weighted arithmetic and harmonic averages are used in most codes that can solve the Nernst-Planck equation. This way of averaging is correct for diffusion coefficients but inappropriate for solute concentrations at interfaces. This averaging approach leads to charge balance problems and thus to numerical instabilities near interfaces separating grid volumes with contrasting properties. We argue that a logarithmic-differential average should be used. Here this result is generalized, and it is demonstrated that it generally leads to improved numerical stability and accuracy of concentrations computed near material interfaces. It is particularly relevant when modeling semipermeable clay membranes or membranes used in water treatment processes.

porousMedia4Foam is a package for solving flow and transport in porous media using OpenFOAM® - a popular open-source numerical toolbox. We introduce and highlight the features of a new generation open-source hydro-geochemical module implemented within porousMedia4Foam, which relies on micro-continuum concept and which makes it possible to investigate hydro-geochemical processes occurring at multiple scales i.e. at the pore-scale, reservoir-scale and at the hybrid-scale. Geochemistry is handled by a third party package (e.g. PHREEQC) that is coupled to the flow and transport solver of OpenFOAM®. We conducted benchmarks across different scales to validate the accuracy of our simulator. We further looked at the evolution of mineral dissolution/precipitation in a fractured porous system. Application fields of this new package include the investigation of hydro-bio-geochemical processes in the critical zone, the modelling of contaminant transport in aquifers, as well as and the assessment of confinement performance for geological barriers.

Molecular dynamics simulations were conducted to evaluate the effect of water desaturation on ion and water density distributions and diffusion properties in a clay mineral mesopore. Simulations were carried out using the polarizable ion model (PIM) force field. Desaturation of the pore was responsible for significant changes in localization and mobility of water molecules and ions, especially in the most unsaturated conditions, for which charge overcompensation by cations adsorbed on the surface through inner-sphere complexes was observed. While cation and anion mobility decreased monotonously with increasing desaturation, the pore-averaged diffusion coefficient of water molecules exhibited a complex evolution with an increased diffusion coefficient value at low desaturation level and then a drop of this value at high desaturation level. Water diffusion was influenced by two antagonistic processes, i.e., acceleration of water molecule diffusion at the water/air interface and slowdown at water/clay mineral interface. These observations can be put in perspective with differences in macroscopic diffusion properties among cationic, anionic, and water tracers observed in experiments conducted on desaturated clay mineral rock samples.

Montmorillonite layer edge surfaces have pH-dependent properties, which arises from the acid-base reactivity of their surface functional groups. Edge surface acidity (with intrinsic reaction equilibrium constant, pK_a) is a chemical property that is affected by crystal structure. While a cis-vacant structure predominates in natural montmorillonites, prior molecular-level studies assume a centrosymmetric trans-vacant configuration, which potentially leads to an incorrect prediction of montmorillonite acid-base surface properties. We computed intrinsic acidity constants of the surface sites of a montmorillonite layer with a cis-vacant structure using the first-principles molecular dynamics-based vertical energy gap method. We evaluated pK_a values for both non-substituted and Mg-substituted layers on common edge surfaces (i.e., surfaces perpendicular to [010], [01̅0], [110], and [1̅1̅0] crystallographic directions). The functional groups ≡Si(OH), ≡Al(OH₂)₂/≡Al(OH)(OH₂), and ≡SiO(OH)Al sites on surfaces perpendicular to [010] and [01̅0] and ≡Si(OH)^U, ≡Si(OH)^L, ≡Al(OH₂), and ≡Al(OH₂)₂ on surfaces perpendicular to [110] and [1̅1̅0] determine the proton reactivity of non-substituted cis-vacant edge surfaces. Moreover, the structural OH sites on edge surfaces had extremely high pK_a values, which do not show reactivity at a common pH. Meanwhile, Mg²⁺ substitution results in an increase in pK_a values at local or adjacent sites, in which the effect is limited by the distance between the sites. A surface complexation model was built with predicted pK_a values, which enabled us to predict surface properties as a function of pH and ionic strength. Edge surface charge of both trans- and cis-vacant models has little dependence on Mg²⁺ substitutions, but the dependence on the crystal plane orientation is strong. In particular, at pH below 7, edge surfaces are positively or negatively charged depending on their orientation. Implications of these findings on contaminant adsorption by smectites are discussed.

Performance assessment in deep geological nuclear waste repository systems necessitates an extended knowledge of the pore water chemical conditions prevailing in host-rock formations. In the last two decades, important progress has been made in the experimental characterization and thermodynamic modeling of pore water speciation, but the influence of experimental artifacts and uncertainties of thermodynamic input parameters are seldom evaluated. In this respect, we conducted an uncertainty propagation study in a reference geochemical model describing the pore water chemistry of the Callovian-Oxfordian clay formation. Nineteen model input parameters were perturbed, including those associated to experimental characterization (leached anions, exchanged cations, cation exchange selectivity coefficients) and those associated to generic thermodynamic databases (solubilities). A set of 13 quantities of interest were studied by the use of polynomial chaos expansions built non-intrusively with a least-squares forward stepwise regression approach. Training and validation sets of simulations were carried out using the geochemical speciation code PHREEQC. The statistical results explored the marginal distribution of each quantity of interest, their bivariate correlations as well as their global sensitivity indices. The influence of the assumed distributions for input parameters uncertainties was evaluated by considering two parametric domain sizes.

Radionuclide transport in smectite clay barrier systems used for nuclear waste disposal is controlled by diffusion, with adsorption significantly retarding transport rates. While a relatively minor component of spent nuclear fuel, 79Se is a major driver of the safety case for spent fuel disposal due to its long half-life (3.3 × 105 yr) and its low adsorption to clay (KD < 10 L/kg), thus a thorough understanding of Se diffusion through clay is critical for understanding the long-term safety of spent fuel disposal systems. Through-diffusion experiments with tritiated water (HTO, conservative tracer) and Se(VI) were conducted with a well-characterized, purified montmorillonite source clay (SWy-2) under a constant ionic strength (0.1 M) and three different electrolyte compositions: Na+, Ca2+, and a Na + -Ca2+ mixture at pH 6.5 in order to probe the effects of electrolyte composition and interlayer cation composition on clay microstructure, Se(VI) aqueous speciation, and ultimately diffusion. The results were modeled using a reactive transport modeling approach to determine values of porosity (ε), De (effective diffusion coefficient), and KD (distribution coefficient for adsorption). HTO diffusive flux was higher in Ca-montmorillonite (De = 1.68 × 10−10 m2 s−1) compared to Na-montmorillonite (De = 7.83 × 10−11 m2 s−1). This increase in flux is likely due to a greater degree of clay layer stacking in the presence of Ca2+ compared to Na+, which leads to larger inter-particle pores. Overall, the Se(VI) flux was much lower than the HTO flux due to anion exclusion, with Se(VI) flux following the order Ca (De = 1.03 × 10−11 m2 s−1) > Na–Ca (De = 2.12 × 10−12 m2 s−1) > Na (De = 1.28 × 10−12 m2 s−1). These differences in Se(VI) flux are due to a combination of factors, including (1) larger accessible porosity in Ca-montmorillonite due to clay layer stacking and smaller electrostatic effects compared to Na-montmorillonite, (2) larger accessible porosity for neutral-charge CaSeO4 species which makes up 32% of aqueous Se(VI) in the pure Ca system, and (3) possibly higher Se(VI) adsorption for Ca-montmorillonite. Through a combination of experimental and modeling work, this study highlights the compounding effects that electrolyte and counterion compositions can have on radionuclide transport through clay. Diffusion models that neglect these effects are not transferable from laboratory experimental conditions to in situ repository conditions.

Quantification of adsorption processes on clay mineral surfaces is often necessary to predict the extent and the evolution of contaminants’ migration in surficial and underground environments. Many studies have been dedicated to retention measurement as a function of pH and ionic strength in relation with the two main identified adsorption processes for clay minerals, that is, cation exchange on their basal surfaces and surface complexation on their edge surfaces. The latest process has been repeatedly assessed as an effective retention mechanism at circumneutral pH conditions, which often prevail in natural environments. This assessment must however be tempered by the lack of information about competitive processes that can take place with the numerous chemical species present in natural settings, compared to simplified systems investigated in laboratory experiments. In this study, we quantified experimentally the competition between Pb2+, Co2+, Zn2+, and Mg2+ for specific adsorption on montmorillonite edge surfaces. Zn2+ was an effective competitor with Pb2+ and Co2+, and our results showed also unambiguously the influence of Mg2+ concentration levels on the specific adsorption of Pb2+ and Co2+. Because of the high ionic strength used in the experiments, cation exchange with Mg2+ was dismissed as a possible reason for such competition process, leaving specific competitive adsorption on edge surfaces as a unique explanation for our observations. Modeling of Pb2+ adsorption data with a state-of-art electrostatic complexation model for montmorillonite edge surfaces, supported by state-of-the-art, made it possible to distinguish two types of possible competition driving forces: Zn2+ competition for adsorption site occupancy but also detrimental changes in surface electrostatic potential following Mg2+ adsorption on sites neighboring those of Pb2+ adsorption. Mg2+ competition observed in our experiment should apply in most of clayey environments. Consequently, adsorption data obtained on pure clay mineral phases, and the associated models that have been built based on these data without considering the geochemical background in competitive species, may overestimate the retention properties of clay minerals when applied to natural settings.

Clay minerals retain or adsorb metal ions in the Earth’s critical zone. Rocks, sediments and soils rich in clay minerals can concentrate rare earth elements (REEs) in ion adsorption-type deposits (IADs) and are similarly effective at metallic contaminant remediation. However, the molecular-scale chemical and physical mechanisms of metal ion retention remain only partly understood. In this Review, we describe the nature, location and energy requirements of metal retention at clay mineral surfaces. Retention originates mainly from electrostatic interactions during cation exchange at low pH and chemical bonding in surface complexation and precipitation at neutral and high pH. Surface complexation can induce surface redox reactions and precipitation mechanisms including neoformation of clay mineral layered structures. In IADs, outer-sphere adsorption is the major retention mechanism of REE ions. By contrast, the use of clay minerals in pollution control relies on various mechanisms that can coexist, including cation exchange, surface complexation and nucleation growth. To more effectively leverage clay mineral–metal interactions in resource recovery and contaminant remediation, complex mechanisms such as surface precipitation and redox reactions must be better understood; for instance, by utilizing advances in quantum mechanical calculations, close combination between synchrotron and simulation techniques, and upscaling of molecular-level information in macroscopic thermokinetic predictive models.