Relevance of charge interactions for contaminant transport in heterogeneous formations: a stochastic analysis

The electrostatic properties of clay (or other charged) mineral surfaces play a significant role in the fate, transport, persistence, and remediation of subsurface contaminant plumes. This study presents a stochastic assessment of the impact and relevance of microscale electrostatic effects on macroscopic, field-scale contaminant transport in heterogeneous groundwater systems involving spatially distributed clay zones. We present Monte Carlo simulations in two-dimensional heterogeneous fields, comprising heterogeneous distributions of physical (i.e., hydraulic conductivity, porosity, tortuosity) and electrostatic (i.e., surface charge) properties, and compare scenarios with different combination and extent of physical and electrostatic processes. The simulations were performed with the multi-continua based reactive transport code, MMIT-Clay, and considering an explicit treatment of the diffuse layer processes. The results reveal that the microscopic electrostatic mechanisms within clay’s diffuse layer can significantly accelerate or retard a particular contaminant depending on its charge, leading to considerably different solute breakthroughs and mass loading/release behaviors in low permeability inclusions. Furthermore, we show that such variations in the macroscale transport behavior, solely driven by charge interactions, are statistically significant over the ensembles of Monte Carlo realizations. The simulations also demonstrate that the omission of electrostatic processes, which is still a common practice in subsurface hydrology, can lead to substantial over- or underestimation of contaminant migration.


Introduction
In the last decades, the investigations of clays, clay rocks, and clay-rich low permeability media have emerged as research areas of growing interest in many scientific disciplines and engineering applications, primarily because of their unique physical, chemical, and mineralogical properties.The low permeability of clay-rich formations makes them ideal for the confinement of radionuclides and/or other contaminants in applications such as long-term nuclear waste disposal in geologic repositories (e.g., Nagra 2002;Posiva 2008;Wersin et al. 2022) or carbon dioxide sequestration in geologic formations (e.g., Bourg et al. 2015;Song and Zhang 2013;Charlet et al. 2017).In the context of contaminant transport in groundwater systems, the presence of such low permeability zones and aquitards can significantly impact the persistence, fate, retention, and migration of contaminant plumes by effectively trapping the contaminants in low-conductive regions and/or subsequently acting as long-term secondary contaminant source (e.g., Liu and Ball 2002;Cherry et al. 2004;Parker et al. 2004;2008;Chapman and Parker 2005;Rasa et al. 2011;Seyedabbasi et al. 2012;Bianchi et al. 2013;Rezaei et al. 2013;Yang et al. 2014Yang et al. , 2016;;Tatti et al. 2016;2019;Muskus and Falta 2018;Sundell et al. 2019;Mosthaf et al. 2021;Muniruzzaman and Rolle 2021;Petrova et al. 2023).Besides hydraulic properties, one of the most distinctive characteristics of clayey media is the abundance of certain minerals containing a higher density of surface charge at the mineral-water interface (e.g., Sposito 1992;Descostes et al. 2008;Gimmi and Kosakowski 2011;Tournassat and Steefel 2015;Tinnacher et al. 2016;Tertre et al. 2018;Soler et al. 2019;Wigger and Van Loon 2017).Electrostatic processes at such charged surfaces lead to the formation of a diffuse ion swarm typically known as the diffuse layer (DL), and give rise to the adsorption of counter-ions and exclusion of co-ions in the proximity of the surface (e.g., Muurinen 2004;Birgersson and Karnland 2009;Tournassat and Appelo 2011).These diffuse layer mechanisms, solely caused by charged minerals, are typically described by the Donnan equilibrium (Donnan and Guggenheim 1932), especially within the framework of continuum-based solute transport models (e.g., Appelo and Wersin 2007;Jougnot et al. 2009;Alt-Epping et al. 2015;Tournassat and Steefel 2019).This approach is based on the assumption that the total pore space is composed of two individual sub-continua: a charge-balanced free water (FW) or bulk water porosity, and a diffuse layer or Donnan porosity containing charged water with a deficit of co-ions and excess of counterions.The diffusive movement of charged solutes in both sub-continua is described by the Nernst-Planck equation, which expresses the flux of a solute as a function of the electrochemical potential gradient, charge density of clay, and porewater solution composition (e.g., Cussler 2009;Appelo 2017).
In subsurface hydrology, low permeability media, aquitards, and clay-rich formations have historically received less attention compared to the permeable porous media, partly because it was traditionally believed that such media are inaccessible for contaminants due to their low hydraulic conductivity.An emerging number of studies have investigated contaminant transport in low-permeability settings, highlighting the importance of diffusioncontrolled mass exchange mechanisms between high and low permeability zones (e.g., Ball et al. 1997;Gusawa and Freyberg 2000;LaBolle and Fogg 2001;Liu et al 2004;Sale et al. 2008;Zhan et al. 2009;Marble et al. 2010;Chapman et al. 2012;Parker and Kim 2015;Yang et al. 2015;2017;Adamson et al. 2016;Falta and Wang 2017;Tatti et al. 2018;Zhang et al. 2020;Rosenberg et al. 2023).Yet, our understanding of chemical transport processes through low permeability media is far from complete because the effects of electrostatic processes, driven by charged surfaces, have rarely been studied for field-scale contaminant transport problems.Most of the published studies dealing with electrostatic effects in clayey and charged materials are exclusively related to the evaluation of the long-term stability of radioactive waste storage systems, where advective flow is typically not considered (e.g., Glaus et al. 2007;Tachi and Yotsuji 2014;Appelo et al. 2010;Wersin et al. 2018).These studies, typically involving lab-scale diffusion experiments, diffusion-dominated field studies, and/or modeling investigations, have established the relevance of charge-induced mechanisms for solute migration in clayey media.Electrostatic effects specifically control the transport of charged species in two ways: (i) via Coulombic interactions between different charged solutes, caused by the diffusion potential derived from the variations in individual ionic diffusivities (e.g., Felmy and Weare 1991;Giambalvo et al. 2002;Boudreau et al. 2004;Appelo and Wersin 2007;Liu et al. 2011;Rolle et al. 2013a;Muniruzzaman et al. 2014;Muniruzzaman and Rolle 2015;2017), and (ii) through surface-solute interactions at the diffuse layer adjacent to the negatively charged surfaces, leading to the adsorption of cations and exclusion of anions (e.g., Muurinen et al. 2004;Leroy et al. 2006;Van Loon et al. 2007;Appelo et al. 2008;2010;Tournassat and Appelo 2011;Glaus et al. 2013;Tertre et al. 2015).The impacts of such electrostatic mechanisms on field-scale contaminant transport in complex, heterogeneous subsurface systems containing both low and high permeability zones are not yet understood, although heterogeneous subsurface sandy-clayey systems are widespread and important in many sedimentary environments (Nichols 2009).
In this study, we perform a stochastic assessment of the key role of charge interactions on field-scale subsurface contaminant transport in complex, heterogeneous sandyclayey domains involving spatially variable distributions of physical and electrostatic properties.The investigation is based on Monte Carlo simulations considering ensembles of two-dimensional binary heterogeneous fields, which are representative of vertical cross-sections of a groundwater aquifer comprising clayey inclusions embedded within a sandy matrix.In each heterogeneous realization, we performed three distinct simulations with different levels of complexity and different combinations of physical and electrostatic processes along with their parameter distributions.The simulations were performed with the multicontinua based reactive transport code, MMIT-Clay (Muniruzzaman and Rolle 2019), and the effects of the diffuse layer processes are thoroughly analyzed by comparing the solute transport and mass-exchange behavior from the analogous scenarios considering and excluding these microscopic electrostatic processes.

