A model for discrete fracture-clay rock interaction incorporating electrostatic effects on transport

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.


Introduction
There is a growing recognition of the importance of diffusion of aqueous species as a fundamental mass transport process in geological and engineered porous media. Diffusion has historically been relegated to a lesser role in the hydrological sciences in part because at larger length scales, advective transport is argued to be more important, a conclusion suggested by a simple Peclet number analysis. However, diffusion can have an outsized impact in heterogeneous media where exchange of chemical species occurs between low and high permeability zones. An excellent example of this is provided by the case of fractures developed in clay-rich shalesthe fractures themselves are the conduits ultimately for most of the chemical transport, but they can derive their chemical signatures via diffusive exchange between the low permeability shale matrix and the fractures. The potentially contrasting fluid chemistry in the fractures and the low permeability shale can result in sharp chemical gradients that drive diffusive mass transport over relatively short length scales (centimeters), whatever the size of the overall fracture flow system. In this regard, it seems clear that to assess the chemical effects of hydraulic fracturing of shales, it will be important to understand these fracture-matrix interactions.
The importance of fracture-matrix interaction has long been recognized in contaminant hydrology . The rock matrix bordering fractures is particularly important in the case where reactions affect either the source of contaminants, which may be derived from the matrix, or where sorption and/or mineral precipitation provides a sink for the contaminants. Even in the absence of reactions, however, the diffusive loss into the rock matrix can result in retardation of a contaminant plume [2]. In such cases, matrix diffusion constitutes a "diffusive loss" from the conductive fracture system, thus slowing the rate at which the plume migrates in the fracture flow system. The magnitude of the loss term is determined by the rate of diffusive transport into the rock matrix. This effect can be considerably amplified in the case where reactions in the rock matrix adjacent to the fracture steepen the gradients and thus enhance the diffusive flux [7,23].
Most if not all analyses of fracture-matrix interaction have assumed that molecular diffusion can be described adequately with Fick's Laws. However, the limitations of Fick's Laws have become increasingly clear as researchers have demonstrated the need to describe diffusion with the Nernst-Planck equation [24][25][26][27][28][29][30][31][32][33][34][35]. This is because charged chemical species diffuse at different rates, thus creating an electric field or diffusion potential that creates a flux that is in addition to the pure Fickian diffusion due to a concentration gradient. For example, the self-diffusion coefficient for the hydrogen ion, H + , is 8.57 × 10 −9 m 2 /s, while the self-diffusion coefficient for chloride is 1.81 × 10 −9 m 2 /s, or about 4.7 times smaller. In the simplest case of a binary salt solution consisting only of hydrogen chloride, for example, electroneutrality requires that the diffusion rate of the hydrogen ion be slowed while the chloride ion is accelerated. In multicomponent mixtures (e.g., seawater or subsurface brines), however, the charge flux can be compensated by a number of different ions, in which case the full Nernst-Planck calculation is required to determine the diffusivity of individual ions.
In the case where the surfaces of the solid or mineral phases in the porous media are charged, the interaction of the ions with these surfaces needs to be accounted for. Clays with their negative charge provide the best example in the subsurface. In this case, the Nernst-Planck equation needs to be supplemented with the Poisson equation to provide an accurate representation of charged chemical species fluxes. Thus, the Poisson-Nernst-Planck (or PNP) equation provides the most complete continuum representation of ion fluxes in the general case of charged porous media, including clay rock and shale [28,30,36]. The Poisson-Boltzmann can be used in the case of static or equilibrium system. A full pore scale continuum treatment may be possible in the future, but is certainly made difficult by the requirement to resolve individual pores within heterogeneous charged porous media. Progress has been made recently in developing approaches for handling discontinuous interfaces for flow and transport problems for pore-scale and other heterogeneous problems [37,38]. In the interim, mean electrostatic potential (MEP) approaches or dual continuum approaches have been proposed that provide an upscaled if simplified approximation of the PNP equations for clay rocks [28].
In this contribution, we extend the dual continuum approach, a modification of the 1D mean electrostatic potential approach that considers bulk water and diffuse layer porosity and has been presented in various studies [25,27,28,[39][40][41][42]. A dual continuum approach that included electrostatic effects in 2D was presented recently by [43]. In this study, the dual continuum approach is applied to a 2D domain problem that contains a single discrete fracture in which single-phase fluid flow occurs. Only diffusion is considered to take place in the remainder of the 2D clay-rich domain. The primary objective of the study is to analyze how contaminant transport within the fracture system is impacted by interactions with the clay-rich matrix within which electrostatic effects are considered, and this analysis is carried out with the code CrunchClay.

Nernst-Planck equation
Although presented previously, the Nernst-Planck equation is so central to the discussion of transport of charged ions that it is useful to derive it from fundamental thermodynamics. The analysis begins with the principles of irreversible thermodynamics that state that the rate of entropy production, σ, is linearly related to the fluxes J i driven by a generalized set of thermodynamic forces X i [44,45]: where T is the absolute temperature. We can also write the generalized thermodynamic forces as gradients of driving forces, F j , (or potentials) and represent the fluxes as where the L ij are the phenomenological or transport coefficients [44,45]. The matrix of transport coefficients, L ij , typically referred to as the Onsager matrix, contains both diagonal entries connecting a generalized force with its conjugate flux, and off-diagonal terms that relate a generalized force to a nonconjugate flux. In our analysis, we will restrict ourselves initially to thermodynamic forces corresponding to the chemical potential ∇μ and electrical potential ∇ψ: We can write the diffusive flux considering only the chemical potential gradient using the diagonal transport coefficient L 11 where u i is the chemical mobility and C i is the species concentration. We can then use the Einstein relation [46]: to relate the flux to the diffusion coefficient normalized to the thermal energy term that includes the gas constant R and the absolute temperature T: To define the L 12 transport coefficient that relates the chemical flux to the electrical potential (i.e., electrophoresis), we make use of the electrophoretic mobility which combined with the definition of the electrophoretic mobility from the Nernst-Einstein equation [36]: where z i is the charge of the ith species and F is the Faraday constant, results in the Nernst-Planck equation Using the definition of the chemical potential where a i and γ i are the activity and activity coefficients, respectively, for the ith chemical species and μ i 0 is the standard state chemical potential defined for the pure phase. Recalling that ∇ ln C = ∇ C/C, we can write the Nernst-Planck equation so that the first term on the right hand side is recognizable as the Fickian diffusion term: The gradients of the activity coefficients in the second term on the right hand side (R.H.S.) of Eq. (11) are important where a strong ionic strength gradient is present; otherwise, this term is typically minor to negligible in normal groundwaters with a relatively constant background chemistry. The gradient in the electrical potential develops where the system is subjected to a current, or in the case of the diffusion of charged ions, when one or more of ions diffuses at different rates. The differing diffusion rates result in an electrical potential (referred to as a diffusion potential) that generates a flux that drives the system towards electroneutrality (zero entropy production). The result is that in a binary system of pure hydrogen chloride, for example, the diffusion of the hydrogen ion is slowed, while the diffusion of the chloride ion is increased, with the result that their effective diffusivities are equal. The diffusion potential effect gives rise to some interesting behavior depending on the background electrolyte-in a natural water like seawater or subsurface brine, the high background concentrations provide anions and cations (typically NaCl) that can compensate the charge of lower concentration charged ions. The result is that lower concentration anions and cations can diffuse at or close to their selfdiffusion rates. In contrast, in relatively dilute natural waters (e.g., lakes and rivers), the absence of a high concentration background electrolyte means that all of the charged ion diffusivities are tightly coupled. In the case where all of the diffusivities are the same (a common assumption in many reactive transport codes [47], but almost never easily justified), the diffusion potential in the third term of the R.H.S of Eq. (11) vanishes.
The Nernst-Planck equation can be extended to include convection due to fluid flow and an electrical field [36].

Poisson-Boltzmann equation
The Nernst-Planck equation discussed above applies to ion transport in the aqueous phase, whether solids are present or not. Another important effect is the result of the presence of charged surfaces on mineral grains, which adds another level of complexity that must be accounted for. Certainly, the most important example of charged mineral surfaces affecting subsurface transport is provided by clay minerals, which are a major component of such rock types as shale and marls. In the case of a mineral surface with fixed charge, a diffuse layer forms in which ions of opposite charge (counter ions) are attracted, while ions of like charge are repelled. This part of the solution, which does not meet electro-neutrality conditions, is referred as the electrical double layer (EDL) or diffuse layer (DL). The EDL is often conceptually subdivided into two regions: the Stern layer located within the first water monolayers of the interface in which ions adsorb as inner and outer sphere surface complexes, and a diffuse layer located beyond the Stern layer in which a diffuse swarm of ions screens the remaining uncompensated surface charge. EDL properties are an important part of all electrostatic sorption models (e.g., surface complexation models as summarized in Davis and Kent [48]). In the case of clays, which have a negative surface charge, the diffuse layer consists of a swarm of mobile ions with a net positive charge that balance the negative charge of the clay surface. Anions, being of like charge to the clay surface, are present at lower concentrations.
The electrical double layer can be described most rigorously with molecular dynamics models [49], but it is also possible to develop quantitative descriptions based on the continuum Poisson-Boltzmann equation if the solvent can be modeled with a mean field theory and ions are assumed to be represented as point charges [36]. The electrical double layer can be represented schematically, proceeding from the mineral out to the aqueous solution, as the fixed immobile charge of the mineral, a thin Stern layer, followed by the diffuse layer that gradually gives way to electrically neutral (bulk) water ( Fig. 1).
Because there is a net charge in the diffuse layer, the electrical potential is not equal to zero as it is in bulk water. We can write the concentration of ions in the diffuse layer with the Boltzmann equation where C i, 0 is the concentration at infinite distance from the charged mineral surface (i.e., the bulk solution) and z i is the charge of the ion. The Boltzmann equation can be combined with the Poisson equation, which relates the local charge density, ρ e , to the electrical potential where ε is the permittivity. Equations (12) and (13) can be combined to yield the Poisson-Boltzmann equation describing the diffuse layer if the solvent can be modeled with a mean field theory and ions are point charges: Various analytical solutions have been provided for the Poisson-Boltzmann equation, although in general, they are applicable to diffuse layers consisting of symmetric electrolytes and low surface charge. For the general case where these conditions do not hold, the Poisson-Boltzmann equation needs to be solved numerically. An important quantity that emerges from these analytical solutions is the Debye length, λ D , which is the characteristic length scale of the diffuse layer potential decay from the charged surface: where I is the ionic strength of the bulk solution.
The Poisson-Boltzmann equation applies where a local equilibrium distribution of ions exists. In general, this is the case in the absence of fluid flow, or in cases of flow down a uniform channel, where a steady equilibrium ion distribution can develop transverse to the charged surfaces. In the case of flow down heterogeneous channels, or flow in channels with inhomogeneous distribution of surface charge, more complex flow patterns and gradients parallel to the direction of flow can develop and the Poisson-Boltzmann equilibrium distribution no longer applies [36,50]. In this case, it is necessary to use the Poisson-Nernst-Planck (or PNP) set of equations that are obtained by combining Equations (11) and (13).

Mean electrostatic potential model
Since the numerical solution of the Poisson-Boltzmann (PB) equation requires a fine discretization close to charged surfaces in order to resolve the diffuse layer, it is not practical for larger scale (>cm) domains. To be usable at the continuum scale, the Poisson-Boltzmann equations need to be upscaled in some fashion. An upscaling approach based on the mean electrostatic potential (or MEP) model has been discussed extensively in the literature [25,27,28,39,51]. The model assumes that the mean concentration in the diffuse layer, C i;DL , can be obtained by integrating over the diffuse layer volume in the Poisson-Boltzmann equation. This allows us to scale to a mean electrical potential, ψ m , that applies to the diffuse layer volume, V DL : The mean electrostatic potential can then be determined from the charge balance between the charged mineral surface (including the Stern layer) and the diffuse layer: where Q DL is the volumetric charge that must be balanced in the diffuse layer [28]. Comparisons between the accuracy of the full Poisson-Boltzmann equation and the mean electrostatic potential model have been presented in Tournassat and Steefel (2019) [28].
The mean electrostatic potential model has often been referred to as a Donnan model in the literature [25,[51][52][53], but as pointed out by Tournassat and Steefel [28], the Donnan model has some assumptions that are not required by the mean electrostatic potential model.

Dual continuum representation of pore space
A more general and flexible representation of the pore space can be provided using a dual continuum approach [28,39]. In this model, the pore space is divided into two compartments, one corresponding to bulk water that is electrically neutral, and a second that is not electrically neutral and that is subject to the mean electrostatic potential required to balance the surface charge in the pores. As we shall see, the volume of the compartment subject to the mean electrostatic potential needs not be strictly the same as the actual diffuse layer volume. A schematic representation comparing the Poisson-Boltzmann equation and the dual continuum model for a case involving a sodium chloride solution is shown in Fig. 2 [28].
Defining the fraction of the pore space subject to the mean electrostatic potential, f DL as: We can transform Equation (17) so that the charge balance is applied over the entire pore [28]: The fraction f DL can be adjusted to be consistent with theoretical Poisson-Boltzmann predictions or with experimental results. As pointed out [28], this refinement of the mean electrostatic potential model does not necessarily imply that electrically neutral bulk water exists in the center of a pore-the Poisson-Boltzmann equation in Fig. 2 shows that in fact, the electrical double layers from either side of the pore overlap in this case, so no electroneutral water is actually present. The objective of the dual continuum is rather to capture the average concentration that is consistent with the Poisson-Boltzmann prediction. The fixed charge of the mineral surface may be modified by adsorption in the Stern layer. In this case, the charge in the diffuse layer balances both the fixed mineral surface charge and the accumulation of ions in the Stern layer. In practice, both the fraction of the pore volume that is considered to be affected by the mean electrostatic potential, f DL , and the diffuse layer charge can be adjusted to achieve the best fit with the available data. Alternatively, the sorption of ions in the Stern layer can be calculated independently using a surface complexation model (e.g., [47]). The volume of the diffuse layer, V DL , can be defined as a multiple of the Debye length and the surface area of the clays where α DL is an empirical multiplying factor. The use of α DL = 2 in the dual continuum model results in an exact fit with the results of the Poisson-Boltzmann equation for the case of a sodium chloride electrolyte when the width of the pore is large compared with 4λ D [40,[54][55][56]. These modifications make it possible to simulate many systems more accurately than is possible with the classical MEP model (see [28] for a more detailed discussion).
Returning to the Nernst-Planck equation while neglecting convection for the moment, we can develop a general form that is applicable to both compartments in the dual continuum (bulk water and diffuse layer) model described above and that accounts for ion mobility in porous media: where ψ e is the electrical potential in the fluid and as a result of a diffusion potential or external electric field, and not the mean electrostatic potential in the diffuse layer [28]. A j is an accumulation factor defined by The value of accumulation factor, A i , is defined for both the diffuse layer and bulk water as, respectively: where it should be noted that A i = 1 for the bulk water because the electrical potential = 0. Defining the chemical and electrophoretic mobilities for porous media to include the tortuosity, τ i , and porosity, ϕ: Equation (21) then becomes: In the case where there is no electrical current and thus no net flux of charge, which leads to [28]: This result makes it possible to write the diffusive flux without the gradient in the electrical potential: In addition, if we neglect gradients in the activity coefficients, the more compact form is: This is a form of the diffusion equation that is included in CrunchClay, although an alternative formulation (not yet implemented) would be to retain the electrical potential as a primary variable rather than eliminating it. It is relatively easy to verify that if the diffusion coefficients for all of the ions are the same, then the second term on the R.H.S of Equation (29) = 0 (the diffusion potential) and we are left with only the Fickian first term.
With the definition of the accumulation factor, the charge balance equation given in Equation (19) becomes Since the sum of charges in the bulk water are assumed equal to 0 (the bulk water is electrically neutral), Equation (30) becomes Equation (31) is solved in CrunchClay together with the mass balance equations for the primary chemical species.

Definition of total concentrations for bulk and diffuse layer porosity
In the mathematical and numerical formulation used by CrunchClay (and CrunchFlow), the rapid complexation reactions (both aqueous and surface) are formally eliminated using mass action equations, thus leading to the definition of a total concentration, Ω i , as [47]: where C l are the concentrations of the N s secondary (equilibrium) complexes that can be written in terms of the primary (or component) species i. In Equation (32), ν i, l are the stoichiometric coefficients used in the mass balance (second term) and mass action expressions, which also use the activity coefficients, γ p , for the primary species and K l , the equilibrium constants for the complexation reactions. The use of total concentrations in Equation (32) turns the original set of partial differential equations (PDE) into a set of differentialalgebraic equations (DAE). The generalization of the total concentration in the diffuse layer is given by combining Equation (32) with Equation (22) to yield:

Numerical approach
The set of equations given above are solved in the code CrunchClay with a fully implicit approach, that is, the transport and reaction terms are both evaluated at the new time level. The accumulation, transport, and reaction terms are discretized at the present and future time step, giving rise to a set of nonlinear ordinary differential equations that are solved numerically. CrunchClay uses a backwards Euler discretization of the time derivative, which for a system without transport yields a system of ordinary differential equations: where R nþ1 i;k are the kinetically controlled aqueous (or homogeneous phase) reactions and R nþ1 i;m are the kinetically controlled mineral reactions.

Discretization of transport terms
The treatment of the transport terms is complicated to some extent by the use of mean concentrations and accumulation factors defined at grid interfaces. Using a finite difference/ volume formulation, the general form of the diffusive flux for a species ik can be written in one dimension as: where N c is the number of primary (component) species, N s is the number of secondary species, D ik A ik is the mean diffusivity defined at the cell interface, and where we define Δx as: where Δx P and Δx + represent the size of the grid cells at the center point and at the neighboring point, respectively. The mean diffusivity is calculated using a logarithmic-differential average discussed in [29], since this approach provides a more accurate representation of concentrations near material interfaces for the Nernst-Planck equation. Note that the summations in the numerator and denominator of Equation (35) are over all of the ions, both primary and secondary.
With the partitioning between primary and secondary species used in CrunchFlow [47], the form given in Equation (35), the diffusive flux of the total concentration i is given by

Cross-fluxes between continua
As discussed in [30], it is necessary to account for cross-fluxes between the two continua, the diffuse layer and bulk water pore space, in order to maintain mass balance. Therefore, two transport terms are not sufficient to describe the system under investigation. We require three flux terms as shown in Fig. 3: & a flux term from bulk 1 to bulk 2 water volumes, & a flux term from DL 1 to DL 2 water volumes, and & a flux term from bulk 1 to DL 2 water volumes.
A fourth term should be considered from bulk 2 to DL 1 water volumes, but only one of the cross-flux terms is non-zero between the two. In the case of the fracture-matrix system considered in this study, transport along the fracture itself occurs only within bulk water porosity. Transport within the matrix is split between bulk-bulk and DL-DL water volumes, with no cross-fluxes present because the fraction of water volumes is the same. However, between the fracture (100% bulk water porosity) and the matrix (consisting of both DL and bulk water porosity), the diffusive cross-flux is important. The full averaging rules for the cross-flux are presented in [27,30].

Newton's method and Jacobian matrix
The nonlinear ordinary differential equations are solved with Newton's method: where f i are the function residuals (mass balance equations written typically in terms of the total concentrations, Ω i (defined below), and ∂ f i = ∂C j are elements of the Jacobian matrix.
Typically, the unknowns are the logarithms of the concentrations, which provide improved numerical stability: Perhaps, the most significant challenge associated with a fully implicit approach based on the Newton method is the calculation of the derivative terms. Although not demonstrated in this study, the use of numerical derivatives results in much longer execution times (mostly due to the recalculation of expensive transcendental functions), but also diminished accuracy and convergence properties. The derivatives of the secondary species concentrations with respect to the primary species are given by which gives a derivative of the total concentration as where δ i, k is the Kronecker delta (= 1 if i = k, otherwise = 0). In the dual continuum system considered here, the concentrations of ions in the diffuse layer also need to be accounted for in the accumulation term. Ignoring surface complexes for simplicity, the derivative of the accumulation term with respect to primary species k is given by: The derivative of the accumulation term with respect to the mean electrostatic potential is: With this partitioning of the chemical system into primary (direct unknowns) and secondary species, the derivative of the diffuse layer charge balance Equation (31) with respect to the kth primary species is given by For the case of surface complex formation in the Stern layer, which modifies the surface charge that needs to be balanced in the diffuse layer, the derivative becomes: where C s is the concentration of the surface complex and ν is is the number of moles of primary species i in the surface complex that is differentiated with respect to primary species k, the stoichiometric coefficient for which is ν k s . Here, we have assumed a nonelectrostatic surface complexation model for the sake of simplicity. The derivative of the charge balance equation given in Equation (31) with respect to the mean electrostatic potential is given by where again we have assumed a non-electrostatic surface complexation model for the sake of simplicity. The derivatives of the diffusive flux terms are given in Appendix 1.

Assembly of Jacobian and Newton solve
In the case of a global implicit treatment of reaction and transport, the size of the Jacobian matrix which must be constructed and solved becomes larger, since each function will include contributions from the concentrations of the grid cell itself and from neighboring grid cells that are used in the discretization of the gradients.
As an example, in the case in this study where we use a simple 2D structured finite difference grid with N c unknown chemical concentrations at each nodal point, the form of the Newton equations to be solved in the interior of the domain is given by: where i refers to the primary species residual (mass balance), k is the primary species number which we are differentiating with respect to, and jx and jy are the nodal points of which there are Nx and Ny, respectively. The logarithms of the concentrations are solved for because of the improved numerical stability this provides. The Jacobian matrix in the case of this 2D system takes a banded form, with the center tridiagonal band corresponding to the discretization in the X direction (i.e., jx-1 and jx + 1), while the right band and left band are separated by a gap of NX-2 and represent the discretization in the Y direction, jy + 1 and jy-1 respectively (Fig. 4). Although this is a general form for a 2D finite difference/volume discretization based on a five point finite difference stencil, here each entry is a sub-matrix of dimension Nc by Nc, the number of primary species in the system [57].  (47)) We solve the system of linear equations within a single Newton iteration using the PETSc libraries [58]. An MPIbased version is under development (normally what the PETSc libraries are used for), but the simulation results shown in this work used only a single processor, that is they are solved with sequential PETSc. The sparse linear system is solved with the GMRES Krylov solver in PETSc, with block Jacobi preconditioning. While the overall system is sparse, the Nc by Nc submatrices are typically dense as a result of the coupling between primary species due to complexation and multicomponent diffusion.

Application to fractures in clay rocks
With the mathematical and numerical formalism developed above, we can develop a model that is applicable to clay rocks, including shales and marls. The focus here is to apply the model for diffusion considering electrostatic effects coupled to flow in a discrete fracture within an idealized 2D system. In this conceptually (if not numerically) simple example, we consider a single fracture in clay rock (Fig. 5). We assume a fracture of 10 m length with a constant aperture of 2 mm as an example calculation. Within the model domain, we consider only ½ the system, so a fracture half-width of 1 mm bordered by 32 mm of clay rock. A constant flow rate of 1000 m per year is assumed for the fracture.
As a representative example of a clay rock, we use approximate properties of the Opalinus Clay. The fracture is assumed to have a tortuosity of 1.0, while the tortuosity of the clay rock is assumed to be 0.35 for the bulk water porosity and 0.1 for cations and uncharged species and 0.0001 for anions within the clay rock. The tortuosity is assumed to be 1000 times lower because of the presence locally of small pore throats where anions are excluded along the various diffusion paths, thus reducing the overall tortuosity of the matrix. However, this value can be set by the user of the software at runtime and is usually calibrated by laboratory or field experiments.
Based on previous studies of the Opalinus Clay, we use a value of 7.8% for the bulk water porosity and 8.2% for the diffuse layer porosity. The clay rock is assumed to contain 34% illite by weight, with a surface charge for the illite of 0.2 C/m −2 corresponding to a cation exchange capacity (CEC) of 0.2 mol c /kg illite (where mol c is a mole of surface charge). Assuming a dry density for the rock of 2.7 g/cm 3 and the combined porosity (bulk and diffuse layer) of 16.0%, this yields 154.6 mol c per m 3 porous medium. The self-diffusion coefficients used in the simulation are given in Table 1. Table 2 gives the initial and boundary conditions for the test problem. The concentration of the background electrolyte (primarily NaCl) is important here because it provides most of the charge compensation for the tracer anion and cations.
For the mineral illite, it is necessary to consider surface complexation and its effects on the net charge that need to be balanced by the diffuse layer (fixed charge + Stern layer charge from the formation of the surface complexes). This is in addition to effects on transport due to sorption of the cations under consideration. For the formation of the surface complexes of Na + (the dominant cation in solution), Ca ++ present in the groundwater, and a trace monovalent cation, Cat + (e.g., 22 Na + ), and trace divalent cation, Cat ++ (e.g., 90 Sr ++ ), the reactions and equilibrium constants are given by Fig. 5 Schematic of a single fracture developed in clay rock, with flow along the fracture in the direction of the arrow. The infiltrating water in the fracture contains trace (but non-zero) anion and cation contaminants, which are transported via advection up the fracture, but also via diffusion into the clay-rich rock matrix. A constant flow rate of 1000 m/yr is assumed in the fracture We also compare the case where no accumulation of cations occurs in the Stern layer so as to quantify the effects of Stern layer adsorption via formation of the surface complexes. The concentration fields for the tracers (HTO, An − , Cat + , and Cat ++ ) are plotted after 10 days in Fig. 6 for the base case that includes surface complexation. From these results, it is apparent that the anion front moves significantly farther through the fracture than does the cation. This is primarily due to the diffusive loss of the ions into the rock matrix, which is much stronger for the cation than it is for the anion. This is conceptually similar to the effect described by Tang et al. (1981), although that study considered a pure Fickian model in which it was not possible to capture the effects for both anions and cations simultaneously. Note that Fig. 6 shows the calculation with the surface complexes formed according to the equilibrium constants in Table 1. Figure 7 shows the behavior of the system when no surface complexation is included in the calculation (the surface complexation reactions in Table 1 are not considered in this case). The absence of surface complexation has two effects: (1) there is no accumulation of charge in the Stern layer, so more charge compensation must occur by cations in the diffuse layer, and (2) there is no sorption of the monovalent and divalent cations. With no Stern layer (Fig. 7), the monovalent cation spreads slightly further into the matrix compared with the case in which a Stern layer is present (Fig. 6), since there is less retardation via sorption this case. For the divalent cation, the effects are subtler. In the case where no Stern layer is present (no surface complexation), the higher surface charge as a result of no Stern layer compensation results in more accumulation of Cat ++ in the matrix, thus reducing its migration rate through the fracture system. The overall behavior is perhaps shown most clearly in Fig. 8, where we plot the breakthrough of the anion and cation at the end of the fracture (i.e., the fracture effluent). The minimum retardation is shown by the trace anion, which is largely excluded from the clay-rich matrix, and also does not sorb. The tritiated water (HTO) shows more loss into the matrix than does the anion because it can diffuse through both the diffuse layer and bulk water porosity. The behavior of the cations is more complex-comparing Fig. 8 a and b (with and without surface complexation) indicates that sorption (surface complexation) in the Stern layer is not a large effect. The cation retardation is accentuated compared with the HTO tracer due to the accumulation in the diffuse layer (compare Cat + and HTO in Fig. 8). Despite the lower diffusivity of the divalent tracer cation, its stronger accumulation in the diffuse layer in Fig. 8b (no complexation) slows its migration compared with the monovalent cation in the fracture. The case with no electrical double layer (no EDL) is shown in Fig. 8c.
The approximate behavior of the anion and cation could be captured with independent calculations, but only by considering differing (and possibly unrealistic) porosity and tortuosity values. More important, the coupled nature of the system, which may be subject to transient changes in a realistic setting like that of a nuclear waste repository, would be lost with the

Conclusions
We have presented a model based on CrunchClay for a discrete fracture-clay matrix system that accounts for electrostatic effects on ion transport. The electrostatic effects affect transport in two ways: (1) through the development of a diffusion potential described by the Nernst-Planck equation and (2) more strongly due to the effects of partial anion exclusion in the charged clay media. A flexible dual continuum formulation has been developed that includes calculation of diffusive fluxes using mean concentrations, diffusion coefficients, and accumulation factors at cell interfaces [29,30]. The resulting set of differential-algebraic equations are discretized in 2D with a five point finite difference/volume formulation and solved with Newton's method using a fully implicit formulation. As an application, we considered a single discrete fracture within which flow and advective transport occurs that is bordered by clay-rich matrix in which the only transport process is diffusion. The results demonstrate the strong retardation of cations due to diffusive loss into the rock matrix, while anions are more weakly retarded in the fracture flow system due to their partial exclusion from the matrix. The diffusive loss of cations is accentuated by their accumulation in the diffuse layer within the clay-rich matrix as ∂ J V;l ∂C i2−X P ν i2;l ≠0 ¼ D l A l Δx ν il C l;X P C i2;X P ! ð66Þ Derivatives of arithmetic mean concentration are given by: ∂C l ∂C i2-X P ¼ Δx P ν i2;l C l;X P C i2;X P Δx þ þ Δx P ∂C l ∂C i2-Xþ ¼ Δx þ ν i2;l C l;X þ C i2;X þ Δx þ þ Δx P Derivatives of logarithmic mean concentration are given by: ∂C l ∂C i2−X P ¼ ν i2;l C l;X p C l;Xþ C l;X P −1−ln C l;Xþ C l;X P ! C i2;X P ln C l;Xþ C l;X P ln C l;Xþ C l;X P ð76Þ ∂C l ∂C i2-X þ ¼ ν i2;l C l;X þ C l;X P C l;X þ −1−ln C l;X P C l;X þ ! C l;X þ ln C l;Xþ C l;X P ln C l;Xþ C l;X P ð77Þ Derivatives of accumulation factor A i with respect to diffuse layer potential ψ m are given by: Derivatives of averaged diffusion coefficient D i A i with respect to diffuse layer potential ψ m are given by: ∂A i;X P ∂ψ m;X P ð80Þ Derivatives of diffusive flux function with respect to diffuse layer potential ψ m are given by: For either X P or X + : ∂ J I;i ∂ψ m;X P ¼ ∂D i A i ∂ψ m;X P C i;X þ −C i;X P Δx ð84Þ ∂ J IV ∂ψ m;X P ¼ ∑ ∂ J V;ik ∂ψ m;X P ¼ ∂D C A ik ∂ψ m;X P C ik;X þ −C ik;X P Δx ð92Þ ∂ J VI;ik ∂ψ m;X P ¼ z ik ∂D ik A ik ∂ψ m;X P C ik ð94Þ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.