Introduction
===========

Radioactive waste disposal systems typically comprise a series of barriers that aim at protecting the environment and human health. The presence of several barriers enhances confidence that the waste will be adequately contained. In deep geological disposal systems, the barriers include the natural geological barrier and the engineered barrier system (EBS).

The efficiency of the EBS and the geological barrier is often tested through complex numerical models. These models serve to evaluate the response of the barrier under different sensitivity scenarios. One important factor to be evaluated is the groundwater flow crossing the barriers under different scenarios, which often consist in changes in the hydraulic properties of the different materials or changes in the waste compartment location. 

The present article shows an example of a COMSOL App set up to allow such sensitivity analyses. The model behind the app contains a hydrogeological groundwater flow model of a specific vault (1BMA) for intermediate level radioactive waste disposal operated by the Swedish nuclear waste management company (SKB).

Model Set-up
===========

The Multiphysics model underlying the application simulates groundwater flow through a rock vault in the Swedish final repository for short-lived radioactive waste [1]. The vault geometry is shown in Figure 1.

<figure className="image">
 <StrapiImage src="/uploads/BMA_app_1_d971a095ec.png" alt=""/>
 <figcaption>Figure 1. On the Geometry page the user can visualize the different structures and material domains available in the application.</figcaption>
</figure>

The Darcy’s law and Solute Transport interfaces have been applied to calculate groundwater flow rates and residence times in the vault. The residence time is defined here as the spatial distribution of the mean groundwater age, taking into account advection, diffusion and dispersion processes. It is calculated by solving a partial differential equation proposed by Goode et al [2].

$$
\frac{\partial A}{\partial{t}} = 1 + \nabla \cdot \left ( \mathbf{D} \cdot \nabla A \right ) - \nabla \cdot \left ( \frac{\mathbf{q}}{\theta} A \right ) + \frac{F}{\theta \rho}
$$

Above, $A$ is the water age ($s$), $\mathbf{q}$ is the Darcy velocity field (m/s), $\mathbf{D}$ is the hydrodynamic dispersion tensor (m2/s), $θ$ is the porosity, $\rho$ is the water density (kg/m3) and $F$ is a source term. 

Application description
===========

Figure 2 shows the ribbon at the top of the application graphical user interface. Clicking the buttons, from left to right, provides a workflow that guides the user through the application.


<figure className="image">
 <StrapiImage src="/uploads/BMA_app_2_19c85c0cd3.png" alt=""/>
 <figcaption>Figure 2. The application ribbon provides the workflow to guide the user.</figcaption>
</figure>

The user can change model inputs on the Setup page. This involves assigning values to the porosity, hydraulic conductivity and effective diffusivity of different materials and engineered structures.

The user can also select among predefined pressure fields to prescribe the boundary conditions for the groundwater flow (Figure 3). As an option, pressure boundary conditions can be imported from file.

<figure className="image">
 <StrapiImage src="/uploads/BMA_app_3_7394511a60.png" alt=""/>
 <figcaption>Figure 3. On the Setup page the user can modify the material transport properties and the boundary conditions (1). Selecting from the list of shoreline positions (2) updates the information graphics and text (3) as well as the displayed pressure field (4).</figcaption>
</figure>

After model setup the user can compute the groundwater flow and the residence time in the rock vault. Results can be analyzed through a set of predefined plots in 2D and 3D (Figure 4). Numerical values can be output in table exports. Formats are standardized such that quantitative simulation results can be used directly in connected modelling activities.

<figure className="image">
 <StrapiImage src="/uploads/BMA_app_4_8a6f9131a5.png" alt=""/>
 <figcaption>Figure 4. On the Plots page the user can select between flow or residence time results (1) displayed in 2D and 3D. Flow results presented in 2D use predefined cutplanes (2). Selecting from the list updates the information graphics and text (3 and 4) as well as the displayed results (5).</figcaption>
</figure>

After simulation, modelling reports can be generated automatically, summarizing the model setup and results.

Perspective
===========
The COMSOL Application Builder has been used to rework a coupled Multiphysics model in 3D into a ready-to-use simulation application. The application provides the opportunity for non-modelling experts and engineers to perform extended simulation trials and production runs. The accessibility of the application format furthermore facilitates review and quality control work associated with the model development process itself.

References
===========

*COMSOL, 2015*. **COMSOL Multiphysics 5.1.** COMSOL AB, 2015.

*Goode D J,* 1996. **Direct simulation of groundwater age.** Water resources research, vol 32. N.12, pages 289-296.

medium_BMA_app_1.png

small_BMA_app_1.png

thumbnail_BMA_app_1.png

BMA_app_1.png

Radioactive waste disposal systems typically comprise a series of barriers that aim at protecting the environment and human health. 

Comsol App for the evaluation of an engineered barrier system

Introduction
=============

Damage mechanics has often been used to study the mechanical behaviour of concrete structures. Damage models can be used to simulate the quasi-brittle behaviour of cementitious materials, since they are capable of representing cracking under various loading conditions. In this post, we present the implementation, regularization and validation of the classical isotropic Mazars’ damage model for concrete in Comsol Multiphysics. We furthermore explore its limitations and present the implementation of an alternative damage model, also proposed by Mazars and co-workers, more suitable for multiaxial loading conditions.

Background
=============

Simulating concrete load bearing capacity and post-peak behaviour is essential to assess the mechanical behaviour of concrete subject to mechanical and environmental loads. The tool presented in this post allows accurate mechanical behaviour characterization through a simple isotropic model. The classical Mazars’ damage model (Mazars, 1986) has been implemented within the Comsol user interface. The regularization of this method to avoid mesh-dependency of the results is carried out by the implicit gradient damage formulation (Peerlings et al., 1996). In addition, we have also implemented a recently developed damage model with updated features (Mazars et al., 2015) to more accurately describe concrete behaviour under biaxial and triaxial loading conditions. The final goal of this work is to develop concrete mechanics models in the Comsol interface that allow their coupling with other processes, such as chemical degradation for durability issues, or moisture and heat transport.

Mazars’ damage model implementation in the Comsol interface
================================================

Classical models following damage mechanics theory (Kachanov, 1958) are based on the definition of a scalar damage variable *d*, varying between 0 for an intact material and 1 for a fully damaged state. This variable affects directly the stiffness tensor, leading to a nonlinear relation between strain (**ε**) and stress (**σ**, MPa) tensors.

One of the most widespread damage models, perhaps due to its simplicity, is the one proposed by Mazars (1986). This model defines an isotropic scalar damage variable described by different evolution laws whether the material is under tensile or compressive conditions. The damage variable *d*, is calculated as a linear combination of compressive and tensile damage contributions based on the evolution of a history variable that stores the maximum value of strains reached. The main model material parameters are At, Bt,Ac and Bc.

In this line, Comsol developed a functionality to access external material models for structural mechanics in a post by Ed Gonzalez (2015). This functionality is intended to enable the programming of any desired constitutive model as an external material model to be used in Comsol applications. In that post, an example is presented with the implementation of Mazars’ damage model. Although this functionality is quite powerful, one of the main drawbacks is the complexity to access the variables defined in the external material code, which complicates the coupling with other physics.

In this post, we present a different implementation of the Mazars’s damage model, in this case using a generic physics (Domain Ordinary Differential Equation) and a specific solver configuration to allow for storing the history variables. The main feature of the implementation is the definition of two history variables needed in the mechanical damage model. Two additional degrees of freedom (DOFs) are added through this new physics.

A compression test with loading/unloading cycles has been modelled with the two above-mentioned approaches. Figure 1 presents the cyclic displacement imposed at the top face of the concrete cylinder and the evolution of the damage variable.

<figure className="image">
<StrapiVideo width="725" height="395" src="/uploads/Mech_Dam_Cyclic_Loading_cce1792b0a.mp4" type="video/mp4"/>
 <figcaption>Figure 1. Imposed displacement in the top face and resulting damage variable evolution with time for the compression test with loading/unloading cycles.</figcaption>
</figure>

Figure 2 shows a comparison of the two implementation methods that confirms that both methods yield the same results. The main advantages of the new implementation are that: (1) it can be fully coupled with other physics and any other constitutive law, (2) all variables are readily available within the Comsol interface, thus facilitating the post-process of results, and (3) the mechanical damage model is easier to adjust or reformulate in the physics tree, without the need of code compilation.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_2_0aec1ae1a0.png" alt=""/>
 <figcaption> Figure 2. Results of the verification test in terms of stress-strain curves using two different damage model implementations.</figcaption>