Problem statement
We perform stochastic reactive transport simulations of differently charged radionuclide solutes in ensembles of two-dimensional heterogeneous domains (20 m long and 2 m thick), representing transects of sandy groundwater aquifer systems containing low-permeability clayey inclusions (Fig. 1).We specifically focus on the surface charge (r Clay ) of the clay minerals and investigate the impact of the charge-driven electrostatic processes in clay nanopores on the overall mass transfer behavior of the transported solutes between the low-permeability clay zones and the background high-permeability sandy matrix.The simulations are performed considering steady-state flow and transient transport conditions, and Fig. 1 illustrates a schematic of the model setup, boundary/initial conditions, and the relevant microscopic mechanisms.We consider four radionuclides (HTO, I -, 22 Na ?, Cs ? ) as model charged contaminants that are injected through the inlet boundary for a period of 500 days (''loading phase''), followed by the injection of a tracer-free solution (NaCl) for another 500 days (''unloading phase'').These radionuclide species are the most extensively studied tracers in the context of clay and nuclear waste repository studies.Investigations based on these tracers have been carried out both in laboratory setups (e.g., Glaus et al. 2007;2013;Tachi and Yotsuji 2014) and in field experiments (e.g., Appelo et al. 2008;2010;Soler et al. 2019).For the initial condition, the entire simulation domain was assumed to be initially in equilibrium with a 1 mM NaCl solution.

