# Resolving experimental biases in the interpretation of diffusion experiments with a user-friendly numerical reactive transport approach

Christophe Tournassat, Carl I Steefel, Patricia M Fox, Ruth M Tinnacher

► To cite this version:

Christophe Tournassat, Carl I Steefel, Patricia M Fox, Ruth M Tinnacher. Resolving experimental biases in the interpretation of diffusion experiments with a user-friendly numerical reactive transport approach. Scientific Reports, 2023, 13 (1), pp.15029. 10.1038/s41598-023-42260-5. insu-04910993

HAL Id: insu-04910993

https://insu.hal.science/insu-04910993v1

Submitted on 24 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Resolving experimental biases in the interpretation of diffusion experiments with a user-friendly numerical reactive transport approach

Christophe Tournassat<fnref n="1" />,²✉, Carl I. Steefel<fnref n="1" />, Patricia M. Fox<fnref n="1" /> & Ruth M. Tinnacher³

The reactive transport code CrunchClay was used to derive effective diffusion coefficients ($D_e$), clay porosities ($\epsilon$), and adsorption distribution coefficients ($K_D$) from through-diffusion data while considering accurately the influence of unavoidable experimental biases on the estimation of these diffusion parameters. These effects include the presence of filters holding the solid sample in place, the variations in concentration gradients across the diffusion cell due to sampling events, the impact of tubing/dead volumes on the estimation of diffusive fluxes and sample porosity, and the effects of O-ring-filter setups on the delivery of solutions to the clay packing. Doing so, the direct modeling of the measurements of (radio)tracer concentrations in reservoirs is more accurate than that of data converted directly into diffusive fluxes. While the above-mentioned effects have already been described individually in the literature, a consistent modeling approach addressing all these issues at the same time has never been described nor made easily available to the community. A graphical user interface, CrunchEase, was created, which supports the user by automating the creation of input files, the running of simulations, and the extraction and comparison of data and simulation results. While a classical model considering an effective diffusion coefficient, a porosity and a solid/ solution distribution coefficient ($D_e-\epsilon-K_D$) may be implemented in any reactive transport code, the development of CrunchEase makes it easy to apply by experimentalists without a background in reactive transport modeling. CrunchEase makes it also possible to transition more easily from a $D_e-\epsilon-K_D$ modeling approach to a state-of-the-art process-based understanding modeling approach using the full capabilities of CrunchClay, which include surface complexation modeling and a multi-porosity description of the clay packing with charged diffuse layers.

A characterization of radionuclide diffusion parameters in clayey materials is essential for the performance assessment of natural and engineered barrier systems in radioactive waste storage concepts<fnref n="1" />,². These parameters include the effective diffusion coefficient ($D_e$), effective porosity ($\epsilon$), and adsorption distribution coefficient ($K_D$), which are characteristic of both clayey materials and the radionuclides under investigation. Values of $D_e$, $\epsilon$, and $K_D$ are empirically determined in the framework of the Fickian diffusion theory. However, these parameters typically become lumped together when a range of physical and chemical processes are present that can be unraveled only by applying multi-scale characterization and modeling techniques³. The simplicity of the $D_e-\epsilon-K_D$ modeling approach, as well as its robustness for cases specific to radioactive waste storage concepts (constant far-field conditions in space and time, trace concentration levels for radionuclides), is in stark contrast to the complexities of the underlying physical and chemical processes. However, the $D_e-\epsilon-K_D$ modeling approach is the current norm in barrier performance evaluation methodology due to its simplicity.

Because the $D_e-\epsilon-K_D$ modeling approach is empirical in nature, parameter estimation is most often carried out based on the fitting of experimental data. For this purpose, extensive diffusion datasets have been acquired over

<footnote n="1">¹Earth and Environmental Sciences Area, Lawrence Berkeley National Laboratory, Berkeley, CA, USA. ²Institut des Sciences de la Terre d'Orléans, Université d'Orléans-CNRS-BRGM, Orléans, France. ³Department of Chemistry and Biochemistry, California State University East Bay, Hayward, CA, USA. ✉email: ctournassat@lbl.gov; christophe.tournassat@univ-orleans.fr</footnote>
the past decades. The interpretation of such data relies on solving Fick's first and second laws of diffusion. Fick's first law states that the diffusive flux of an aqueous solute ($J$ in mol m⁻² s⁻¹) is proportional to its concentration gradient (with $c$ in mol m⁻³ water), and its effective diffusion coefficient $D_e$ (in m² s⁻¹)⁴:

$$J = -D_e \frac{\partial c}{\partial x} \qquad (1)$$

The effective diffusion coefficient can be further expressed as a function of the effective porosity $\varepsilon$ (in m³ water m⁻³ medium), a geometrical factor $\tau$ (dimensionless), hereafter referred to as tortuosity, and the self-diffusion coefficient of the aqueous solute $D_0$ (in m² s⁻¹):

$$D_e = \varepsilon \tau D_0 \qquad (2)$$

Fick's second law is derived from the mass conservation law:

$$\frac{\partial C_{tot}}{\partial t} = -\frac{\partial J}{\partial x} \qquad (3)$$

where $C_{tot}$ is the total concentration of the element of interest in the porous medium (in mol m⁻³ medium), including aqueous species, surface species, and species in solid phases. The rock capacity factor $\alpha$ (in units of m³ water m⁻³ medium) relates the total concentration of the (radio)tracer in the porous medium to its aqueous concentration only:

$$\alpha = \frac{C_{tot}}{c} \qquad (4)$$

Rearranging Eq. (4) to $\alpha c = C_{tot}$, and combining this expression with Eqs. (1) and (3) yields:

$$\frac{\partial \alpha c}{\partial t} = \frac{\partial}{\partial x} \left( D_e \frac{\partial c}{\partial x} \right) \qquad (5)$$

We will now make two assumptions to relate the rock capacity factor $\alpha$ (in m³ water m⁻³ medium) to the effective porosity in water-saturated condition $\varepsilon$ (in m³ water m⁻³ medium), the bulk dry density $\rho_d$ (in kg solid m⁻³ medium) and the adsorption distribution coefficient $K_D$ (in m³ water kg⁻¹ solid). First, we assume that the tracer is present only in solution and on surfaces (with $C_{surf}$ in mol kg⁻¹ solid) so that $C_{tot}$ (in mol m⁻³ medium) can be expressed as:

