Articles | Volume 58
Adv. Geosci., 58, 109–119, 2023
Adv. Geosci., 58, 109–119, 2023
18 Jan 2023
18 Jan 2023

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

TransPyREnd: a code for modelling the transport of radionuclides on geological timescales
Christoph Behrens1, Elco Luijendijk1,2, Phillip Kreye1, Florian Panitz1, Merle Bjorge1, Marlene Gelleszun1, Alexander Renz1, Shorash Miro1, and Wolfram Rühaak1 Christoph Behrens et al.
  • 1Bundesgesellschaft für Endlagerung (BGE) mbH, Eschenstraße 55, 31224 Peine, Germany
  • 2now at: Department of Earth Science, University of Bergen, 5020 Bergen, Norway

Correspondence: Wolfram Rühaak ( and Christoph Behrens (


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 finite-differences code for modeling the transport of radionuclides in the subsurface 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.

1 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; StandAG2017), 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 (BGE2020). 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 Loon2014, 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 third-party 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 (EndlSiUntV2020), 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 Dentz2019), and differing accessible porosities (Van Loon2014) 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 (Freundlich1907), 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 (Larue et al.2013; 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.

Table 1The nuclides included in the simplified nuclide scheme of Larue et al. (2013), Fischer-Appelt et al. (2013).

Download Print Version | Download XLSX

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.

Figure 1Schematic depiction of the model domain. Note that the direction of transport is not necessarily vertical as shown here. Note that while the zoom-in only shows the host rock, all geological layer are discretized.


2 Mathematical approach

2.1 Governing equations

The model computes the evolution of concentration ci 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 ns radionuclides considered. The 1D transport equation for the nuclide with the index i reads (compare i.e. Bear1972; de Marsily1986; Kinzelbach1992; Clauser2003):

(1) j ϕ i R i c i t Concentration change = x ( j D e , i c i x Diffusion - j q c i Advection ) + j j R j ϕ i c j λ i , j - R i ϕ i c i Λ i Decay + j W i Source term

Here, ci is the amount of nuclide i per fluid volume, measured in mol m−3. Ri is the retardation factor, defined as

(2) R i = 1 + K d , i ρ b ϕ i ,

and is the factor by which sorption slows down transport due to advection and diffusion.

Kd,i is the linear equilibrium sorption coefficient, ρb is the bulk density of the rock, and ϕi is the accessible porosity. De,i is the effective diffusion coefficient,

(3) D e , i = D i ϕ i ,

with Di being the diffusion coefficient for nuclide i. If an additional dispersion term is to be considered, De,i should be replaced by

(4) D e , i = D e , i + D disp

with Ddisp 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 t1/2 by Λi=log2t1/2-1. The individual rates λi,j are connected to this total decay rate by the branching fraction of the relevant decay channel. Finally, Wi 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. Ri 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 Ri, effective diffusion coefficient Kd,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.

2.2 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 Harleman1967; Marsal1989; Clauser2003; Rühaak et al.2008; Luijendijk2012). We divide our domain into n nodes at location xk. The distance between grid points, hk=xk-xk-1 is allowed to vary in space. An illustration of the grid structure is shown in Fig. 2.

Figure 2Illustration of the grid structure.


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.

2.2.1 Discretization of derivatives

We begin by writing down discretized approximations of the right hand side of Eq. (1). As mentioned before we assume De 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. LeVeque2007):

(5) x D e , i c x 2 h k + 1 + h k 1 h k D e , i , k - 1 2 c i , k - 1 - h k + 1 + h k h k h k + 1 D e , i , k - 1 2 + D e , i , k + 1 2 c i , k + 1 h k + 1 D e , i , k + 1 2 c i , k + 1

Here, the index k±12 denotes the location half-way between node k and node k±1, xk±1/2. The diffusion coefficient at these interfaces is evaluated using a harmonic mean (Huysmans and Dassargues2007). Clearly, in the case of a homogeneous diffusion coefficient and a constant grid spacing hk=h, Eq. (5) reduces to the familiar second order central difference.

The advection term is also approximated using central differences:

(6) q c x q c i , k + 1 - c i , k - 1 x k + 1 - x k - 1

Evaluation of the decay and source term in Eq. (1) is straight-forward, as they do not include derivatives.

2.2.2 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 Nicolson1947). Schematically, for a given differential equation of the form

(7) y t = f ( x , y , t )

the θ-scheme yields

(8) y m + 1 - y m Δ t = θ F m + 1 + ( 1 - θ ) F m

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 Euler-forwards 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 ci on potentially all other cj with ji that could decay into ci. 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


2.2.3 Matrix equation

Identifying the right hand side with a vector b, this can be represented as a matrix equation, Micim+1=bi, for each nuclide i. Mi is a tridiagonal matrix with the structure:

(21) M i = B l , 0 B l , 1 0 θ ( - u 1 - v 1 - ) 1 + θ v 1 θ ( u 1 - v 1 + ) 0 0 θ ( - u 2 - v 2 - ) 1 + θ v 2 θ ( u 2 - v 2 + ) 0 . 0 θ ( - u n - 2 - v n - 2 - ) 1 + θ v n - 2 θ ( u n - 2 - v n - 2 + ) 0 B r , n - 2 B r , n - 1

The terms Bl,0, Bl,1, Br,n-2, Br,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.

2.2.4 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. ci,0m=ci,0m+1=Ci,left,ci,n-1m=ci,n-1m+1=Ci,right yielding:

(22) 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, b0=Ci,left or bn-1=Ci,right.

Neumann boundary conditions can be considered as well. For the left boundary at k=0, they are defined by

(23) c i x | x = 0 = j l c i , 1 m - c i , - 1 m 2 h 0 ,

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 h0=h1. Furthermore, we assume Dk-12=Dk+12.

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:

(25) B l , 0 = 1 + 2 θ v 0 + , B l , 1 = - 2 θ v 0 + b 0 = 1 - 2 θ v 0 + c i , 0 m + 2 θ v 0 + c i , 1 m - 2 h 0 j l u 0 + v 0 + +

A similar derivation can be made for Neumann boundaries on the right boundary of the model.

2.2.5 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.

2.2.6 Implementation

The matrix Eq. (21) can be solved using standard methods. Our implementation is based on the PyBasin model (Luijendijk2020).

Given the sparse structure of the matrix, we employ the spsolve solver in the Python package scipy (Virtanen et al.2020), which in our case uses the SuperLU (Li2005) 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).

2.2.7 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 Δtmax. There are three timescales determining Δtmax, namely the typical timescales of advection ta=Δxvmax, diffusion tdiff=Δx2Dmax, and decay td; vmax is the maximum groundwater velocity in the setup, and Dmax 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 Δtmax=αmin(ta,td) ensures stability and physical correctness, with α<1.

2.2.8 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. Geiser2004) 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. LeVeque2007) 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 c0=C0, whereas as cr=0 mol m−3 for all other r. Then, the Bateman solution reads (Pressyanov2002):

(26) c r ( t ) = C 0 i = 0 i = r - 1 λ i × i = 0 r e - λ i t j = 0 ; j i r ( λ j - λ i )

To include the more general case in which all cr 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 (C0,0,0), (C1,0), (C2) 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 as c̃im+1=fBateman(cim,Δt), where cim is the concentration at time step m and Δt is the size of the time step. We put a tilde on the result, c̃im+1, as a reminder that the outcome concerns only the decay-part of the problem. Note that in this case, Δtmax becomes independent of the half-life of the nuclides, resulting in larger allowed time steps.

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:

(27) θ - u k - v k - c i , k - 1 m + 1 + 1 + θ v k c i , k m + 1 + θ u k - v k + c i , k + 1 m + 1 = θ u k + v k - c i , k - 1 m + 1 - θ v k c i , k m + θ - u k + v k + c i , k + 1 m + s f Bateman ( R i , k ϕ i , k c i , k m , Δ t ) - R i , k ϕ i , k c i , k m Δ t + W i , k m

Currently, the user has to decide if the operator splitting is to be used or not before starting a program run.