Generation of heterogeneous conductivity fields and groundwater flow equation
To create the desired heterogeneous sandy-clayey domains, we first generate random fields of an auxiliary variable (a = lnK) using a Gaussian covariance model and by following the spectral approach of Dykaar and Kitanidis (1992).We consider an ensemble of 100 Monte Carlo realizations for the auxiliary variable, a, with zero mean, variance of 1, and length scale parameters of l x = 2 m and l z = 0.1 m in the horizontal and transverse directions, respectively.As a next step, we define a cutoff value for the variable a in order to generate the binary sandy-clayey conductivity fields (e.g., Werth et al. 2006; Chiogna et al. ).The locations with the a-values lower than the cutoff value are assigned to the clayey inclusions, whereas the remaining points containing a-values higher than the cutoff are assigned to the background sandy matrix (Fig. 2).In this thresholding process, the clay inclusions are designed to specifically occupy * 20% of the total domain area in all Monte Carlo realizations by following the approach of Lee et al. (2018).Finally, the hydraulic conductivity value is set to 6.14 9 10 -4 m/s in the sandy matrix, whereas the conductivity that is several orders of magnitude smaller (6.14 9 10 -9 m/s) is assigned in the clayey inclusions.The porosity of the matrix and inclusions are kept as 0.41 and 0.71, respectively.
Based on the generated hydraulic conductivity fields, we solve groundwater flow in two-dimensional sandy-clayey domains.Under steady-state flow conditions, the governing equations for hydraulic head, h [m], and stream function, w [m 2 /s], in such domains read as (e.g., Cirpka et al. 1999a): where K [m/s] is the tensor for hydraulic conductivity.Equation ( 1) and ( 2) are solved with a mixed finite element method in rectangular grids (with Dx = 0.4 m, Dz = 0.02 m) to obtain the velocity field required in the solute transport problem.The hydraulic head is solved by applying Dirichlet boundary conditions at the left and right boundaries and no-flow boundaries at the top and bottom boundaries.For the stream function, we consider fixed values at the bottom and top boundaries.The difference between these two boundary values is equal to the total volumetric water flux through the domain.
Figure 2 illustrates the complex patterns of the simulated streamlines (black lines) for two selected realizations.The results show that the streamlines clearly diverge around the low permeability clay zones (blue-colored zones), where the flow velocity is orders of magnitude smaller (blue color) than in the sandy regions where the streamlines are focused (red colors) (Fig. 2c, d).In all simulations, we considered an average flow velocity of * 0.1 m/d, which lies in the typical range of groundwater values.Based on the computed distribution of hydraulic head and stream function, we construct streamline-oriented grids for each realization, following the approach of Cirpka et at. (1999a).As detailed in the next section, the governing multicomponent transport equations are resolved in the streamline-oriented grids.This is advantageous since it allows for minimizing numerical dispersion by eliminating the off-diagonal entries of the hydrodynamic dispersion tensor (Cirpka et al. 1999b).

Multicomponent transport equation including electrostatic effects
In charged porous media, the consideration of the chargeinduced electrostatic processes, occurring at the micro-to nanopores located at the vicinity of the charged surfaces (i.e., diffuse layer), is of utmost importance in assessing migration of charged solutes.Therefore, the description of solute transport in such systems requires explicit treatment of these microscopic diffuse layer processes in connection with the classical advective-dispersive transport mechanisms.In the continuum description of reactive transport processes, these electrostatic effects are treated by where h [-] is the total porosity, f FW and f DL [-] represent the fractions of the total porosity occupied by the free water and diffuse layer, respectively (f ] is the specific discharge vector in the free water porosity, R r [mol/m 2 /s] is the reactive source/sink term, t ir [-] is the stoichiometric coefficient of species i for r-th reaction, and c FW i and c DL i [mol/L] are the concentrations in the free water and diffuse layer porosities, respectively.These two sub-domains are assumed to be in chemical equilibrium and their concentrations are related via Boltzmann equation (e.g., Donnan and Guggenheim 1932;Appelo et al. 2010): where C i [-] represents the ratio of the activity coefficients in the free and diffuse layer water is the charge number of species i, u DL [V] is the mean electrostatic potential in the diffuse layer water (Fig. 1), F [J/V/eq] is Faraday's constant, R [J/mol/K] is the ideal gas constant, and T [K] is the temperature.While the free water porosity is electroneutral, the net charge balance in the diffuse layer water is described by the Donnan polynomial (e.g., Appelo and Wersin 2007;Gimmi and Alt-Epping 2018): where V DL [L] is the volume of water in the diffuse layer, and Su [mol] is the total charge at the mineral surface.In Eq. ( 3), J FW i and J DL i [m 2 /s] are the vectors for the diffusive/dispersive fluxes in the free water and diffuse layer porosities, respectively.Instead of the classical Fickian formulation, we consider the Nernst-Planck equation, which allows a more rigorous description of charged solutes' fluxes by explicitly considering ion-ion Coulombic interactions and local electroneutrality constraints (e.g., Boudreau et al. 2004;Rasouli et al. 2015;Muniruzzaman and Rolle 2016;Rolle et al. 2018;Huang et al. 2022;Trinchero et al. 2022): where, J FW L;i , J FW T;i and J DL L;i , J DL T;i represent the longitudinal and transverse components of the diffusive/dispersive fluxes in the free water and diffuse layer, respectively, / FW d and / DL d [V] are the electrical potentials in the two subcontinua, and /s] are the tensors for the local hydrodynamic self-dispersion coefficients.In a two-dimensional local coordinate system oriented along the principal flow directions (x L , x T ), as in the streamline-oriented grids considered in this study, the local dispersion tensors have only nonzero diagonal entries (e.g., Cirpka et al. 1999a): where D FW L;i , D FW T;i and D DL L;i ,D DL T;i [m 2 /s] are the longitudinal and transverse components of the hydrodynamic self-dispersion coefficients (i.e., when a particular charged solute is transported independently, i.e., ''liberated state'', from the other ions in a multicomponent solute environment, e.g., Boudreau et al. 2004;Muniruzzaman et al. 2014) in the corresponding domains.In the free water, the longitudinal and transverse components of the dispersion coefficients are parameterized by the following linear (Guedes de Carvalho and Delgado 2005; Kurotori et al. 2019) and nonlinear compound-specific (Chiogna et al. 2010;Rolle et al. 2012) relations: where s] is the seepage velocity in the free water domain, d [m] is the average grain size diameter, d [-] is the ratio between the length of a pore channel and its hydraulic radius, b [-] is an empirical exponent that accounts for the effects of incomplete mixing in the pore channels, and Pe i [-] ¼ vd=D aq;i À Á is the grain Pe ´clet number with D aq,i [m 2 /s] being the aqueous diffusion coefficient.We consider d = 5.37 and b = 0.5, which were provided by a comprehensive analysis of transverse dispersion datasets in twodimensional and three-dimensional setups (Ye et al. 2015a).In both expressions, the grain size, d is required as an input parameter at each location in the domain, and we calculate this quantity from the local hydraulic conductivity values by adopting the Hazen approximation (Hazen 1892): , (where c = 0.01 m 0.5 s 0.5 is an empirical proportionality constant), which was also adopted in previous studies on high-resolution transport in heterogeneous domains (e.g., Eckert et al. 2012;Rolle et al. 2013b).In this study, we consider advective transport only in the free water porosity, hence, the hydrodynamic self-dispersion coefficients in the diffuse layer water reduce to pore-diffusion coefficients (i.e., D DL L;i ¼ D DL P;i and D DL T;i ¼ D DL P;i ).In both sub-porosities, the pore-diffusion terms are parameterized as the ratio between the aqueous diffusion coefficient and the tortuosity, s [-] (D FW P;i ¼ D FW aq;i =s FW and D DL P;i ¼ D DL aq;i =s DL ).In Eq. ( 6) and ( 7), the electrical potentials arise due to Coulombic interactions driven by the variation in the diffusive mobility of charged solutes or because of the application of an external electric potential.However, in the absence of any external electric field, there is no electrical current , and P N i¼1 z i J DL T;i ¼ 0 , hence, the longitudinal and transverse components of the individual electrical gradient terms in two sub-continua can be expressed as: After substituting these expressions, Eq. ( 6) and ( 7) can be recast as: These equations can be further rearranged to a more compact notation as: where D FW ij and D DL ij [m 2 /s] are now the tensors for the cross-coupled diffusive/dispersive terms in the free water and diffuse layer, respectively (e.g., Muniruzzaman and Rolle 2016;2019): where d ij is the Kronecker delta function, which is equal to 1 when i = j and equal to 0 if i = j.
Compared to classical single continuum description of solute transport and/or common Fick's law, this multicontinua formulation ensures a greatly improved representation of the conservative and reactive multicomponent transport of charged solutes in physically and electrostatically heterogeneous charged porous media by explicitly considering the i) electrostatic processes in the diffuse layer, ii) inter-ionic coupling between diffusive/dispersive fluxes, concentration and activity gradients, and iii) spatially variable description of the local hydrodynamic self-dispersion coefficients.We solve the reactive transport problem with a finite volume method in streamline-oriented grids, as explained in Sect.2.1.To solve the flow and transport equations presented above, we use a recently developed multi-continua based reactive transport simulator, MMIT-Clay (Muniruzzaman and Rolle 2019).MMIT-Clay uses PHREEQC (Parkhurst and Appelo 2013) as a reaction engine, with the PhreeqcRM module (Parkhurst and Wissmeier 2015) allowing accessing all the chemistry and thermodynamic capabilities of PHREEQC (e.g., Jara et al. 2017;Healy et al. 2018;Rolle et al. 2018;Sprocati et al. 2019;Muniruzzaman et al. 2020).

