Modeling diffusion processes in the presence of a diffuse layer at charged mineral surfaces: a benchmark exercise

The electrostatic properties of clay mineral surfaces play a significant role in their diffusion properties. The negative electrostatic potential field at clay mineral surfaces results in the presence of a diffuse layer that balances the mineral surface charge. The diffusion properties of the porosity fraction that is affected by this phenomenon are different from the diffusion properties of electroneutral bulk water. These properties have attracted growing interest from diverse communities in the past years, especially in the field of study of radioactive waste disposal. The influence of the diffuse layer can be described at the continuum scale by a set of equations that are formulated in terms of the Nernst-Planck equation. The number of codes that can handle the coupling between transport properties in clay affected by the presence of a diffuse layer in the porosity and chemical reactions is very limited, and no benchmark exercises have been published yet that make it possible to validate the numerical implementation of these equations in reactive transport codes. The present study proposes a set of benchmark exercises of increasing complexity that highlight caveats related to the finite difference (volume) treatment of the Nernst-Planck equation in the presence of a diffuse layer in heterogeneous systems. Once these problems are identified and solved, the codes PHREEQC, CrunchClay, and a new Fortran routine written for this study gave results in very good agreement for most of the benchmark exercises. When present, the differences in results were directly traceable to the differences in averaging methods at grid cell boundaries, and to the consideration or the omission of the activity gradient term in the Nernst-Planck equation.