</figure>

Regularization by the implicit gradient damage formulation
============================================

The use of damage models for simulation of nonlinear stress-strain response including the softening regime has the disadvantage that the results depend on the size of the finite elements, i.e. they are mesh dependent. To remedy this, several methods have been proposed in the literature to regularize the solution so that the results converge to a single solution upon mesh refinement. One of such methods is the so-called gradient-enhanced formulation, which we have implemented in Comsol following the implicit gradient formulation proposed by Peerlings et al. (1996) and Simone (2007). This model consists of a Helmholtz type differential equation in terms of an equivalent strain, . The formulation includes the use of a characteristic length l (m) to compute the nonlocal variable . Thus, in the regularized model, the history variable κ is calculated with this nonlocal variable instead of the equivalent strain. The regularization of Mazars’s damage model has been implemented in Comsol using the built-in Helmholtz equation physics.

To validate this regularized implementation, a classical 3-point bending beam test has been simulated. This kind of test has been widely studied experimentally and computationally by several authors (Kormeling and Reinhardt, 1983; Jirásek, 2011).

Results are presented with and without regularization in Figure 3. Without regularization, the damage model results in different peak loads depending on the finite element size. The model convergence is also sensitive to the element size, resulting in stability problems in the post-peak regions for finer meshes (similar to Jirásek, 2011). On the other hand, with the implicit gradient formulation, the results are much more objective and little mesh dependence is observed. In this case, the characteristic length is 1 mm. The results correspond very well with the regularized test reported by Jirásek (2011) when comparing load-displacement curves. Therefore, it is concluded that the regularized mechanical damage model implementation works properly.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_3_1f0a09af90.png" alt=""/>
 <figcaption> Figure 3. Load (N) vs. displacement (mm) curve of the unnotched three-point beam test with no regularization method (left), and with the regularized damage model (right). Results are shown for different mesh refinements.</figcaption>
</figure>

Mu model
========

Based on Mazars’ classical damage model (Mazars, 1986) a new formulation named “μ model” (Mazars et al., 2015) has been recently proposed. The improvements of this formulation include among others a more accurate representation of concrete behaviour against biaxial and triaxial compression. This is a useful development when modelling materials under multi-axial loading regimes, as is the case of underground structures. In these cases, high compression stresses are normally acting in the three principal directions of the stress tensor.

It is known that the Mazars’ classical damage model under plane stress underestimates the strength of concrete under biaxial compression (e.g. Jirásek, 2011). We have implemented the μ model in Comsol following a similar approach than for the Mazars’ damage model.

The results of the Mazars’ damage model implemented in Comsol under plane stress conditions for different biaxial loads are compared to classical experimental data from Kupfer et al. (1969) in Figure 4. It is clearly observed that under biaxial compression and plane stress conditions the model predicts a significantly lower concrete strength compared to experimental data. On the other hand, strength increase under biaxial compression is obtained with the μ model, as shown in Figure 4 from the current implementation in Comsol.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_4_9a77569630.png" alt=""/>
 <figcaption> Figure 4. Comparison of the results of several biaxial loading tests from Kupfer et al. (1969) and the model results using the Comsol implementation of (1) the Mazars’ damage model and (2) the μ model.</figcaption>
</figure>

To quantify the improvements of the μ model under confined conditions, a set of simulations representing triaxial compression tests were implemented. The triaxial test simulated here were carried out by Ramtani (1990) and were extracted from Mazars et al. (2015). The experiments are carried out at several confinement levels (S= 0, 20, 35, 50 and 100 MPa). The results of the μ model obtained by Mazars et al. (2015) and those extracted from the Comsol implementation are presented in Figure 5. The model implementation is again validated, as results from both models are perfectly matching. Although the implemented scheme gives accurate results in terms of strength, some convergence problems were found after the peak stress with increasing strains.

<figure className="image">
 <StrapiImage src="/uploads/comsol_mechanical_damage_5_b4e74e3f30.png" alt=""/>
 <figcaption> Figure 5. Comparison of triaxial test modelled results from Mazars et al. (2015) in blue and the implemented μ model results in Comsol Multiphysics in red.</figcaption>
</figure>

Conclusions
==========

In this post, we presented an alternative implementation of damage models for concrete mechanical behaviour to the external materials functionality of Comsol. We used a generic physics (Domain Ordinary Differential Equation) and a specific solver configuration to allow for storing the history variables of the constitutive model. The presented models can represent concrete mechanical behaviour under diverse loading paths. Models are implemented in Comsol taking special care to history variables description. Classical mesh-dependency problems are overcome by the regularization of the damage variable through a gradient-enhanced scheme. The models were validated by comparison with results from different experimental tests. The implemented models are able to predict concrete mechanical behaviour under confined conditions and can be readily couple with other mechanical processes, as well as chemical degradation or hydraulic and thermal processes.


Acknowledgements
==============

This work has received funding from the European Union's Horizon 2020 Research and Training Programme of the European Atomic Energy Community (H2020-NFRP-2014/2015) under grant agreement n° 662147 (CEBAMA) and from the Swedish Nuclear Fuel and Waste Management Company (SKB), which are gratefully acknowledged. Fruitful discussions with Henrik von Schenck and Per Mårtensson from SKB are also gratefully acknowledged.

References
==============

**González E, 2015.** Accessing External Material Models for Structural Mechanics. COMSOL Blog December 2015, accessible at: https://www.comsol.com/blogs/accessing-external-material-models-for-structural-mechanics/

**Jirásek M, 2011.** Damage and smeared crack models. In Numerical Modeling of Concrete Cracking. Editors: Günter Hofstetter, Günther Meschke, p. 1-49.

**Kachanov L M, 1958.** Time of the rupture process under creep conditions. Isv. Akad. Nauk. SSR, 8, 26–31.

**Kormeling H A, Reinhardt H W, 1983.** Determination of the fracture energy of normal concrete and epoxy modified concrete. Stevin Laboratory, Delft University of Technology, Report No. 5-83-18, 1983.

**Kupfer H, Hilsdorf H K, Rusch H, 1969.** Behavior of concrete under biaxial stresses. Journal of the American Concrete Institute, 66, 656–666.

**Mazars J, 1986.** A description of micro- and macroscale damage of concrete structures. Engineering Fracture Mechanics, 25(5–6), 729-737.

**Mazars J, Hamon F, Grange S, 2015.** A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings. Materials and Structures, 48, 3779–3793.

**Peerlings R H J, de Borst R, Brekelmans W A M, de Vree J H P, 1996.** Gradient-enhanced damage for quasi-brittle materials. Int. J. for Numerical Methods in Engng., 39, 3391-3403.

**Ramtani S, 1990.**Contribution à la modélisation du comportement multiaxial du béton endommagé avec description du caractère unilateral. PhD thesis. University Paris 6, 1990.

**Simone A, 2007. **Explicit and implicit gradient-enhanced damage models. Revue Européenne de Génie Civil, 11(7-8), 1023-1044.

small_comsol_mechanical_damage_2.png

thumbnail_comsol_mechanical_damage_2.png

comsol_mechanical_damage_2.png

Damage mechanics has often been used to study the mechanical behaviour of concrete structures

Comsol implementation of mechanical damage models for concrete

The safety of Deep Geological Repositories (DGR) for spent nuclear fuel is a highly important issue, and should be assessed using the Best Available Technology (BAT). With this motivation in mind, Amphos 21 has recently developed a numerical tool, denoted as iDP (interface between DarcyTools and PFLOTRAN; Molinero et al., 2016), which provides aJuqueen unique High Performance Computing(HPC)-based framework for the simulation of multicomponent reactive transport problems in fractured crystalline media. iDP combines two powerful state-of-the-art software: DarcyTools (Svensson and Ferry, 2014), which is used to compute groundwater flow in sparsely fractured rocks using a stochastic continuum approach, and the massively parallel reactive transport code PFLOTRAN (Hammond and Lichtner, 2010). iDP has been tested during the JUQUEEN Extreme Scaling Workshop 2016 (Brömmel et al., 2016) showing excellent scalability up to 131,072 processes (8 racks).

<figure className="image">
 <StrapiImage src="/uploads/HPC_Blog_Post_1_c4a1155a9d.png" alt=""/>
 <figcaption> </figcaption>
</figure>

