A longer introduction

TransPyREnd is short for “Transportmodell in Python für den Austrag von Radionukliden aus einem Endlager”. This code is being used to estimate radionuclide transport in the subsurface as part of the safety assessments that are currently on the way in the site selection process for a nuclear waste repository in Germany. The code is accessible at https://gitlab.opengeosys.org/owf/transpyrend.

The motivation for developing TransPyREnd was the desire to use an open-source code, which is beneficial for a transparent workflow in the safety assessments of proposed nuclear waste repository sites. In addition, in the current stage of the site selection in Germany there is a need for a relatively simple and fast code. The code that is presented here simulates transport in one dimension only and does not contain any overhead such as a graphical user interface. This means that its increased numerical performance allows for a large number of model experiments. This is helpful for the evaluation of a large number of potential repositories and to quantify the effects of uncertainty of model parameters and geological developments over the lifetime of the repository. The code was written in Python, which enables a modular object oriented structure and allows the code to benefit from the large number of Python libraries for input and output, visualization, multi-processing and statistics.

in the following, the basic equations that represent the main transport mechanisms of radionuclides in the subsurface, diffusion, advection and sorption, as well as radioactive decay (section 2) are described. Detail of the implementation can be found in Behrens et al. 2023.

Mathematical representation of transport processes

The transport of radionuclides in the subsurface can be described by the following equation:

(1)\[\begin{split}%\underbrace{ R_i}_ %{\text{Retardation}} \underbrace{ \phi R_i \frac{\partial c_i}{\partial t} }_ {\substack{\text{Concentration}\\\text{change}}} {=} \underbrace{ -\frac{\partial j_i}{\partial x}}_ {\text{Diffusion}} {-} \underbrace{ \frac{\partial (q c_i)}{\partial x}}_ {\substack{\text{Advection}}} {+} \underbrace{ W}_ {\text{Decay}}\end{split}\]

where \(\phi\) is porosity (dimensionless), \(R_i\) is the retardation factor of nuclide \(i\) (dimensionless), \(c_i\) is the concentration of nuclide \(i\) in the pore fluid (\(mol \; m^{-3}\)), \(t\) is time (\(s\)), \(j_i\) is the diffusive flux of radionuclide \(i\) (\(mol\; m^{-2} \; s^{-1}\)) and \(q\) is the Darcy velocity (\(m \; s^{-1}\)).

The term \(W\) represents the source/sink term, which is defined as:

\[W = \sum_{j} R_j c_j \lambda_{j,i} - R_i c_i \Lambda_i + S\]

The source/sink term describes the change rate of nuclide \(i\) from decays of other species \(j\) (first term), and from loss of nuclide \(i\) due to decay (second term). \(\lambda_{j,i}\) is the matrix of decay constants for the decay of species \(j\) to species \(i\). \(\Lambda_i\) is the total decay rate of species \(i\) (note the flipped indices):

\[\Lambda_i = \sum_j \lambda_{i,j}\]

\(S\) represents an additional source or sink (\(kg\;m^{-3}\;s^{-1}\)) that can for instance represent the release of a nuclide over time.

The diffusive flux of radionuclide \(i\) is defined as:

\[j_i = -D_{e,i} \dfrac{\partial c_i}{\partial x}\]

where \(D_{e,i}\) is the effective diffusion coefficient of nuclide \(i\) (\(m^2 \; s^{-1}\)). Combining this with equation (1), we arrive at:

(2)\[\begin{split}%\underbrace{ R_i}_ %{\text{Retardation}} \underbrace{ \phi R_i \frac{\partial c_i}{\partial t} }_ {\substack{\text{Concentration}\\\text{change}}} {=} \underbrace{ \frac{\partial}{\partial x} \left(D_{e,i} \frac{\partial c_i}{\partial x} \right)}_ {\text{Diffusion}} {-} \underbrace{ \frac{\partial (q c_i)}{\partial x}}_ {\substack{\text{Advection}}} {+} \underbrace{ W }_ {\text{Decay}}\end{split}\]

This equation represents the physical processes diffusion, advection and the decay of radionuclides. The left-hand side of this equation and the diffusion and advection terms follow standard solute transport formulation (Diersch and Kolditz 2002; Ingebritsen, Sanford, and Neuzil 2006). The full equation, including the decay terms, is similar to the mathematical formulations in other model codes that are frequently used to simulate radioactive transport, such as TOUGH2 (Oldenburg and Pruess 1995) and PFLOTRAN (Hammond et al. 2012).

The sorption of nuclides to the solid rock matrix is described by the retardation factor \(R_i\), which is defined as:

\[R_i = 1 + \frac{K_{d,i}\rho}{\phi}\]

where \(K_{d,i}\) the sorption coefficient (\(m^3 \; kg^{-1}\)) and \(\rho\) is the bulk density (\(kg \; m^{-3}\)). If we divide equation (2) by \(R_i\), we can easily see the effective transport is indeed slowed down by a factor of 1/\(R_i\), hence the name.

Numerical implementation