Simulation scenarios
We perform a series of numerical experiments involving stochastic reactive transport simulations to explore the impact and relevance of charge-induced mechanisms within the clay diffuse layer during solute transport in complex, physically and electrostatically heterogeneous domains.In each Monte Carlo realization, we consider three distinct simulation scenarios specifically focusing on the different combinations of spatially variable physical and electrostatic properties in the model domain.These scenarios are summarized in Table 1.They are designed to include levels of increasing complexity in terms of incorporating different mechanisms and distributions of physical and electrostatic properties.
Physical heterogeneity (P): This scenario reflects typical stochastic numerical analysis of flow and solute transport in heterogeneous groundwater systems (e.g., Rubin et al. 1994, Dagan 2004;Fiori 2003;Jankovic et al. 2003).The focus is on the physical solute transport mechanisms of advection and dispersion, without considering the electrostatic processes generated by the charge of solutes and/or clay surface.These simulations are run with a single-continuum model description, including Fick's law for diffusive/dispersive fluxes.The stochastic realizations of the heterogeneous domains consist of spatially variable distributions of physical properties, such as hydraulic conductivity (K), porosity (h) and tortuosity (s), resulting in spatially variable flow velocity (v) (Fig. 2) and dispersion coefficients (D L and D T , Eqs. 10 and 11).
Physical and electrostatic heterogeneity 1 (P/E 1): This scenario explicitly considers both the physical transport processes and the electrostatic processes within the diffuse layer water in binary sandy-clayey domains.We use the multi-continua formulation including full Nernst-Planck description for solute fluxes in these simulations.In addition to the heterogeneous physical parameters, as considered in scenario P, this scenario also resolves a spatially variable distribution of the electrostatic properties in the domain, homogeneous within each clay inclusion.In this study, we specifically refer to the surface charge of the porous matrix as the key electrostatic property.In each realization, we set a surface charge of 0.11 eq/kg in the clay inclusions, which is representative of the charge density of Opalinus clay (e.g., Appelo et al. 2008;2010).In contrast, the background matrix is assumed to be charge neutral for simplicity, although sandy media can also contain charged minerals (e.g., McNeece and Hesse 2017; Stolze et al. 2020).
Physical and electrostatic heterogeneity 2 (P/E 2): This scenario is essentially an extension of the previous case (P/ E 1), and the ultimate difference is associated with the distribution of surface charge density within the clay inclusions.Unlike a fixed value of surface charge as assigned in P/E 1, we generate a random distribution of this electrostatic property within each clay inclusion.In this step, we consider a range of values between 0 and 2 eq/kg, which is representative of a wide variety of clay minerals (e.g., Appelo and Postma 2005).For specific surface area (A s = 37 m 2 /g) and fractions of free water and diffuse layer porosities (f FW = 20% and f DL = 80%) of clay, we use the values reported in Opalinus clay studies (e.g., Appelo and Wersin 2007;Appelo et al. 2008;2010) in P/E 1 and P/E 2 scenarios.
In all three scenarios, a tortuosity of 6.25 and 2.44 is assigned in the clayey and sandy zones, respectively.Among these cases, Scenario P can be considered representative of conventional high-resolution subsurface solute transport simulations, whereas the electrostatic heterogeneity scenarios (P/E 1, P/E 2) represent a novel reactive transport description that requires a multi-continua based simulator, such as MMIT-Clay.All simulations are performed with a suite of radionuclide tracers, and Table 2 reports the composition of the initial and boundary solutions along with the self-diffusion coefficients of the solute species.

