TransPyREnd: a code for modelling the transport of radionuclides on geological timescales

. The German site selection procedure for a high-level nuclear waste repository is entering a stage in which preliminary safety assessments have to be conducted and the release of radionuclides has to be estimated for a large number of potential sites. Here, we present TransPyREnd, a 1D ﬁnite-differences code for modeling the transport of radionuclides in the sub-surface at geological timescales. The code simulates the processes advection, diffusion, equilibrium sorption, decay of radionuclides, and the build-up of daughter nuclides. We summarize the modeled physical processes, their mathematical description and our numerical approach to solve the governing equations. Finally, two simple tests are shown, one considering diffusion, sorption, and radioactive decay, the other involving diffusion and a radioactive decay chain. In both tests, the code shows good agreement with the reference solutions. Caveats of the model and future additions are discussed.


Introduction
In the German site selection procedure for a high-level nuclear waste repository preliminary safety assessments have to be evaluated in three consecutive phases. In the first phase (starting 2020) representative preliminary safety assessments have to be carried out for a large number of potential sites (Hoyer et al., 2021;StandAG, 2017), with three different types of host rock considered (clay rock, rock salt, crystalline rock). The legal requirements include strict criteria on both the maximum mass and the maximum amount of radionuclides that are allowed to be released from the repository or its surroundings over one million years (10 −4 of the total mass and amount initially disposed), and a maximum rate of release (10 −9 yr −1 of the total mass and amount initially disposed). These criteria need to be evaluated individually for each potential site. The detailed legal requirements can be found in EndlSiAnfV (2020) (Disposal Safety Requirements Ordinance) and StandAG (2017) (Repository Site Selection Act).
In general, the large number of radionuclides to consider (several hundreds up to a thousand) make such evaluations difficult. Apart from the complexity of considering three different types of host rock, the particular challenge in the German site selection procedure is two-fold: First, the number (and area) of potential sites is large, with 54 % of the area of Germany being a potential site after the initial stage of the site selection procedure (BGE, 2020). Secondly, data for parametrizing transport/reaction models is both heterogeneous and scarce. Some geological information and basic parameters such as porosity may be available site-specific from different sources. However, other parameters governing the details of the transport process are widely unavailable from local measurements and need to be estimated in a systematic way, as on-site exploration/measurements will only begin in subsequent stages of the site selection procedure. Given the sparse data, it is beneficial to employ a computationally efficient 1D code that can be used to run large sets of models in order to capture the uncertainties in the governing parameters.
At present, a workflow for estimating the transport parameters is set up, which is inspired by the one developed by Nagra for the Swiss site selection procedure (Van Loon, 2014, note that these authors focus on clay formations only). Details can be found in BGE (2022a).
In the literature, several models and codes for assessing the migration of radionuclides out of a deposit exist, e.g. Norman and Kjellbert (1990), Reiche (2016), Nagra (2013), Garibay-Rodriguez et al. (2022), Trinchero et al. (2020). For the safety assessments in the German site selection procedure, one constraint was the desire to be able to publish the code as open-source. Additional requirements for the code include simplicity, ease of use, and performance. The code should be light-weight and should not rely on large thirdparty tools to reduce the complexity. In the version presented here, the code focuses on porous media and does not include a treatment of fractured media, which could be important for transport in crystalline rocks.
Our transport model for the representative preliminary safety assessments and its implementation is named TransPyREnd ("Transportmodell in Python für Radionuklide aus einem Endlager", English translation "Transport model in Python for radionuclides in a nuclear waste repository"). It will be released as an open-source project in early 2023, and will provide a unique opportunity for the community and public to assess the model, test and/or modify it. This will both help in assuring quality of the tools used in the site selection procedure, and aid in trust-building for the general public.

