Reactive transport simulations of uranium migration in the Opalinus Clay depend on ion speciation governed by underlying thermodynamic data

Safety assessments must demonstrate that radionuclides in potential disposal sites are retained within the containment providing rock zone. The impact of thermodynamic data on calculated migration lengths resulting from reactive transport simulations is quantified for the example of uranium in the hydrogeological system of the Opalinus Clay at Mont Terri. In this geochemical system, speciation is controlled by the calcite-carbonate-ion system. Aqueous uranium is mainly present as U(VI) as ternary complexes with calcium or magnesium together with carbonate. Previous simulations using the first NEA update of thermodynamic data for uranium indicated that the anionic complex ${\mathrm{CaUO}}_{\mathrm{2}}({\mathrm{CO}}_{\mathrm{3}}{)}_{\mathrm{3}}^{\mathrm{2}-}$ is the predominant species with a maximum migration distance of 50 m after one million years. The NEA published an update of the thermodynamic data for uranium, what, in turn, changes the predominant species from anionic to almost only the neutral ternary complex Ca₂UO₂(CO₃)₃. With identical simulations, except for the application of the second NEA update, a maximum distance of 80 m was obtained. This can be attributed to a decrease in sorption capacity due to a stronger complexation of uranium with calcium and carbonate. Therefore, the impact of the change in the underlying thermodynamic data can be quantified with +30 m. Our work clearly shows how sensitive migration lengths resulting from reactive transport simulations are to the model conceptualisation and selection of underlying data. Consequently, the compilation and further development of data sets and a site specific investigation are indispensable for reliable outcomes of transport simulations, and thus of performance assessments.

1 Introduction

Isolation of high-level radioactive nuclear wastes from the environment is to be ensured by the disposal in deep geological formations (IAEA, 2003). For safety assessments in Germany, it is of particular relevance to demonstrate that the thickness of the containment providing rock zone (CPRZ) is sufficient to retain the radionuclides for time periods of up to one million years. A thickness of at least 100 m is prescribed by law for this purpose (§ 23 Article 5 StandAG¹). To cover the temporal and spatial scales required in the context of safety assessments, numerical simulations are indispensable to quantify radionuclide migration lengths.

Argillaceous formations are among favoured host rocks due to their low permeability only allowing diffusive transport of radionuclides that can be retarded by sorption processes taking place on the surfaces of the inherent clay minerals. Both processes depend on pore water geochemistry and mineralogy. Reactive transport simulations numerically solved with geochemical codes, such as PHREEQC (Parkhurst and Appelo, 2013), enable a process-based quantification of diffusion lengths and sorption effects as a function of both, pore water geochemistry and mineralogy, by the application of advanced approaches like surface complexation modelling. For this, comprehensive thermodynamic databases are required including the relevant species, mineral phases, surface complexes, exchangeable species and associated stability constants necessary to describe the investigated geochemical system and radionuclide. Therefore, changes in the underlying thermodynamic data affect simulated migration lengths. A comprehensive and up-to-date thermodynamic database is essential to quantify radionuclide migration lengths by the application of advanced approaches with reactive transport simulations.

Migration of uranium, the main component of spent fuel, in the potential host rock Opalinus Clay has been quantified as a function of mineralogical heterogeneities between three litho-facies (shaly, sandy and carbonate-rich, Fig. 1) and for varying partial pressures of carbon dioxide (pCO₂, Hennig et al., 2020; Hennig and Kühn, 2021a) within the hydrogeological system at Mont Terri (Hennig and Kühn, 2021b). It has been shown that the application of the process-based multi-component diffusion approach (Appelo and Wersin, 2007) provided no additional profit on the host rock scale. Calculated migration lengths did not differ significantly (±5 m) to simulations following Fick's law, and thus experimentally determined effective diffusion coefficients can be used instead of the computationally intensive multi-component diffusion approach. Sorption must be quantified as a function of the geochemical conditions. Further, the results showed that present-day pore water profiles can be explained by diffusive exchange with the embedding aquifers (Mazurek et al., 2011; Hennig and Kühn, 2021b). Previous simulations were conducted using the PSI/Nagra thermodynamic database (Thoenen et al., 2014) including the Nuclear Energy Agency (NEA) data for uranium published by Guillaumont et al. (2003). With this, the anionic ternary uranyl complex ${\mathrm{CaUO}}_{\mathrm{2}}({\mathrm{CO}}_{\mathrm{3}}{)}_{\mathrm{3}}^{\mathrm{2}}-$ is the predominant species in the geochemical system. In 2020, the NEA published an update of the thermodynamic data for uranium (Grenthe et al., 2020) including recent data from the literature. In particular, the values for the stability constants of the ternary complexes with calcium and carbonate were revised. As shown by Hennig and Kühn (2021a), this might lead to significant changes in the speciation, and thus resulting migration lengths within the various facies of the Opalinus Clay, as uranium speciation is dominated by ternary complexes with calcium and carbonate. For instance, a difference of 1.33 log units in the value for the stability constant of the neutral complex Ca₂UO₂(CO₃)₃ changes the migration lengths by 5 to 7 m for the sandy and shaly facies, respectively. The next step, presented here, is the investigation what such a change in underlying thermodynamic data mean in metres difference on the host rock scale after one million years quantified with one-dimensional reactive transport models based on simulations of the hydrogeological system of the Opalinus Clay at Mont Terri (Hennig and Kühn, 2021b). For this, resulting migration lengths are compared between simulations with varying thermodynamic data sets of the NEA for uranium (Guillaumont et al., 2003; Grenthe et al., 2020).