<div className="grid grid-cols-5 gap-4">
 <div className="col-span-2">
 <StrapiImage className="" src="/uploads/HPC_Blog_Post_2_9212e50e5d.jpg"/>
</div>

 <div className="col-span-3">
In cooperation with the Forschungszentrum Jülich and the Swedish Nuclear Fuel and Waste Management Company (SKB), Amphos 21 has used iDP to carry out two challenging reactive transport simulations. The first set of calculations has focused on evaluating the potential ingress of oxygenated glacial-meltwater into initially anoxic fractured media. This is an issue of concern for the safety of deep geological Pocos de Caldasrepositories for spent nuclear fuel as oxidizing conditions would increase the solubility and mobility of many radionuclides and O2-rich waters would corrode canister materials containing the nuclear waste.
</div>
</div>

<div className="grid grid-cols-7 gap-4">

 <div className="col-span-5">
The reactive transport models have pointed out the strong buffering capacity of a typical crystalline medium (e.g. Forsmark in Sweden), with the infiltration of oxygen being limited to few highly transmissive deformation zones. These findings are in agreement with data from the natural analogue site of Poços de Caldas (Brazil) (Nordstrom et al., 1992). For more details about these reactive transport simulations, the reader is referred to a paper recently published in Journal of Contaminant Hydrology (Trinchero et al, 2017a).
</div>
 <div className="col-span-2">
 <StrapiImage className="" src="/uploads/HPC_Blog_Post_4_7e2c2e6e88.jpg"/>
</div>
</div>



Another important issue for the safety of a DGR is related to the potential of the host rock to retain those radionuclides that are released due to a possible failure of the engineered barriers.Total CS This potential, in turn, depends on two coupled processes: diffusion into the rock matrix and sorption onto the mineral surfaces of the available reactive minerals. A typical modelling assumption is to use an equivalent (homogeneous) description of sorption sites. This assumption is here assessed using a very detailed micro-scale model of the rock matrix, which is built upon a numerical grid consisting of more than 50 million elements. The results of the calculations, in which the distribution of sorption sites is described at the grain scale, show that mineralogical heterogeneity has a significant impact on the retention of radionuclides, namely cesium. Not considering that sorption sites are sparse and unevenly distributed results in an important overestimation of the rock retention capacity. Further details can be found in a paper recently published in Transport in Porous Media (Trinchero et al., 2017b).

<figure className="image">
 <StrapiImage src="/uploads/HPC_Blog_Post_5_7a7da1df76.png" alt=""/>
 <figcaption> </figcaption>
</figure>


Acknowledgement
==============
This work has been funded by the Swedish Nuclear Fuel and Waste Management Company (SKB). The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC Partition part of the supercomputer JUQUEEN at Forschungszentrum Jülich. The authors also thank Ignasi Puigdomenech, Björn Gylling, Guido Deissmann, Urban Svensson, Glenn Hammond, Peter Lichtner and the PFLOTRAN development group for their help and valuable suggestions during the project.

References
=========
*Brömmel, D., Wylie, B. J., Frings, W.*, 2016.** JUQUEEN extreme scaling workshop 2016. **Tech. Rep. FZJ-JSC-IB-2016-01, Jülich Supercomputing Centre.

*Hammond, G. E., Lichtner, P. C.*, 2010. **Field-scale model for the natural attenuation of uranium at the Hanford 300 area using high-performance computing. **Water Resources Research 46 (9), W09527, 131.

*Molinero, J., Trinchero, P., Ebrahimi, H., de Vries, L., Luna, M., Svensson, U., Lichtner, P. *2016. **The BRIDGE project: development, testing and application of a High Performance Computing framework for reactive transport modelling in crystalline roks (iDP)**. R-15-17, Svensk Kärnbränslehantering AB, Stockholm, Sweden.

*Nordstrom, D., McNutt, R., Puigdomènech, I., Smellie, J. A., Wolf, M. *1992. **Ground water chemistry and geochemical modeling of water-rock interactions at the Osamu Utsumi mine and the Morro do Ferro analogue study sites,** Poços de Caldas, Minas Gerais, Brazil. Journal of Geochemical Exploration 45 (1-3), 249-287

*Svensson, U., Ferry, M., *2014. **Darcytools: A computer code for hydrogeological analysis of nuclear waste repositories in fractured rock. **Journal of Applied Mathematics and Physics 2 (06), 365.

*Trinchero, P., Molinero, J., Deissmann, G., Svensson, U., Gylling, B., Ebrahimi, H., Hammond, G., Bosbach, D., Puigdomenech, I.*, 2017b.** Implications of grain-scale mineralogical heterogeneity for radionuclide transport in fractured media.** Transport in Porous Media 116: 73. doi:10.1007/s11242-016-0765-0

*Trinchero, P., Puigdomenech, I., Molinero, J., Ebrahimi, H., Gylling, B., Svensson, U., Bosbach, D., Deissmann, G.*, 2017a. **Continuum-based DFN-consistent numerical framework for the simulation of oxygen infiltration into fractured crystalline rocks. **Journal of Contaminant Hydrology DOI: 10.1016/j.jconhyd.2017.04.001

large_HPC_BlogPost_1.png

medium_HPC_BlogPost_1.png

small_HPC_BlogPost_1.png

thumbnail_HPC_BlogPost_1.png

HPC_BlogPost_1.png

The safety of Deep Geological Repositories (DGR) for spent nuclear fuel is a highly important issue, and should be assessed using the Best Available Technology (BAT)

High Performance Computing as a prognostic tool for Deep Geological Repositories for spent nuclear fuel

The problem
================
Disposal in a deep underground repository is currently considered by several countries as the preferred option for long-term management of spent nuclear fuel. For regulatory reasons the safety of such a disposal facility must be demonstrated quantitatively through Performance Assessment calculations. One of the factors that have the potential to significantly influence the results of Performance Assessment is the radionuclide release rate from the waste immediately upon contact with groundwater (the so called **Instant Release Fraction**). Performance Assessment calculations indicate that instant release is important for the overall safety assessment of the facility, but there are significant uncertainties about the nature of this phenomenon and how it should be represented in Performance Assessment.

 With the above in mind, the **FIRST-Nuclides** Collaborative Project was funded by the European Commission under the 7th Framework Programme of the European Atomic Energy Community (EURATOM). The project has been initiated to provide better (conceptual and quantitative) understanding of the processes involved in instant release of radionuclides from spent nuclear fuel.

In the context of instant release, the knowledge of the rate of spent fuel saturation with water is key as the release of some radionuclides is contingent upon wetting of the fuel surfaces by water. Therefore, the instant release rate of these radionuclides depends on the rate of spent fuel saturation. An improved understanding of the rate of spent fuel saturation with water will help elucidate radionuclide instant release both under repository conditions and with regard to laboratory experiments conducted under near-room conditions.

Approach
=======

Under repository conditions the exposure of spent fuel to groundwater would be conditional upon the failure of multiple interconnected engineered barriers. The processes involved in such a phenomenon are multiple and coupled, making the system highly complex.

As a point of departure therefore, the saturation of a single idealised spent fuel pellet was considered (henceforth referred to as the *Reference Pellet*). This situation can be thought of as representative of batch laboratory leaching experiments performed on complete or fragmented fuel pellets, which are performed in order to measure radionuclide release rates from spent fuel upon contact with water. The results of modelling of saturation with water of a single spent fuel pellet are directly applicable to the interpretation of laboratory leaching results, but are also the necessary first step towards more complex multiphysics models capable of representing the fate of spent fuel upon wetting under repository conditions.


The conceptual model assumed that the pellet was composed of an impervious to water UO2 matrix and a network of (thermal) cracks, which can be filled with water (Figure 1). By saturation with water of the pellet is therefore understood the process of filling up of these cracks with water.

 

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modelling_saturation_figure_1_9d08583700.png" alt=""/>
 <figcaption>Figure 1: The conceptualised Reference Pellet (left – lines indicate intersection of cracks with the pellet’s surface) and the finite element mesh (right) of the idealised crack system. Note that water flow is considered to occur through the cracks only (the fuel matrix is impervious to water). Dimensions shown are in metres.</figcaption>
</figure>

Implementation in Comsol
======================
Water penetration into the interconnected system of cracks was simulated by solving the Richard’s equation (in the pressure-based form). This was implemented in Comsol utilising the Fracture Flow interface (calculation of flow in discrete fractures), which was modified by manually implementing a storage-retention term and a set of constitutive equations (the water retention model) so as to account for variably-saturated flow within the cracks. More details on the conceptual and numerical details of the model can be found in [1].