TransPyREnd solves the transport equations with a finite-difference approach. Time evolution is done using the second-order accurate Crank-Nicolson scheme. Additionally, higher-order schemes from the scipy library can be used, e.g. DOP853, an 8th-order accurate Runge-Kutta scheme. In principle, the semi-implicit Crank-Nicolson method has a higher degree of numerical stability at the cost of numerical accuracy due to numerical dispersion. The explicit method is a higher-order method, which ensures a relatively high degree of accuracy. However this method has a lower degree of numerical stability, i.e. it is only stable at relatively small timesteps, which makes this method computationally more demanding.

The implicit method was based on the existing model code PyBasin (Luijendijk et al. 2011; Luijendijk 2012, 2020) with modifications to the advection component and the addition of sorption and decay.

Crank-Nicolson method

We first derive a discretized form of the transport equation (2). Our domain consists of point \(x_k\), with the index index \(_k\) and distance \(h_k = x_k - x_{k-1}\) between the nodes. The diffusive flux can be discretized as

\[\frac{\partial}{\partial x} D_{e,i} \frac{\partial c}{\partial x} \approx \frac{2}{h_{k+1}+h_{k}}\left( \frac{1}{h_{k}}D_{e,i,k-\frac{1}{2}} c_{i,k-1} - \frac{h_{k+1}+h_{k}}{h_{k}h_{k+1}} (D_{e,i,k-\frac{1}{2}}+D_{e,i,k+\frac{1}{2}}) c_{i,k} + \frac{1}{h_{k+1}} D_{e,i,k+\frac{1}{2}} c_{i,k+1} \right)\]

where subscript i denotes radionuclide i, subscript \(_k\) denotes the grid node for which the equation is evaluated, \(k-1\) denotes the node one position to the left and \(k+1\) denotes the node one position to the right, and an index \(k\pm \frac{1}{2}\) denotes the interface between the respective nodes.

Similarly, the advective flux reads

\[q \frac{\partial c}{\partial x} \approx q \frac{c_{i,k+1} - c_{i,k-1}}{x_{k+1} - x_{k-1}}\]

The Crank-Nicolson scheme can be written as

\[\frac{y^{m+1}-y^m}{\Delta t} = \frac{1}{2} F^{m+1} + \frac{1}{2} F^{m}\]

with \(F\) an appropriate discretization of a given differential equation

\[\frac{\partial y}{\partial t} = f(x, y, t)\]

and \(m\) the index of the considered time step.

The PDE in equation (2) can therefore be approximated as

\[\begin{split}R_{i, k} \phi_{i, k} \frac{c_{i,k}^{m+1}-c_{i,k}^m}{\Delta t} &= \theta H_k \left( H_k^- D_{e,i,k-\frac{1}{2}} c_{i,k-1}^{m+1} - G_k (D_{e,i,k-\frac{1}{2}}+D_{e,i,k+\frac{1}{2}}) c_{i,k}^{m+1} + H_k^+ D_{e,i,k+\frac{1}{2}} c_{i,k+1}^{m+1} \right) \\\\ &- \theta q \frac{c_{i,k+1}^{m+1} - c_{i,k-1}^{m+1}}{x_{k+1} - x_{k-1}} \\\\ &+ (1-\theta) H_k \left(H_k^- D_{e,i,k-\frac{1}{2}} c_{i,k-1}^{m} - G_k (D_{e,i,k-\frac{1}{2}}+D_{e,i,k+\frac{1}{2}}) c_{i,k}^{m} + H_k^+ D_{e,i,k+\frac{1}{2}} c_{i,k+1}^{m} \right) \\\\ &- (1-\theta) q \frac{c_{i,k+1}^{m} - c_{i,k-1}^{m}}{x_{k+1} - x_{k-1}} \\\\ &+ \sum_j \left(R_{j, k} \phi_{j, k} c^m_{j,k} \lambda_{i,j}\right) - R_{i, k} \phi_{i, k} c^m_{i,k} \Lambda_i + W^m_{i,k}\end{split}\]

with the following definitions

\[\begin{split}H_k &= \frac{2}{h_{k+1}+h_{k}}\\\\ H_k^- &= \frac{1}{h_{k}}\\\\ H_k^+ &= \frac{1}{h_{k+1}}\\\\ G_k &= \frac{h_{k+1}+h_{k}}{h_{k}h_{k+1}}\end{split}\]

and

\[\begin{split}s &= \dfrac{\Delta t}{\phi_{i, k} R_{i, k}}\\\\ u_k &= \dfrac{qs}{x_{k+1} - x_{k-1}} \\\\ v_{k}^{-} &= H_k H_k^- D_{e,i,k-1/2} s\\\\ v_k &= H_k G_k (D_{e,i,{k-1/2}} + D_{e,i,{k+1/2}}) s\\\\ v_{k}^{+} &= H_k H_k^+ D_{e,i,k+1/2}s\end{split}\]

We can simplify the equation to