Results and discussion
In this section, we present the results illuminating the electrostatic controls on the transport behavior of different solute species in heterogeneous fields containing different combinations of physical and electrostatic properties (Table 1).The results obtained from the three defined scenarios (P, P/E 1, and P/E 2, Sect.2.3) are compared across the entire ensemble of Monte Carlo realizations.We interpret and discuss the specific impacts of diffuse layer processes on solute breakthrough, mass-exchange between the low-and high-permeability zones, and contaminant storage within (and release from) the low-permeability clay inclusions.

Concentration distribution
Figure 3 shows an example of the results obtained from the multi-continua transport simulations in scenario P/E 1 for a single Monte Carlo realization.The top row panels illustrate maps of the free water concentrations, whereas the solute concentrations in the charged diffuse layer water are shown in the bottom row panels.Upon injection from the left-hand side boundary, the solute species travel along the streamlines, as illustrated in Fig. 2, and gradually enter the low permeability clay inclusions via sandy-clayey interfaces.The effects of physical heterogeneity (hydraulic conductivity and flow velocity, Fig. 2) are evident because the plumes of different solutes show clearly irregular shapes in the free water porosity.Particularly, all the tracer species in the free water domain tend to diverge around the clayey inclusions due to their considerably low hydraulic conductivity (Fig. 3a-d).This implies that the masstransfer between the clayey inclusions and the background sandy matrix is predominantly controlled by diffusive/ dispersive mechanisms.Directly at the inclusions, the electrostatic processes occurring in the clay diffuse layer dominate the transport of different radionuclide species, leading to considerably different plume shapes for the uncharged, negatively charged, and positively charged tracers.For instance, the charge-neutral species, HTO, shows similar concentration distribution and magnitude in both the free water and diffuse layer porosities in the clay zones (Fig. 3a, e).In contrast, the anionic species (I -) tends to show relatively higher concentrations in the free water (dark red colors, Fig. 3b) and a clear depletion (e.g., by almost 3 orders of magnitude compared to the inlet concentration) in the diffuse layer concentrations (orange colors, Fig. 3f) compared to the other species.These phenomena can be explained by the fact that HTO is not affected by the diffuse layer mechanisms, as it can be present and travel similarly in both sub-porosity domains, whereas the electrostatic repulsion in the charged clay domain leads to the anion exclusion effect for I -(e.g., Gvirtzman and Gorelick 1991).The situation is opposite for the cationic tracers ( 22 Na ? and Cs ?), which show significant enrichment (up to almost ten-fold compared to HTO) in the diffuse layer (Fig. 3g, h) but relatively smaller concentration in the free water (Fig. 3c, d).This behavior is due to the electrostatic attraction of these solutes with the negatively charged clay surface.