3 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.

3.1 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 103 and 106 years in Fig. 3.

Table 2The parameters used in the test (Van Genuchten1981), and the additional model parameters for the numerical treatment (bottom).

Download Print Version | Download XLSX

Figure 3Comparison between the analytic solution by Van Genuchten (1981) (lines) and TransPyREnd (circles) for various times.


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.

Figure 4Convergence of the Van Genuchten (1981) test. (a) RMSE as a function of grid spacing. (b) RMSE as a function of time step size.


3.2 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 106 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 m2 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.

Table 3The simplified Neptunium chain from Larue et al. (2013).

Download Print Version | Download XLSX

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 Lemoine2022). 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.

Figure 5Comparison of a simplified nuclide chain (Actinium chain) between radioactivedecay and TransPyREnd. While radioactivedecay (Malins and Lemoine2022) 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.


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.

4 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 (Larue et al.2013) 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 (Bilke et al.2022). 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.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 (Thomas1949). 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, 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.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “European Geosciences Union General Assembly 2022, EGU Division Energy, Resources & Environment (ERE)”. It is a result of the EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022.


We would like to thank Olaf Kolditz, Haibing Shao and Renchao Lu for fruitful discussions of the validation strategy. TransPyREnd makes use of Python (Van Rossum and Drake1995), IPython (Pérez and Granger2007), (Virtanen et al.2020), Numpy (Harris et al.2020), Numba (Lam et al.2015), Matplotlib (Hunter2007), Pandas (Pandas development team2020), Pathos (McKerns et al.2012), and radioactivedecay (Malins and Lemoine2022).

Review statement

This paper was edited by Michael Kühn and reviewed by two anonymous referees.


Bateman, H.: The solution of a system of differential equations occurring in the theory of radioactive transformations, Proc. Cambridge Philos. Soc., 15, 423–427, 1910. a

Bear, J.: Dynamics of Fluids in Porous Media, Dover Publications, Inc., New York, ISBN13 9780486656755, 1972. a

BGE: Zwischenbericht Teilgebiete gemäß § 13 StandAG, (last access: 13 January 2023), 2020. a

BGE: Konzept zur Durchführung der repräsentativen vorläufigen Sicherheitsuntersuchungen gemäß Endlagersicherheitsuntersuchungsverordnung, Tech. rep., BGE mbH, (last access: 13 January 2023), 2022a. a

BGE: Methodenbeschreibung zur Durchführung der repräsentativen vorläufigen Sicherheitsuntersuchungen gemäß Endlagersicherheitsuntersuchungsverordnung, Tech. rep., BGE mbH, (last access: 13 January 2023), 2022b. a

Bilke, L., Fischer, T., Naumov, D., Lehmann, C., Wang, W., Lu, R., Meng, B., Rink, K., Grunwald, N., Buchwald, J., Silbermann, C., Habel, R., Günther, L., Mollaali, M., Meisel, T., Randow, J., Einspänner, S., Shao, H., Kurgyis, K., Kolditz, O., and Garibay, J.: OpenGeoSys, Zenodo [code],, 2022. a

Bjorge, M., Kreye, P., Heim, E., Wellmann, F., and Rühaak, W.: The role of geological models and uncertainties in safety assessments, Environ. Earth Sci., 81, 190,, 2022. a

Clauser, C.: SHEMAT and Processing SHEMAT – Numerical Simulation of Reactive Flow in Hot Aquifers, Springer, Heidelberg-Berlin, Berlin, Heidelberg,, 2003. a, b

Crank, J. and Nicolson, P.: A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, Math. Proc. Cambridge, 43, 50–67,, 1947. a

de Marsily, G.: Quantitative Hydrogeology, Academic Press, San Diego, ISBN 9780122089169, 1986. a

EndlSiAnfV: Endlagersicherheitsanforderungsverordnung vom 6. Oktober 2020 (BGBl. I S. 2094), Bundesanzeiger Verl.-Ges., ISSN 0341-1095, 2020. a

