MULTICOMPONENT DIFFUSION IN CLAYS

Diffusion in clays depends on a number of factors. Solutes have different diffusion coefficients. Anions are excluded from the diffuse double layer (DDL), and thus diffuse over a smaller surface area than cations and neutral species. Cations accumulate in the DDL or are sorbed, which endows them with higher concentration gradients than other species.
Let's do model experiments to illustrate the effects following Appelo and Wersin, 2007. Two cores are taken from Opalinus Clay and contacted for 1 year with a boundary solution made of porewater to which tritium (HTO), I- and 22Na+ are added as tracers. The concentrations of the tracers in the boundary solution decrease as they diffuse into the core, and we want to see how that depends on the properties of the tracers and the clay. PHREEQC's transport option -stagnant is used to compute the concentrations in the boundary solution.

First, 'classic' PHREEQC and multicomponent diffusion (MCD) are compared, with all solutes having the same tracer diffusion coefficient of HTO (2.24*10-9 m2 / s). The concentrations diminish by diffusion in the sediment. The results are identical for the 3 tracers with the 2 model options (PHREEQC input file opa_col0.phrq):

Next, multicomponent diffusion (MCD) is calculated, with each solute having its own tracer diffusion coefficient.
The first core contains no clay and has neither DDL nor sorption capacity (for this example). The figure shows the decrease of the tracer concentrations with time in the outer solution (PHREEQC input file opa_col1.phrq). HTO declines faster than I-, which in turn decreases quicker than 22Na+. This is in agreement with the order of the tracer diffusion coefficients, HTO > I- > 22Na+ (cf. phreeqd.dat).

The 2nd core has the same porosity, but it contains montmorillonite with an exchange capacity of 1.6 eq/L, compensated by counter ions in the DDL. Half of the pores is occupied by the DDL from which I - is excluded completely with option '-only_counter true'. As a result, the accessible surface for I- is smaller, and this anion disappears slower than in the previous case. 22Na+ accumulates in the DDL, therefore the concentration in the free pore solution is initially smaller than of anions and neutral species. Consequently, the concentration gradient is higher and 22Na+ diffuses faster than in the former case (PHREEQC also calculates the diffusion in the DDL with option -Donnan of keyword SURFACE). The decrease of HTO is the same in all cases. (PHREEQC input file opa_col2.phrq)


Recent diffusion experiments have shown that equal-charged tracers are enriched (or depleted) differently in the DDL, probably because of increased complexation as a result of the decreased dielectric permittivity of water close to a charged surface (Appelo, Van Loon and Wersin, 2010). Furthermore, surface or interlayer diffusion can be significant for very strongly sorbed cations such as Cs+. For modeling these effects, PHREEQC 2.16 has been extended with option '-erm_ddl' (for enrichment in the DDL), and option '-interlayer_d' in TRANSPORT (for diffusion of exchangeable cations).
As an illustration, the 2nd core is modeled by distributing the exchange capacity over interlayer sites X - = 0.96 eq/L and counter ions in the DDL = 0.64 eq/L (40% of the CEC). Furthermore, Cl - and Cs+ are added as tracers (PHREEQC input file opa_col3.phrq).
The DDL enrichment for Cl- and I- is set to 0.6 and 0.3, respectively (SOLUTION_SPECIES). Anion exclusion from the DDL shows up in an increased tortuosity factor: when DDL's overlap in pore constrictions, anions are forced to circumnavigate the obstruction, while HTO and cations can still pass. The higher tortuosity can be modeled by introducing additional diffusional paths that are inaccessible for anions or by lowering the free water diffusion coefficient (-dw in SOLUTION_SPECIES) as in this example.
The different effects of surface (or interlayer) diffusion for Na+ and Cs+ are illustrated in the graph: compare the Na and Cs lines indicated as such (with interlayer diffusion) and with (no_IL_d) (without interlayer diffusion).
Note again that the concentration trend of HTO remains the same.

In the files, option -stagnant of keyword TRANSPORT is used to calculate diffusion over the interfaces of the cells. This enables to model concentration changes in the boundary solution. Mixing factors are calculated with Eqn (128) of the PHREEQC-2 manual (p. 52) for the physical dimensions given in the table below. ε is the porosity, θ2 is the tortuosity-factor, VH2O is the volume of water in a cell, and V0 is the volume of the external solution with tracers in contact with the clay.

-- For classic PHREEQC, the mixing factors are calculated with the cell volumes, as explained with Eqn (128) in the manual.
-- For multicomponent diffusion, the mixing factors must be calculated with a volume of 0.001 m3, and the actual volume of water in a cell is entered when the SOLUTION is defined. The diffusion coefficient and porosity that were used for the mixing factors must be entered with option -multi_D in the input file. PHREEQC adapts the mixing factor for individual solutes by the ratio of the tracer diffusion coefficient and the default diffusion coefficient given with -multi_D. For each cell, only mixing factors with higher numbered cells need be entered. If PHREEQC warns that concentrations turn negative, the mixing factor must be reduced, which is generally possible by letting PHREEQC take time-substeps:
-time 5e5 3 # 5e5 seconds are substepped into 5e5 / 3 = 1.666e5 sec.

References
Appelo, C.A.J. and Wersin, P., 2007. Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide and sodium in Opalinus Clay. Env. Sci. Technol., 41, 5002-5007.
Appelo, C.A.J., Van Loon, L.R. and Wersin, P., 2010. Multicomponent diffusion of a suite of tracers (HTO, Cl, Br, I, Na, Sr, Cs) in a single sample of Opalinus Clay. Geochim. Cosmochim. Acta 74, 1201–1219.

«—