# A longer introduction

TransPyREnd is short for “*Trans*portmodell in *Py*thon für den Austrag von *R*adionukliden aus einem *End*lager”. 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:

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:

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):

\(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:

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

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:

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

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

The Crank-Nicolson scheme can be written as

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

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

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

with the following definitions

and

We can simplify the equation to

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.

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.