The primary variable solved for is water pressure, which, through the water retention model, links directly with water saturation. In order to represent the physics of water penetration a negative water pressure (positive capillary suction pressure) was prescribed homogenously within the cracks, while a positive pressure (here corresponding to 5 cm of water column) was applied to all external crack edges. The driving force for water penetration into the cracks is water pressure (capillary suction) gradients (Figure 2), which disappear after the saturation process is complete.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modelling_saturation_figure_2_070a9e6cf2.png" alt=""/>
 <figcaption>Figure 2: Successive stages of saturation with water (from left to right) of an interconnected system of cracks in a single spent fuel pellet. The colour scale represents effective saturation (-), which varies from the initial value of 0.1 ‰ (deep blue) to unity (red) at full saturation.</figcaption>
</figure>


Insights Gained
=============
The model enables quantification of the time needed to completely saturate a given system of cracks, which may vary depending on the statistical properties of the cracks (e.g. their number, geometry, interconnectivity and aperture), and depending on the initial and boundary conditions used (e.g. representing the presence of partial cladding or depth at which water penetration occurs). This is useful information as it makes it possible to compare the crack system surface wetting rate with the rate of radionuclide release observed in laboratory leaching experiments. Preliminary results of this on-going work indicate that, indeed, the time scales for cracks saturation with water are comparable with those operationally defined for instant/fast release of some radionuclides (hours to weeks).

We are currently working on an implementation of tracer release from the crack surface and transport by diffusion through the crack system towards an external reservoir (experimental batch solution or rock formation water in case of repository conditions). Also, the saturation model is being extended to include two distinct, but interconnected, populations of cracks (of different mean apertures). Ultimately, we intend to include relevant chemical reactions (e.g. dissolution of the fuel matrix at later stages) and upscale the model to the rod assembly scale. Work progress and results obtained in this work will continue to be presented and discussed at this blog.

Acknowledgements
================
The research leading to these results has received funding from the European Union's European Atomic Energy Community's (Euratom) Seventh Framework Programme FP7/2007-2011 under grant agreement n° 295722 (FIRST-Nuclides project). 

This work was carried out with important contributions from Amphos 21 colleagues: Andrés Idiart, Lara Duro and Olga Riba.

References
================
[1] *Pękala et al.* 2012. **Modelling of spent fuel saturation with water – approach, preliminajry results and potential implications. **In: KIT Scientific Reports 7639, 1st Annual Workshop Proceedings of the Collaborative Project “Fast/Instant Release of Safety Relevant Radionuclides from Spent Nuclear Fuel (7th EC FP CP FIRST-Nuclides).

medium_modelling_saturation_figure_2.png

small_modelling_saturation_figure_2.png

thumbnail_modelling_saturation_figure_2.png

modelling_saturation_figure_2.png

In this post, Marek Pękala describes a recent application of Comsol at Amphos 21 to study the saturation with water of spent nuclear fuel.

Modelling saturation with water of spent nuclear fuel in the context of its geological disposal

Scope and objectives
==================
In the framework of safety assessment studies for depth geological disposal of nuclear waste, important scientific efforts are needed in order to understand the feasibility of the planned repositories. This is precisely the case of Olkiluoto in Finland, where the land uplift, related to the withdrawal of ice from the last glacial Age, plays a pivotal role in the hydrogeological patterns of the site. It is thus very important to improve the knowledge about the past evolution of the hydrogeological and geochemical conditions. In fact, as a result of land uplift, the infiltration of meteoric water could have triggered weathering processes, which, in turn, have impact on key parameters controlling the long-term stability of the planned repository (pH and redox conditions, availability of mineral buffers, etc). Furthermore, the presence of sulphate-rich waters along with bacterial communities and substrate might activate microbial mediated sulphate reduction processes with the related generation of, which is a detrimental agent for the copper canisters where the waste is placed.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/coupled_palaehydrogeological_microbial_figure_1_4d26cbef8c.png" alt=""/>
 <figcaption>Geographical location of Olkiluoto Island and related fractured system object of simulations. It is shown the initial density field in depth.</figcaption>
</figure>

Numerical model
==============
In order to improve the understanding of the effect of the aforementioned events on the current hydrochemical conditions, we have developed a complex three-dimensional palaehydrogeological reactive transport model of the Olkiluoto site. The model consist of four hydrogeological zones of the Well Characterized Area (WCA), namely BFZ100, HZ19A, HZ19B and HZ19C. The Palaeohydrogelogical Reactive Transport (PRT) simulations, which spans 2000 years (from 0AD to present), aim at capturing the mutual interplay of different processes (i.e. land uplift, intrusion of meteoric water, weathering, mixing of waters etc.). The hydrogeological zones are simulated as two-dimensional planes within a three-dimensional domain. 

 <figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/coupled_palaehydrogeological_microbial_figure_2_d60154d447.png" alt=""/>
 <figcaption>Detail of the hydrogeological zones and FE grid</figcaption>
</figure>

The numerical model has been implemented using the Interface COMSOL-PHREEQC (iCP), which has been developed and validated by Amphos 21. To increase the reliability of the PRT simulations, groundwater flow parameters (i.e. residual pressures and density fields) are taken from an existing palaeohydrogeological flow model (Löfman and Karvonen, 2012)

Some Results
===========
After some preliminary runs, different sets of three-dimensional PRT simulations have been carried out. In each of these simulations the geochemical complexity of the system has been progressively increased by including new processes (e.g. mineral equilibrium → mineral equilibrium and kinetics → mineral equilibrium, mineral kinetics and cation exchange). 

 
 <figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/coupled_palaehydrogeological_microbial_figure_3_0ff817dfc0.png" alt=""/>
 <figcaption>Detail of sulphide distribution in mol/L (x10-5) for 2000 AD</figcaption>
</figure>


 <figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/coupled_palaehydrogeological_microbial_figure_4_ba61628795.png" alt=""/>
 <figcaption>Detail of calcite distribution in mol/Kg for 2000 AD</figcaption>
</figure>

The results of the PRT simulations have been compared with the available field data and the general good agreement of the results has highlighted both the reliability of the model and its potential for the prediction of the future evolution of the geochemical conditions at Olkiluoto.

 <figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/coupled_palaehydrogeological_microbial_figure_5_b154cae116.png" alt=""/>
 <figcaption>Sulphate distribution in depth for 2000 AD</figcaption>
</figure>

Acknowledgements
=================
This work has been funded by POSIVA. Amphos 21 wants to thank Petteri Pitkänen and Lasse Koskinen for their support and very helpful comments. Miguel Luna also wants to thank the coauthors of this work: Paolo Trinchero and Jorge Molinero.

medium_coupled_palaehydrogeological_microbial_figure_2.png

small_coupled_palaehydrogeological_microbial_figure_2.png

thumbnail_coupled_palaehydrogeological_microbial_figure_2.png

coupled_palaehydrogeological_microbial_figure_2.png

A large-scale, long-term simulation of the hydrogeochemical behaviour of fractured rocks is shown in this post. 

Coupled Palaehydrogeological Microbial and Geochemical Reactive Transport Model of the Olkiluoto Site (Finland)

The basics
=============

Helium originates mainly from three sources: the atmosphere, the mantle and as a result of radioactive decay processes taking place in the crust. Each source of helium is characterized by a particular 3He/4He ratio. The isotopic composition of helium in the atmosphere is characterized by a 3He/4He ratio of about 1.4x10-6. Radiogenic helium is primarily 4He, with a 3He/4He ratio in the terrestrial rocks closer to 10-8, and primordial helium in the deep earth is enriched in 3He with a variable isotopic composition (from 1.1x10-5 in the upper mantle to 8.7x10-5 in the lower one). Helium in the crust originates from the alpha decay of natural uranium and thorium series, whereas 3He is formed by thermal neutron capture by 6Li.

Helium production in the crust is dominated by the alpha decay of the uranium and thorium decay chains, and therefore directly proportional to the concentration of 235,238U and 232Th in the rocks.

The Helium production rate in the rock can be evaluated as:

$$
G=\frac{\rho}{V_m}\left( 1.19 \cdot 10 ^{-16} [U] +2.88 \cdot 10 ^{-17} [Th] \right)
$$