Model approach
A number of choices have to be made to develop a model for the transport of radionuclides. Here only one spatial dimension (1D) is considered, that is, we model transport along a line, e.g. in the vertical, horizontal, or any other meaningful direction. Given the limited available data and the large number of potential sites, this is a justified approximation also applied in other safety assessments, e.g. Nagra (2002). A 1D model also has advantages in terms of the low computational effort needed. Hence, uncertainties can easily be incorporated in the parametrization by evaluating ensembles of model runs, i.e. sets of models in which uncertainties are accounted for by varying the parameters within a distinct range. The general concept of dealing with uncertainties in the preliminary safety assessments can be found in BGE (2022b). In short, the workflows used for generating the model parameters are also used to estimate the uncertainties in these parameters. Sensitivity analyses are used to assess which parameters are important: the combined uncertainties in sensitive parameters are then evaluated using Monte-Carlo simulations.
Furthermore, the scope of the model is limited to the transport through the geological barriers far from the emplacement area (also termed as the "far-field" in contrast to the "near-field" that would include the technical barriers near the emplacement area). The transport along e.g. mining shafts or other technical elements is currently not included in the model.
Transport is included in the host rock layer and in geological layers above and below (if transport is vertical) or next (in case of horizontal transport) to the host rock layer. A sketch of an example model domain for vertical transport is shown in Fig. 1. As the current application is limited to the question of how much of the disposed radionuclides escape from the host rock, following the legal and regulatory demands (Endl-SiUntV, 2020), what we extract from the completed model is the amount and mass of radionuclides released from the host rock layer. Nevertheless, the model contains all relevant geological layers, as these have an effect on the release from the host rock.
Regarding the modeled processes, diffusion, advection, mechanical dispersion, and sorption are considered. Diffusion is assumed to be Fick'ian (Tartakovsky and Dentz, 2019), and differing accessible porosities (Van Loon, 2014) for individual nuclides (depending on whether they are present as anions or cations) are taken into account, assuming a fixed speciation. We assume Darcy's law holds for advection, as flow velocities in the host rock are expected to be very small due to the low hydraulic permeability expected in the host rocks (in particular in rock salt and clay). Consequently, the model is not applicable where turbulent flows are expected. All pore space is assumed to be saturated, following e.g. Nagra (2014). Ignoring the unsaturated zone introduces only a small error in most cases. In arid environments, the unsaturated zone might be very large, but even in this case, the assumption will still be conservative: It will lead to an overestimation of the radionuclide release, not to an underestimation. Dispersion can be included as an effectively increased diffusion. Sorption is modeled as linear equilibrium sorption (Freundlich, 1907), i.e. we do not consider competition for sorption sites.
Additionally, we include radioactive decay and the formation of daughter nuclides, as required by the regulatory stipulations. The inventory of nuclear waste that has to be stored contains hundreds of different radionuclides: however, many of them can be shown to be not relevant for the effective release of radionclides. For example, a radionuclide with a half-life of hours will not be transported significantly in the designated host rock before it decays, even in a pessimistic scenario. However, its daughter nuclide might be relevant, if it has a long half-life. At present, a simplified nuclide scheme following this line of argument is employed, taking into account 47 radionuclides Fischer-Appelt et al., 2013). It neglects nuclides with small half-lifes, but includes long-lived daughters of short-lived radionuclides. The 47 nuclides in the scheme are shown in Table 1.
Several additional physical processes and their coupling are not taken into account which is discussed together with other simplifications at the end of Sect. 4. The model computes the evolution of concentration c i for each nuclide i over time.
The aforementioned processes can be included into a set of coupled differential equations, namely one equation for each of the n s radionuclides considered. The 1D transport equation for the nuclide with the index i reads (compare i.e. Bear, 1972;de Marsily, 1986;Kinzelbach, 1992;Clauser, 2003): Here, c i is the amount of nuclide i per fluid volume, measured in mol m −3 . R i is the retardation factor, defined as and is the factor by which sorption slows down transport due to advection and diffusion. K d,i is the linear equilibrium sorption coefficient, ρ b is the bulk density of the rock, and φ i is the accessible porosity. D e,i is the effective diffusion coefficient, with D i being the diffusion coefficient for nuclide i. If an additional dispersion term is to be considered, D e,i should be replaced by with D disp being the dispersion. λ i,j is the reaction rate for the decay of nuclide j to nuclide i, and i is the total decay rate of nuclide i. The total decay rate i is related to the nuclide's half-life t 1/2 by i = log 2t −1 1/2 . The individual rates λ i,j are connected to this total decay rate by the branching fraction of the relevant decay channel. Finally, W i is a generic source term for radionuclides, that is, it has units of amount of substance per volume and time and measures the inflow of radionuclides from an external source. This can be used, for example, to connect the model to a dedicated near field or canister release model. R i appears on both sides of the equation, but not in front of the advective and diffusive terms, reflecting that it expresses the factor by which both advection and diffusion are reduced due to sorption. The decay term consists of a term with a positive sign, which measures the rate at which other nuclides decay into nuclide i, and a term with negative sign, measuring the rate at which nuclide i itself decays into other nuclides, including those that are not explicitly tracked in the model (e.g. individual stable nuclides).
Note that we consider transport in various geological layers, so the retardation factor R i , effective diffusion coefficient K d,i , and porosity φ i are in general functions of location x. q is the Darcy velocity.
In the following, we do not consider steady state solutions to the system in Eq. (1). The physical problem we want to solve involves a finite total amout of radioactive waste, meaning that the steady-state solution is the trivial one in which all radionuclides have decayed to stable nuclides. Thus, we focus on the transient solution of the system of equations.