$$C_{tot} = c\varepsilon + C_{surf}\rho_d \qquad (6)$$

Secondly, we assume that the tracer surface concentration $C_{surf}$ is linearly related to the aqueous concentration through $C_{surf} = cK_D$ (with $K_D$ in m³ water kg⁻¹ solid). Based on this assumption,

$$C_{tot} = c\varepsilon + cK_D\rho_d, \qquad (7)$$

and based on Eq. (4), $\alpha$ can be expressed as:

$$\alpha = \varepsilon + \rho_d K_D \qquad (8)$$

When interpreting diffusion data, the adsorption distribution coefficient is commonly assumed to be representative of an instantaneous and reversible adsorption process. If it is further assumed that the medium is homogeneous, and hence $D_e$, $\varepsilon$, $\rho_d$ and $K_D$ are independent of $x$, then Eq. (5) reduces to:

$$\frac{\partial c}{\partial t} = \frac{D_e}{(\varepsilon + \rho_d K_D)} \frac{\partial^2 c}{\partial x^2} \qquad (9)$$

In practice, predictions based on Eq. (9) for a specific set of geometry, initial and boundary conditions are compared to experimental data obtained from a setup representative of the same conditions. Through-diffusion experiments with constant boundary conditions (Fig. 1A), i.e., a constant concentration gradient across the diffusion cell over the course of the experiment, allow for an independent characterization of $D_e$ and $\alpha$ in a single experiment. This is accomplished by using the stationary ($D_e$) and transient ($\alpha$) regime of measured diffusive flux data⁵.

Under ideal conditions, the boundary and initial conditions for $c(x, t)$ in the through-diffusion problem are given by:

$$\begin{aligned} c(x, 0) &= 0 \quad \text{for } 0 < x < L_s \\ c(0, t) &= c_0 \quad (\text{high-concentration reservoir}) \\ c(L_s, t) &= 0 \quad (\text{low-concentration reservoir}) \end{aligned} \qquad (10)$$

where $L_s$ is the clay sample thickness in the diffusion cell (m), $c_0$ (in mol m⁻³) is the aqueous concentration of the element of interest (the tracer, whether radioactive or not) in the high and constant concentration reservoir, and $t$ is time (in s). The instantaneous flux $F_i$ (in mol m⁻² s⁻¹) of a diffusive species entering the low-concentration reservoir can then be calculated with an infinite series as a function of time according to⁴:
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic of an ideal (A) and actual (B) through-diffusion experimental setup with the goal of achieving constant boundary conditions, i.e., with a constant concentration gradient across the diffusion cell at steady state. (C) 1D diffusion geometry. (D,E) 2D-axisymmetric cylindrical geometries. Tracer cannot diffuse through the plain O-ring of the filter support (D), making the diffusional path 2D. The cross-sectional area available for diffusion in the filters is equal to $(\frac{d_{filter}}{2})^2 / (\frac{d_{sample}}{2})^2$ of the cross-sectional area of the sample.</figcaption></figure>

$$F_t = \frac{c_0 \cdot D_e}{L_s} + \frac{2 \cdot \alpha \cdot c_0 \cdot L_s}{\pi^2} \sum_{j=1}^{\infty} \frac{(-1)^j D_e \cdot j^2 \cdot \pi^2}{j^2 L_s^2 \cdot \alpha} \exp\left(-\frac{D_e \cdot j^2 \cdot \pi^2 \cdot t}{L_s^2 \cdot \alpha}\right) \quad (11)$$

The numerical implementation of Eq. (11) is a truncated version of the infinite series, in which the maximum j value can be set high enough to limit the truncation error value to machine precision.

In real experiments (Fig. 1B), the concentration in the high-concentration reservoir is kept as constant as possible by using a large fluid volume. In contrast, the concentration in the low-concentration reservoir is kept as close to zero as analytically/practically possible by replacing the low-concentration reservoir periodically. The experimental tracer diffusive flux ($F_{exp}$) is then evaluated from the tracer accumulation in the low-concentration reservoir. The cumulated total amount of tracer in the low-concentration reservoir at time $t_n$, $Q(t_n)$ (in mol), is obtained with:

$$Q(t_n) = \sum_n c_L(t_n) V_L(t_n) \quad (12)$$

where $c_L(t_n)$ and $V_L(t_n)$ are respectively the tracer concentration (in mol m⁻³), and the volume (in m³) of the low-concentration reservoir at sampling time $t_n$ (in s). The experimental tracer diffusive flux ($F_{exp}$) is then evaluated by using one of the following backward or central numerical derivative expressions:

$$F_{exp}(t_n) = \frac{Q(t_n) - Q(t_{n-1})}{(t_n - t_{n-1})A} \quad (13)$$

or

$$F_{exp}(t_n) = \frac{Q(t_{n+1}) - Q(t_{n-1})}{(t_{n+1} - t_{n-1})A} \quad (14)$$

where A is the cross-sectional area of the sample packed into the diffusion cell (in m²). If tubing is used to connect the reservoirs and filters holding the solid sample in place, then $V_L$ also includes the corresponding tubing volumes.
In principle, the least square fitting of $D_e$ and $\alpha$ by matching the prediction of $F_t$ based on Eq. (11) with experimental results for $F_{exp}(t_n)$ calculated from Eq. (13) or (14) is a very convenient way to derive diffusion parameters. The simultaneous fitting of diffusion parameters of non-sorbing tracers for which $K_D=0$ and $\alpha=\varepsilon$ (e.g. tritiated water HTO, or $^{36}$Cl$^-$), and sorbing tracers, for which $\alpha=\varepsilon+\rho_d K_D$, makes it possible to obtain a comprehensive and operational set of $D_e-\varepsilon-K_D$ values.

In practice, however, experimental conditions do not often meet the idealized conditions. First, filters at both ends of the sample and tubing connections between the filters and the reservoirs are usually required (Fig. 1B). This complicates the calculations since filters have diffusion properties that differ from those of the solid samples under investigation. Yet, analytical solutions to the diffusion equation can still be obtained under these conditions⁶:

$$F_t = \frac{c_0 \cdot D_e}{L_s(1 + \beta k)} + 2 \cdot \alpha \cdot c_0 \cdot L_s \sum_{j=1}^{\infty} \frac{D_e \cdot \omega_j^2}{L_s^2 \cdot \alpha} \exp \left( - \frac{D_e \cdot \omega_j^2 \cdot t}{L_s^2 \cdot \alpha} \right) \quad (15)$$

where $\omega_j$ is the jth root of the equation (with $\omega_j \neq 0$):

$$a\sin(\omega_j(1+k)) + b\sin(\omega_j(1-k)) - c\sin(\omega_j) = 0 \quad (16)$$

with:

$$\begin{align}
\beta &= \sqrt{\frac{\varepsilon D_e}{\varepsilon_f D_f}} \nonumber \\
k &= 2r \frac{L_f}{L_s} \nonumber \\
r &= \sqrt{\frac{\varepsilon_f D_e}{\varepsilon D_f}} \tag{17} \\
a &= \frac{1}{4}(\beta + 1)^2 \nonumber \\
b &= \frac{1}{4}(\beta - 1)^2 \nonumber \\
a &= \frac{1}{2}(\beta^2 - 1) \nonumber
\end{align}$$

where $\varepsilon_f$ is the porosity of the filters, $D_f$ is the effective diffusion coefficient through the filter, and $L_f$ is the thickness of the filter. Second, the concentration in the high-concentration reservoir can decrease significantly over time in the case of a sorbing tracer. Maintaining a constant concentration level is technically difficult to implement since it requires a very large reservoir volume, potentially a problem since many radiotracers are expensive. Hence, this ideal condition is seldom achieved in experimental practice⁷, and analytical expressions of the diffusive flux must be adapted in conditions where concentrations change significantly in the reservoirs⁸. Third, diffusion equations should, in principle, be solved in 2D- axisymmetric cylindrical coordinates if filters with solid O-rings are used (Fig. 1D,E). Fourth and last, the periodic replacement of low-concentration reservoir solutions generates changes in the concentration gradient across the diffusion cell over the course of the experiment. This further influences the overall diffusive flux as a function of sampling events⁹. No analytical solution has been proposed yet to handle these last two experimental features.

An accurate evaluation of diffusion parameters requires that these experimental biases are specifically taken into account in the fitting procedure, or are mitigated with special experimental procedures¹⁰. Experimental mitigation strategies include the development of flushed filters in which the reservoir solution is circulated inside of the filters (as opposed to filters at which the solution is circulated at the surface opposite to the sample), thus enabling a more homogeneous solute concentration in the filter and the connected tubing and reservoir⁹. However, all biases and their combinations cannot always be addressed using experimental mitigation strategies or analytical solution development. Thus, ultimately, the estimation of diffusion parameters often relies on approximate analytical solutions¹¹ or fully numerical solution strategies⁹. A literature survey reveals that many numerical tools have been developed to interpret diffusion experiments⁹,¹²⁻¹⁴. However, these tools lack either generality concerning the series of identified biases, easy accessibility by the scientific community, or both.

Reactive transport codes have now been developed for over four decades to model complex interactions in natural systems, while coupling advection, dispersion, diffusion, and sorption reactions¹⁵. These codes can numerically and iteratively consider an arbitrary number of kinetic or equilibrium reactions together with flow and diffusion equations. Their robustness and accuracy have been tested repeatedly and successfully in benchmark exercises covering a range of applications including diffusion problems¹⁶,¹⁷. Most of these codes are freely available to the scientific community, and in many cases open-source as well¹⁵. It is thus very surprising that reactive transport models have not been used more extensively to interpret experimental diffusion data. From our own experience, we can identify at least five assumptions that may explain this observation: (1) reactive transport codes are perceived as too difficult to handle by experimentalists; (2) experimentalists and reactive transport modelers are two separate scientific communities with limited communication; (3) reactive transport calculations are too slow to efficiently handle the fitting of diffusion parameters; (4) numerical algorithms in
reactive transport codes are not sufficiently accurate compared to analytical solutions; (5) reactive transport codes are not adapted to handle a $D_e-\epsilon-K_D$ modeling approach because they are focused on process understanding.

With the present contribution, we would like to refute assumptions (3) to (5) that we consider are incorrect and overcome assumptions (1) and (2) by providing the community with a user-friendly, fast, generic, reliable, and free tool, namely the CrunchEase user interface to be used with the CrunchClay code for use interpreting diffusion data with $D_e-\epsilon-K_D$ conceptual models. We will also show that CrunchEase in conjunction with CrunchClay can be used as a first modeling step that can help to set up a more complex modeling approach that takes into account state-of-the-art theoretical concepts for diffusion and surface processes of contaminants in clayey materials.

## Reactive transport model and investigated systems

**Reactive transport simulator and user interface.** Reactive transport calculations of this study were carried out with the code CrunchClay. CrunchClay code capabilities, equations, and the solver scheme have been described in detail in previous publications$^{15,17-20}$.

A graphical user interface, CrunchEase, has been designed to help users build CrunchClay input files and to extract radionuclide (or another tracer) diffusion results. This further allows the readers to re-run the simulations shown in the present study, and to easily use the code to fit their experimental diffusion data. The CrunchEase package is freely available. Installation instructions and snapshots describing the interface are available in the supplementary information. Data and modeling results are plotted in the interface as normalized tracer flux, $F_{norm}$ (in m s$^{-1}$) as a function of time:

$$F_{norm}(t_n) = \frac{F_{exp}(t_n)}{c_{H,init}} \qquad (18)$$

where $c_{H,init}$ is the initial concentration of tracer in the high-concentration reservoir. In addition, a tabulated text file is produced as an output file that contains modeling results together with corresponding experimental results.

**Model systems.** A series of through-diffusion experimental setups with increasing complexity were used to demonstrate the advantages and reliability of reactive transport simulations in a step-wise manner (Table 1). First, experimental setups without filters were modeled with varied $K_D$ values and sample lengths. Second, the same systems were modeled in the presence of filters. Third, sampling events were additionally included. Last, stagnant tubing volumes were added to the simulation as well. For each of these conditions, we characterized how the calculated diffusive fluxes would change by including a particular feature or not. This approach enabled us to quantify the specific influence of each feature on the estimation of diffusion parameters, and thus to quantify the errors induced in the estimation of diffusion parameters by overlooking these features. Where possible, CrunchEase results were compared to analytical solutions to demonstrate the numerical accuracy of our reactive transport modeling approach.