Breakthrough curves and solute mass evolution
Figure 4 shows the flux-averaged breakthrough curves of the different tracer species at the outlet of the domain.The results are obtained in the considered heterogeneous scenarios for the entire ensemble of Monte Carlo realizations.The blue lines represent the outlet concentration profiles for the scenario considering only physical processes (P), whereas the red and green lines correspond to the scenarios including the effects of the electrostatic processes (P/E 1 and P/E 2).The difference between the two electrostatic Cs ? 10 -3 0 2.07 9 10 -9 I - 10 -3 0 2.00 9 10 -9 Na ?0 1 0 -3 1.33 9 10 -9 Cl - 0 1 0 -3 2.03 9 10 -9 heterogeneity scenarios stems from the fact that P/E 1 considers a uniform charge density distribution in clayey inclusions (red lines), whereas a random distribution of charge density within each clay inclusion is considered in the P/E 2 case (green lines).In each plot, the thin lightcolored lines denote the results from the individual realizations of the heterogeneous fields, whereas the thick dark-colored lines in each simulation scenario represent the median of the computed breakthrough curves over the ensemble of realizations.All four solutes show trends of smoothening of the invading and receding fronts of the temporal concentration  profiles with signatures of long-term tailing, which is indicative of mass-transfer limitations in a dual-domain type system such as sandy-clayey domains.The computed breakthrough curves exemplify a consistent behavior, clearly indicating the influence of the charge driven processes on each tracer species as observed in the 2-D concentration maps (Fig. 3).For instance, HTO breakthroughs do not show any visible difference between the specific scenarios considering (P/E 1 and P/E 2) and neglecting (P) electrostatic heterogeneity, and exhibit overlapping profiles in each realization (Fig. 4a).Conversely, the effects of electrostatic interactions at the charged clay surfaces and the resulting distinct behavior between the three simulation types are clearly observed for the charged species over the entire ensemble (Fig. 4b-d).I -shows distinct patterns with relatively sharper fronts and less pronounced tailings in P/E 1 and P/E 2 compared to scenario P (Fig. 4b).The breakthroughs of the cationic species have remarkably different shapes with significantly lower peak concentrations in the simulation scenarios where charge interactions were considered (Fig. 4c, d).However, the cationic breakthroughs still exhibit less pronounced tailings behavior with much smaller concentration in P/E 1 and P/E 2 scenarios.Such tailings behavior of the ionic tracers is controlled by the electrostatic mechanisms because I -mass storage in clayey inclusions is inhibited by the anion exclusion effect, leading to lower loading and thus lower available mass to release during elution.In contrast, the electrostatic attraction in the clay's diffuse layer enhances cationic adsorption in the inclusions and also hinders the release of the trapped 22 Na ? and Cs ?from clay inclusions when concentration gradient is reversed.This leads to much smaller tailing concentrations for these cationic solutes in the electrostatic scenarios.The observed breakthrough trends can be further explained by visualizing the temporal evolution of the total solute mass within clay inclusions in the domain (Fig. 5).This quantity is calculated by integrating the solute concentrations in the free water and diffuse layer sub-continua over the entire volume of clay inclusions (V clay ): In most cases, the stored mass in the low permeability inclusions generally shows an increasing trend, which peaks at * 500 days when the injection of the tracer species was turned off at the left-hand side boundary of the domain (Fig. 5).Afterwards, the solute mass gradually decreases because of the reversal of the concentration gradients due to flushing the domain with a tracer-free solution and the induced mass-transfer from the clayey inclusions to the background sandy matrix.However, the evolution of the tracer masses for different species (Fig. 5) is consistent with the observed trends in the solute breakthroughs (Fig. 4), and confirms that the diffuse layer mechanisms do not influence the mass transfer of the uncharged solute (HTO) at sandy-clayey interfaces, as reflected in the identical mass profiles of HTO in P, P/E 1 and P/E 2 scenarios (Fig. 5a).Conversely, the anion exclusion effect significantly affects the mass evolution of I -because the electrostatic repulsion within the clay's diffuse layer hinders mass loading in low permeability zones.This is reflected in the lower extent (i.e., less than half) of the stored I -masses in P/E 1 and P/E 2 scenarios compared to scenario P (Fig. 5b).
For positively charged species ( 22 Na ? and Cs ?), the mass storage in the clay inclusions is significantly enhanced (i.e., * three-fold in P/E 2 and * five-fold in P/E 1 compared to P) when electrostatic processes were considered (Fig. 5c, d).In addition to the magnitudes, the shapes of the cationic mass profiles also considerably vary in the electrostatic scenarios (P/E 1 and P/E 2) compared to the analogous physical heterogeneity case (P) and/or other tracer species in different scenarios.Instead of a gradually decreasing pattern after 500 days, 22 Na ? and Cs ?profiles appear to reach a stable plateau, which does not seem to decrease within the simulation time considered in this study (red and green lines, Fig. 5c, d).Such behavior is driven by the electrostatic forces in the charged sub-continuum, which attract positively charged ions to counterbalance the negatively charged surface of clay.Consequently, 22 Na ? and Cs ?mass storage in the clayey zones is significantly enriched during the loading phase (t = 0-500 days).The predominant electrostatic attraction by the charged clay surface greatly inhibits the release of these loaded cations from the clayey inclusions upon reversal of the concentration gradient during the unloading phase (t = 500-1000 days).In scenario P, these chargedriven mechanisms were not considered, and this allowed the cations to behave as uncharged species, leading to breakthroughs and mass profiles similar to those of HTO.These electrostatic mechanisms, inducing enhanced mass storage for cations and diminished for anions, explain the observed trends in the higher peak concentrations in the cationic and relatively sharper profiles for the anionic breakthrough curves in the electrostatic scenarios (Fig. 4).Among the electrostatic heterogeneity scenarios, the effects of the heterogeneous distribution of charge density in P/E 2 are also persistently observed throughout the entire ensemble because both the solute breakthroughs (Fig. 4) and the mass temporal profiles (Fig. 5) deviate in P/E 2 compared to P/E 1. Specifically, the mass storage is comparatively lower for the cations but higher for the anion in P/E 2 compared to P/E 1, which may be due to slightly higher charge capacity, and thus enhanced cation retention, and anion exclusion effects in P/E 1.

