This example illustrates how PHREEQC version 3 can simulate a diffusion experiment, as is now often performed for assessing the properties of a repository for nuclear waste in a clay formation. A sample is cut from a core of clay, enveloped in filters, and placed in a diffusion cell (see Van Loon and others, 2004, for details). Solutions with tracers are circulated at the surfaces of the filters, the tracers diffuse into and out of the clay, and the solutions are sampled and analyzed regularly in time. The concentration changes are interpreted with Fick’s diffusion equations to obtain transport parameters for modeling the rates of migration of elements away from a waste repository. Transport in clays is mainly diffusive because of the low hydraulic conductivity, and solutes are further retarded by sorption (cations) and by exclusion from part of the pore space (anions).

For calculating diffusion, one needs to account for the different diffusion coefficients of the tracers, the effect of the filters, and the properties of the clay. Figure 22 presents a schematic diagram of the processes in the clay (Appelo and others, 2010). The pores in the clay are lined with clay minerals that have a negative surface charge. The charge is partly neutralized by cations that are bound to the surface and partly by the electrostatic double layer that extends some distance into the pore. The double layer contains an excess of cations (counter-ions, in general) and a deficit of anions (co-ions, in general). In swelling clay minerals, such as montmorillonite, cations in the interlayer space also neutralize part of the negative charge. In comparison with concentration gradients that drive diffusion in free (uncharged) pore water, the gradient for a cation is magnified in the double layer, whereas the gradient for an anion is diminished, and the gradient for a neutral species remains the same. The charge in the double layer lowers the dielectric permittivity of the water, which enhances the ion-association of cations and anions into neutral species. Also, the viscosity of water may be higher than in free pore water. Where double layers overlap in pore constrictions, anions are impeded and forced to take longer paths, whereas cations and neutral species pass through unhindered.

PHREEQC has capabilities to model many of these diffusion processes. An averaged double layer composition can be calculated by using the identifier -Donnan in keyword SURFACE, which, in essence, neutralizes the surface charge. Solute species can be assigned an enrichment factor in the Donnan pore space to emulate the additional concentration change as a result of different ion association. For diffusion, the viscosity can be set differently with respect to free pore water. In this example, all these properties are adjusted for the tracers tritium (HTO), Na+, Cs+ and Cl- to simulate experimental results.

The experiments to be modeled were done in a radial diffusion cell, shown in figure 23. The radial cell enables measurement of diffusion parallel to the bedding plane of the clay (Van Loon and others, 2004). A solution with tracers is circulated at the surface of the inner filter, and another solution with the same major ions, but without tracers, contacts the surface of the outer filter. The latter solution is removed regularly and analyzed for the tracer that has diffused through the filters and the clay. Meanwhile, the outer solution is replaced by fresh solution, thus, maintaining tracer concentrations near zero in the outer solution. In the simulations, low solubility phases were defined to maintain small concentrations of tracers in the outer solution. The moles of precipitates of these phases give the accumulated mass that has diffused through the column.

For a linear column, PHREEQC calculates diffusion automatically by using the parameters entered with identifier -multi_D in keyword TRANSPORT. However, diffusion in the experiment is radial, and the filters have different properties than the clay. Also, the boundary solutions for the default (linear) column have a constant composition, whereas we want to know the concentration changes in these solutions during the diffusion experiment. All these experimental details can be simulated for a stagnant column by defining mixing factors among the individual cells with keyword MIX.

The mixing factors can be derived from Fick’s diffusion equations, and , which transform to finite differences for an arbitrarily shaped cell j:

where
is the concentration in cell *
j*
at the current time (mol/m3),
is the concentration in cell *
j*
after the time step, Dw is the tracer diffusion coefficient (m2/s),
is the time step (s), i is an adjacent cell,
is the porosity over the interface of cells i and j (unitless), Gij is the geometrical factor that corrects for tortuosity of the porous medium (unitless), hij is the distance between midpoints of the cells (m), Aij is the shared surface area of cells i and j (m2), Vj is the volume of cell j (m3), and fbc is a factor for boundary cells (unitless). The summation is for all cells (up to n) adjacent to j. When hij is equal for all cells, a central difference algorithm is obtained that has second-order accuracy [O(h)2] (Appelo and others, 2008). It is therefore advantageous to make the grid regular. However, the same accuracy is achievable for a heterogeneous domain, even with widely variable gridsize, if the harmonic mean of the parameters in
is used. These parameters together, translate the tracer diffusion coefficient Dw into the effective diffusion coefficient De.