where r is the rock bulk density (kg/m3), Vm the molar volume of the gas (m3/mol) and [U] and [Th] are the uranium and thorium contents of the rock matrix (ppm), respectively. Then, if we know (measure) the concentration of dissolved helium in a groundwater sample, we can estimate how much time this water had (theoretically) been in contact with the rock to reach such Helium content. But this is assuming all the Helium comes from radiogenic origin without mantelic and/or atmospheric contribution. And also assuming a closed 0-dimensional system, i.e. no diffusive transport of Helium occurs in the rock porewater.

New analytical solution for interpretation of He data in deep groundwater
============================================================

In our paper we present the details of a new analytical solution for modeling Helium production and diffusion in a 1D finite profile. It is important to remark the finite nature of the model domain because this is one of the key innovations. So far the only available analytical solution, known as the “Andrews model” (Andrews 1985), was for an infinite thickness of the crust. The other key innovation with respect to “Andrew’s model” is the possibility of discretizing the crust in different domains with heterogeneous [U] and [Th] compositions (i.e. to include different rocks in the model). The analytical solution looks like this:

$$
c\left( z, t_r \right) = \frac{2d}{D_p} \sum^{\infty}_{k=1} \frac{ sin \left(\alpha_k \frac{z}{d} \right)}{\alpha_k^2} \left[ 1 - exp \left( \alpha_k^2 \frac{D_p}{d^2} t_r \right) \right]
\times \\
\left\{ \frac{G_1 d}{\alpha_k \phi} \left[ cos \left( \alpha_k \frac{d_l}{d} \right) - cos \left( \alpha_k \frac{d_u}{d} \right) \right] + \left( -1 \right)^k j_0 \right\}
$$

Where $D_p$ is the pore diffusion coefficient, tr is the residence time of groundwater, $G$ is the Helium production, f is the rock porosity, $J_0$ is the mass flux of Helium coming from the mantle (bottom boundary condition) and d is the total depth of the model (i.e. the deep boundary of the finite domain where $J_0$is applied as a boundary condition). In the example we have published, the discretization of the domain is done by 3 intervals like:

$$
 G \left( z \right) =
 \begin{cases}
 0, & \ 0>z>-d_u \\
 G_ 1, & -d_u \geq z>-d_l \\
 0, & -d_l \geq z >-d
 \end{cases}
$$

where $d_u$ is the depth of the “shallow” aquifer with rapid groundwater renewal (and we assume the production of he is negligible), and dl is the maximum depth of the rock unit containing [U] and [Th]. It is worth noting that we adopt this “discretization” to mimic the application example shown below, but it can be changed or extended to any other settings.

 

Applying the new methodology at Olkiluoto (Finland)
=======================================

Olkiluoto is the site where POSIVA (the finish company for nuclear spent fuel management) is conducting field investigations to evaluate the suitability of constructing an underground facility to dispose high-level radioactive waste. Olkiluoto is located on the coast of Gulf of Bothnia in south-western Finland. The two main lithological units of the Olkiluoto site are (1) high grade metamorphic rocks and (2) igneous rocks. Hydrogeological and geochemical investigations shows that the site has a “shallow” part (from surface to about 250 meters) with active hydrogeological circulation, and a deeper zone with groundwater movement is much lower and limited to some water conducting fracture zones.

Radiogenic helium production at Olkiluoto is related with the presence of granites in the bedrock. Recent seismic studies have suggested that in the Olkiluoto area, different granites extend to a depth of at least 4 km. Thus, we have considered two models, one where a constant production rate is applied from the upper boundary (which is located at 250 m depth, after the zone with more active hydrogeological circulation) down to 4000 m depth, and the second where in situ production is present only in the first 1000 m of bedrock. A deep source of helium ($J_0$), which mimics deep-seated granite production, is located at the bottom of the domain.

The Uranium and Thorium contents of the rock are well known, as well as diffusion coefficients and porosities. Then, we did a fitting exercise by calibrating the groundwater residence time, tr and the deep source of helium, $J_0$. Both models were able to fit very well the Helium concentrations measured in the field, as shown in the next figure (dots are field measurements and lines are the calibrated models).

 
<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/00_Groundwater_Age_Helium_2_0694b5b652.png" alt=""/>
 <figcaption>Figure 1. Comparison among measured Helium concentrations and the best fit calibration with the new analytical solution. The fitting parameters are the residence time of groundwater and the Helium flux at the bottom boundary.</figcaption>
</figure>

The resulting fitted parameters are the same for 2 models, which indicates that extending the Helium production deeper than 1 kilometer is not too sensitive. The best fit was obtained for a tr (residence time) of 3.2x107y, and $J_0$(Helium flux from the mantle) of 3.1x10-9 mol/m2y.

So far we do not have data about Helium isotopic ratios in the groundwater. Future determinations of such isotopes could help us to quantify the mantelic contribution and then test the validity (or at least the degree of confidence) about our estimation of 30 million years of residence time for the groundwater in the granitic rocks of Olkiluoto.

ACKNOWLEDGEMENTS
=====================

I would like to acknowledge to all my coauthors of this paper. It was a great project and I think we all enjoy very much working together. Paolo Trinchero was the leader of the project at that time, together with Anne Delos who worked very hard in the first stages. Marco Dentz was the one really behind the maths for the analytical solution and Petteri Pitkanen provided the best expertise for Olkiluoto geological settings and the behavior of dissolved gasses in the groundwater. Ignasi Puigdomenech (SKB) was also a key person for this project, and provided high level comments and constructive suggestions to the first versions of the report and paper.

References
=========

*Trinchero, P.; Delos, A.; Molinero, J.; Dentz, M. & Pitkanen, P. *(2014). **Understanding and modelling dissolved gas transport in the bedrockof three Fennoscandian sites.** Journal of Hydrology, 512, 506–517.

*Andrews, J. 1*(985). **The isotopic composition of radiogenic helium and its use to study groundwater movement in confined aquifers.** Chemical Geology 49, 339–351.

medium_00_Groundwater_Age_Helium_1.png

small_00_Groundwater_Age_Helium_1.png

thumbnail_00_Groundwater_Age_Helium_1.png

00_Groundwater_Age_Helium_1.png

Dating groundwater is something that hydrogeologists like to do for different practical purposes

How old is groundwater? Estimation of residence time with Helium

Background
=========
Groundwater inflows into tunnels may constitute a potential hazard as well as an important factor controlling the rate of advancement in the excavation. In fact, one of the major practical difficulty during tunnel construction is often associated with the presence of groundwater. Predictions of groundwater inflows are also needed to design suitable drainage systems. Furthermore, drawdowns produced by tunnel construction may induce land subsidence, water table decline and environmental impacts on rivers and wetlands.

In the literature, there are several analytical expressions to calculate water discharges into tunnels but a major development was presented recently by Perrochet and Dematteis (2007). The program TunFlow makes use of the above-mentioned analytical expression to compute groundwater infiltration into tunnels and galleries. The program allows an unlimited number of geological zones (heterogeneous system) to be accounted for. Furthermore, fully transient hydrogeological conditions (accounting for the velocity of tunnel advancement) are considered. TunFlow is easy to use and is especially suited for the evaluation of groundwater filtrations in scoping studies for tunnels construction.

 

Input parameters
=============
The analytical expression developed by Perrochet and DeMatteis (2007) to calculate the groundwater tunnel inflows and used in TunFlow is shown below:

$$
Q(t) = 2 \pi \sum^{N}_{i=1} H(t-t_i) \int^{v_i\left( t-t_i \right)}_{0}
\frac{K_i s_i H \left( L_i - x \right)}
{ln \left( 1 + \sqrt{ \frac{\pi K_i}{S_i r_i^2} \left( t-t_i - \frac{x}{v_i} \right) } \right) } dx
$$

The parameters of this formula are:

- $t$ [d]: time
- $L$ [m]: length of the section of a given geological layer
- $v$ [m/d]: drilling velocity
- $K$ [m/d]: permeability of the considered layer
- $S$ [1/m] Specific storage of the considered layer
- $r$ [m]: tunnel radius
- $s$ [m] level of water table (distance from the tunnel roof)
 