2 Methods

One-dimensional diffusion simulations following Fick's laws were conducted with the geochemical code PHREEQC Version 3.5.0 (Parkhurst and Appelo, 2013). The thermodynamic data is based on the PSI/Nagra database version 12/07 (Thoenen et al., 2014). This database was updated with the thermodynamic data for uranium recently published by the NEA (Grenthe et al., 2020) and supplemented with the sorption data for the hydroxo-complexes of uranium on the inherent clay minerals. This sorption data set is based on published K values and was compiled in the context of an uranium sorption experiment in Opalinus Clay (Joseph et al., 2013a). All surface parameters and reactions can be found in Table 1 of the Supplement of Hennig et al. (2020). Aqueous uranium speciation in the modelled system is dominated by ternary complexes with carbonate and calcium (Hennig et al., 2020; Hennig and Kühn, 2021a; Joseph et al., 2011) that are not considered for sorption (Joseph et al., 2013a). Diffusion and sorption experiments of uranium on clay minerals could be modelled without the requirement to include uranium surface complexes with carbonate as they do not sorb or only weakly (Joseph et al., 2013a; Hennig et al., 2020; Hennig and Kühn, 2021b; Stockmann et al., 2022; Tournassat et al., 2018). As ionic strengths in the hydrogeological system of the Opalinus Clay at Mont Terri do not exceed 0.5 mol L⁻¹ (Pearson et al., 2003; Hennig and Kühn, 2021b), the Davies approach (Davies and Shedlovsky, 1964) is used for the database to represent ion-ion interactions (Stockmann et al., 2017; Noseck et al., 2018). After closing of the repository, elevated temperatures are expected due to heat generated by the waste packages. However, this does not impact the migration behaviour of uranium (Joseph et al., 2013b). Based on diffusion experiments for 25 and 60 ^∘C (Joseph et al., 2013b), it has been shown that diffusion as well as sorption increase with temperature, and thus compensate each other. Furthermore, changes in pore water composition for higher temperatures (45 ^∘C) are <5 % (Wersin et al., 2016). Therefore, all simulations were conducted isothermal for a reference temperature of 25 ^∘C according to the database.

The hydrogeological system at Mont Terri is characterised by diffusive exchange over millions of years between the embedding aquifers in the Dogger and Lias and the 210 m thick CPRZ consisting of Opalinus Clay and Liassic shales (Pearson et al., 2003; Mazurek et al., 2011). The Mont Terri anticline was uplifted during the Jura folding about 10 Ma ago and subsequent erosion enabled freshwater infiltration into both aquifers, and thus activated the current hydrogeological boundaries. Only after erosion of the overlying stratigraphic layers, the footwall aquifer was activated leading to a substantial time lag in freshwater infiltration, and thus to an asymmetric shape of present-day pore water profiles (Pearson et al., 2003; Mazurek et al., 2011). Based on purely diffusion-driven simulations (Mazurek et al., 2011, 2009; Hennig and Kühn, 2021b), modelled profiles match with pore water data measured at the underground laboratory Mont Terri (Pearson et al., 2003).