EndlSiUntV: Endlagersicherheitsuntersuchungsverordnung vom 6. Oktober 2020 (BGBl. I S. 2094, 2130), Bundesanzeiger Verl.-Ges., ISSN 0341-1095, 2020. a

Fischer-Appelt, K., Baltes, B., Buhmann, D., Larue, J., and Mönig, J.: Synthesebericht für die VSG: Bericht zum Arbeitspaket 13. Vorläufige Sicherheitsanalyse für den Standort Gorleben, GRS-290, Tech. rep., Gesellschaft für Anlagen-und Reaktorsicherheit (GRS) mbH, Köln, ISBN 978-3-939355-66-3, 2013. a, b

Freundlich, H.: Über die Adsorption in Lösungen, Z. Phys. Chem., 57U, 385–470,, 1907. a

Garibay-Rodriguez, J., Chen, C., Shao, H., Bilke, L., Kolditz, O., Montoya, V., and Lu, R.: Computational Framework for Radionuclide Migration Assessment in Clay Rocks, Front. Nucl. Eng., 1,, 2022. a

Geiser, J.: Diskretisierungsverfahren für Systeme von Konvektions-Diffusions-Dispersions-Reaktions-Gleichungen und Anwendungen, PhD Thesis, Naturwissenschaftlichen-Mathematischen Gesamtfakultät der Ruprecht-Karls-Universität Heidelberg, Heidelberg, (last access: 13 January 2023), 2004. a

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array programming with NumPy, Nature, 585, 357–362,, 2020. a

Hoyer, E.-M., Luijendijk, E., Müller, P., Kreye, P., Panitz, F., Gawletta, D., and Rühaak, W.: Preliminary safety analyses in the high-level radioactive waste site selection procedure in Germany, Adv. Geosci., 56, 67–75,, 2021. a

Hunter, J. D.: Matplotlib: A 2D graphics environment, Comput. Sci. Eng., 9, 90–95,, 2007. a

Huysmans, M. and Dassargues, A.: Equivalent diffusion coefficient and equivalent diffusion accessible porosity of a stratified porous medium, Transport Porous Med., 66, 421–438,, 2007. a

Javandel, I., Doughty, C., and Tsang, C. F.: Groundwater Transport: Handbook of Mathematical Models, vol. 10 of Water Resources Monograph, AGU, ISBN 9780875903132, 1984. a

Kinzelbach, W.: Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser, 2nd Edn., Oldenbourg, ISBN 9783486263473, 1992. a

Lam, S. K., Pitrou, A., and Seibert, S.: Numba: A LLVM-Based Python JIT Compiler, in: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in: HPC, LLVM '15, Association for Computing Machinery, New York, NY, USA,, 2015. a, b

Larue, J., Baltes, B., Fischer, H., Frieling, G., Kock, I., Navarro, M., and Seher, H.: Radiologische Konsequenzenanalyse, GRS-289, Tech. rep., Gesellschaft für Anlagen-und Reaktorsicherheit (GRS) mbH, Köln, ISBN 9783939355656, 2013. a, b, c, d, e

LeVeque, R. J.: Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems, SIAM, ISBN 9780898716290 ,2007. a, b

Li, X. S.: An Overview of SuperLU: Algorithms, Implementation, and User Interface, ACM T. Math. Software, 31, 302–325, 2005. a

Luijendijk, E.: 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, PhD thesis, Vrije Universiteit Amsterdam, 2012. a

Luijendijk, E.: PyBasin: Numerical model of basin history, heat flow and thermochronology, Zenodo [code],, 2020. a

Malins, A. and Lemoine, T.: radioactivedecay: A Python package for radioactive decay calculations, J. Open Source Softw., 7, 3318,, 2022. a, b, c

Marsal, D.: Finite Differenzen und Elemente, Numerische Lösung von Variationsproblemen und partiellen Differentialgleichungen, Springer-Verlag, Berlin, ISBN 978-3540501923, 1989. a

McKerns, M. M., Strand, L., Sullivan, T., Fang, A., and Aivazis, M. A. G.: Building a Framework for Predictive Science, arXiv [preprint],, 2012. a, b