The parameter values chosen were representative of values reported for actual diffusion experiments in the literature$^{10,21-24}$. The self-diffusion coefficient of the tracer was set at $D_0 = 2.1 \times 10^{-9}$ m$^2$ s$^{-1}$, which corresponds to an effective diffusion coefficient of $D_{e,sample} = 5.25 \times 10^{-11}$ m$^2$ s$^{-1}$ [see Eq. (2) and parameters in Table 1]. The sample grain density was $\rho_{grain} = 2840$ kg m$^{325}$, which corresponds to a bulk dry density of $\rho_d = \rho_{grain} \times (1 - \varepsilon_{sample}) = 1420$ kg m$^{-3}$.

**Application to real experimental data.** As a proof-of-concept, we also used CrunchEase to simulate published through-diffusion data for HTO, Br$^-$ and Ca$^{2+}$ in montmorillonite. Raw data (collected volumes and concentrations as a function of time) are available in the published Supplementary Information from Tinnacher et al.$^{23}$, as well as input parameters summarized in Table S1.

## Results and discussion

**Verification of CrunchEase/CrunchClay calculation accuracy.** For simulations 1–6 (Table 1, no filter), diffusion results calculated with CrunchEase/CrunchClay are in almost perfect agreement with the results obtained with analytical solutions, independent of $K_D$ and $L_s$ values (Fig. 2; the graphical user interface CrunchEase makes it possible to visualize this comparison directly).

**The influence of filters.** The presence of filters has a marked influence on the outcomes for simulations 7–12 (with filters) compared to simulations 1–6 (no filters) (Fig. 2 and S10), as reported in the literature. The error ($\Delta D_e$, $\Delta \alpha$) associated with neglecting filter effects can be illustrated by comparing the actual sample diffusion parameters ($D_{e,ref}$ and $\alpha_{ref}$), which can be calculated from parameters in Table 1, with the values that can be fitted ($D_{e,fit}$ and $\alpha_{fit}$) with Eq. (11) to match the results obtained in the presence of filters (dashed lines on Fig. S10):

$$\begin{aligned} \Delta D_e &= \frac{D_{e,ref} - D_{e,fit}}{D_{e,ref}} \\ \Delta \alpha &= \frac{\alpha_{ref} - \alpha_{fit}}{\alpha_{ref}} \end{aligned} \qquad (19)$$

If the presence of filters is neglected in the fitting procedure, effective diffusion coefficients are underestimated and rock capacity factors are overestimated. Furthermore, the influence of filter effects on the absolute $D_e$ value
<table><caption>**Table 1.** Summary of input parameters of reactive transport simulations 1–24. *2D-axisymmetric simulations.</caption><thead><tr><th>Sim</th><th>Sample diameter (mm)</th><th>Sample length (mm)</th><th>Sample porosity (-)</th><th>Sample tortuosity (-)</th><th>Filters diameter (mm)</th><th>Filters length (mm)</th><th>Filters porosity (-)</th><th>Filters tortuosity (-)</th><th>K<sub>D</sub> (L kg<sup>-1</sup>)</th><th>High reservoir volume (mL)</th><th>Low reservoir volume (mL)</th><th>High and low reservoir conc</th><th>Sampling events</th><th>Volume of tubing (mL)</th><th>Figure(s)</th></tr></thead><tbody><tr><td>1</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>2</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>3</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>4</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>5</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>6</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3</td></tr><tr><td>7</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4</td></tr><tr><td>8</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4, 6</td></tr><tr><td>9</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>1</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4</td></tr><tr><td>10</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>1</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4</td></tr><tr><td>11</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4</td></tr><tr><td>12</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4</td></tr><tr><td>13</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>7*</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4, 5</td></tr><tr><td>14</td><td>10</td><td>10</td><td>0.5</td><td>0.05</td><td>7*</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4, 5</td></tr><tr><td>15</td><td>10</td><td>10</td><td>0.2</td><td>0.025</td><td>7*</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4, 5</td></tr><tr><td>16</td><td>10</td><td>10</td><td>0.5</td><td>0.8</td><td>7*</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Constant</td><td>None</td><td>-</td><td>3, 4, 5</td></tr><tr><td>17</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0.25 day</td><td>-</td><td>6</td></tr><tr><td>18</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0.5 day</td><td>-</td><td>6</td></tr><tr><td>19</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Variable</td><td>Every 1 day</td><td>-</td><td>6, 7</td></tr><tr><td>20</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0.25 day</td><td>-</td><td>6</td></tr><tr><td>21</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0.5 day</td><td>-</td><td>6</td></tr><tr><td>22</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>2</td><td>1000</td><td>20</td><td>Variable</td><td>Every 1 day</td><td>-</td><td>6</td></tr><tr><td>23</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0 to 1 day (random)</td><td>-</td><td>8</td></tr><tr><td>24</td><td>10</td><td>5</td><td>0.5</td><td>0.05</td><td>10</td><td>1</td><td>0.3</td><td>0.3</td><td>0</td><td>1000</td><td>20</td><td>Variable</td><td>Every 0 to 1 day (random)</td><td>1.5</td><td>8</td></tr></tbody></table>

increases with increasing $K_D$ values and $L_{filter}/L_{sample}$ ratios (Fig. S10 and Table S2), as previously reported¹⁰. Since it is so easy to include filters in the proposed new modeling package, one should do so as a good standard practice, even when the filters are not expected to affect the results substantially compared to other sources of bias¹⁰.