The conceptual model applied here is equal to the one presented in Hennig and Kühn (2021b) despite the application of the updated thermodynamic data for uranium published by the NEA (Grenthe et al., 2020). Therefore, the simulations of the hydrogeological system at Mont Terri presented by Hennig and Kühn (2021b) were used to represent the initial geochemical conditions in the investigated system as the NEA update only affects the speciation of uranium. Calcium, magnesium, pH and pCO₂ are the governing parameters for uranium speciation (Hennig et al., 2020; Hennig and Kühn, 2021a), and therefore only these modelled profiles are shown as an example for the other pore water components (Fig. 1). With differences <15 %, all modelled profiles match with the measured data. Within the CPRZ, pH and pCO₂ are coupled via the calcite-carbonate-ion system (Pearson et al., 2003; Hennig et al., 2020). The scattering of measured values can be explained by the high sensitivity of both parameters, and thus associated uncertainties with measurements (Wersin et al., 2022). Therefore, the profiles modelled by Hennig and Kühn (2021b) with differences <12 % to the measured values were considered as sufficient to represent the initial conditions within the system. The model boundaries were defined by the water composition of the aquifers as Dirichlet-conditions at the model outlets. Numerical stability is ensured by the Neumann criteria. The uranium source term to represent the failed high-level waste canisters is integrated via a mineral equilibrium with uraninite (UO_2(am,hyd)) as well as pyrite and siderite to control the redox conditions. The saturation index of uraninite is chosen so that a concentration of 1 µmol L⁻¹ is applied in the centre of the model according to Keesmann et al. (2005) and Joseph et al. (2013b). For further details on the conceptual model as well as initial and boundary conditions for the simulation of the hydrogeological system, we refer to Hennig and Kühn (2021b).

Figure 1 Concentrations of pore water components decrease within the CPRZ (Opalinus Clay and Liassic Shales) towards embedding aquifers in the Liassic and Dogger. Initial geochemical conditions in the investigated system are represented by simulations of the hydrogeological system of the Opalinus Clay at Mont Terri conducted by Hennig and Kühn (2021b). Calcium, magnesium, pH and pCO₂ are the governing parameters for uranium speciation, and thus exemplary shown for the other pore water components. Deviations between modelled and measured data are given by the relative Root Mean Square Errors (rRMSE). Modified from Hennig and Kühn (2021b).

Sorption is integrated via a bottom-up approach (Marques Fernandes et al., 2015; Stockmann et al., 2017) using surface complexation models based on the two-layer model of Dzombak and Morel (1990) as well as cation exchange (Kim, 2001). An average clay mineral composition of 66 wt. % is used for the entire CPRZ as the geochemistry of the pore water is more decisive for uranium sorption processes than the amount of clay minerals (Hennig et al., 2020). The individual clay mineral quantities of illite, montmorillonite and kaolinite as well as the values for the hydro-physical parameters are equal to Hennig and Kühn (2021b). The distribution coefficient K_d (m³ kg⁻¹) is calculated from the PHREEQC results and used to evaluate the differences in the sorption capacity. The K_d is defined as ratio between the species concentration adsorbed on the solid phase and present in the liquid phase. Diffusion processes are quantified following Fick's law using the D_e with a value of $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{12}}$ m² s⁻¹ as determined by Joseph et al. (2013b) in an uranium diffusion experiment with Opalinus Clay. In previous studies (Hennig et al., 2020; Hennig and Kühn, 2021b), the present conceptual model was calibrated against this experiment. The neutral ternary complex was identified as predominant species in the diffusion experiment (Joseph et al., 2013b) based on speciation calculations using stability constants similar to the NEA data of Grenthe et al. (2020). The D_e of Joseph et al. (2013b) was used for all simulations. Accordingly, anion exclusion effects that would be reflected by a reduction of D_p are not considered. In the literature, no data on the anion-accessible porosity ratios for the different ternary uranyl complexes are available. Therefore, the reduction of D_p can only be estimated based on the porosity ratio for other species, like chloride. Furthermore, process-based simulations using the multi-component diffusion approach, which takes the interaction of charged species with the diffuse double layer into account (Appelo and Wersin, 2007), showed for the hydrogeological system at Mont Terri that anion exclusion effects only have a minor impact on the resulting migration lengths of uranium (Hennig and Kühn, 2021b).

3 Results

3.1 Speciation of uranium in the Opalinus Clay system

