Radionuclides migration in Opalinus clay

The diffusion of Cs-135 and U-238 through a 1-D column of fully water-saturated Opalinus clay for 1 million years is simulated using a mechanistic approach via the coupling of OGS to PHREEQC (OGS#iPHREEQC). Radioactive decay has not been considered due to the long half-lives (Nagra, 2002) of Cs-135 (2.3e6 years) and U-238 (4.5e9 years). The diffusion of radionuclides into the clay porewater is modeled with a Dirichlet boundary condition applied at the left side of the domain, whereas the right-side is modeled as a closed boundary (Neumann). The 1-D domain length is 40 m, which is sufficient to ensure that Cs does not reach the right-side boundary.

The following Opalinus clay properties are used:

The domain is discretized in 200 linear elements (avg. element size = 0.2 m) with a non-uniform mesh refinement at the inlet boundary (to account for the high concentration gradient produced due to the constant concentration). The minimum and maximum element sizes are 0.04 m and 0.35 m, respectively (element length multiplier = 1.00625). Mesh-independence has been verified by running the simulations with a mesh with 400 elements (0.1 m element size). The results did not show a significant difference. The temporal discretization follows the maximum time-step size given by the von Neumann stability critetion for diffusive transport:

therefore, a uniform time-stepping scheme is adopted with $\Delta t = \text{2e8}~s~(\approx 6.34~\text{years})$. This means that the simulation needs a total of 157680 time-steps to estimate the migration of 1 million years. This causes that large CPU times are needed when using OGS#iPHREEQC.

Finally, the total concentrations in Opalinus clay of the main species, including background concentrations of Cs and U, are as follows:

Cs-135 migration

A Cs concentration at the left boundary of 1.019e-7 mol/kgw is set constant throughout the simulation. The sorption model is taken from Águila et al., 2021. Cs-135 migration is simulated by using OGS#iPHREEQC via the <surface> keyword with the following parameters:

The corresponding reactions of the sorption model are added to the latest version (2020) of the PSI/Nagra thermodynamic database (PSINagra2020v1-1_davies.dat).

In the Cs.prj file, the surface input under the <chemical_system> keyword is:

<chemical_system chemical_solver="Phreeqc">
    ...
        <surface site_unit="mole">
            <site>
                <name>Su_</name>
            </site>
            <site>
                <name>Su_fes</name>
            </site>
            <site>
                <name>Su_ii</name>
            </site>
        </surface>
    ...

The amounts of each site are given through the <media> keyword as follows:

                <phase>
                    <type>Solid</type>
                    <components>
                        <component>
                            <name>Su_</name>
                            <properties>
                                <property>
                                    <name>molality</name>
                                    <type>Constant</type>
                                    <value>1.4250e+00</value>
                                </property>
                            </properties>
                        </component>
                        <component>
                            <name>Su_fes</name>
                            <properties>
                                <property>
                                    <name>molality</name>
                                    <type>Constant</type>
                                    <value>1.8000e-03</value>
                                </property>
                            </properties>
                        </component>
                        <component>
                            <name>Su_ii</name>
                            <properties>
                                <property>
                                    <name>molality</name>
                                    <type>Constant</type>
                                    <value>1.3800e-01</value>
                                </property>
                            </properties>
                        </component>
                    </components>
                </phase>

Where the site names should be an exact match in both <media> and <chemical_system> keywords. This means that it is possible to setup complex simulations comprised of multiple physical domains with distinct chemical systems (e.g., different sorption sites, diffusion, porosity) such as the multiple barriers in a nuclear waste repository.

1) Solve the numerical model

The simulation takes approximately 23 hrs. to complete when t = 1e6 years. Therefore, only the first 50 time-steps are simulated in this notebook. To run the full simulation, the time loop parameters can be easily adapted.

import os
import time

import matplotlib.pyplot as plt
import numpy as np
import vtuIO
from matplotlib.pyplot import cm

from pathlib import Path

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)

prj_name = "Cs.prj"

t0 = time.time()
print(">>> OGS started execution ... <<<")

! ogs {prj_name} -o {out_dir} > {out_dir}/outCs.txt

tf = time.time()
print(">>> OGS terminated execution  <<< Elapsed time: ", round(tf - t0, 2), " s.")

>>> OGS started execution ... <<<

>>> OGS terminated execution  <<< Elapsed time:  36.55  s.

2) Plot the results

Results can be conveniently visualized with VTUinterface. Note that the results of the full simulation are plotted (provided in the folder ./Cs_full_simulation). The results of the first 50 time steps can be plotted by replacing pvdfile with the commented line below.

length = 40
xaxis = [(i, 0, 0) for i in np.linspace(start=0.0, stop=length, num=201)]
x = np.array(xaxis)[:, 0]
fieldname = "Cs"