Numerical methods
In order to solve Eq. (1), we resort to a numerical treatment via a finite difference scheme. The finite difference method is a well established method for solving ordinary and partial differential equations (i.e. Shamir and Harleman, 1967;Marsal, 1989;Clauser, 2003;Rühaak et al., 2008;Luijendijk, 2012). We divide our domain into n nodes at location x k . The distance between grid points, h k = x k − x k−1 is allowed to vary in space. An illustration of the grid structure is shown in Fig. 2.
Accordingly, the transient problem in the time domain is split into specific time-steps.
The solution of transient partial differential equations requires initial and boundary conditions.

Discretization of derivatives
We begin by writing down discretized approximations of the right hand side of Eq. (1). As mentioned before we assume D e to vary with location, whereas q is assumed to be constant. In the following, we use the index k to denote the node at which we evaluate a quantity. As expected, special treatment of the boundary nodes (k = 0, k = n − 1) is necessary.
For the diffusive term, we use the following central difference scheme (e.g. LeVeque, 2007): Here, the index k ± 1 2 denotes the location half-way between node k and node k ± 1, x k±1/2 . The diffusion coefficient at these interfaces is evaluated using a harmonic mean (Huysmans and Dassargues, 2007). Clearly, in the case of a homogeneous diffusion coefficient and a constant grid spacing h k = h, Eq. (5) reduces to the familiar second order central difference.
The advection term is also approximated using central differences: Evaluation of the decay and source term in Eq. (1) is straight-forward, as they do not include derivatives.

Time-evolution scheme
To solve the equation, we need to choose a time-evolution scheme. We use a semi-implicit scheme governed by a meta parameter 0 ≤ θ ≤ 1 sometimes called the θ -method as it is robust and relatively fast.
The θ -method is a one-step scheme and a combination of the Euler forwards and backwards schemes. For θ = 0.5 it becomes the second-order accurate Crank-Nicolson scheme (Crank and Nicolson, 1947). Schematically, for a given differential equation of the form the θ -scheme yields where F is an appropriate discretization of f , and where we use the index m to denote the discrete time step m at which a quantity is evaluated. For θ = 0, this becomes the Eulerforwards scheme which is prone to numerical instability. For θ = 1, it becomes the Euler-backwards scheme which is unconditionally stable, and for θ = 0.5 we recover the Crank-Nicolson scheme which is second order accurate and unconditionally stable as well.
Inserting our discretizations from Eqs. (5) and (6), as well as the decay and source term into Eq. (8), we obtain: Where we have defined As can be seen from Eq. (9), the decay terms are evaluated at time step m and not m+1, like the additional external source term. This is necessary to deal with the dependency of c i on potentially all other c j with j = i that could decay into c i . As the currently used nuclide chains are relatively simple, it would be possible to sort the equations in a way that allows including the decay term in an implicit fashion. We leave this to future work. By defining and rearranging the terms in Eq. (9) such that unknown terms are on the left, we find