Before running the model, it is preferable to create an Excel sheet with all the information requested by the code. The numbers of layers is specified in the Graphical User Interface (GUI) of TunFlow (1 and 2, **Figure 1**).The number of time steps must be also selected (3, **Figure1**) . By default, 1000 time steps are used. It is worthwhile noting that increasing the number of time steps improves the accuracy of the solution but also increases the computational costs of the simulations (i.e. simulation time).

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/tunflow_1_03_72cf2245f5.png" alt=""/>
 <figcaption>Figure 1. Screenshot of the Graphical User Interface (GUI) of TunFlow.</figcaption>
</figure>

Once all the parameters are imported in TunFlow, the program can be executed (4, **Figure 1**).

 

Results and output options 
=============
 

Tunflow provides diferent outputs. Calculated data can be exported in txt and xlsx file formats (1 and 2, **Figure 2**).On the other hand, the program also plots the results as the tunnel inflows Q (l/s) versus time (d) (3, **Figure 2**).

 

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/tunflow_1_04_d4f5cb40c1.png" alt=""/>
 <figcaption>Figure 2. Results and output options of TunFlow.
</figcaption>
</figure>


The input parameters can be also exported as an xlsx file (1 and 2 in **Figure 3**) and used for new simulation runs.

 


<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/tunflow_1_05_55b81820f6.png" alt=""/>
 <figcaption>Figure 3. Preparation of the input parameters in an Excel spreadsheet.
</figcaption>
</figure>

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/tunflow_1_06_20968f3d42.png" alt=""/>
 <figcaption>Figure 4. Import of the previously created Excel file into Tun Flow.
</figcaption>
</figure>

 
Acknowledgements
===============
This work has been conducted at Amphos 21 in collaboration with Santos Jiménez, Paolo Trinchero and Jorge Molinero.

References
=========
*Perrochet, P., Dematteis, A.*, 2007. **Modeling Transient Discharge into a Tunnel Drilled in a Heterogeneous Formation.** Ground Water 45 (6), 786–790.

medium_tunflow_1_01.png

small_tunflow_1_01.png

thumbnail_tunflow_1_01.png

tunflow_1_01.png

In this post, Jordi Font, Senior Consultant and Project Manager at Amphos 21, blogs about a program (developed at AMPHOS 21)

Tunflow 1.0 - How to calculate groundwater inflow into tunnels and galleries

The challenge
=============

The Swedish Organization for Radioactive Waste Management, SKB, plans the disposal of spent nuclear fuel in a deep geological repository. Bentonite clay is planned to be used around the cylindrical waste packages as buffer material due to its low permeability and high swelling properties. To assess alternative clay materials, SKB started a field experiment in 2006 [1] where packages of different types of compacted bentonite hollow cylinders were placed and heated (> 100 ºC) from the inside to target values during ~2.5yr (Figure 1). Post-mortem analysis indicated that initially contrasting cation-exchanger compositions between different bentonite blocks significantly homogenized already after 880 days (Figure 2). The question arose of whether these observations can be reasonably well described by our conventional tools and understanding of the underlying processes. 




<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modeling_bentonite_figure_01_902f919667.jpg" alt=""/>
 <figcaption>Installation of the package and conceptualization in Comsol model. </figcaption>
</figure>

 

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modeling_bentonite_figure_02_12c1723b4f.png" alt=""/>
 <figcaption>Initial (before) and final (after the experiment) exchanger composition in terms of unit fraction Na, Ca, Mg and K occupancy (empty: no data available). </figcaption>
</figure>



 
Comsol: an ideal platform for simple tailored models
=======================================
 

To verify whether a simple reactive transport model coupling diffusion, cation-exchange and temperature evolution can explain the observed homogenisation of cation-exchanger composition, a simplistic geochemical model was implemented in Comsol and coupled to physical processes. The model considers different porewater compositions for each clay material assuming chemical equilibrium and only 5 ionic species: Na+, Ca2+, Mg2+, K+ and Cl-. Fickian diffusion (coupled to temperature evolution based on experimental data) is the only transport process while cation-exchange is the only type of chemical reactions implemented. The implementation of diffusion uses the “Solute Transport” Interface, while cation-exchange is set through a system of coupled non-linear algebraic equations. Transport-Chemistry coupling uses the approach proposed in [2] and solved using Operator Splitting. More details can be found in [3].

 

Results: when modeling meets real experiments
=====================================

In conjunction with a composite boundary type (Figure 3) based on fracture logging measurements, the simple model is capable of reproducing the experimental data with reasonably good accuracy, especially for the dominant cations (Ca and Na). A relatively poorer fit for Mg is observed in the top (although still unclear, it could be due to additional chemical reactions or altered groundwater composition). Thus, we see how a very simple model (Figure 4) can be very useful to gain insight in the experiment evolution and serve as a basis for discussion between experimentalists. The model independently showed that part of the homogenization is most probably due to the higher water renewal in the fractured zone of the deposition hole.

 


<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modeling_bentonite_figure_03_d0d83660ea.png" alt=""/>
 <figcaption>Left – schematic representation of the model domain and its boundaries (blue: no-flux, red: constant concentration). Right – final experimental and simulated exchanger compositions (after 880 days) across the package.</figcaption>
</figure>

 

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/modeling_bentonite_figure_04_9c89975751.png" alt=""/>
 <figcaption>Time evolution of Ca unit fraction on the cation exchanger. From initial 
</figcaption>
</figure>
 

Acknowledgments
==============
This work was made possible by the financial support from SKB and has been conducted at Amphos 21 by Marek Pekala, Andrés Idiart, Albert Nardi, and David Arcos. Their contributions are gratefully acknowledged.

References
=========
[1] *Svensson D. et al.,* 2011. Alternative Buffer Material – **Status of the ongoing laboratory investigation of reference materials and test package **1. SKB TR-11-06.

[2] *Saaltink M. et al., *1998.** A mathematical formulation for reactive transport that eliminates mineral concentrations. **Water Resources Research, 34, 1649-1656.

[3] *Idiart AE, Pękala M, Nardi A, Arcos D,* 2012. **Reactive transport processes in compacted bentonite. Application to a prototype experiment of underground repository for nuclear waste. **Int. Meeting Clays in Natural & Engineered Barriers for Radioactive Waste Confinement, Montpellier, France.

large_modeling_bentonite_figure_04.png

medium_modeling_bentonite_figure_04.png

small_modeling_bentonite_figure_04.png

thumbnail_modeling_bentonite_figure_04.png

modeling_bentonite_figure_04.png

This post is about a recent application case of Comsol in the field of thermo-hydro-geochemical behavior of clayey materials 

Modeling reactive transport processes of an in-situ test in compacted bentonite

Background and motivation
======================

In performance assessment (PA) studies of nuclear waste repositories, numerical simulations are used to study the potential impact of the release and transport of radionuclides in the near and far field. These simulations are typically carried out using a lumped approach: the different and highly non-linear underlying radionuclide retention processes are represented by a single parameter, the distribution coefficient (Kd), which relates the radionuclide mass retained in the solid phase to its aqueous concentration.

Kd-based simulations, which are particularly appealing due to their numerical efficiency, are mathematically robust but usually rely on an important assumption: the Kd value of each radionuclide is constant in time. This assumption could be critical when geochemical conditions are changing dramatically (e.g. in the transition from temperate to glacial conditions), as distribution coefficients depend on the pH, redox conditions and major chemistry of the system.

We have recently developed and tested a framework (interface to Marfa - iFM) that constitutes a significant improvement over the traditional constant Kd-based models. The basic premise of this approach is that the “background” geochemistry (excluding radionuclides) is first be pre-calculated using a reactive transport code in order to generate a dataset of Kd values in space and time that are representative of the actual geochemical conditions in the considered spatial and temporal domain. This information is then used by the computer code Migration Analysis of Radionuclides in the Far Field (MARFA) [1,2] to calculate radionuclide transport based on these “intelligent” Kd values.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/smart_Kd_figure_1_bb6a5d7027.png" alt=""/>
 <figcaption>Sketch showing a typical radionuclide transport problem. Radionuclides, which are released at repository depth, move along a number of transport pathways and are discharged into the biosphere.</figcaption>
</figure>

The tool: iFM version 2.0
=====================
The figure below shows a workflow of iFM v2.0.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/smart_Kd_figure_2_2eeb82baaa.png" alt=""/>
 <figcaption>Functionality overview of iFM v2.0. </figcaption>
</figure>

The numerical tool consists of two sets of numerical calculations:

- FASTREACT [3] simulates transport and geochemical reactions along the same set of transport pathwasy used by MARFA. This simulation only accounts for the background geochemistry (without radionuclides).
- PHREEQC batch calculations. At each time step, the background geochemical composition computed along the pathways (previous step) and an arbitrarily low concentration of the considered radionuclide(s) are used to perform zero-dimensional calculations. These batch calculations use a mechanistic description of the relevant sorption processes (e.g. surface complexation, cation exchange, etc) and provide as output quantitative information about the distribution coefficients, Kd, of each radionuclide at each location and at the considered simulation times.
 

Testing the tool – a synthetic example
======================

To test the robustness of the tool, we have solved a “typical" radionuclide waste performance assessment (PA) problem: the simulation of the migration of radionuclides from a deep geological repository to the surface.

A snapshot of a groundwater velocity field has been taken from an existing regional-scale three-dimensional hydrogeological model. Using particle tracking simulations, a number of transport pathways have been delineated.

Radionuclides have been released during a small time interval compared to the time frame of the problem. Initially, glacial conditions are assumed to prevail in the whole domain. After around 20,000 years, a different and less diluted water, representative of new temperate conditions, enters the model domain from the release boundary and propagates undisturbed along the streamlines.

The numerical exercises have been solved using (1) iFM+MARFA (see Figure below) and (2) FASTREACT.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/smart_Kd_figure_3_11bb3789da.png" alt=""/>
 <figcaption>Sketch that shows the approach and the different processes included in the iFM+MARFA simulations.</figcaption>
</figure>



The strontium breakthrough curves computed using iFM+MARFA and FASTREACT are shown in the Figure below. When strontium is overtaken by the geochemical perturbation, it is mobilized and moves nearly unretarded. Both simulations capture this effect and their results closely match each other. Thus, the verification exercise can be considered successful.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/smart_Kd_figure_4_671d85dc4e.png" alt=""/>
 <figcaption>Calculated breakthrough curve for strontium (iFM+MARFA: blue solid line; FASTREACT: orange dashed line).</figcaption>
</figure>

Acknowledgements
===============
This work is a collaborative project between SKB and POSIVA, the organizations in charge of managing nuclear fuel and radioactive waste in Sweden and Finland, respectively. Amphos 21 and Scott Painter want to thank Lasse Koskinen (POSIVA) and Jan-Olof Selroos (SKB) for their support and very helpful comments.

Paolo Trinchero also wants to thank the coauthors of this work: Luis Manuel de Vries, Jorge Molinero and Scott Painter.

References
=========
[1] *Painter, S., V. Cvetkovic, J.-O. Selroos, and O. Pensado.* (2007). **Time domain particle tracking methods for simulating transport with retention and first-order transformation. **Water Resources Research 43(9), W01406.

[2] *Painter, S., and Mancillas, J. *(2013). **MARFA User’s manual: migration analysis of radionuclides in the far field, **Working Report 2013-01, Posiva Oy.

[3]* Trinchero, P., Molinero, J., Román-Ross, G.* (2013) **A streamline-based approach for the solution of multicomponent reactive transport problems,** SKB report R-10-45.

medium_smart_Kd_figure_1.png

small_smart_Kd_figure_1.png

thumbnail_smart_Kd_figure_1.png

smart_Kd_figure_1.png

Distribution coefficients, also known as Kd values, are lumped parameters that contain information about the radionuclude-soil/rock interaction processes

An “intelligent” Kd approach for the simulation of radionuclide transport in the geosphere - iFM version 2.0

The safety assessment of deep geological repositories (DGR) for used spent fuel (SF) requires a fundamental understanding of the processes controlling fuel alteration and the release of radionuclides into the geosphere.

The corrosion of the containers in which the spent fuel is stored will eventually lead to container failure and groundwater inflow. Water entering the failed container will be devoid of oxygen due to the corrosion of steels and cast iron, ensuring a low redox potential by the production of H2 and Fe(II). Subsequently, the oxygen depleted water will get in contact with the spent fuel experimenting radiolysis by the effect of alpha radiation (emitting gamma and beta fission products are expected to be decayed at repository times greater than 10.000 years, when failure of the container is expected). Water radiolysis generates both oxidizing and reducing species in solution: H+, OH-, O2, H2O2, H2, HO2-, HO2·, O·, O-, O2-, H·, ·OH and e- [1]. Taking into account that the reducing effect of H2 is inhibited by its slow kinetics, radiolysis could potentially lead to relatively fast oxidative dissolution despite the bulk reducing environment. However, some fission products contained in the irradiated fuel have shown to be able to activate H2 through their catalytic effect. We refer to the epsilon particles (Rh, Pd, Ru), whose catalytic effect can result in low H2 concentrations, suppressing the oxidation of UO2 [2-4]. The high complexity of the described system is schematically represented in **Figure 1**. 

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_1_0931a293fb.png" className="mx-auto" alt=""
 />
 <figcaption>Figure 1. Interrelated processes involved in the spent fuel matrix alteration.</figcaption>
</figure>

In the framework of EURATOM Project “DisCo” (modern spent fuel Dissolution and chemistry in failed Container conditions), carried out within the Horizon2020 framework program, Amphos 21 has integrated the processes occurring at the spent fuel surface (see **Figure 1**) with a 1D reactive transport model. The model has been implemented in iCP [5], using the two-way coupling approach shown in **Figure 2**. 

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_2_0a79e24c73.png" className="mx-auto" alt=""
 />
 <figcaption>Figure 2. Integration in iCP of the different processes involved in the SF alteration.</figcaption>
</figure>

The complete water radiolysis scheme, integrating more than 35 kinetic reactions together with the yields of primary products [1, 6] (included in Appendix 1) has been implemented in COMSOL though a set of ordinary differential equations (ODE) and were coupled with the solute transport equation (Eq. 1) through a source term in Riba et al. (2020) [7].

$$
\phi\frac{\partial c_i}{\partial t}+\mathrm{\nabla}\left(-D_e\mathrm{\nabla}c_i\right)=S_i+\sum_{j} R_{ij} 							$$

The SF matrix has been modelled as amorphous UO2 with a homogeneous chemical composition and containing 1 atom % of Pd, as representative element of the epsilon particles. The kinetic reactions related to the alteration of SF and Fe(s) corrosion are implemented in PhreeqC (see **Figure 2**).

The interrelated processes described above have been implemented in a 1D model with the geometry presented in Figure 3 to simulate a DGR concept where the iron insert of the copper canister is at 10 mm from spent fuel rods [8]. A uniform alpha dose rate of 2.86·10-2 Gy·s-1 affecting the first 13 &mu;m in water adjacent to the SF surface was assumed in the 1D simulation. This uniform alpha dose rate is equivalent to a non-uniform exponential distribution of dose rate affecting 35 μm in water adjacent to the SF surface, in the sense that both produce the same total dose rate [9].

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_3_5ca2e9d13f.png" className="mx-auto" alt=""
 />
 <figcaption>Figure 3. Geometry (not to scale) considered in the 1D models that consider a geometry comparable to that of the spent fuel in the DGR.</figcaption>
</figure>

**Figure 4** shows the temporal evolution of the total U, U(IV) and Fe2+ concentration profiles. The total U concentration is in accordance with the Eh in the different domains and at different simulation times. Notice the sharp increase of both [U] and Eh after 20 days of simulation time. 

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_4_affde9483b.gif" className="mx-auto" alt=""
 />
 <figcaption>Figure 4. Modelling results considering the geometry detailed in Figure 3 and porous matrix of SF of 4mm.</figcaption>
</figure>

This modelling approach has solved the challenge of coupling the complete water radiolysis system with chemical complexation reactions and dissolution/precipitation processes which occur at very different time scales, often with rates differing by more than 6 orders of magnitude.

The model not only integrates the radiolysis of water but also these two important systems: uranium (as the main element of the fuel) and the iron of the steel-container and metallic insert, having, both of them, a very complex chemistry. Also, it has the flexibility for future inclusions of other chemical processes to account for the complex system of SF, steel container and water chemistry. Therefore, it represents a step forward towards the development of the models needed to assess the performance and the safety of a DGR for spent nuclear fuel.

References
========
[1] *Kelm, M., Bohnert, E. *(2004) **A kinetic model for the radiolysis of chloride brine, its sensitivity against model parameters and a comparison with experiments**, Forschungszentreum Karlsruhe, FZKA 6977.