# Read simulation results
pvdfile = vtuIO.PVDIO("./Cs_full_simulation/out.pvd", dim=1)
# pvdfile = vtuIO.PVDIO(f"{out_dir}/Cs.pvd", dim=1)

color_map = iter(cm.rainbow(np.linspace(0, 1, len(pvdfile.timesteps))))

# Plot the results
fig, ax = plt.subplots()

for t, c in zip(pvdfile.timesteps, color_map):
    y = pvdfile.read_set_data(t, fieldname, pointsetarray=xaxis)
    ax.plot(
        x,
        y,
        color=c,
        marker=".",
        label=f"{round(t / 3600 / 24 / 365, 0):.0e} yr",
    )

ax.set_xlim(0, 5)
ax.set_title("Cs migration in Opalinus clay")
ax.set_ylabel("Cs concentration [mol/kgw]")
ax.set_xlabel("x [m]")
ax.set_yscale("log")
ax.legend()
plt.tight_layout()

The results show that Cs-135 does not migrate for more than 3 m even after 1 million years.

U-238 migration

Similarly, a U concentration at the left boundary of 1.060e-6 mol/kgw is applied to the domain. In this case, the sorption model is taken from Hennig et al., 2020. U-238 migration is simulated by using OGS#iPHREEQC via the <surface> keyword definition with the following parameters:

The Opalinus clay parameters remain the same as in the Cs case, as well as the space-time discretization. In the same way, the reactions of the sorption model are added to the themodynamic database PSINagra2020v1-1_davies.dat.

1) Solve the numerical model

The model is solved for the first 50 time steps only to minimize the CPU time of this notebook. The time loop parameters in U.prj can be adapted accordingly to cover the full simulation time of 1 million years.

prj_name = "U.prj"

t0 = time.time()
print(">>> OGS started execution ... <<<")

! ogs {prj_name} -o {out_dir} > {out_dir}/outU.txt

tf = time.time()
print(">>> OGS terminated execution  <<< Elapsed time: ", round(tf - t0, 2), " s.")

>>> OGS started execution ... <<<

>>> OGS terminated execution  <<< Elapsed time:  79.16  s.

2) Plot the results

Again, the full simulation results are provided in /U_full_simulation. They are plotted to show the migration of U after 1 million years.

length = 40
xaxis = [(i, 0, 0) for i in np.linspace(start=0.0, stop=length, num=201)]
x = np.array(xaxis)[:, 0]
fieldname = "U"

# Read simulation results
pvdfile = vtuIO.PVDIO("./U_full_simulation/out.pvd", dim=1)
# pvdfile = vtuIO.PVDIO(f"{out_dir}/U.pvd", dim=1)

color_map = iter(cm.rainbow(np.linspace(0, 1, len(pvdfile.timesteps))))

# Plot the results
fig, ax = plt.subplots()

for t, c in zip(pvdfile.timesteps, color_map):
    y = pvdfile.read_set_data(t, fieldname, pointsetarray=xaxis)
    ax.plot(
        x,
        y,
        color=c,
        marker=".",
        markevery=5,
        label=f"{round(t / 3600 / 24 / 365, 0):.0e} yr",
    )

ax.set_xlim(0, 40)
ax.set_title("U migration in Opalinus clay")
ax.set_ylabel("U concentration [mol/kgw]")
ax.set_xlabel("x [m]")
ax.set_yscale("log")
ax.legend()
plt.tight_layout()

In this case, it is clear that the migration of U-238 is much more pronounced than Cs-135. After 1 million years, U-238 shows a breakthrough to around 20 m.

References

Águila, J. F., Montoya, V., Samper, J., Montenegro, L., Kosakowski, G., Krejci, P., & Pfingsten, W. (2021). Modeling cesium migration through Opalinus clay: a benchmark for single-and multi-species sorption-diffusion models. Computational Geosciences, 25(4), 1405-1436.

Garibay-Rodriguez, J., Chen, C., Shao, H., Bilke, L., Kolditz, O., Montoya, V., & Lu, R. (2022). Computational framework for radionuclide migration assessment in clay rocks. Frontiers in Nuclear Engineering (under review).

Hennig, T., Stockmann, M., & Kühn, M. (2020). Simulation of diffusive uranium transport and sorption processes in the Opalinus Clay. Applied Geochemistry, 123, 104777.

Hummel W. Thoenen T. (2021): The PSI Chemical Thermodynamic Database 2020, Nagra NTB 21-03.

Nagra, 2002. Project Opalinus Clay: Models, Codes and Data for Safety Assessment. Technical Report NTB 02–06. Nagra, Switzerland.

Credits:

Jaime Garibay-Rodriguez, Chaofan Chen, Haibing Shao, Lars Bilke, Olaf Kolditz, Vanessa Montoya, Renchao Lu