The harmonic mean can be derived in general (omitting the activity coefficient to simplify the formulas) as follows. The fluxes inside cells i and j and over the interface of the two cells must be the same and are given by:

where h is the cell length (m), and cij is the concentration at the interface. Substituting = gi, in equation 46 and similarly in equation 47, the two equations can be combined while eliminating the concentration cij to give:

Multiplying the flux by the surface area, the time step, and the boundary factor, and dividing by the volume of the cell for which the concentration change is calculated (here, cell j), the mixing factor is obtained:

where fbc is 2 for cells in contact with a constant concentration cell, and 1 otherwise. When calculating mixf, Vj is set to 10^{
-3}
m^{
3}
, and for Dw the default diffusion coefficient is used, entered with identifier -multi_D. In multicomponent-diffusion mode, PHREEQC adapts the volume Vj (entered as 10^{
-3}
m^{
3}
in equation 49), to the actual volume of water in cell j, and multiplies the mixing factor for each solute species a by the ratio of the tracer diffusion coefficient for a and the default diffusion coefficient (Dw, a / Dw).

To avoid numerical oscillations, it is necessary that mixf < 0.5, which can be realized by limiting the time step. This maximum permissible time step is usually determined by the cell with the smallest volume (in the model) and the fluxes of the proton, which normally has the highest tracer diffusion coefficient. However, if the proton concentration is sufficiently buffered by alkalinity or other species, and the solutions are uniform except for the tracers, the tracer with the highest Dw may be selected for calculating the maximum permissible time step. If, nevertheless, time steps are too large, PHREEQC warns that negative concentrations are calculated (and the program stops because the system reached an infeasible state), the time step can be subdivided. For example, -time_step 5e2 3 will subdivide the time step of 500 s in three equal ones of 166.7 s.

The input file in table 60 defines the physical and chemical properties of the clay pore space and writes the mixing factors for diffusional transport in the filters and the clay. The filters used by Van Loon have a geometrical factor of 4 for all the tracers, whether charged or not (Glaus and others, 2008). In the clay, the geometrical factor for HTO is 6.2. For cations, the geometrical factor appears to be 2 to 4 times smaller than for tritium. However, in this example (example 21), the same geometrical factor is used for both HTO and the cations; the smaller apparent geometrical factors are modeled by subdividing the pore space into free pore water and double-layer water. The concentrations of cations are higher in double-layer water than in free pore water, and hence, the concentration gradients of cations are also higher, which, as noted previously enhances their diffusion. In models that do not account for this physical aspect of the clay pore space, the faster diffusion is treated by diminishing the geometrical factor.

For anions, the geometrical factor is about 1.5 times larger than for tritium. This larger factor is related to narrowing of the pores, where overlapping double layers decrease the accessible porosity for anions and obstruct their passage (as noted above). The model accounts for the observed, smaller accessible porosity for anions relative to tritium and cations by the calculation of anion exclusion in the double layer.

The file starts with the tracer species, where, for the monovalent tracers 22Na+ (=“Na_tr+”) and Cs+, an enrichment factor for the double layer is entered with -erm_ddl in the SOLUTION_SPECIES data block. The enrichment is related to increased complexation of the polyvalent cations in the low dielectric permittivity of the double layer, and thus, more charge must be counterbalanced by monovalent cations. Sorption of Cs+ is much stronger than that of Na+, which is modeled by two surface complexes and one exchange reaction with large constants. The constants are based on the measured adsorption isotherm for Opalinus Clay, but may be generally applicable because they are associated primarily with strong sorption sites on illite. Next, the file writes SOLUTION 0-2 for a regular column, followed by SOLUTION 3,

which forms the start of the stagnant column and circulates at the inner filter of the diffusion cell. This solution contains the tracers and the amounts of water used in the experiment. The lines can be uncommented one by one to run the file successively with the different tracers and amounts of water. When SOLUTION 3 is calculated, USER_PUNCH is processed to write a SELECTED_OUTPUT file named radial, which contains the SOLUTION definitions for the cells in the filters and the clay, the mixing factors, and the TRANSPORT and USER_GRAPH data blocks. The Basic lines in USER_PUNCH do the following tasks:

- Lines 1-4 define variables that facilitate printing of special symbols used in printing the output--end of line, comment sign, and semicolon--and define π, which is used to write definitions for the radial configuration of the experimental cell.
- Lines 10-150 define the dimensions of the experimental cell and properties of the filters and the clay that have been measured and, thus, should be considered as constant.
- Lines 160-310 fill arrays with the names of the tracers, the experimental times, and the scale factors for the graphs.
- Lines 350-480 give model parameters that may be varied to simulate the experiments, and can be changed to check the diffusion model. Typically, for checking the numerics, the number of cells for the filters (variables nfilt1 and nfilt2) and for the clay (nclay), and the time step (punch_time) can be altered without affecting the calculated results. It also is interesting to set nfilt1 and/or nfilt2 to zero and inspect the effects that the filters have on the fluxes. (The program will fail if nclay is set to zero.) Values of the other model parameters were derived from the geometrical factors obtained by fitting the measured tracer diffusion curves (Appelo and others, 2010). By adjusting these parameters, the Basic program and the functioning of PHREEQC can be verified. For example, f_free sets the fraction of free pore water, partitioning the pore space in free and double layer water. Values may range from 0.01 to 1. The variable has little effect on tritium which, as an uncharged species, diffuses equally quickly in free and double-layer water (if the latter is given the viscosity of free pore water). However, the variable has a major effect on the diffusion of 36Cl-, which diffuses more quickly in the free pore water. The parameter f_DL_charge partitions the Cation Exchange Capacity (CEC) over the double layer and exchange sites. Increasing its value does not affect the diffusion of tritium, but decreases the flux of Cl- and increases the flux of Na+.
- Line 490 checks which tracer is present to set the experimental times and scale factors for the graphs.
- Lines 520-540 define clay pore-water composition.
- Lines 550-630 define the solid phases in which the tracers precipitate in the outer solution. The moles of this phase will record the amounts that have diffused. The phases have such a low solubility that the tracer concentration is essentially zero.
- Lines 650-1100 write the solutions for the filter cells and the clay, with radially increasing amounts of water. For the clay, keyword SURFACE is used to define the moles of the surface sites of Su_, Su_ii and Su_fes, for the double layer, and the sites on illite that sorb the alkaline cations with intermediate and high strength, respectively. The fixed sites of the Cation Exchange Capacity (CEC) are defined with EXCHANGE.
- Lines 1200-1320 write the external SOLUTION and the EQUILIBRIUM_PHASES in which the tracers are captured.
- Lines 1400-1900 calculate and write the mixing factors as explained above in equations 44-49. First, the maximum time step is derived from one of the following: the innermost filter cell, at the transition of the inner filter to the clay, the innermost clay cell, or the outer filter cell. The time step is decreased when it is larger than desired by punch_time. With this time step, the mixing factors are calculated for each cell and written to the file, taking care of the heterogeneities at physical boundaries and the radial outline of the field.
- Lines 2000-2520 write data blocks for TRANSPORT and USER_GRAPH. The experimental data (L.R. Van Loon, Paul Scherrer Institute, written commun., 2011) are plotted (“-plot_tsv_file file_name”) together with the calculated accumulated mass in the outer solution and the flux, obtained by dividing the mass that has accumulated by the time interval and the outer surface area of the clay.
- Lines 3000-3390 contain the subroutine to write TRANSPORT and USER_GRAPH data blocks for plotting concentration profiles in the filters and the clay.

The experiments with the tracers 22Na and Cs can be calculated together, because the tracer solutions contained identical amounts of water. Thus, when the 22Na profile has been plotted after 45 days, another TRANSPORT block is written to calculate diffusion of Cs during 1,000 days in total.

Following the END after USER_PUNCH, the file radial is loaded in the input file with INCLUDE$ and then processed. Before the file is included, writing to file radial is suspended (PRINT; -selected_output false), thus avoiding rewriting the file over and over again, each time a solution is calculated.

Three transport calculations with different tracers can be run with the input file (1) HTO, (2) 36Cl, and (3) 22Na and Cs together. To perform one of these three calculations, uncomment the appropriate line defining the tracer(s) in SOLUTION 3. Note that the length of the Cs experiment is about 1,000 days compared to less than 50 days for the other experiments; consequently, the Cs calculation is 20 times longer than the others.