Matrix equation
Identifying the right hand side with a vector b, this can be represented as a matrix equation, M i c m+1 i = b i , for each nuclide i. M i is a tridiagonal matrix with the structure: The terms B l,0 , B l,1 , B r,n−2 , B r,n−1 control how the boundary nodes -the nodes on the left and right edge of the model domain -evolve and are set by the boundary conditions. Again, it becomes evident here why we treat the decay terms as an external source term evaluated at time step m: otherwise, we could not write the equation as system of linear equations, because the decay term couples the concentration of different species with each other.

Boundary conditions
Dirichlet boundary conditions are used to fix the concentration at the outer edges of the domain (here termed left and right) to a constant value, i.e. c m i,0 = c m+1 i,0 = C i,left , c m i,n−1 = c m+1 i,n−1 = C i,right yielding: B l,0 = 1, B l,1 = 0, B r,n−2 = 0, B r,n−1 = 1 and the vector entries are set to the specified concentrations, b 0 = C i,left or b n−1 = C i,right . Neumann boundary conditions can be considered as well. For the left boundary at k = 0, they are defined by where we have applied a central differencing scheme on the right hand side. The node at k = −1 is a fictitious grid point outside of the domain, and we set h 0 = h 1 . Furthermore, we assume D k − 1 2 = D k + 1 2 . We can use Eq. (23) to eliminate the fictitious grid point from the transport equation for the boundary, which results in where we have omitted the decay and source terms for brevity. Consequently, the boundary terms for Neumann boundary conditions on the left boundary read: A similar derivation can be made for Neumann boundaries on the right boundary of the model.

Initial conditions
We assume instantaneous release of the radionuclides from canisters and the emplacement area and neglect solubility limits. Thus, the initial conditions are simply given by the total amount of each radionuclide in the repository and the volume of the repository. Future development will be focused on overcoming this simplification, see Sect. 4.

Implementation
The matrix Eq. (21) can be solved using standard methods. Our implementation is based on the PyBasin model (Luijendijk, 2020). Given the sparse structure of the matrix, we employ the spsolve solver in the Python package scipy , which in our case uses the SuperLU (Li, 2005) library to perform an LU decomposition followed by Gaussian elimination. For testing purposes and benchmarking, we also implemented an interface for external solvers, namely scipy's solve_ivp, that operate directly on the right hand side of Eq. (9).
Since the nuclide scheme with its 47 nuclides contains a few linear decay chains and many activation or fission products that directly decay to a stable nuclide, many of the nuclides can be solved for in parallel. For example, all of the activation and fission products (see Table 1) in our model decay directly to stable nuclides. As in our model, they do not interact with each other nor with other nuclides, we can as well run the model for each nuclide individually on a different CPU, parallelizing the work. For the nuclide chains, this is not possible, as the different nuclides' concentrations are coupled: however, one model can be run for each of the (four) chains in parallel.
To exploit this, we make use of the pathos (McKerns et al., 2012) package and run the model in parallel for appropriate chunks of the nuclide scheme. To further speed up the calculation, we use just-in-time compilation through Numba (Lam et al., 2015).

Choice of time step size
Predefined time steps can be used with an initial "ramping up" of the time step size from a small initial value to a specified size t max . There are three timescales determining t max , namely the typical timescales of advection t a = x v max , diffusion t diff = x 2 D max , and decay t d ; v max is the maximum groundwater velocity in the setup, and D max is the maximum diffusion coefficient in the setup.
In the given geological context, the timescale associated with radioactive decay or advection is usually the shortest. The selection of t max = αmin(t a , t d ) ensures stability and physical correctness, with α < 1.