**1D versus 2D axisymmetric cylindrical geometry.** 2D- axisymmetric cylindrical simulations 13–16 were carried out with a filter inner diameter (Fig. 1D) of 7 mm, which corresponds to a filter diffusive central sectional area equal to 49% of the total sample cross-sectional area (Table 1; simulations 13–16; Fig. 3). The result of a 1D simulation with a filter porosity scaled down to 49% of its initial value (0.147 instead of 0.3) and a tortuosity value decreased from 0.3 to 0.11 after trial and error iterations proved to be in good agreement with the result of the 2D axisymmetric cylindrical geometry simulation for simulations 13 and 14 (Fig. 3A,B). However, the filter tortuosity value had to be adjusted to other values for simulation 15 and 16 (Fig. 3C,D). If the sample diffusivity is very low compared to that of the filter (Sim. 15, Fig. 3 C), the corrected filter tortuosity should be set at a value far lower than the actual value (0.032 instead of 0.3 in Sim. 15). In contrast, if the sample diffusivity is high compared to that of the filter (Sim. 16, Fig. 3 D), the corrected filter tortuosity should be set to a value similar to that of the actual value (0.27 instead of 0.3 in Sim. 16). In this last case, the lateral diffusion in the sample from the filter diffusive central sectional area is rapid compared to the longitudinal diffusion through the filter and the sample. Thus, the effect of considering a 1D versus 2D axisymmetric cylindrical geometry is almost entirely accounted in the reduction of the porosity that corresponds to a re-scaling of the diffusion sectional area.

The consideration of a 2D axisymmetric cylindrical geometry instead of a 1D geometry (Fig. 1D,E) increases the calculation time significantly. For example, the 1D calculation for simulation 12 (see Table 1) took 23 s, while the equivalent simulation for the 2D axisymmetric cylindrical geometry (Fig. 1D) took 13 min (Intel Core i7, 10th Gen; CrunchClay calculations ran on 4 cores). It may be possible to improve the speed of the 2D axisymmetric cylindrical geometry simulation e.g., by using a coarser grid resolution, or larger time steps. However, hardwired parameters in CrunchEase were chosen to maximize the accuracy of simulations while minimizing users' concerns in terms of their need to intervene in the parameter optimization process. Instead, the fitting
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Verification of the accuracy of reactive transport modeling calculation by comparison with analytical equations results for simple model systems, without (left) or with (right) the presence of filters. Simulations 1–12 (Table 1): Comparison of CrunchEase/CrunchClay results (open circles) with analytical (solid lines) Eq. (11) (no filters) or Eq. (15) (with filters) results as a function of sample length ($L_{sample}$ = 5 mm or 10 mm) and $K_D$ value (from 0 to 2 L kg$^{-1}$).</figcaption></figure>

procedure workflow may be greatly accelerated by simulating a representative 1D geometry, before explicitly taking into account the 2D axisymmetric cylindrical geometry, for further model verification and parameter fine-tuning. CrunchEase makes it easy to implement this approach.

**Influence of sampling events.** Regular sampling time intervals. The periodical replacement of low-concentration reservoir solutions with fresh, tracer-free solutions makes it possible to calculate a diffusive flux based on Eq. (13) or (14). However, this flux calculated at time $t_n$ corresponds to an average (or time integrated) flux from time $t_{n-1}$ to $t_n$ or $t_{n-1}$ to $t_{n+1}$ respectively. If Eq. (13) is used, then time-integrated fluxes are lower than instantaneous fluxes in the early, transient regime of the flux data before diffusion steady-state is achieved (simulations 13–18, Fig. 4A,B). This shifts the breakthrough curves to later times. This bias could be incorrectly interpreted as a retardation factor, and thus to an overestimation of the rock capacity factor. A simple approach for a correction consists in plotting the flux value with $t_{plot} = (t_n + t_{n-1})/2$, which cancels most but not all of the bias (Fig. 4C). As a matter of example, the residual bias in Fig. 4C would correspond to an incorrect estimation of the $K_D$ value at 0.04 L kg$^{-1}$ instead of the actual value at 0 L kg$^{-1}$.

The central difference scheme of Eq. (14) does not produce the same bias as the backward difference scheme of Eq. (13), especially if the sampling time intervals are constant. If they are not constant, then Eq. (14) leads to either a backward or a forward shift in time depending on the position of $t_n$ with respect to $t_{n-1}$ and $t_{n+1}$. The resulting uncertainty of the flux values may be higher for the central difference than the backward difference because uncertainties on measured concentrations are propagated over two intervals of time with Eq. (14) instead of one with Eq. (13).

*Irregular sampling time intervals and influence of tubing volumes.* In practice, sampling time intervals vary over the course of a diffusion experiment. These variations have a very significant influence on the time-integrated
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of 2D axisymmetric cylindrical geometry simulations (plain lines) with 1D geometry simulations (symbols) in which filters porosity was scaled down to 49% of its true value (0.147 instead of 0.3) and its tortuosity was decreased from 0.3 to 0.11 for all 1D simulations. Depending on the contrast of diffusional properties between the sample and the filters, 1D geometry simulations with adapted porosity and tortuosity values for the sample can be in good agreement with the results obtained with 2D axisymmetric cylindrical geometry simulations (A sim.13 and 14; B sim. 15), or not (B sim. 16).</figcaption></figure>

flux values if tubing volumes are present, which is often the case in through-diffusion experiments (Fig. 4D). First, computed time-integrated fluxes with tubing are always higher than without tubing. This can be explained by the “memory” effect as follows⁹. The solution contained in the tubing volume is not replaced with a fresh, tracer-free solution at the sampling event. This means that this solution provides a significant amount of tracer to the next zero-concentration solution as soon as the new reservoir is connected and circulated. Because the time-integrated flux is calculated based on the accumulation of the tracer in this reservoir and tubing, this initial “flash” increases the tracer concentration results, leading to an overestimation of the flux compared to the real instantaneous flux. Second, variations in sampling time intervals can lead to time-integrated flux variations greater than 100% (Fig. 4D). Again, tubing volume plays a major role here. If a long sampling time interval is followed by a short one, then the tracer concentration increases in the tubing volume proportionally with the length of the time interval. Then, the high-concentration volume contributes significantly more to the total concentration measured in the low-concentration reservoir in the next short sampling time interval. Thus, this leads to an apparent increase in the diffusive flux. This sampling artefact can be corrected by subtracting the contribution of the solution in the tubing to the total activity measured in the low-concentration reservoir. However, this correction is not always carried out in published studies. In most cases, no indications about the corrections applied to the activity accumulation and flux calculation are provided.