Speciation of uranium in the Opalinus Clay is controlled by the calcite-carbonate-ion equilibrium (Hennig et al., 2020; Hennig and Kühn, 2021a). Table 1 shows the speciation of uranium in the pore water of the shaly facies, i.e. at a depth of around 140 m corresponding to the highest calcium and magnesium concentrations (Fig. 1), using the second NEA update on thermodynamic data for uranium (Grenthe et al., 2020). The species distribution is compared with similar simulations (Hennig and Kühn, 2021a) using the previous published NEA data (Guillaumont et al., 2003).

Guillaumont et al. (2003)Grenthe et al. (2020)

Table 1 The aqueous species distribution changes from predominantly anionic to neutral ternary complexes with the application of the second NEA update for thermodynamic data of uranium (Grenthe et al., 2020) compared to the values published in Hennig and Kühn (2021a) using the data of the first NEA update (Guillaumont et al., 2003). Mol fraction (%) of main uranium species are exemplarily given for the pore water composition of the shaly facies and the predominant species is printed bold.

With the application of the second NEA update on thermodynamic data for uranium the species distribution changes from predominantly anionic to almost only the neutral ternary uranyl complex Ca₂UO₂(CO₃)₃ (Table 1). The anionic complexes with calcium, magnesium and strontium decreased from roughly 90 % to ≤12 %, and therefore only represent a minor proportion of the total. U(IV) species account for <2 %. Accordingly, uranium is mainly present as U(VI) with the neutral ternary complex as predominant species.

3.2 Uranium migration in the hydrogeochemical Opalinus Clay system at Mont Terri

Uranium migration in a geochemically heterogeneous system is quantified using the second NEA update on thermodynamic data for uranium (Grenthe et al., 2020) for the example of the hydrogeological system of the Opalinus Clay at Mont Terri (Hennig and Kühn, 2021b). Resulting migration lengths are shown in Fig. 2a (green line) and compared with simulations using the NEA data published by Guillaumont et al. (2003) to evaluate the impact of the second NEA update on the migration lengths (blue line, Hennig and Kühn, 2021b). Corresponding distribution coefficients K_d (m³ kg⁻¹) are given by the dots (Fig. 2b).

Figure 2 Changes in the speciation alter the sorption behaviour, and thus migration lengths. Uranium migrates farther through the Opalinus Clay (a) using the thermodynamic NEA data according to Grenthe et al. (2020) with the neutral ternary uranyl complex as predominant species (green lines) compared to the simulations of Hennig and Kühn (2021b), that are dominated by the anionic complex (blue lines), following the data of Guillaumont et al. (2003). The change in speciation is associated with a decrease in sorption capacity (b) represented by the distribution coefficients K_d (m³ kg⁻¹) calculated from PHREEQC.

In the simulations using the second NEA update on thermodynamic data for uranium (green lines, Fig. 2a) maximum migration distances were 80 and 70 m towards the Dogger and Liassic, respectively. Accordingly, uranium migrates about 30 m farther compared to the simulations presented in Hennig and Kühn (2021b) based on the first NEA update (blue lines, Fig. 2a) with maximum distances of 50 and 30 m towards Dogger and Lias, respectively. With the change in the predominant species due to the application of the second NEA update for uranium data (Grenthe et al., 2020), the sorption capacity represented by the K_d decreases, e.g. at a depth of 150 m from $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m³ kg⁻¹ (Fig. 2b).

4 Discussion

The impact of the second NEA update for the thermodynamic data for uranium (Grenthe et al., 2020) on the migration behaviour is quantified for the example of the geochemically heterogeneous Opalinus Clay system at Mont Terri. This is done by comparing resulting migration lengths with equal simulations except for the application of the first NEA update (Guillaumont et al., 2003).

Speciation calculations for the pore water composition of the shaly facies, i.e. at a depth of 140 m corresponding to the highest ionic strength (Fig. 1), showed that the predominant species changes from anionic to almost only the neutral ternary complex (Table 1). This is in line with speciation measurements done in the frame of an uranium diffusion experiment with Opalinus Clay (Joseph et al., 2013b). The change in the predominant species (Table 1) is associated with a farther migration of uranium of up to 30 m (Fig. 2a) due to a decrease in the sorption capacity (Fig. 2b). This can only be explained by the speciation, in particular the stronger complexation of uranium with calcium and carbonate (Table 1) with the application of the second NEA update (Grenthe et al., 2020), as the same D_e was used for all simulations. Anion exclusion effects are not taken into account as the composition of the diffuse double layer is only used in the model concept to counterbalance the surface charge. With the first NEA update (Guillaumont et al., 2003), about 94 % of uranium is present as ternary complexes as can be read from the percentages given in Table 1. However, this percentage proportion increases to almost 100 % with the second update (Grenthe et al., 2020). This, in turn, means that uranium increasingly forms ternary complexes, and thus less uranyl-ions are in total available to be retarded by sorption processes on the clay minerals (Table 1). The sorption capacity or K_d value, respectively, decreases and migration lengths increase (Fig. 2), because only uranyl ions are considered for sorption in the database. Consequently, the uranyl concentration in the modelled system has a huge impact on the sorption effects and ultimately on the calculated migration lengths.