Statistics of breakthrough concentrations and solute masses
The solute transport behavior in the different scenarios can be further analyzed by illustrating the distribution of the breakthrough concentration and solute mass populations over the entire ensemble of heterogeneous realizations.
Figure 6 shows the histograms of the probability density function (pdf) for the breakthrough peak concentrations (i.e., concentration at t = 500 days, left column panels) and the final stored solute masses (i.e., loaded mass at t = 1000 days, right column panels).Similar to the breakthrough and mass temporal profiles, the pdfs of HTO tend to show an identical distribution in the three different scenarios (P, P/E 1 and P/E 2) as reflected by their overlapping histograms (Fig. 6a, b).In contrast, the pdfs of the charged species concentrations and masses are significantly different in the electrostatic heterogeneity cases compared to the analogous P scenario.For the anion I -, there is a partial overlap between the concentration and mass populations obtained from P/E 1 and P/E 2 (red and green histograms), but these distributions from the multi-continua simulations are clearly different from the one obtained from the single-continuum based simulations in P (blue histogram, Fig. 6c, d).Furthermore, the distributions for the mass and concentration datasets for the cations are clearly distinct in three different scenarios (Fig. 6e-h).Table 3 summarizes the statistics for the concentration and mass populations.It is interesting to note that the electrostatic processes lead to a relatively wider distribution for the cations and a narrower distribution for the anion in both datasets.For instance, the computed standard deviation (SD) for the breakthrough concentrations (left panel, Fig. 6) reveals that HTO has an identical value in all cases, whereas such values are relatively smaller for I -and much higher for 22 Na ? and Cs ? in the electrostatic scenarios (P/E 1 and P/E 2).
The statistical significance of these deviations in the concentration and solute mass distributions for different solutes in different scenarios can also be assessed.To this end, we perform the Wilcoxon rank sum test (5% significance level in all tests) to analyze these datasets.For HTO, the rank sum test does not show enough evidence to reject the null hypothesis for both the concentration and mass datasets, and suggests that the datasets in P vs P/E 1 (pvalue = 0.93 and 0.95 for concentration and mass, respectively) and P versus P/E 2 (p-value = 0.97 and 0.96 for concentration and mass, respectively) cases indeed come from an identical continuous distribution.In contrast, the test result confirms that the corresponding datasets for I -, 22 Na ?, and Cs ?, representing the breakthrough (scenario P/E 1), green lines-physical and electrostatic heterogeneity (scenario P/E 2).Bold lines: median of the ensemble of realizations concentrations and solute masses from the two types of simulations neglecting (P) and including (P/E 1 and P/E 2) electrostatic processes, reject the null hypothesis and, hence, show that they are not samples from the same continuous distribution.We conclude that an explicit treatment of the charge-driven processes in clay nanopores is required to correctly characterize this electrostatically driven macroscopic transport behavior.Between the two electrostatic scenarios (P/E 1 and P/E 2), the rank sum test confirms that the datasets for 22 Na ? and Cs ?are not samples from identical distributions, which is somewhat expected.Interestingly, I -datasets also suggest a clear shift between P/E 1 and P/E 2 cases by rejecting the null hypothesis, despite the partial overlaps of their mass and concentration populations.
Figure 7 summarizes the datasets of the solute mass distributions for the different heterogeneous scenarios.These box plots further confirm that HTO shows the same mass statistics with an identical median, 25%-75% bounds, and min-max limits for P, P/E 1, and P/E 2 scenarios (Fig. 7a).Conversely, the distributions of I -, 22 Na ?, and Cs ?masses have considerably different ranges when diffuse layer processes are active compared to scenario P, which neglects electrostatic interactions.Specifically, the Blue-only physical heterogeneity (scenario P), red-physical and electrostatic heterogeneity (scenario P/E 1), green-physical and electrostatic heterogeneity (scenario P/E 2) anionic mass populations have relatively smaller medians with narrower distribution in both P/E 1 and P/E 2 scenarios (Fig. 7b).In fact, the standard deviation (SD) for HTO and I -reveals that HTO has an identical value in all scenarios, whereas the anion shows smaller values when charge effects are considered compared to the simulations with only physical processes (Table 3).Furthermore, the cationic masses show an opposite trend with much higher median values in P/E 1 and P/E 2 compared to scenario P (Fig. 7c, d).Interestingly, the electrostatic heterogeneity leads to an enhancement in the variations of cationic mass distributions as reflected in their wider distributions and much higher standard deviation in P/E 1 and P/E 2 scenarios compared to the ones of the P case (Table 3).Among the electrostatic scenarios, the corresponding values of SD and the observed trends in Fig. 7 for the charged species are also consistent with the breakthrough and mass evolution behaviors reported in Figs. 4 and 5.