**Modeling actual experimental diffusion data with CrunchEase.** Diffusion data published in the literature sometimes exhibit large variations in diffusive flux calculated at steady-state²⁶⁻²⁸. These variations are most often interpreted in terms of uncertainties, e.g., related to analytical tracer concentration measurements or experimental temperature fluctuations. The present analysis supports previous findings⁹ that a significant fraction of these apparent variations might instead be due to inherent biases resulting from the experimental setups. These biases can be fully taken into account with our reactive transport modeling approach in CrunchEase/ CrunchClay, which makes it possible to further decrease the uncertainties related to diffusion parameter estimation. This statement is exemplified with a simulation case based on real experimental data in the next section.
<figure><img src="image_6.png" /></figure>

**Figure 4.** (A,B) Comparison of diffusion breakthrough curves, computed with CrunchEase/CrunchClay, for different sampling time intervals in simulations 17–19 (A $K_D=0$ L kg$^{-1}$), and simulations 20 to 22 (B $K_D=2$ L kg$^{-1}$). Full lines: instantaneous fluxes [calculated with very small time steps, or equivalently with the analytical solution shown in Eq. (15)]. Symbols: fluxes integrated using Eq. (13) over 0.25-day (open circles; simulations 17 and 20), 0.5-day (open squares; simulations 18 and 21) and 1-day (open triangles; simulations 19 and 22) sampling intervals. The dotted lines are eye guides only. (C) Instantaneous flux (plain line) and integrated flux with sampling events (triangles; dotted line is an eye-guide) in simulation 19 (1-day sampling interval) with plotting time set at $t_{plot} = (t_n + t_{n-1})/2$. (D) Instantaneous flux (full line) and integrated flux (triangles and circles) for a $K_D$ value of 0 and a sample length $L_{sample}=5$ mm, simulated with CrunchEase/CrunchClay. Time intervals were generated randomly. Triangles show results in the absence of tubing volumes (simulation 23), while circles are results for simulations in the presence of tubing (dead) volumes (simulation 24). The dotted lines are eye guides for time-integrated flux only. The plotting time for integrated flux was set at $t_{plot} = (t_n + t_{n-1})/2$.

Through-diffusion data are most often published (plotted) in the form of flux and/or total diffused concentrations as a function of time. Unfortunately, none of these representations convey the original raw data in terms of tracer concentrations and exact volumes of reservoir solutions and tubings.

However, Tinnacher et al.²³ have provided original raw and tabulated data, which can be used as described in the following. Experimental data and modeling results are plotted as normalized flux as a function of time, as well as low-concentration reservoir concentrations as a function of time (Fig. 5). We believe that the latter representation is preferred for evaluation of the agreement between experimental and modeling results because it does not include any propagated error from the data interpretation and the consideration of two measurements at time $t_{n-1}$ and $t_n$ or $t_{n+1}$. The former representation was included nonetheless to allow for comparisons with previously published data.

The experiment was modeled with a 1D geometry using the same dimensions as those provided by Tinnacher et al.²³. HTO data were modeled by adjusting only the tortuosity of the clay materials ($\tau_{clay,HTO}=0.049$, corresponding to $D_{e,HTO}=7.51 \times 10^{-11}$ m² s$^{-1}$), other parameters (sample porosity $\varepsilon_{clay,HTO}=0.72$, filter porosity $\varepsilon_{filter,HTO}=0.25$, and filter tortuosity $\tau_{filter,HTO}=0.43$) having been determined independently from the diffusion experiments²³.

Modeling results, which almost perfectly match the measured concentrations (Fig. 5), clearly demonstrate the power of the proposed reactive transport modeling approach in terms of reproducing the effects of biases that are unavoidable due to the experimental setup, including the presence of filters, tubing ‘dead’ volumes, and variable duration sampling events. The same was true for experimental Br⁻ diffusion data for which two parameters
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Modeling with CrunchEase/CrunchClay of HTO, Br, and Ca diffusion data from Tinnacher et al.²³. Open squares: experimental data. Closed and open circles: modeling results. HTO model parameters: $\tau_{clay}$=0.049; $\varepsilon_{clay}$=0.72. Br model parameters: $\tau_{clay}$=0.04; $\varepsilon_{clay}$=0.55. Ca model parameters: $\tau_{clay}$=1.1–2.3; $\varepsilon_{clay}$=0.72; $K_D$=79 L kg⁻¹. Closed circles: model with $\tau_{clay}$=1.6; Open circles: model with $\tau_{clay}$=2.3 (upper curve) or $\tau_{clay}$=1.1 (lower curve). The dotted lines are eye guides for modeling results only. Other (fixed) parameters are given in Table S1 in the supplementary materials.</figcaption></figure>

were fitted: $\tau_{clay,Br}$ = 0.040; $\varepsilon_{clay,Br}$ = 0.55, corresponding to $D_{e,Br}$ = 4.4 × 10⁻¹¹ m² s⁻¹ and $\alpha_{Br}$ = 0.55 (Figs. 5), as well as for Ca²⁺ diffusion data for which three parameters were fitted: $\tau_{clay,Ca}$ = 1.1–2.3 (Figs. 5 shows that it is not possible to give a better estimated than this range of values); $\varepsilon_{clay,Ca}$ = 0.72; and $K_{D,Ca}$ = 79 L kg⁻¹, corresponding to $D_{e,Ca}$ = 6.28 × 10⁻¹⁰ m² s⁻¹ to 1.26 × 10⁻⁹ m² s⁻¹ and $\alpha_{Ca}$ = 63.5 (Fig. 5). The $\tau_{clay,Ca}$ value fitted here differs significantly from the value reported in Tinnacher et al.²³. The difference is due to improper handling of filter tortuosity in PHREEQC calculations carried out by Tinnacher et al.²³ for the Ca²⁺ diffusion simulation.

**Beyond the $D_e-\varepsilon-K_D$ modeling approach.** Although CrunchEase is not (yet) intended to handle the creation of CrunchClay files with complex solution chemistry and porosity distributions, the code can be used to create input files with a specific geometry of interest (1D or 2D rotational, presence or absence of filters, tubings, actual reservoirs volumes, etc.), and can also handle sampling events of the low-concentration reservoir. These files (text files with free format) can then be modified directly by the user to include additional chemical reactions and porosity characteristics, such as surface charge and the diffuse layer pore distribution of the clay (a model input file is provided in the Supplementary Materials as a template for interested users).