\[\begin{split}&\theta\left( - u_k - v_{k}^{-} \right) c_{i,k-1}^{m+1} + \left( 1 + \theta v_k \right) c_{i,k}^{m+1} + \theta\left( u_k - v_{k}^{+} \right) c_{i,k+1}^{m+1} =\nonumber\\\\ &\theta'\left( u_k + v_{k}^{-} \right) c_{i,k-1}^{m} + \left( 1 - \theta' v_{k} \right) c_{i,k}^{m} + \theta'\left( -u_k + v_{k}^{+} \right) c_{i,k+1}^{m} + s \left( \sum_j \left(R_j \phi_i c^m_{j,k} \lambda_{i,j}\right) - R_i \phi_i c^m_{i,k} \Lambda_i + W^m_{i,k}\right)\end{split}\]

This linear system of equations can be solved by standard methods; TransPyREnd uses a sparse matrix solver from scipy.

References

Baeyens, B, T Thoenen, MH Bradbury, and M Marques Fernandes. 2014. “Sorption Data Bases for Argillaceous Rocks and Bentonite for the Provisional Safety Analyses for SGT-E2.” Paul Scherrer Institute (PSI).

Bateman, Henry. 1910. “The Solution of a System of Differential Equations Occurring in the Theory of Radioactive Transformations.” In Proc. Cambridge Philos. Soc., 15:423–27.

Behrens, C. and Luijendijk, E. and Kreye, P. and Panitz, F. and Bjorge, M. and Gelleszun, M. and Renz, A. and Miro, S. and Rühaak, W.. 2023. “TransPyREnd: a code for modelling the transport of radionuclides on geological timescales” Advances in Geosciences 58 (109-119)

Diersch, H. Jg, and Olaf Kolditz. 2002. “Variable-Density Flow and Transport in Porous Media: Approaches and Challenges.” Advances in Water Resources 25 (8-12): 899–944.

Hammond, G. E., P. C. Lichtner, Chuan Lu, and Richard T. Mills. 2012. “PFLOTRAN: Reactive Flow & Transport Code for Use on Laptops to Leadership-Class Supercomputers.” Groundwater Reactive Transport Models 5: 141–59.

Huysmans, Marijke, and Alain Dassargues. 2007. “Equivalent Diffusion Coefficient and Equivalent Diffusion Accessible Porosity of a Stratified Porous Medium.” Transport in Porous Media 66 (3): 421–38. https://doi.org/10.1007/s11242-006-0028-6.

Ingebritsen, Steven E, Ward E Sanford, and Christopher E Neuzil. 2006. Groundwater in Geologic Processes. Cambridge University Press.

Larue, Jürgen, Bruno Baltes, Heidermarie Fischer, Gerd Frieling, Ingo Kock, Martı́n Navarro, and Holger Seher. 2013. “Radiologische Konsequenzenanalyse.” Gesellschaft für Anlagen-und Reaktorsicherheit.

Li, Xiaoye S. 2005. “An Overview of SuperLU: Algorithms, Implementation, and User Interface.” ACM Trans. Math. Softw. 31 (3): 302–25. https://doi.org/10.1145/1089014.1089017.

Luijendijk, Elco. 2012. “The role of fluid flow in the thermal history of sedimentary basins: Inferences from thermochronology and numerical modeling in the Roer Valley Graben, southern Netherlands [Ph.D. thesis].” PhD, Vrije Universiteit Amsterdam.

Luijendijk, Elco. 2020. PyBasin: Numerical model of basin history, heat flow and thermochronology (version v0.91-alpha). Zenodo. https://doi.org/10.5281/zenodo.4263427.

Luijendijk, Elco, R. T. Van Balen, Marlies Ter Voorde, and P. A. M. Andriessen. 2011. “Reconstructing the Late Cretaceous inversion of the Roer Valley Graben (southern Netherlands) using a new model that integrates burial and provenance history with fission track thermochronology.” Journal of Geophysical Research 116 (B06402): 1–19. https://doi.org/10.1029/2010JB008071.

Navarro, Martin, Torben Weyand, Jens Eckel, and Heidemarie Fischer. 2019. “Indikatoren Zur Bewertung Des Einschlusses Und Der Isolation Mit Exemplarischer Anwendung Auf Ein Generisches Endlagersystem Mit Dem Wirtsgestein Tongestein.” Gesellschaft für Anlagen und Reaktorsicherheit (GRS) gGmbH.

Norsett, E Hairer SP, and G Wanner. 1987. “Solving Ordinary Differential Equations i: Nonstiff Problems.” Springer-Verlag.

Oldenburg, Curtis M, and Karsten Pruess. 1995. “Eos7r: Radionuclide Transport for Tough2.” Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States).

Van Loon, L. R. 2014. “Effective Diffusion Coefficients and Porosity Values for Argillaceous Rocks and Bentonite: Measured and Estimated Values for the Provisional Safety Analyses for SGT-E2.” *Report 1015-2636. http://inis.iaea.org/search/search.aspx?orig_q=RN:48088310.

Virtanen, Pauli, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, et al. 2020. “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.” Nature Methods 17: 261–72. https://doi.org/10.1038/s41592-019-0686-2.