Introduction
The mineralogical and chemical properties of clays have been the subject of longstanding study for the long-term disposal of nuclear wastes in geological repositories [1,2]. The low permeability of clay materials, including shales, ion swarm, or diffuse layer, as opposed to the bulk porosity where electroneutrality prevails [13].
Numerical methods for modeling macroscopic properties of clay media with the consideration of the presence of a diffuse ion swarm have attracted growing interest from diverse communities in the past years [9,[14][15][16][17]. Yet, the number of codes that can handle the coupling between transport properties in clay affected by the presence of a diffuse layer in the porosity and chemical reactions is very limited. To our knowledge, only PHREEQC [18] and CrunchFlowMC (now CrunchClay) [19] have been released with these full capabilities. A recently modified version of FLOTRAN [20,21] includes some of the capabilities with a different treatment on the basis of an alternative conceptual model. With these codes, the influence of the diffuse layer is described by a specific set of equations that are formulated in terms of the Nernst-Planck equation. These equations are currently limited to applications where the advective flow is restricted to the bulk porosity. This limitation makes these codes suitable for simulating systems where advection is indeed negligible and diffusive processes dominate instead, such as clay barriers. The present benchmark aims at providing a set of exercises that may be used by code developers to verify their implementation of the Nernst-Planck equation together with the presence of diffuse layer capabilities in their software.

Nernst-Planck equation
Tournassat and Steefel [12] reviewed the basics of the application of the Nernst-Planck equation applied to diffusive processes in the diffuse layer bordering charged surfaces in reactive transport codes. The fundamental hypothesis underlying such an application of the Nernst-Planck equation lies in the achievement of rapid equilibrium of the diffuse layer composition with the bulk water composition in a representative elementary volume (grid cell). With this assumption, it is possible to describe the diffusive flux within the diffuse layer with an equation having the same form as the Nernst-Planck equation for the bulk porosity, given here in one dimension: where J i is the diffusive flux of the species i, which has a diffusion coefficient D i , a charge z i , a bulk concentration b c i , and an activity coefficient b γ i , and with diffusion potential defined here as the second term of Eq. 1. The superscript b/DL indicates that the parameter must be defined for the bulk water or the diffuse layer water, depending on the diffusive flux of interest. It should be noted that the concentration gradient terms always apply to the bulk concentrations present in the cell, but the DL-DL and bulk-DL fluxes are dependent on the term A i that corresponds to an accumulation/depletion factor of cations and anions in the diffuse layer compared to bulk water [12,22]. A i = 1 in the bulk water porosity and A i = e −z i Fψ DL RT in the diffuse layer porosity where ψ DL is the mean electrostatic potential in the diffuse layer porosity, F is the Faraday constant, R is the ideal gas constant, and T is the temperature. The mean electrostatic potential model corresponds to an approximation of the Poisson-Boltzmann model that makes it possible to calculate effectively average concentrations in the diffuse layer without solving the full Poisson-Boltzmann equation. One might view the overall treatment here as representing an averaging of the Poisson-Nernst-Planck equation once diffusive transport is included. More details on this model can be found in previous publications [9,12,22,23].

Diffusion equation
The 1D Cartesian diffusion equation in a saturated porous medium takes the form where φ is the porosity. When applied to a system where diffuse-layer water and bulk water are present, i.e., a dual continuum model, Eq. 2 becomes where b J i and DL J i are fluxes in bulk and diffuse layer water as defined in Eq 1. Expanding (1) for bulk and diffuse layer water, we obtain where D 0,i is the self-diffusion coefficient of the species i, and τ i is a tortuosity term that is specific to each species and to each type of porosity and type of medium. The flux terms b J i and DL J i are defined at the interface between two adjacent numerical cells. The concentration, porosity, tortuosity, and other terms must be evaluated at the interface based on the values defined at the center of each of the adjacent cells so as to calculate fluxes within and between the different types of water present in the porosity, i.e. bulk water and DL water. The problems related to these calculations, and particularly those associated with cross fluxes (between bulk and DL water) are highlighted in the following paragraphs together with suggestions on how to solve them.

Properties at the interface between two cells
When only bulk water diffusion is considered, it is common to define the value of the product (φτ ) int at an interface between two cells 1 and 2 as the harmonic mean of the product values in the two cells (here the space-grid is assumed to be regular): Whenever the porosity must be divided into bulk and diffuse-layer contributions, it is necessary to include additional rules. A first idea would be to consider separately diffusion from bulk to bulk and from diffuse-layer to diffuse-layer water compartments. The harmonic mean of b φ b τ and DL φ DL τ would thus be written as In this case, if the diffuse layer porosity of one cell approaches zero, the harmonic mean given in Eq. 8 tends to a value of zero, meaning that there is no flux associated with the diffuse layer. However, this result would be unsatisfactory in the special situation depicted on Fig. 1.
The left and right systems can be considered to be strictly equivalent from a conceptual point of view, if the diffuse layers DL 1 and DL 2 have the same properties as the bulk water bulk 1 and bulk 2 respectively (meaning that there is no surface charge associated to the diffuse layer). In the system depicted on the left part of the figure, there would be no flux from bulk 1 to DL 2 because DL (φτ ) int = 0. On the right system, there would be a flux from EDL 1 to EDL 2 that would be, considering our hypotheses, equal to the flux from bulk 1 to bulk 2. The total diffusive flux at the interface between cell 1 and cell 2 in the left system would thus be half of this at the interface in the right system. In this special situation, the straightforward use of Eqs. 7 and 8 would lead to a problem of continuity. One should note that this problem illustrates conditions similar to a boundary condition between a filter (without a diffuse layer) and a clay plug (with a diffuse layer), so it is far from an academic scenario.
In the case in which there are only two flux terms defined in the transport equation (one for the DL, the other for the bulk), there is no "good" answer to this "filter paradox" problem. Depending on how a modeler considers that the interface between a domain without a DL and another domain with a DL is treated, the associated numerical scheme will not be the same.

Rules for flux summation
Two transport terms are not sufficient to describe the systems under investigation if we want to avoid the filter paradox shown in the previous section. A model with three flux terms would be far more accurate as exemplified in Fig. 2:  • A flux term, J a , from bulk 1 to bulk 2 water volumes • A flux term, J b , from DL 1 to DL 2 water volumes, and • A flux term, J c , from bulk 1 to DL 2 water volumes A fourth term should be considered, J d , from bulk 2 to DL 1 water volumes, but only one of the cross-flux terms, J c or J d , is non-zero between cells 1 and 2.
The relative importance of these fluxes must be scaled to the relevant porosity volumes with scaling factors S a , S b , S c , and S d .
We describe here first the system depicted on Fig. 2 where the total porosities of the two adjacent cells are the same: where min and max are the minimum and maximum functions respectively. Note that S a + S b + S c + S d = 1 If the porosities of the two adjacent cells are not the same, φ 1 = φ 2 , it is possible to define scaled porosity terms: Note that, again, S a + S b + S c + S d = 1. Note also that Eqs. 9-12 are a special case of Eqs. 13-17 These equations are those implemented in CrunchClay and in 3Diff. At the interface and for all species i, the flux equation must be discretized. After neglecting the activity coefficient gradient term, as is done in most of the codes, we obtain where "mean" refers to an averaging method at the interface, which depends on the specific code. As noted above, it is common to apply a harmonic mean formula to calculate average porosity and tortuosity, and an arithmetic mean formula to calculate concentrations as in PHREEQC. In systems with irregular grid cells, it is also common to weigh the values by the lengths of the cells on each side of the interface. With the reminder that a product of averages is not numerically equivalent to an average of products, the averaging method chosen may have a non-negligible effect on the results. The flux summation rules implemented in the last version of PHREEQC, however, are different. In version 3.4.5 of the code, the total flux is computed with the following equation:

Code capabilities
PHREEQC and CrunchClay are currently the only two codes that can handle diffusion in the diffuse layer following Eq. 18, although a modified version of FLOTRAN was recently adapted to deal with the diffuse layer on the basis of an alternative conceptual model [20,21]. The numerical resolution scheme of CrunchClay follows a global implicit approach while that of PHREEQC follows a sequential noniterative approach [19]. Averaging methods for interfaces between grid cells with contrasting properties may also be different for the two codes. For these reasons, solving the same problem with the two codes may give slightly different results, and it is important to quantify these differences within the framework of a benchmark exercise. For further insights into the comparison of these to codes, the diffusion equations have been implemented into a new set of Fortran routines, which we refer to here as 3Diff, that can handle test cases in which there are no chemical reactions other than non-specific adsorption in the diffuse layer. The resolution scheme of the diffusion equations in this code is either forward in time and central in space as in PHREEQC (3Diff-Expl), or backward in time and central in space as in CrunchClay (3Diff-Impl). The calculation of the Jacobian matrix is fully analytical in CrunchClay, while it is calculated numerically in 3Diff-Impl. In 3Diff, the properties at the boundary between two cells are averaged with a harmonic mean applied to the products of all parameters except concentration for which a log average was applied. Additional tests were also conducted using a range of combinations of products of harmonic means with 3Diff-Expl. The diffusion equation in PHREEQC includes an additional activity coefficient gradient term, which influences the diffusive flux in the presence of gradient of ionic strength (see Eqs. 1 and 19).

"Simple" 1D Cartesian systems
The first benchmark exercises were carried out on the apparently simple 1D Cartesian geometry shown on Fig. 2 with the following conditions: • The size of each domain was 1.5 mm.
• Only diffusion and reactions took place (no advection).
• The system was closed (no-flux) at the left and right boundaries. • The domains i with i = 1 and i = 2 were homogeneous and contained respectively a solid i having a surface area, ssa i (in m 2 m −3 of porous medium) and a surface charge σ i (in mol m −2 ), and a porosity φ i that was subdivided into bulk and DL porosity (φ b,i and φ DL,i respectively). • Two monovalent tracers, one cation (Nat + ) and one anion (Clt − ), were added with concentrations (c Nat,i , c Clt,i ) that were different in the two domains. Their selfdiffusion coefficients in bulk water were set at D 0,Nat = 1.3 · 10 −9 m 2 s −1 , and D 0,Clt = 2.1 · 10 −9 m 2 s −1 . • The NaCl background electrolyte concentration was set at 0.1 mol kg −1 w , and self-diffusion coefficients of Na + and Cl − in bulk water were set at D 0,Na = 1.3 · 10 −9 m 2 s −1 and D 0,Cl = 2.1 · 10 −9 m 2 s −1 .
• Na + , Cl − , Nat + , and Clt − self-diffusion coefficients in DL water (D 0,DL,Na , D 0,DL,Cl , D 0,DL,Nat , and D 0,DL,Clt respectively) were varied as a function of the calculation case. • No reactivity other than the cation enrichment and anion depletion in the diffuse layer was considered. • In PHREEQC, the water density in the diffuse layer has a hard coded value of 1 kg dm −3 , while in CrunchClay it has the same value as in bulk water, i.e., 0.997 kg dm −3 for a dilute electrolyte at 25 • C. In 3Diff, the water density has a hard coded value of 1 kg dm −3 in the diffuse layer and bulk water.
These systems, the parameters for which are given in Table 1, are labeled "case 1.x". Case 1.0 was a simple test of the diffusion equation. No salinity gradient was present, and the concentration of the NaCl electrolyte background was very large compared to the concentration of the tracers. Under these conditions, the diffusion potential in the Nernst-Planck equation has a negligible effect on the diffusive flux: The results were thus compared to the analytical solution of the Fickian diffusion equation [24]: where l is the total length of the system (l = 0.003 m) and h is the position of the interface at which the concentration of the tracers dropped from C 0 to zero at t = 0 (h = 0.0015 mm). The accuracy of the calculation was checked with the calculation of C/C analytical with C = C numerical − C analytical , where C numerical was obtained from the different codes.
*(Case 1.4): a value of bulk porosity of 10 −3 was the minimum that could be handled by 3Diff-Impl. CrunchClay was able to handle a value of 10 −4 , PHREEQC a value of 10 −6 without issuing warnings, and 3Diff-Expl a value of 10 −8

Complex 1D Cartesian systems
In a second step, we considered systems similar to those given in Table 1, but where tortuosity values depended on the domain, on the type of water (bulk, τ b,1(2) , vs. DL, τ DL,1(2) ), and/or on the type of ions (Table 2). In CrunchClay, the consideration of spatially dependent tortuosity values was straightforward using the keywords Tortuosity and Tortuosity MP for bulk and DL water respectively. The diffusion coefficient of each anion and cation can be defined independently for the bulk and the DL waters as well. In PHREEQC, the tortuosity values were defined cell by cell using a porosity value that was raised to the power 1 defined in the exponent value in the "multi d" input line (see PHREEQC help for more details).
In PHREEQC, the tortuosity term that applies to the DL water can be different from that of the bulk water with the use of a "viscosity" term that can be defined for each cell. However, the diffusion coefficients, which can be defined for each anion and cation, apply to both bulk and DL waters. Strictly speaking, case 2.5 cannot be calculated with PHREEQC, because to simulate this case the "viscosity" of the DL should be set at different values for each of the species.
In the simulation cases 1.1 to 2.5, no salinity gradient was present, and the concentration of the NaCl electrolyte background was again very high compared to the concentration of the tracers. Under these conditions, the diffusion potential in the Nernst-Planck equation has a negligible effect on the diffusive flux. These cases can thus be considered as reference cases to test the correctness of the implementation of the left part of diffusion Eq. (1) over a wide range of geometrical and physical conditions with respect to the bulk and diffuse layer porosity properties. These simulation cases can be run at low computational cost, making them ideal for first benchmarking exercises of implementation of flux summation rules in reactive transport codes.

Diffusion experimental systems
In a third step, we modeled a 1D system representative of a through diffusion experiment setup in which we introduced chemical processes with increasing complexity as well as variations in the representation of the bulk and DL water  1.3 · 10 −9 1.3 · 10 −9 1.3 · 10 −9 1.3 · 10 −9 1.3 · 10 −9 D 0,DL,Cl (m 2 s −1 ) 2.1 · 10 −9 2.1 · 10 −9 2.1 · 10 −9 2.1 · 10 −9 2.1 · 10 −9 D 0,DL,Nat (m 2 s −1 ) 1.3 · 10 −9 1.3 · 10 −9 1.3 · 10 −9 1.3 · 10 −9 0.65 · 10 −9 D 0,DL,Clt (m 2 s −1 ) 2.1 · 10 −9 2.1 · 10 −9 2.1 · 10 −9 2.1 · 10 −9 0.21 · 10 −9 properties. A setup similar to that described in the literature was considered [15,25,26]. We considered a diffusion cell of thickness l = 5 mm and diameter d = 10 mm, filled with montmorillonite with a total porosity of φ tot = 0.72, corresponding to a dry density of ρ dry = 0.8 kg dm −3 . The specific surface area was taken as ssa = 750 m 2 g −1 clay , i.e., 6 · 10 8 m 2 /m 3 , and the surface charge at σ clay = −0.9 mol kg −1 clay , i.e., −1.2 · 10 −6 mol m −2 clay . We considered that 90% of the total porosity was made up of diffuse layer porosity (R DL = φ DL /φ tot = 90%). The tortuosity of the diffuse layer and bulk porosities were set at τ DL = τ b = 0.05. Each end of the diffusion cell was connected to a 10-L reservoir (100% bulk water; τ b = 1), in order to simulate an "infinite" volume. The presence of a filter between the reservoir and the clay material was not considered. In the simulation, a cocktail of tracers diffused out of the reservoirs, including a charge neutral tracer, HTO, with D 0 = 2.1 · 10 −9 m 2 s −1 ; a single-charged cation, Cat + , with D 0 = 1.3 · 10 −9 m 2 s −1 ; a single-charged anion, An − , with D 0 = 2.1 · 10 −9 m 2 s −1 ; and a double-charged cation, Ca 2+ , with D 0 = 0.793 · 10 −9 m 2 s −1 , all taking place in a 0.1 mol L −1 NaCl electrolyte. The flux of the tracers was calculated as a function of tracer accumulation in the down-gradient reservoir. In a first simulation, specific reactions between the cationic tracers and the charged surface were not considered. In a second simulation, the following reactions were considered: 2 > S − + Ca + => S 2 Ca log K S2Ca = 0.5 The same spatial discretization was used for all codes with x = 0.0001 m, but the numerical representation of the system geometry varied depending on the codes. Each reservoir was represented by one grid cell. In PHREEQC, the volume of one grid cell is always 1 dm 3 , and the total surface area available for diffusion from one grid cell to another was thus S diff.tot = 10 m 2 (1 dm 3 / x). In the system studied here, the clay plug had a surface in contact with the reservoir of S clay/res = d 2 × π = 7.854 · 10 −5 m 2 , of which S poro/res = ε tot × S clay/res = 5.655 · 10 −5 m 2 was available for diffusion. The porosity entered in the -water keyword of the SOLUTION block in PHREEQC was then ε PHREEQC = (1 − R DL ) × S poro/res /S diff.tot = (1 − R DL ) × 5.655· 10 −6 for the grid cells in the clay domain, while the porosity (−water) entered in the two reservoirs cells was 10 in order to have an effective reservoir volume of 10 L. The corresponding mass of clay was m clay,cell = S clay/res × x× ρ dry = 6.283 · 10 −3 g. The corresponding surface charge was m clay,cell × σ clay = 5.655 · 10 −6 mol of negative charge and the total clay surface area was s clay,cell = m clay,cell × ssa = 4.7124 m 2 . The volume of the diffuse layer was set at ε tot × R DL by setting the diffuse layer thickness (-donnan keyword) to R DL × S poro/res /S diff.tot /s clay,cell = 1.2 · 10 −9 R DL . In CrunchClay, and 3Diff, the dimensions in the Y and Z dimension were set to 100 m × 1 m, so that the porosity values in the clay plug were defined as one tenth of the value used in PHREEQC. The porosity in the reservoirs was set to a value of 1. In CrunchClay, the surface charge was linked to the presence of a very insoluble fictional phase (bogusite) having an infinitesimally slow dissolution rate, and whose fictional constituent (bogus) had a negligible diffusion coefficient. Its molar volume was set at 50 cm 3 mol −1 and its molar mass at 1000 g mol −1 , i.e., a density ρ bogusite = 20 000 g dm −3 . The specific surface area of bogusite was set at that of the clay, i.e., ssa = 750 m 2 g −1 .

Gradient of ionic strength
In a fourth step, 1D Cartesian systems were re-run in the presence of a gradient of ionic strength created by a difference of NaCl concentration from 0.1 to 0.01 M. This difference of concentration was applied either from the right to the left, or from the left to the right of the systems.

Grid cell size and time step values
In a first step, calculations were carried out with a CrunchClay maximum time step equal to the fixed time step used with the explicit scheme resolution codes. This time step was set at 1 s for a 30-cell system (0.0001 m The results were very similar with the two types of resolution schemes and a very good agreement with the analytical solution was obtained (Fig. 3). The discrepancy between the numerical and the analytical solutions was larger at the down-gradient side of the system than at the up-gradient side. This result was related to the low value of the tracer concentration on the down-gradient side: similar errors in the absolute value of the concentration gave larger relative errors for a low concentration than for a high concentration. In the conditions tested, the relative error obtained with the implicit scheme was larger than with the explicit one. Sensitivity analysis on the maximum time step showed that increasing the maximum time step of calculations with the implicit scheme had a limited effect on the precision of the calculation (Fig. 4). At 2000 s, the relative error in concentration compared to the analytical solution increases towards the right-hand side of the system, i.e., in the down-gradient direction, thus showing a slight underestimation of the diffusive flux obtained with the implicit resolution scheme compared to that obtained with the analytical solution. This effect is a well-known consequence of the use of an implicit scheme that calculates the fluxes on the basis of the concentration gradients present at the next time step. The reverse is true for the explicit scheme, which overestimates the diffusive fluxes. The overestimation obtained with the explicit scheme is however one order of magnitude lower than the underestimation obtained with the implicit scheme. The maximum time step used with the implicit resolution scheme of CrunchClay was set to a value larger than the simulated time, i.e., only the convergence criteria limited the actual time steps in the following simulation cases. For 3Diff-Impl, the maximum time step was set at 10 s. The effect of the slightly faster diffusion coefficient of the anion tracer in case 1.0 is evident from the faster concentration drop in the left part, as seen in Fig. 3 after 1000 s. The charge imbalance between the anion tracer and the cation tracer is achieved through the counter diffusion of background electrolyte ions. However, the difference of concentration between tracer and background electrolyte is very large (11 orders of magnitude difference), thus leading to an absence of significant diffusion of the electrolyte ions.

Flux summations at the interface
All codes tested calculated the same values for the mean electrostatic potential and the ionic concentrations in the diffuse layer at the beginning of the calculations (Table  A1). The diffusion simulation results obtained for case 1.1 with the different codes were also in remarkable agreement (Fig. 5). The fluxes calculated with CrunchClay and 3Diff-Impl were slightly smaller than the fluxes calculated with PHREEQC and 3Diff-Expl, in agreement with the underestimation of the flux with an implicit scheme compared to its overestimation with an explicit scheme.
In case 1.2, the overall fluxes calculated with PHREEQC were slightly larger than the fluxes calculated with 3Diff-Expl, 3Diff-Impl, and CrunchClay (see Fig. 6 where the concentration profiles are flatter with PHREEQC than with the other codes). This discrepancy was not due to the difference of implicit versus explicit resolution scheme, since 3Diff-Expl, 3Diff-Impl, and CrunchClay gave nearly identical results. Hence, the discrepancy must be due to differences in the flux summation rules at the interface between two grid cells. Since the diffusion potential term is negligible, the discrepancy must originate from a difference between the value of mean( b φ + DL φA i ) (PHREEQC) and that of ab,c,d S (ab,c,d) mean(φA i,(ab,c,d) ) (CrunchClay and 3Diff). Indeed, the direct calculation of these averaging terms at the interface between the domains 1 and 2, using the parameters given in Table 1 and a harmonic mean, results in a value 33% higher for the cations with the calculation method of PHREEQC than with the calculation method of CrunchClay and 3Diff. For anions, the difference is 9%, in agreement with the observation of a lesser discrepancy of results with anions than with cations (Fig. 6). By increasing the magnitude of the surface charges σ 1 and σ 2 by a factor 10 (and thus the magnitude of the A i value), it was possible to exacerbate the difference of results given by PHREEQC and the other codes (Fig. 7). This result is in agreement with the direct evaluation of the averaging terms, with a value 145% higher for cations with PHREEQC than with the other codes (43% larger for anions).
The results obtained with the four codes were in very good agreement for cases 1.3 and 1.4 (Figs. S-1 and S-2). In these cases, the averaging terms calculated with the different methods gave values that were similar within 0.6%. This was also true for cases 2.1 to 2.4 (Figs. S-3, S-4 and S-5 and Fig. 8). Case 2.5 was specific to CrunchClay and 3Diff, since PHREEQC does not offer the capabilities to implement it (Fig. 9).

Diffusion setup
The diffusion experimental setup provided a system with large contrasts of porosity, tortuosity, and accumulation factors (A i ) values at the boundary between the clay and the reservoir cells. The results of the codes can be compared to an analytical solution by considering that the reservoirs have an infinite volume. In that case, the instantaneous flux F t (in mol m −2 s −1 ) of a diffusive species entering the lowconcentration reservoir can be calculated as a function of time according to [24,27] where C 0 (in mol m −3 ) is the concentration of the species in the high-and constant-concentration reservoir, D e (in m 2 s −1 ) is the effective coefficient of the species (D e = φτ D 0 ), L (in m) is the thickness of the sample, t is the time (in s), and α (dimensionless) is the rock capacity factor with α = φ + ρ·K D , where ρ (kg dm −3 ) is the dry density of the clay and K D (kg dm −3 ) is a distribution/partition coefficient that is representative of an instantaneous and reversible adsorption process, which is linearly related to the equilibrium concentration. The flux F t can be normalized, using F Nt = F t /C, to make the comparison of tracers easier for those with concentrations that are different in the high concentration reservoir.
In the case of the diffusion of a tracer with a neutral charge (HTO), there was excellent agreement between the results of 3Diff-Expl, PHREEQC and the analytical solution (Fig. 10). The HTO steadystate flux predicted with CrunchClay was larger than the flux calculated with the analytical solution, pointing out a difference in the evaluation of the average term mean(D HTO,0 τ HTO φA HTO ) in Eq. 18. It was possible to reproduce the results of CrunchClay with 3Diff-Expl by changing the averaging method from mean(D HTO,0 τ HTO φA HTO ) = harm(D HTO,0 τ HTO φA HTO ) (referred to as "harm(prod)" method) to mean D HTO,0 τ HTO φA HTO = harm(D HTO,0 τ HTO A HTO ) × harm(φ) (referred to as "prod(harm)" method), where harm represents a harmonic mean. Both The analytical solution for the Nat + flux can be obtained by considering a K D value of zero and an accessible porosity value of φ Nat + = φ b + φ DL A Nat + = 7.32 The discrepancy between the 3Diff-Expl/PHREEQC/analytical solution and CrunchClay was less for Na + than for HTO or Br − (Fig. 11), because the values of D Nat + ,0 τ Nat + A Nat + were similar in the clay domain and in the reservoirs (τ Nat + A Nat + = 0.05 × 11.2 = 0.56 in the clay domain versus τ Nat + A Nat + = 1 in the reservoirs).
The analytical solution for the Ca 2+ flux in the absence of surface complexation reactions can be obtained by considering a K D value of zero and an accessible porosity value of φ Ca 2+ = φ b + φ DL A Ca 2+ = 81.2. CrunchClay predicted a lower flux than the analytical solution, because the value of harm(D Ca 2+ ,0 τ Ca 2+ A Ca 2+ ) × harm(φ) was smaller than the value of harm(D Ca 2+ ,0 τ Ca 2+ φA Ca 2+ ), due to a large accumulation of Ca 2+ in the diffuse layer present in the clay domain. This explanation was verified with Alternative averaging methods were tested with 3Diff-Expl and CrunchClay for the calculation of mean D i,0 τ i φA i in Eq. 18 (prod(harm)); see text for details the results of 3Diff-Expl using the "prod(harm)" averaging method (Fig. 11).
Following these results, an option was added to CrunchClay (MultiplyPorosityTortuosity with flag false or true in the CrunchClay POROSITY block) in order to use the "harm(prod)" method for the product D i,0 τ i φA i (MultiplyPorosityTortuosity in the CrunchClay POROSITY block). With this option set to true, CrunchClay gave results in very good agreement with the other codes and with the analytical solution.
In the presence of surface complexation reactions, the analytical solution for Nat + , Br − , and Ca 2+ fluxes can be calculated by considering the same equation as in the absence of surface complexation, with a non-zero K D value for Nat + and Ca 2+ . In the simulation conditions, these species were present at trace concentration compared to the concentration of surface sites, and compared to the concentration of competing Na + for adsorption. Hence, Nat + and Ca 2+ surface concentrations were linearly correlated with Nat + and Ca 2+ concentration in solution.  The corresponding K D values were calculated according to a batch calculation and resulted in K D (Nat + ) = 4.39 L kg −1 and K D (Ca 2+ ) = 21.5 L kg −1 . Because of the surface complexation reactions, the charge compensated in the diffuse layer was less than in the absence of surface complexation reactions, and the accessible porosity for Br − , Nat + , and Ca 2+ in the analytical solution were 0.187, 3.87, and 22.3 respectively. A very good agreement was found between PHREEQC, CrunchClay, and the analytical solution (Fig. 12). Contrary to the results of the simulation in the absence of surface complexation, a good agreement was found between the Ca 2+ flux results of PHREEQC and of CrunchClay with the "prod(harm)" averaging method. In the presence of surface complexation, the value of τ Ca 2+ A Ca 2+ in the clay domain was 0.05 × 34.4 = 1.7, i.e., a value sufficiently close to the value in the reservoir (= 1), so that harm(D Ca 2+ ,0 τ Ca 2 + A Ca 2 + )× harm(φ) ∼ harm(D Ca 2+ ,0 τ Ca 2 + φA Ca 2 + ).

Influence of a gradient of ionic strength
The presence of a NaCl concentration gradient had a very significant influence on the diffusion of the anion and cation tracers, even in the absence of a diffuse layer. In the simulation case 1.0 with a NaCl concentration gradient from left to right, CrunchClay and 3Diff gave very similar results, showing a slowdown of the anion tracer diffusion, and an acceleration of the cation tracer diffusion (compare Figs. 13 and 3). This effect was due to the diffusion potential term in Eq. 1. Although PHREEQC results showed the same tendency, they were significantly different from those of CrunchClay and 3Diff. In contrast to CrunchClay and 3Diff, PHREEQC does not neglect the gradient of activity coefficient, and this difference likely explains the observed discrepancy. To further test this hypothesis, the gradient of activity coefficient in the PHREEQC calculation was set to zero by forcing the activity coefficient to take a single value of 1. This can be achieved by redefining the Debye-Huckel parameters to a value of zero in a LLNL AQUEOUS MODEL PARAMETERS block. After this correction, PHREEQC gave the same results as CrunchClay and 3Diff (Fig. 13). The same conclusion was reached in the presence of a NaCl concentration gradient from the right to the left (Fig. 14).
In the presence of a diffuse layer, the results obtained with CrunchClay and 3Diff were very similar, while they differed from the results obtained with PHREEQC, whether the activity coefficient gradient term was omitted or not (Fig. 15). In the presence of a diffuse layer, the diffusion potential term calculated using Eq. 18 was not the same as that calculated with Eq. 19. In Eq. 19, the concentrations in the bulk and DL porosity are lumped together, while the contributions of the bulk and DL porosities are calculated separately in Eq. 18. In previous versions of PHREEQC, an equation similar to Eq. 18 was used. Indeed, running the same simulation case with PHREEQC version 3.2.2 gave results in far better agreement with CrunchClay and 3Diff than with version 3.4.5 (Fig. 16).

Conclusions
The benchmark exercises presented in this study provide useful tests to verify the correctness of the implementation of a dual-continuum averaged method to represent the Poisson-Nernst-Planck equation in reactive transport codes for modeling diffusion in nanoporous media. The treatment of the Nernst-Planck equation with a finite difference (volume) numerical method requires the evaluation of parameters averaged at the boundary between grid cells. The benchmarks were specifically designed to be sensitive to the averaging methods used to evaluate these parameters. Only two reactive transport codes, CrunchClay and PHREEQC, are currently available that have the full capability of simulating multicomponent and dual continuum diffusion in bulk and diffuse layer water, electrochemical migration, and reactions including surface complexation. To increase confidence in the comparisons, a third code, 3Diff, was developed with no reaction capability, but with flexibility on the choice of the averaging methods at grid cell boundaries. In addition, analytical solutions of the problems investigated were provided when available. Overall, very good agreement between the simulation results was obtained with the different codes and the analytical solutions. When present, the differences in results were directly traceable to the consideration or omission of the activity gradient term in the Nernst-Planck equation, and to differences in averaging methods at grid cell boundaries. While the proposed benchmark exercises made it possible to highlight these differences, they did not provide any information on the best method that must be implemented in the codes. Indeed, it may be possible that the adequacy of a method is dependent on the characteristics of the system under investigation. More information and modeling at small scales are necessary to decide on the best method(s) at the continuum scale investigated here. At last, the "diffusion set-up" exercise, which simulated the geometry through a diffusion experiment, provided the opportunity to validate the use of reactive transport modeling codes to evaluate diffusion parameters of clay media. We hope that it will encourage an increased use of reactive transport modeling to interpret laboratory diffusion experiments results in a mechanistic and quantitative manner.