Nagra: Project Opalinus Clay, Technical Report 02-05, ISSN 1015-2636, 2002. a

Nagra: PICNIC-TD, User Guide for Version 1.4, Arbeitsbericht NAB 13-21, 2013. a

Nagra: Provisional Safety Analyses for SGT Stage 2: Models, Codes and General Modelling Approach, Technical Report 14-09, ISSN 1015-2636, 2014. a

Norman, S. and Kjellbert, N.: FARF31 – a far field radionuclide migration code for use with the PROPER package, Tech. rep., Report Number SKB-TR-90-01, Svensk Kärnbränslehantering AB, (last access: 13 January 2023), 1990. a

Pandas development team: pandas-dev/pandas: Pandas, Zenodo [code],, 2020. a

Pérez, F. and Granger, B. E.: IPython: a System for Interactive Scientific Computing, Comput. Sci. Eng., 9, 21–29,, 2007. a

Pressyanov, D. S.: Short solution of the radioactive decay chain equations, Am. J. Phys., 70, 444–445,, 2002. a

Reiche, T.: RepoTREND Das Programmpaket zur integrierten Langzeitsicherheitsanalyse von Endlagersystemen, GRS-413, Tech. rep., Gesellschaft für Anlagen-und Reaktorsicherheit (GRS) mbH, ISBN 9783944161952, 2016. a

Rühaak, W., Rath, V., Wolf, A., and Clauser, C.: 3D finite volume groundwater and heat transport modeling with non-orthogonal grids, using a coordinate transformation method, Adv. Water Resour., 31, 513–524,, 2008. a

Shamir, U. Y. and Harleman, D. R. F.: Numerical solutions for dispersion in porous mediums, Water Resour. Res., 3, 557–581,, 1967. a

StandAG: Standortauswahlgesetz vom 5. Mai 2017 (BGBl. I S. 1074) (das zuletzt durch Artikel 1 des Gesetzes vom 7. Dezember 2020 (BGBl. I S. 2760) geändert worden ist), Bundesanzeiger Verl.-Ges., ISSN 0341-1095, 2017 (updated 2020). a, b

Tartakovsky, D. M. and Dentz, M.: Diffusion in Porous Media: Phenomena and Mechanisms, Transport Porous Med., 130, 105–127,, 2019. a

Thomas, L. H.: Elliptic problems in linear difference equations over a network, Watson Sci. Comput. Lab. Rept., Columbia University, New York, 1, 71 pp., 1949. a

Trinchero, P., Cvetkovic, V., Selroos, J.-O., Bosbach, D., and Deissmann, G.: Upscaling of radionuclide transport and retention in crystalline rocks exhibiting micro-scale heterogeneity of the rock matrix, Adv. Water Resour., 142, 103644,, 2020. a

Van Genuchten, M.: Analytical solutions for chemical transport with simultaneous, zero-order production and first-order decay, J. Hydrol., 49, 213–233,, 1981. a, b, c, d

Van Genuchten, M. T.: One-dimensional analytical transport modeling, in: Symposium on Unsaturated How and Transport Modeling, Rep. PNL-SA-10325, Pacific Northwest Lab., Richland, Wash., also available as U.S. Nuc./Regul. Comm. Rep. NUREG/CP-0030, 231–248, 1982.  a, b

Van Loon, L.: Effective diffusion coefficients and porosity values for argillaceous rocks and bentonite: measured and estimated values for the provisional safety analyses for SGT-E2, Tech. Rep. 1015-2636, ISSN 1015-2636, (last access: 13 January 2023), 2014. a, b

Van Rossum, G. and Drake Jr., F. L.: Python tutorial, Centrum voor Wiskunde en Informatica Amsterdam, Rep. No. CS-R9526, ISSN 0169-118X, Amsterdam, Netherlands (last access: 13 January 2023), 1995. a

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 17, 261–272,, 2020. a, b

Short summary
The mathematical basics of a numerical code developed specifically for the search of a site for high-level radioactive waste in Germany is presented. The code is developed in accordance to the specific regulations. First tests of the code are shown.