The strength of this approach is illustrated with a dual-porosity model to simulate Ca²⁺, Br⁻ and HTO diffusion data from Tinnacher et al.²³ (Fig. 6). Model input files were modified from the input files built by CrunchEase and used to simulate and plot Fig. 5. The porosity of the clay sample was split into two contributions: a bulk porosity ($\varepsilon_{bulk}$) and a diffuse layer porosity ($\varepsilon_{DL}$) in which the surface charge (0.9 mol kg⁻¹ of montmorillonite was compensated (see¹⁸,²³,²⁹⁻³¹ for full details about theory and parameters). An almost perfect data fit was obtained with the following parameters (Fig. 6): $\varepsilon_{bulk}$ = 0.19; $\varepsilon_{DL}$ = 0.53; $\tau_{bulk}$ = 0.1; $\tau_{DL}$ = 0.032; log $K_{Na}$ = 0 for the surface complexation reaction:

$$> \text{Surf}^{-} + \text{Na}^{+} \rightleftharpoons > \text{SurfNa} \qquad (20)$$

and log $K_{Ca}$ = 0.7 for the surface complexation reaction:

$$2 > \text{Surf}^{-} + \text{Ca}^{2+} \rightleftharpoons > \text{Surf}_2\text{Ca} \qquad (21)$$

While the model prediction using a dual-porosity model (Fig. 6) is not better than the model prediction made with a simple $D_e-\varepsilon-K_D$ modeling approach (Figs. 5), it has the advantage to produce model parameters that are more realistic with regards to sample porosity properties. In this respect, CrunchEase can help simplify the use of CrunchClay to achieve process level understanding.

## Conclusions

We have demonstrated that a reactive transport modeling approach that includes the full geometrical and operational complexity of an experiment, such as heterogeneities and time-dependent sampling effects is a very powerful tool for the interpretation of diffusion data. While a classical $D_e-\varepsilon-K_D$ model may be implemented
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Dual porosity (bulk and diffuse layer) model (black circles; dotted lines are eye guidelines) applied to diffusion data (open squares) from Tinnacher et al.²³. Model input files were modified from files built by CrunchEase and used to simulate and plot Fig. 5. Ca²⁺, Br⁻, and HTO diffusion were modeled simultaneously. Ca²⁺ and Br⁻ input files were the same because of identical sampling times (same experiment), but the HTO input file was handled separately (different experiment, sampling times, and sampled volumes). Modeling parameters were otherwise the same (see text).</figcaption></figure>

in any reactive transport code, the development of the new graphical user interface CrunchEase makes it easy to apply by experimentalists without background in reactive transport modeling. CrunchEase can also be used to fit through-diffusion data from experiments with changing boundary conditions (no reservoir replacement; not shown in this paper). In the future, the code may be easily extended to other types of diffusion experiments, such as in-diffusion experiments.

Since experimental biases can be modeled explicitly, the accuracy of calculated diffusion parameters with a reactive transport approach is expected to be better than for the application of analytical solutions, which cannot include the full complexities of real experiments. The numerical efficiency of CrunchClay makes our approach suitable to fit model parameters very quickly. For instance, for 1D problems simulating a diffusion experiment with an experimental time frame of approximately 50 days, the model running time is most often less than 1 min. This approach could also be used easily with parameter estimation softwares such as PEST³².

On a final note, CrunchEase makes it easier to transition from a $D_c-\epsilon-K_D$ modeling approach to a reactive transport modeling approach taking into account surface-enhanced diffusion processes and advanced chemical reactivity models such as surface complexation models. Very few reactive transport codes can handle all of these features³¹, which currently makes the CrunchClay-CrunchEase combination unique. This will hopefully help to
facilitate the dialog between experimentalists, modelers, and radioactive waste agency decision-makers regarding the most recent concepts applied to diffusion in clayey materials. This may ultimately lead to an improved decision-making process by radioactive waste agencies.

## Data availability

The CrunchEase interface and data are available on an open Github repository: https://github.com/Tournassat/ CrunchEaseForAll/releases/tag/v1.0.0.

Received: 23 March 2023; Accepted: 6 September 2023

Published online: 12 September 2023

## References

1. Altmann, S. et al. Diffusion-driven transport in clayrock formations. *Appl. Geochem.* **27**, 463–478 (2012).

2. Bourg, I. C. & Tournassat, C. Chapter 6—Self-diffusion of water and ions in clay barriers. In *Natural and Engineered Clay Barriers* Vol. 6 (eds Tournassat, C. et al.) 71–100 (Elsevier, 2015).

3. Liu, X. et al. Molecular-level understanding of metal ion retention in clay-rich materials. *Nat. Rev. Earth Environ.* **3**, 461–476 (2022).

4. Crank, J. *The Mathematics of Diffusion* (Oxford University Press, 1975).

5. De Cannière, P. et al. Diffusion and sorption of ³²Si-labelled silica in the boom clay. *Radiochim. Acta* **82**, 191–196 (1998).

6. Melkior, T. *Etude méthodologique de la diffusion de cations interagissant dans les argiles. Application: mise en œuvre expérimentale et modélisation du couplage chimie-diffusion d'alcalins dans une bentonite synthétique* (CEA, 2000).

7. Yaroshchuk, A. E., Glaus, M. A. & Van Loon, L. R. Diffusion through confined media at variable concentrations in reservoirs. *J. Membr. Sci.* **319**, 133–140 (2008).

8. Tachi, Y. & Yotsuji, K. Diffusion and sorption of $Cs^+$, $Na^+$, $I^-$ and HTO in compacted sodium montmorillonite as a function of porewater salinity: Integrated sorption and diffusion model. *Geochim. Cosmochim. Acta* **132**, 75–93 (2014).

9. Glaus, M. A., Aertsens, M., Maes, N., Van Laer, L. & Van Loon, L. R. Treatment of boundary conditions in through-diffusion: A case study of $^{85}$Sr$^{2+}$ diffusion in compacted illite. *J. Contam. Hydrol.* **177**, 239–248 (2015).

10. Glaus, M. A., Rosse, R., Van Loon, L. R. & Yaroshchuk, A. E. Tracer diffusion in sintered stainless steel silters: Measurement of effective diffusion coefficients and implications for diffusion studies with compacted lays. *Clays Clay Miner.* **56**, 677–685 (2008).

11. Takeda, M., Nakajima, H., Zhang, M. & Hiratsuka, T. Laboratory longitudinal diffusion tests: 1. Dimensionless formulations and validity of simplified solutions. *J. Contam. Hydrol.* **97**, 117–134 (2008).

12. Moridis, G. J. Semianalytical solutions for parameter estimation in diffusion cell experiments. *Water Resour. Res.* **35**, 1729–1740 (1999).

13. Savoye, S., Beaucaire, C., Fayette, A., Herbette, M. & Coelho, D. Mobility of cesium through the callovo-oxfordian claystones under partially saturated conditions. *Environ. Sci. Technol.* **46**, 2633–2641 (2012).

14. Chen, C.-L. & Wang, T.-H. An innovative method for determining the diffusion coefficient of product nuclide. *Nucl. Eng. Technol.* **49**, 1019–1030 (2017).

15. Steefel, C. I. et al. Reactive transport codes for subsurface environmental simulation. *Comput. Geosci.* **19**, 445–478 (2015).

16. Alt-Epping, P. et al. Benchmark reactive transport simulations of a column experiment in compacted bentonite with multispecies diffusion and explicit treatment of electrostatic effects. *Comput. Geosci.* **19**, 535–550 (2015).

17. Tournassat, C. & Steefel, C. I. Modeling diffusion processes in the presence of a diffuse layer at charged mineral surfaces. A benchmark exercise. *Comput. Geosci.* https://doi.org/10.1007/s10596-019-09845-4 (2020).

18. Tournassat, C. & Steefel, C. I. Ionic transport in nano-porous clays with consideration of electrostatic effects. *Rev. Mineral. Geochem.* **80**, 287–330 (2015).

19. Soler, J. M., Steefel, C. I., Gimmi, T., Leupin, O. X. & Cloet, V. Modeling the ionic strength effect on diffusion in clay. The DR-A experiment at Mont Terri. *ACS Earth Space Chem.* **3**, 442–451 (2019).

20. Steefel, C. I. & Tournassat, C. A model for discrete fracture–clay rock interaction incorporating electrostatic effects on transport. *Comput. Geosci.* **25**, 395–410 (2021).

21. Glaus, M. A., Frick, S., Rosse, R. & Van Loon, L. R. Comparative study of tracer diffusion of HTO, Na-22(+) and Cl-36(-) in compacted kaolinite, illite and montmorillonite. *Geochim. Cosmochim. Acta* **74**, 1999–2010 (2010).

22. Motellier, S., Devol-Brown, I., Savoye, S., Thoby, D. & Alberto, J.-C. Evaluation of tritiated water diffusion through the Toarcian clayey formation of the Tournemire experimental site (France). *J. Contam. Hydrol.* **94**, 99–108 (2007).

23. Tinnacher, R. M., Holmboe, M., Tournassat, C., Bourg, I. C. & Davis, J. A. Ion adsorption and diffusion in smectite: Molecular, pore, and continuum scale views. *Geochim. Cosmochim. Acta* **177**, 130–149 (2016).

24. Van Loon, L. R., Soler, J. M. & Bradbury, M. H. Diffusion of HTO, ³⁶Cl⁻ and ¹²⁵I⁻ in opalinus clay samples from Mont Terri: Effect of confining pressure. *J. Contam. Hydrol.* **61**, 73–83 (2003).

25. Tournassat, C. & Appelo, C. A. J. Modelling approaches for anion-exclusion in compacted Na-bentonite. *Geochim. Cosmochim. Acta* **75**, 3698–3710 (2011).

26. Glaus, M. A., Frick, S., Rossé, R. & Van Loon, L. R. Consistent interpretation of the results of through-, out-diffusion and tracer profile analysis for trace anion diffusion in compacted montmorillonite. *J. Contam. Hydrol.* **123**, 1–10 (2011).

27. Melkior, T. et al. Na⁺ and HTO diffusion in compacted bentonite: Effect of surface chemistry and related texture. *J. Hydrol.* **370**, 9–20 (2009).

28. Savoye, S. et al. Effect of temperature on the containment properties of argillaceous rocks: The case study of Callovo-Oxfordian claystones. *J. Contam. Hydrol.* **125**, 102–112 (2011).

29. Appelo, C. A. J., Van Loon, L. R. & Wersin, P. 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 (2010).

30. Appelo, C. A. J. & Wersin, P. Multicomponent diffusion modeling in clay systems with application to the diffusion of tritium, iodide, and sodium in Opalinus clay. *Environ. Sci. Technol.* **41**, 5002–5007 (2007).

31. Tournassat, C. & Steefel, C. I. Reactive transport modeling of coupled processes in nanoporous media. *Rev. Mineral. Geochem.* **85**, 75–110 (2019).

32. Doherty, J. *Calibration and Uncertainty Analysis for Complex Environmental Models* (Watermark Numerical Computing Brisbane, 2015).

## Acknowledgements

The authors would like to acknowledge the funding support by the Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, through its Geoscience program at LBNL under Contract DE-AC02-05CH11231, the U.S. Department of Energy under Contract DE-AC02-05CH11231 under
the auspices of the Spent Fuel and Waster Disposition program (Office of Nuclear Energy) and the Geosciences program (Office of Science, Office of Basic Energy Sciences), the DOE Office of Nuclear Energy's Nuclear Energy University Programs (NEUP, Federal grants DE-NE0008683 and DE-NE0008938), the Agence Nationale de la Recherche (ANR) under Contract No. ANR-18-CE05-0035-01, a grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" Programme LabEx VOLTAIRE, 10-LABX-0100, and EC Horizon 2020 project EURAD under Grant Agreement 847593 (WP DONUT).

**Author contributions**

C.T. and C.I.S. conceived the modeling approaches. C.T. developed the CrunchEase interface, modeled the data and wrote the initial version of the manuscript. C.T. and C.I.S. developed CrunchClay. P.M.F. and R.M.T. provided experimental data and tested the code. All authors reviewed the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-42260-5.

**Correspondence** and requests for materials should be addressed to C.T.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**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/.

© The Author(s) 2023