Radioactive decay and convergence
As the time step is usually determined by the radioactive decay of the shortest-lived radionuclide in question, a large number of time steps might be required to reach the designated accuracy. For some nuclides, this can become very expensive in terms of computational effort. It is possible to increase both accuracy and performance by applying operator splitting (i.e. Geiser, 2004) between the transport problem on the one hand (advection, diffusion, sorption) and radioactive decay on the other hand. The idea here is to evaluate the decay with smaller time steps than the employed global time step, and feed the result back into the transport scheme at every global time step. Generally, this could be done e.g. by employing a higher-order Runge-Kutta (e.g. LeVeque, 2007) solver for the sub-time steps. As we currently only consider linear decay chains with few convergent branches, we instead use the analytic solution given by Bateman (1910).
Consider an ordered, linear decay chain with R members in which each species r directly decays into species r + 1. Assume that initially, only species r = 0 is present at a concentration c 0 = C 0 , whereas as c r = 0 mol m −3 for all other r. Then, the Bateman solution reads (Pressyanov, 2002): To include the more general case in which all c r are non-zero initially, the above equation is employed multiple times on the subchains and the results are summed.
For illustration, consider a chain with three members, R = 3. It can be decomposed into three subchains, containing the species (0, 1, 2), (1, 2), and (2). We can employ Eq. (26) on the three subchains, with the initial conditions (C 0 , 0, 0), (C 1 , 0), (C 2 ) and the indices shifted appropriately. Branched decay chains can be implemented in the same fashion, as long as the chain can be cut into linear subchains. We write the resulting Bateman solution asc m+1 i = f Bateman (c m i , t), where c m i is the concentration at time step m and t is the size of the time step. We put a tilde on the result,c m+1 i , as a reminder that the outcome concerns only the decay-part of the problem. Note that in this case, t max becomes independent of the half-life of the nuclides, resulting in larger allowed time steps.  (Van Genuchten, 1981), and the additional model parameters for the numerical treatment (bottom). To employ this operator splitting scheme in the transport model, we drop the original decay terms in Eq. (20), and instead introduce the Bateman solution averaged over the time step: Currently, the user has to decide if the operator splitting is to be used or not before starting a program run.

First tests
Testing TransPyREnd is an on-going effort that will be continuously documented. Here, we briefly show first tests on analytic or semi-numeric solutions.