[2] *Carbol, P. Fors, T.Gouder, K. Spahiu *(2009). **Hydrogen suppresses UO2 corrosion**. Geochimica et Cosmochimica Acta 73 (2009) 4366–4375.

[3] *Trummer, T., M. Jonsson M, J.* (2010). **Resolving the H2 effect on radiation induced dissolution of UO2-based spent nuclear fuel **Nucl. Mater., 396, 163–169.

[4]* Ollila K* (2011). **Influence of Radiolysis on UO2 Fuel Matrix Dissolution Under Disposal Conditions - Literature Study**. Posiva Working Report 2011-27, 64 pp. Posiva Oy, Olkiluoto, Finland.

[5] *Nardi, A., Idiart, A., Trinchero, P., de Vries, L. M., and J. Molinero *(2014). **Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry. **Computers & Geosciences 69: 10-21.

[6] *Eriksen, T. E., Jonsson, M., Merino, J.* (2008). Modelling of time resolved and long contact time dissolution studies of spent nuclear fuel in 10 mM carbonate solution–a comparison between two different models and experimental data. Journal of Nuclear Materials, 375(3), 331-339.

[7] *Riba, O., Coene, E., Silva, O., Duro, L.* (2020). **Spent fuel alteration model integrating processes of different time-scales.** MRS Advances, 5(3), 159-166.

[8] *Raiko *(2012). **Canister design 2012**. Posiva Report 2012-13.

[9] *Wu, L., Beauregard, Y., Qin, Z., Rohani, S., Shoesmith, D. W.* (2012). **A model for the influence of steel corrosion products on nuclear fuel corrosion under permanent disposal conditions**. Corrosion Science, 61, 83-91.

[10] *Christensen, H., Sunder, S., Shoesmith, D. W.* (1994).** Oxidation of nuclear fuel (UO2) by the products of water radiolysis: development of a kinetic model. Journal of alloys and compounds**, 213, 93-99.

[11] *Cera, E., Bruno, J., Duro, L., Eriksen, T. *(2006). **Experimental determination and chemical modelling of radiolytic processes at the spent fuel/wáter interface**. SKB TR 06-07, Svensk Kärnbränslehantering AB.

Appendix A
=========

Kinetic reaction scheme for water radiolysis [1]

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_appendix_1_8b4339d20d.png" className="mx-auto" alt=""
 />
</figure>

Yields of primary products (mol·l-1 x 107) formed on radiolysis of water

<figure className="image">
 <StrapiImage
 src="/uploads/euratom_safety_assessment_appendix_2_875e8e8d0d.png" className="mx-auto" alt=""
 />
</figure>

large_euratom-safety-assessment_1.png

medium_euratom-safety-assessment_1.png

small_euratom-safety-assessment_1.png

thumbnail_euratom-safety-assessment_1.png

euratom-safety-assessment_1.png

The safety assessment of deep geological repositories (DGR) for used spent fuel (SF) requires a fundamental understanding of the processes controlling fuel alteration 

Spent fuel matrix alteration model under Deep Geological Repository conditions

<div className="">
 <a className="inline-flex items-center text-gray-500 bg-transparent border border-solid border-gray-500 hover:bg-gray-100 hover:text-white active:bg-gray-600 font-bold uppercase text-sm px-6 py-3 rounded-full outline-none focus:outline-none mr-1 mb-1" href="https://dose.amphos21.com/" target="_blank">Launch Experiment</a>
</div>

How often have you needed to convert mass to activity units or vice versa? This is something that we constantly need to do at Amphos 21 in:

- **Dose calculations**: When converting radionuclide concentration in the environment (soil, groundwater, air…) to activity.
- **Radionuclide inventory**: activity units are commonly used to express radionuclide concentration and they should be converted to mass units to perform geochemical and transport modelling.
- **Geochemical modelling**: this type of calculations is mainly performed in molality but concentration limits in groundwater could be given in mg/L.
- **Bibliographic search**: Activity concentration in old publications is given in Curies but currently, the international unit for activity is Becquerel.
- **Experimental work**: In the lab, conversion from mass to activity, and vice versa, is needed day to day. There are several units (Becquerel, mol, gram…) used to express the amount of a radionuclide in any media (soil, waste, solid or liquid matrix, groundwater, air, etc.) and the conversion from one unit to another is not always direct.

The Experiment we propose in this post is a fast and easy-to-use tool to help on this type of calculations. Amphos 21 presents here a RadioUnits Converter that allows the user to perform instant conversion between the following units: mol, g, Bq, Ci.

large_units_conversor.png

medium_units_conversor.png

small_units_conversor.png

thumbnail_units_conversor.png

units_conversor.png

Amphos 21 presents here a RadioUnits Converter that allows the user to perform instant units conversion

RadioUnits Converter

Knowledge management has increased its importance during the last years. It is crucial to transfer knowledge to newcomers at the radioactive waste management sector, students and young generations. Games are an excellent tool to engage them to learn more on concepts related with PREDIS.

 

The project:
=========
 

PREDIS (Pre-disposal management of radioactive wastes) – 2021-2024 – Funding received from the Euratom research and training programme 2019-2020 under grant agreement No 945098.

 

The PREDIS project targets the development and implementation of activities for pre-disposal treatment of radioactive waste streams other than nuclear fuel and high-level radioactive waste. The PREDIS project develops and increases the Technological Readiness Level of treatment and conditioning methodologies for wastes for which no adequate or industrially mature solutions are currently available, including metallic materials, liquid organic waste and solid organic waste. PREDIS project also develops innovations in cemented waste handling and pre-disposal storage by testing and evaluating.

 

Games:
=========

 

**Release 1** – 25-04-2022:

- Wordle [link](https://mywordle.strivemath.com/?word=lfvgtw)
- Words Soup [link](https://nepazing.com/WordSearch/Play/88C431A9)
- Answers & more information  [link](https://back.amphos21.com/uploads/Answers_predis_1_22_04_25_d26b67f65c.pdf)

**Release 2** – 16-05-2022:

- Wordle: [link](https://mywordle.strivemath.com/?word=ncrgxel)
- Quiz: [link](http://take.quiz-maker.com/Q6J10N8DD)
- Puzzle: [link](https://puzzel.org/es/jigsaw/play?p=-N13zKupUCJn7c0BNhry)
- Answers & more information  [link](https://back.amphos21.com/uploads/Answers2_22_05_16_f58fe324fb.pdf)
 

large_predis_logo_green-2048x613.png

medium_predis_logo_green-2048x613.png

small_predis_logo_green-2048x613.png

thumbnail_predis_logo_green-2048x613.png

predis_logo_green-2048x613.png

Knowledge management has increased its importance during the last years. It is crucial to transfer knowledge to newcomers at the radioactive waste management sector, students and young generations.

PREDIS knowledge management: Learning by playing

Since 2000, Amphos 21 has been deeply involved in the development and maintenance of ThermoChimie, the database managed by Andra (France), Nuclear Waste Services (UK) and Ondraf/Niras (Belgium).

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/thermochimie_11_1_088c1c1dd8.png" alt=""/>
</figure>
The ThermoChimie database is under constant revision and development. After its public release in 2014 and further updates in 2015, 2018 and 2021, a new database version (11a) is now <a href="http://www.thermochimie-tdb.com"> available on-line</a> .

ThermoChimie version 11 contains significant updates on U, Np, Pu, Am/Cm, and Tc following the NEA-TDB update (Grenthe et al. 2020) and beyond. In addition, an exhaustive and accurate review of Fe thermodynamic properties has been carried out, using the NEA-TDB reviews (Lemire et al. 2013, 2020) and additional literature publications. Furthermore, PFlotran extraction has been newly added to the list of geochemical compatible codes.

<figure className="image">
 <StrapiImage className="mx-auto" src="/uploads/thermochimie_11_2_11c25fd05a.png" alt=""/>
</figure>
A new XCheck tool to consult, filter, retrieve and export the information available in the latest database version is also now available for the user’s at <a href="https://register.thermochimie-ctool.com">https://register.thermochimie-ctool.com</a>.

Please stay tuned to the website http://www.thermochimie-tdb.com/ or contact consortium-tc3-tech@thermochimie-tdb.com for additional information.

large_thermochimie_11_1.png

medium_thermochimie_11_1.png

small_thermochimie_11_1.png

thumbnail_thermochimie_11_1.png

thermochimie_11_1.png

After its public release in 2014 and further updates in 2015, 2018 and 2021, a new database version (11a) is now available on-line.