The results of the three calculations are shown in figures See Mass outflow (red) and corresponding flux (green) by diffusion through the radial cell for (A) HTO, (B) 36Cl-, (C) 22Na+, and (D) Cs+. Lines are modeled; symbols indicate measured data. and See Concentration profiles in the filters and free pore water in the clay (blue), and total moles per kilogram water in the double-layer water in the clay (red) at the end of the calculations in the radial cell for (A) HTO, (B) 36Cl-, (C) 22Na+, and (D) Cs+. Vertical lines indicate the positions of the filters.. As shown, the arrival times of the tracers that accumulate in the outer solution are delayed by the storage in pore water and the sorption on minerals in the clay. The delay increases in the order 36Cl- < HTO < 22Na+ < Cs+ (fig. 24). The total storage can be obtained from the graphs by extrapolating the straight-line segment of the accumulated mass to time zero (fig. 24) and reading the value from the secondary Y-axis (a negative number, because mass is lost). The flux, the derivative of the mass with time, shows that the accumulation of HTO, and of Cs+ in particular, decreases during the experiment (fig. 24) because the concentration is diminishing in the tracer solution. The volume of the solution with HTO, the first experiment performed, was relatively small (0.2 L) and insufficient to maintain a constant concentration for the inner solution. Sorption of Cs+ is so strong (more than 99.5 percent of Cs+ resides in the solid phase) that 1 L of inner solution simply contains insufficient mass to fill all the sorption sites on the clay. Figure 25 shows on the primary Y axis the concentration profile at the end of the diffusion periods, including concentrations in the inner filter, the free pore water of the clay, and the outer filter. The figure shows on the secondary Y axis the total moles, which includes moles sorbed and moles dissolved in free and double layer water, expressed per kilogram of total water. The two concentrations are the same for HTO (the concentration of HTO is the same in free and in double layer water), the free pore-water concentration of Cl- is twice the total concentration (the concentration of Cl- in the double layer is smaller than in free pore water), and the free pore-water concentrations of Na+, and particularly of Cs+, are smaller than the total concentrations because the cations have a higher concentration in the double layer than in free pore water, and because cations are sorbed on the surface and exchange sites.

The model simulates the experimental results very well (fig. 24), except for Cs+. The calculated arrival time of Cs+ is almost 100 days later than observed and then the mass accumulates too slowly. This behavior of Cs+ has been found in many similar experiments. It has been modeled by increasing the diffusion coefficient and decreasing the sorption capacity for Cs+ relative to batch experiments; these adjustments can be easily incorporated in this example.

The diffusion of Cs+ can be increased by setting interlayer diffusion to true in Line 430:

430 interlayer_D$ = 'true'

(It is of interest to see the different effects of interlayer diffusion on 22Na+ and on Cs+.)

The results, shown in figure 26, illustrate that the arrival time of Cs+ can be matched by increasing diffusion, but that the mass accumulates too quickly. Another option for reducing the delay of the tracer arrival is to decrease the sorption of Cs+, either by decreasing the moles of surface and exchange sites or by decreasing the complexation constants, or both. By adjusting both the sorption and the diffusion coefficient, it may be possible to simulate the experimental data for Cs+. However, it remains difficult to explain why the sorption capacity is different between batch and diffusion experiments for Cs+, but not for the other cations.

Alternatively, and in line with the heterogeneous distribution of Cs+ in the clay after the experiment, the relatively fast arrival time can be modeled with a dual-porosity structure in which the pore space is subdivided into continuous and stagnant pores that can exchange by diffusion. The continuous pores allow Cs+ to move more rapidly through the clay than is calculated for a homogeneous medium, depending on the proportion, the flow velocity, and the diffusional exchange with the stagnant pores. With equation 49, and some effort, the dual-porosity structure can be introduced in the Basic program. Otherwise, the C program that is given as supplementary information in Appelo and others (2010) can be used. Similarly to the Basic program, the C program writes a complete PHREEQC input file for diffusion of Cs+, but in a dual porosity clay; in addition, the program permits distributing the surface and exchange sites differently between the stagnant and continuous pores.