Comparison with Van Genuchten (1982)
First, we show a comparison of TransPyREnd with an analytic solution for the transport of nuclides in porous media given by Van Genuchten (1981, 1982, Javandel et al. (1984). The solution includes decay, diffusion, sorption, and advection in a homogeneous geological layer. We choose the parameter values stated in Table 2 and use operator splitting for the radioactive decay as detailed in Sect. 2.2.8. The result is shown for various output times between 10 3 and 10 6 years in Fig. 3. We find a good agreement with the analytic solution: The root mean square error (RMSE) is 6×10 −5 mol m −3 . The total run time is about 17 s on a normal desktop machine. Note that the analytical solution is derived under the assumption of a semi-infinite domain, with the concentration gradient vanishing at infinity. Since we cannot exactly reproduce this  boundary condition in a finite domain, we instead chose a no-flow boundary condition on the right edge and made the domain considerably larger than the region where we compare with the analytic solution.
In Fig. 4, we show how the RMSE changes when varying both the time step size t and the grid spacing x. In the left panel, we vary the grid spacing and keep the time step size fixed, while in the right panel we vary the time step size and keep the grid spacing constant. In general, both methods perform similarly well for reasonable parameter choices. The advantage of the operator splitting, namely allowing for time step sizes larger than the half-life of the relevant nuclide, is not in play in this test, as the half-life is very large compared to the timescale for advection and diffusion. As can be seen in the right panel, the error made by treating the processes of diffusion/advection on the one hand and decay on the other hand quickly rises with growing time step size.

Testing Decay Chains
Secondly, we show a test involving both the radioactive decay in a decay chain and diffusion to test the mass conservation of our scheme. We arbitrarily place 1 mol of Cm-245 in the center of a domain with a total length of 1000 m and let the system evolve under decay and diffusion for 10 6 years. The boundaries are held fixed at a concentration of 0 mol m −3 . All species in the chain have the same effective diffusion coefficient of 10 −12 m 2 s −1 , the porosity is 0.05. The spatial resolution for this test is 1 m. We use the simplified Neptunium chain from Larue et al. (2013) here; the chain and the half-lifes are shown in Table 3, and show results both with and without operator splitting. Figure 5 shows the total amount of each of the chain members, that is, the concentration summed over the spatial domain, compared with the solution to the pure decay problem as calculated with the Python package radioactivedecay (Malins and Lemoine, 2022). Again, we are in good agreement, with an RMSE of 2×10 −6 mol m −3 both for the run with operator splitting and 5×10 −4 mol m −3 for the run without. The execution times differ as well: with operator splitting, the test takes about 4 s on a Desktop machine. Without operator splitting, the test takes about 10 s, due to the stronger constraint on the timestep. In this specific case, we use a timestep of 1000 years for the operator splitting case, and 216 years for the case without operator splitting. We also show a run without the operator splitting scheme for which we enforced a time step of 1000 years. In this case, the RMSE increases to 10 −2 mol m −3 . The deviations in the short-lived radionuclides become visible towards the end of the simulation time.
Naturally, the performance gain grows as the minimum half-life in the problem goes down. For example, for the Actinium decay chain, the allowed time step is of the order of 10 years, resulting in a runtime of 36 s without operator splitting, and 4 s with operator splitting. Enforcing the same time step as before will render the system unstable in this case. Thus, depending on the details of the decay chains in question, the operator splitting scheme can provide a relevant improvement in performance. Figure 5. Comparison of a simplified nuclide chain (Actinium chain) between radioactivedecay and TransPyREnd. While radioactivedecay (Malins and Lemoine, 2022) solves the pure decay problem, TransPyREnd was run with a simple diffusion+decay setup. In order to compare the two in terms of the conservation of amount of substance, the TransPyREnd concentrations were integrated over the domain to yield the total amount per species. Results are shown with and without operator splitting.

Conclusions and Outlook
In this paper, we have presented the radionuclide transport code TransPyREnd, developed specifically for application in the representative preliminary safety assessments of the German site selection procedure for high-level radioactive waste. The 1D model includes the processes diffusion, advection, sorption, and radioactive decay in the geological barrier. For the model we use a nuclide scheme  from the literature, pending updates from current BGE research projects. We have summarized the mathematical description of the model, the discretization steps required and our strategy to solve the equations numerically. The code is undergoing rigorous testing and comparison e.g. to OpenGeoSys . The benchmarks will be part of a forthcoming paper. TransPyREnd will be made public under an Open-Source license for the community and the interested public in early 2023.
As noted in Sect. 1, we have reduced the complexity of the problem by neglecting, in particular, hydraulic, geochemical, mechanical, and thermal processes as well as their coupling. As an example for the possible effect of these additional coupled processes, consider the heating from high-level radioactive waste in the first few hundred years, peaking at up to a few kW per fuel element. The heating has a potential effect on the mechanical state of the deposit and its surroundings (due to thermal expansion) which in turn can affect e.g. the fluid pressure and the porosity of the host rock layer. Both parameters affect the migration of radionuclides.
Our current approach of neglecting these effects is justified by the degree of knowledge about the parameters governing these processes at the current stage of the site selection procedure. However, it is clear that in later stages, these processes need to be carefully addressed to assess their importance to radionuclide migration.
Also, the one-dimensional approach might become less appropriate at later stages, in particular to assess effects arising from complex geometries. While the latter is generally not expected to lead to an increased discharge of radionuclides in the model, the former might have effects in both directions, lowering or enhancing the release of radionuclides in the geosphere/biosphere. At the current stage, the arising uncertainties (i.e. Bjorge et al., 2022) could be considered by systematically varying the model parameters, in particular the diffusion and sorption coefficients. Apart from these considerations, it should be noted that decay processes could also be solved in a semi-implicit fashion by taking advantage of the simple, linear decay chains in our nuclide scheme. The performance of the code could be be imprived by the implementation of a more specialized solving method like the Thomas algorithm (Thomas, 1949). We leave this to future work. Additionally, a near-field model which could consider for instance retention by backfill material is planned for the future.
In the next phases of the site-selection procedure more complex 3D numerical models will be computed. For this aim the OpenWorkFlow project was initiated by BGE (e.g. see https://www.ufz.de/index.php?en=48378, last access: 13 January 2023).
Code and data availability. TransPyREnd will be made public under an Open-Source license.
Author contributions. CB and EL developed the model code, contributed to writing, conceptualizing the manuscript, ran the calculations and created the plots. EL reviewed the manuscript. SM contributed to writing, calculations, and review. MG contributed to writing, review, and created Fig. 1. FP and PK contributed to writing and review. MB and AR reviewed the paper. WR contributed to conceptualizing, writing, and reviewing the manuscript, and did supervision.