Conclusions
We have presented a stochastic analysis of subsurface solute transport simulations in two-dimensional heterogeneous sandy-clayey formations, including different combinations of spatially distributed physical and electrostatic properties.In particular, we systematically analyzed the relevance of the electrostatic effects occurring in the diffuse layer of the clay in the context of field-scale contaminant transport through a series of numerical experiments based on Monte Carlo simulations in different heterogeneous domains.The results obtained from the simulated scenarios highlighted the following key points: • The electrostatic interactions, resulting from the clay negatively charged surface and mainly occurring at the length scales of nano-to micrometers adjacent to the surface, are still relevant at the field scale.These microscopic Coulombic interactions have macroscale impacts on subsurface transport of charged contaminants in geologic formations with low permeability zones and/or charged porous media.Depending on the contaminant charge, its migration is strongly affected by the electrostatic attraction or repulsion within the charged diffuse layer water, thus considerably retarding or accelerating its transport, persistence, and mass storage/release behavior in clay inclusions.• Significant differences have been observed in breakthrough curves (up to 10% increase for the anion and 80% decrease for the cations in peak concentrations) and solute mass evolution (up to 14-fold lower and 15-fold higher stored mass in the clay inclusions for the anion and cations, respectively) when electrostatic mechanisms were explicitly considered.Moreover, the charged solutes consistently showed distinct behavior compared to the uncharged species, with a maximum relative difference of 10% and 93% for the anionic breakthrough concentrations and masses, and 80% and 1481% for the cationic breakthrough concentrations and masses, respectively, in the electrostatically heterogeneous scenarios.• The observed transport behavior and the associated impacts of the charge-driven mechanisms are significant not only in each individual realization, but also in a statistical average sense.The mean behavior over the ensemble of Monte Carlo realizations still retains clear signatures of the control by microscopic diffuse layer processes on macroscopic solute transport.Distinct distributions of the electrostatic properties (charge density) also led to considerably different statistical distributions of breakthrough concentrations and solute masses.
This study highlights the need for an explicit treatment of the diffuse layer processes in subsurface solute transport studies at field scales, at least when charged media such as clay are in focus.Neglecting such mechanisms can extremely under-or overpredict the behavior of charged contaminants.Multi-continua based reactive transport formulation, as used in this study, are convenient approaches to rigorously capture coupled physical and electrostatic processes, as well as the effects of the spatially variable distribution of the associated parameters.In this view, classical subsurface solute transport codes and their internal formulations may need to be further developed in order to include important microscale effects such as diffuse layer processes.
The findings of this study are relevant for a wide variety of contaminant transport problems, especially including charged solutes/contaminants and/or porous/fractured media containing charged mineral surfaces.However, the key controls of surface-charge driven mechanisms and the specific effects on contaminant transport should be thoroughly explored beyond the assumptions and limitations of the current study, such as steady state flow, dimensionality, and the lack of reactive processes.Therefore, it would be of interest to investigate the effects of dynamic water flow conditions (e.g., Haberer et al. 2012;Qi et al. 2020), fully 3-D domains involving anisotropy and complex flow topology (e.g., Chiogna et al. 2014;2015;Ye et al. 2015b), presence of other charged minerals (e.g., Prigiobbe and Bryant 2014;McNeece and Hesse 2017;Stolze et al. 2019;Cogorno et al. 2022), biogeochemical (e.g., Bauer et al. 2009;Purkamo et al. 2022) and mineral dissolution/precipitation reactions (e.g., Tartakovsky et al. 2008;Poonoosamy et al. 2016;Soltanian et al. 2015;Fakhreddine et al. 2016;Battistel et al. 2019;2020;Pieretti et al. 2022), multiphase flow and transport (e.g., Ahmadi et al. 2022;Muniruzzaman et al. 2021).Further extension of this approach is also envisioned in synergy with hydrogeophysical exploration (e.g., ERT, DCIP), as well as in the context of engineered remediation involving electrokinetic transport (e.g., Wu et al. 2012;Martens et al. 2021;Alizadeh et al. 2019;Sprocati et al. 2020;Lo ´pez-Vizcaı ´no et al. 2022).In addition to the developments of theoretical concepts, modeling tools and approaches, and numerical studies, we think that further high-resolution experimental investigations, specifically targeting sandy-clayey domains in laboratory setups and/or field environments, will contribute to significantly advance our capability to understand and describe mass transfer and diffuse layer processes in complex subsurface environments.

Fig. 1
Fig. 1 Schematic of the model setup, Monte Carlo approach, boundary and initial conditions in the physical (hydraulic conductivity, porosity, and tortuosity) and electrostatic (surface charge density) heterogeneous fields.The bottom insets show details of electrostatic processes in the clay inclusions and illustrate the pore-scale

Fig. 2
Fig. 2 Examples of the binary hydraulic conductivity fields studied (a, b), and the simulated streamlines (a, b) and velocity fields (c, d) in each realization

Fig. 3 .
Fig. 3. 2-D concentration maps of HTO (a, e), I -(b, f), 22 Na ?(c, g), and Cs ?(d, h) in free water (a-d) and diffuse layer (e-h) porosities during P/E 1 scenario at t = 500 days for a single realization.Flow direction is from left to right

Fig. 4
Fig. 4 Flux-averaged breakthrough curves of HTO (a), I -(b), 22 Na ?(c), and Cs ?(d) at the right-hand side boundary of the domain.Blue lines-only physical heterogeneity (scenario P), red lines-physical

Fig. 5
Fig. 5 Evolution of total mass of HTO (a), I -(b), 22 Na ?(c), and Cs ?(d) in the clay inclusions.Blue lines-only physical heterogeneity (scenario P), red lines-physical and electrostatic heterogeneity

Table 1
Description of the simulation scenarios used in each Monte Carlo realization

Table 3
Statistics of the breakthrough concentrations and loaded solute masses in the clayey inclusions