Maximum migration distance of uranium is 80 m after one million years. The retention capacity of the Opalinus Clay seems to be sufficient to retain uranium, since concentrations close to the aquifers do not exceed the natural background concentration (Fig. 2a). However, the total thickness of the CPRZ is >200 m in the investigated scenario. In Germany, one of the minimum requirements defined by law (§ 23 Article 5 StandAG¹) is a thickness of at least 100 m for the CPRZ, regardless of the host rock type: argillaceous formation, crystalline or salt rock. This does not imply that such a thickness is sufficient in any case. Assuming that a potential repository is constructed in the centre, this would mean for the investigated case that uranium reaches adjacent aquifers as the retention capacity of the CPRZ is not sufficient. In this case, the site would not be suitable, because the safety requirements are not fulfilled. Should a thickness of at least 200 m therefore be favoured to fulfill safety requirements? The underground rock laboratory Mont Terri solely serves for the investigation of the Opalinus Clay as host formation. The site was and will never be considered as a potential disposal site due to its specific geological structure. For instance, the geochemical gradients at the potential disposal sites in Switzerland are less pronounced compared to Mont Terri (Mazurek et al., 2009; Nagra, 2014; Wersin et al., 2016, 2018). Sorption would decrease less towards the aquifers leading to reduced uranium migration. Hence, the impact of the hydrogeological system at the potential disposal sites on the migration lengths is diminished compared to Mont Terri. Furthermore, our simulations have to be considered as maximum scenarios as the impact of the engineered barriers is not taken into account (Hennig and Kühn, 2021b). They are supposed to minimize the source term concentration. Accordingly, the concentration gradient in Fick's laws is smaller and so diffusive transport. Site-specific analyses of the hydrogeological system as well as a more realistic integration of the source term concentration is required in order to calculate migration lengths of uranium at a potential disposal site. Nevertheless, our work shows how sensitive migration lengths resulting from reactive transport simulations are to the model conceptualization as well as selection of underlying thermodynamic data. Consequently, compilation and further development of comprehensive data sets is required and a site specific analyses of the hydrogeological system indispensable.

5 Conclusions

The impact of thermodynamic data on migration lengths resulting from reactive transport simulations is quantified for the example of uranium in the hydrogeological system of the Opalinus Clay at Mont Terri. In the geochemical system, speciation is controlled by the calcite-carbonate-ion system. In the pore water, uranium is mainly present as U(VI) as ternary complexes with calcium or magnesium together with carbonate. Previous simulations using the first NEA update of thermodynamic data for uranium indicated that the anionic complex ${\mathrm{CaUO}}_{\mathrm{2}}({\mathrm{CO}}_{\mathrm{3}}{)}_{\mathrm{3}}^{\mathrm{2}-}$ is the predominant species with a maximum distance of 50 m after one million years. The NEA published an update of the thermodynamic data for uranium, what, in turn, changes the predominant species from anionic to almost only the neutral ternary complex Ca₂UO₂(CO₃)₃. With identical simulations, except for the application of the second NEA update, a maximum distance of 80 m was obtained. Consequently, the impact of the change in the underlying thermodnamic data can be quantified with +30 m.

A stronger complexation of uranium with calcium and carbonate decreases sorption capacity. Sorption is quantified for the uranyl-hydroxo complexes. Due to the stronger complexation of uranium with calcium and carbonate with the application of the second NEA update, less uranyl-ions are in total available for sorption. This, in turn, decreases the sorption capacity, as shown by calculated K_d values, and leads to a farther migration through the formation.

The presented simulations show, how sensitive migration lengths of uranium resulting from reactive transport simulations are to the model concept, underlying thermodynamic data and applied approaches. Consequently, a site specific analysis of the hydrogeological system at potential disposal sites in combination with a comprehensive and up-to-date database are required to conduct reliable and powerful process-based simulations to quantify radionuclide migration lengths.