<figure><img src="image_1.png" /></figure>

HAL
open science

# Role of reactive transport in the alteration of vitrified waste packages: the MOS model

Pierre Frugier, Nicole Godon, Yves Minet

► To cite this version:

Pierre Frugier, Nicole Godon, Yves Minet. Role of reactive transport in the alteration of vitrified waste packages: the MOS model. npj Materials Degradation, 2024, 8 (1), pp.91. 10.1038/s41529-024-00496-0. cea-04716779

HAL Id: cea-04716779

https://cea.hal.science/cea-04716779v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Role of reactive transport in the alteration of vitrified waste packages: the MOS model

<figure><img src="image_4.png" /></figure>

Check for updates

Pierre Frugier✉, Nicole Godon & Yves Minet

The MOS model (acronym coming from the French MOdèle Simplifié) was born from the desire to have a simple tool that can quantify the contribution of the diffusive reactive environment to the alteration of a vitrified nuclear waste package in deep geological disposal conditions. In the model, this environmental contribution consists partly of the ability of iron, metallic casing corrosion products, and argillite to consume silicon, and partly of the brake on diffusive transport provided by silicon through the successive layers of environmental material. It is a modeling tool serving as an intermediary between operational modeling for the calculation of the source term from the glass, mathematically more simple and giving higher upper margins, and models that use geochemistry and transport, giving greater accuracy for the interactions between glass and its environment. The goal of the MOS model is to calculate the possible impact of silicon reactive diffusion on the alteration rate within the different layers of material surrounding nuclear glass. This article lists the simplifying hypotheses on which the MOS is based, presents the digital resolution method for an environment consisting of several successive layers with different reactivity and transport properties, and explains the model's implementation.

The non-recyclable residues from spent fuel reprocessing that are produced by the electro-nuclear industry account for the bulk of French High Level waste. To contain this type of waste, France uses a vitrification process. The country has also carried out research studies into ways of managing the resulting vitrified waste. These studies led to the 2006 decision to implement a deep geological repository, and the French National Nuclear Waste Management Agency (ANDRA) was entrusted with further investigations into this subject.

During the preparation of the file accompanying the Creation Authorization Request for the repository site¹,², preliminary work was carried out by ANDRA in order to establish the release models used in the performance calculations. The release model decided on for R7T7 type High-Level glasses was based on the Operational Model (MOP) $V_0 \rightarrow V_r$ developed by the CEA in 2004³,⁴. It is a very simplified representation where the glass source term dominates in repository conditions. It gives a global value for the release of radionuclides from the vitreous matrix through time ("glass term source"), and from which the "package lifetime" can be deduced, i.e., the time necessary to release the entire radioactive inventory contained in the package. In the MOP $V_0 \rightarrow V_r$, the alteration rate is assumed to be equal to the initial alteration rate (maximum rate measured for a given glass composition, temperature, and solution composition),

while the surrounding medium is reactive and does not allow the development of conditions appropriate for the rate to decline. It then drops immediately to the residual rate value (several degrees of magnitude lower than the initial rate) when the medium's reactivity is exhausted. Two fracturing ratios associated with each of the rate regimes enable the associated reactive surface to be taken into account⁵.

In an approach even closer to the real situation, the CEA developed a geochemical model for glass alteration, called GRAAL⁶⁻⁸. It was implemented in a reactive transport code and coupled with behavior models describing the environmental materials (iron, argillite). This enabled the model to describe the main physical and chemical mechanisms that govern glass alteration at a macroscopic scale. Compared to short-term laboratory experiments, and used to identify the dominant long-term mechanisms, the GRAAL model has been useful in confirming simplifications and in giving the best margin options in operational release models⁷⁻¹⁶.

This article presents a new modeling tool with a simplification level between those of the MOP $V_0 \rightarrow V_r$ and of the GRAAL model. The modeling tool is called MOS, for (in French) MOdélisation Simplifiée (Simplified Modeling). The MOS was designed to meet the need for an easy to use tool able to quantify the contribution from a diffusive reactive environment to glass alteration. The focus was on:

<figure><img src="image_5.png" /></figure>
- the ability of iron, corrosion products, and argillite to consume silicon, whose concentration in the vicinity of the glass controls the glass alteration rate.

- the diffusive transport of silicon through the successive layers of material in the environment and the way the resulting flux decreases over time. This slowing of matter exchanges decreases glass alteration.

The brake on transport means that only materials present in the immediate vicinity of the glass will be able to maintain a high alteration rate. The current ANDRA design for the disposal of High-Level waste foresees the existence of different materials surrounding the containment glass (Fig. 1). Moving outward from the glass, we find (1) a stainless steel container into which the glass was poured, forming the primary package, (2) the storage container in non-alloy steel, which constitutes the storage overpack, (3) a tunnel liner, a long non-alloy steel tube inserted into a storage tunnel dug through the rock, (4) an alkaline filling material injected between the liner and the host rock, called Annular Space Filling Material (in French, Matériau de Remplissage de l'Espace Annulaire or MREA), (5) the host rock (Callovo-Oxfordian clay). During the repository's operation period, air oxidation of iron sulfate in the clay closest to the liner will free sulfuric acid. The MREA will neutralize this acidic transition to avoid it accelerating steel corrosion.

When water in liquid or gaseous form eventually arrives in contact with the glass, the materials surrounding it will have evolved and interacted, in particular forming corrosion products, oxides, carbonates, sulfurs, and iron silicates. Each successive layer of transformed and/or newly formed residual materials will have its own ability to consume the silicon coming from the glass and slow its transfer into the environment by diffusion. These two properties, consumption in a broad sense and diffusion, are at the heart of the MOS.

The “Theory and Equations” section of this article is divided into two sub-sections. The first part gives the hypotheses regarding the dominant alteration mechanisms in the long term within a reactive environment. These simplifying hypotheses in chemistry and physics are intended to contribute to modeling that does not require the use of a calculation code by finite elements to solve the system of equations, as is the case in implementing the GRAAL model. Thus the MOS hypotheses lead to a simple analytical equation in the case of a semi-infinite homogeneous diffusive environment. This is explained in the second part.

In the “Results” section, the model’s functioning is shown. The section illustrates the behavior of the equations and gives the magnitudes of alteration rates on the vitrified waste package and its lifetime, depending on the selection of the parameters to be used.

The “Discussion” gives a summary of the main conclusions and perspectives concerning the contributions and the uses of such a model.

The last part includes i) the method for equation resolution in the general case of several successive layers with different reactivity and transport properties, ii) the input parameters, iii) the implementation of the calculations. The article is followed by two appendices, one concerning the

analytical demonstration of the resolution method and the other describing the entire digital resolution.

## Theory and equations

### MOS hypotheses

In a repository situation, and when water arrives in contact with the metallic materials in the near vicinity, the corrosion of iron maintains a dihydrogen pressure that slows saturation of the medium. The time necessary for complete saturation may well reach several hundreds of thousands of years²,¹⁷.

When the block is submerged in water, the hypothesis made in the MOS is that the glass alteration can be controlled by the diffusion of silicon in the package environment and by its consumption in contact with reactive minerals. The central role of silicon in the silicate glass alteration process is generally agreed upon: it is the main element in the glass and in the more or less passivating gel that develops on its surface during alteration. It is also one of the most soluble elements among all the poorly soluble elements that precipitate to form the gel. It must be present in solution in sufficient quantities for a passivating gel to form, and for the alteration rate to decrease by several magnitudes¹⁸⁻²¹. Silicon's major role is exploited in the MOS to simplify process modeling. The objective of the MOS is to quantify the consumption of silicon by the environment and, from this to deduce, in a first approximation, the effect it can have on glass alteration.

A matter balance for silicon is written for the interface between the package and its environment. This involves writing the flux equality between the silicon from the glass and the silicon consumed by its environment (Fig. 2). The simplified model shows the link between the alteration rate and the distances between the package and the reactive materials. It highlights the reactivity limitation of materials in the vicinity depending on their nature and on their distance from the glass package.

The hypotheses of the MOS calculation tool had to be validated to enable this tool to be used. They came from previous operational modeling approaches³,⁴ and from simplifications necessary for the calculation resolution, which is the subject of this document. The hypotheses enabled the listing of the mechanisms explicitly not to be taken into account before justifying this choice and the way in which the calculation is intentionally somewhat overestimated. They are summarized in Table 1.

A number of these hypotheses could be proved and validated with a more sophisticated chemistry-transport modeling. Table 2 summarizes the mechanisms simplified in the MOS compared to the mechanisms taken into account via the couple formed by a geochemical model for glass alteration, i.e., the GRAAL model⁸, and a chemistry-transport calculation code²².

## The equation for a homogeneous environment

A spherical geometry is used in the MOS for equation resolution (Hypothesis G1, Table 2). The use of a cylindrical 2D geometry would lead to an underestimation of the ability of the environment to remove matter, as one removal direction would be forgotten. Furthermore, 3D modeling would require the use of a reactive transport code, given the complexity of the equation resolution. For a spherical medium, the equation resolutions are simple.

For a spherical, homogeneous, and semi-infinite medium, the equality of the flux between (i) the quantity of matter coming from the glass dissolution (taking into account the retention of silicon in the alteration film)

<figure><img src="image_6.png" /><figcaption>Fig. 1 | Cross-section view of the materials surrounding vitrified waste in the current storage design²⁷.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 2 | In the MOS, a matter balance for silicon is written for the interface between the package and its environment. The diagonal dash-dotted line represents a linear concentration profile in the diffusive barrier.</figcaption></figure>
<table><caption>**Table 1 | Summary of the base hypotheses for the MOS**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td colspan="2"><strong>Intrinsic behavior (without the presence of water)</strong></td></tr><tr><td>I1</td><td>The long-term evolution of the glass in terms of crystallization is ignored.</td></tr><tr><td>I2</td><td>The long-term evolution of the glass in terms of self-irradiation is ignored.</td></tr><tr><td>I3</td><td>The long-term evolution of the glass in terms of fracturing is ignored.</td></tr><tr><td>I4</td><td>The volatility of radionuclides is ignored.</td></tr><tr><td colspan="2"><strong>Transport in solution</strong></td></tr><tr><td>T1</td><td>Ion transport in the solution within the pores of the package can slow glass alteration. This is not taken into account. Pore water in the block saturates in silicon at the maximum rate.</td></tr><tr><td>T2</td><td>Convection (solvent movement) in the package environment is ignored.</td></tr><tr><td>T3</td><td>Silicon is the only element whose diffusion in the package environment is taken into account.</td></tr><tr><td>T4</td><td>The porosity of each surrounding material is input data for the calculation. Its evolution through time is not taken into account.</td></tr><tr><td>T5</td><td>The diffusion coefficient within the pores of each surrounding material is input data for the calculation. Its evolution through time is not taken into account.</td></tr><tr><td>T6</td><td>The reactive diffusion calculation is intentionally an overestimation, using simplifying hypotheses presented in this article.</td></tr><tr><td colspan="2"><strong>Physico-chemistry in solution</strong></td></tr><tr><td>C1</td><td>The chemistry of the dissolved ions is not described: the acid-base, complexation, and oxidoreduction reactions are not taken into account.</td></tr><tr><td>C2</td><td>Temperature is constant in time and space.</td></tr><tr><td colspan="2"><strong>Geometry</strong></td></tr><tr><td>G1</td><td>The diffusive fluxes are calculated in spherical 1D geometry outward from the package radius.</td></tr><tr><td>G2</td><td>The exterior surface geometry of the package is input data.</td></tr><tr><td>G3</td><td>The positions and thicknesses of the different layers of surrounding material are input data. Their evolution through time isn't taken into account.</td></tr><tr><td colspan="2"><strong>Chemical reactivity of materials in the environment</strong></td></tr><tr><td>R1</td><td>The chemical reactivity of the surrounding materials is taken into account via the silicon partition coefficients. These coefficients, specific to each material, are the ratio between the maximal concentration in the solid after chemical reaction and the maximal concentration in solution that can be found in contact with the glass. This mathematical formalism is very useful to simplify the resolution of the diffusive transport equations and is relevant if the silicon consumption reaction is rapid compared to its diffusive transport.</td></tr><tr><td colspan="2"><strong>Glass reactivity in the vapor phase</strong></td></tr><tr><td>V1</td><td>The alteration rate during the vapor phase is input data, and its evolution through time isn't taken into account.</td></tr><tr><td>V2</td><td>The fraction of the package altering in the vapor phase at a given moment is input data.</td></tr><tr><td colspan="2"><strong>Glass reactivity in aqueous phase</strong></td></tr><tr><td>A1</td><td>The total and maximal concentration of silicon in solution in contact with the glass is input data.</td></tr><tr><td>A2</td><td>The maximal glass alteration rate is input data.</td></tr><tr><td>A3</td><td>The minimal glass alteration rate is input data.</td></tr><tr><td>A4</td><td>The glass block fracturing ratio is input data. The "effective" fracturing ratio of the glass block under initial rate condition is input data.</td></tr><tr><td>A5</td><td>The silicon retention ratio in the glass alteration products under initial rate condition is input data. The silicon retention ratio in the glass alteration products applicable in all other rate regimes is input data.</td></tr></tbody></table>

and (ii) the quantity of silicon diffusing into the environment leads to Eq. 1. This equation includes:

* $t(s)$ is the time

* $V(t)(g_{\text{glass}} m^{-2}s^{-1})$, is defined as the glass alteration rate per surface unit of the reactive glass/water interface. This reactive surface is assumed to be the developed surface of the glass package accessible to renewed water under initial-rate conditions ($S_{\text{acc}} = S_0 T_{x0}$). This

surface area is greater than the geometric surface area of the package ($S_0$) but less than the total surface area of the cracks ($S_0 T_{xf}$). The majority of these cracks, namely ($S_0(T_{xf} - T_{x0})$), are assumed to be already altered under residual rate condition ($V_f$) at the initial time of alteration.

* $x_{\text{Si}}(g_{\text{Si}} \cdot g_{\text{glass}}^{-1})$, the mass fraction of silicon in the unaltered glass,

* $T_{x\text{Si}}^{\text{ret}}$, the silicon retention ratio within the alteration film,

* $M_{\text{Si}}$ (28.1 g · mol$^{-1}$), its molar mass,

* $D(m^2 \cdot s^{-1})$, its effective diffusion coefficient in the environment,

* $\omega(\text{SU})$, the porosity of the medium,

* $C_{\text{sat}}(\text{mol} \cdot \text{m}^{-3}_{\text{solution}})$, the silicon concentration in contact with the package,

* $C_{\infty}(\text{mol} \cdot \text{m}^{-3}_{\text{solution}})$, the silicon concentration in the medium,

* $R_0(m)$, the package radius,

* And $T_{x0}(\text{SU})$, the effective fracturing ratio under initial rate condition, a ratio between the package glass surface accessible to renewed water, based on the geometric surface area $S_0 = 4\pi R_0^2$ of the interface between the package and its environment, and on which the matter balance is written.

Writing the equation in spherical geometry leads to the addition of a one-dimension equation of the term in $1/R_0$ (Eq. 1). The details of Eq. 1 are provided in Supplementary Information Note 1.

$$V(t) = D(C_{\text{sat}} - C_{\infty}) \left( \sqrt{\frac{\omega(1 + \Re)}{\pi Dt}} + \frac{1}{R_0} \right) \frac{M_{\text{Si}}}{x_{\text{Si}}(1 - T_{x\text{Si}}^{\text{ret}}) T_{x0}} \quad (1)$$

The reactivity coefficient $\Re(\text{SU})$ expresses the way the medium reactivity maintains the difference in the chemical potential of silicon, accelerates its transport, and slows the progression of the silicon front dissolved in the medium. $\Re$ is a partition coefficient for silicon between solid and liquid. This formalism is necessary for the equation to preserve a mathematical form that is simple to solve. $\Re$ is the ratio between the maximal concentration in the solid after the chemical reaction and the maximal concentration in the solution that can be found in contact with the glass (Cf. 5.2).

Equation 1 is a function of time that describes reactive transport in the semi-infinite medium. There is no analytical resolution in the general case of several successive layers with different properties (different reactivity, porosity, and diffusion coefficients). Meshing must therefore be used. Digital resolution for a heterogeneous medium is the main difficulty solved by the MOS. This is presented in detail in the “Method” section. This section also examines the question of water arrival, in water vapor or liquid form.

## Example of the results

The objective here is to show the model's functioning, to illustrate the behavior of the equations, and to give orders of magnitude. The values of the parameters presented in the “Method” section constitute a basic case used as a reference in this article. Certain values are known with high precision, others involve uncertainties. Thus in spite of the model's potential, work remains to be done before this tool can be used to quantitatively predict the lifetime of vitrified waste packages: (i) Validate the MOS hypotheses with the help of modelings dedicated to each hypothesis. For example, the GRAAL model coupled to a reactive transport code would enable the qualification of numerous effects linked to chemistry and transport, (ii) enter relevant, substantiated values for the input parameters and give their uncertainties, and (iii) qualify the propagation of these uncertainties throughout the package lifetime.

The potentials of the model are illustrated by the sensitivity of three different responses to two parameters: on the one hand, the silicon diffusion coefficient $D_0$ through the metallic envelopes, and on the other, the reactivity $\Re$ of these metallic envelopes with silicon. These two parameters have been chosen among those least easy to estimate. As the uncertainty of their values is high, they are, therefore, likely to have a high impact on the calculation
**Fig. 4** | Example of silicon concentration gradients in a HYTEC calculation, in spherical geometry (axial symmetry). The concentration decrease is more pronounced as one moves farther away from the glass block (colored in dark red, with the lowest concentrations illustrated in dark blue). The environments surrounding the block can be schematically represented by adapting the mesh at the interfaces.

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Fig. 5** | Time before the alteration rate drops by one order of magnitude. Sensitivity to the diffusion coefficient D and to the reactivity of the medium are here at 20 °C. The MOS automatically recalculates the D value to the calculation temperature, which is 50 °C.

## The time necessary for the glass alteration rate to reach $V_f$

Figure 6 shows the number of years necessary for the glass alteration rate to decrease to the final rate $V_f$. This corresponds to the low rate value retained for the glass and applies in operational modelings once the reactivity of the surrounding materials is exhausted. It comes from rate measurements made in laboratories over relatively long periods, during glass alteration in non-renewed solution and in the absence of any interaction with surrounding materials<sup>23,24</sup>.

The waiting time before reaching this drop in the rate is much longer than that in the previous calculation for $V_0/10$. This illustrates the fact that while it is easy to guarantee that the alteration rate does not exceed $V_0/10$, guaranteeing lower rates down to $V_0/10^4$, the order of magnitude for the ratio $V_0/V_f$ requires a better mastery of the uncertainties associated with the modeling parameters.

The greater the reactivity of the materials to silicon, the longer the alteration rate takes to equal the final rate for a given diffusion coefficient. For $\mathfrak{R} < 10^3$ regardless of the value of D, and for $D \le 10^{-11} \text{ m}^2 \text{ s}^{-1}$ regardless of the value of $\mathfrak{R}$ (case 1: ascending left side of the curve), the waiting time to reach $V_f$ increases with the diffusion coefficient. For $\mathfrak{R} \ge 10^3$ and $D \ge 10^{-11} \text{ m}^2 \text{ s}^{-1}$ (case 2: decreasing right side of the curve), the model gives an illustration of a phenomenology that is difficult to foresee. The waiting time to reach $V_f$ decreases when the diffusion coefficient increases in the metallic envelopes. This is actually due to fact that in case 1 the diffusive barrier consisting of the metallic envelopes alone is enough to slow the transfer of matter and impose $V_f$. In this situation, the diffusion front goes no further than the thickness of the metallic envelopes. In case 2, the

diffusion front extends beyond the metallic envelopes, and it is the diffusive medium outside the envelopes that enables $V_f$ to be obtained. The right side of the curve, therefore, indicates that the brake on transport in the zone outside the metallic envelopes is established all the more rapidly when the diffusion coefficient in the metallic envelope zone is high. By only slightly slowing the arrival of the matter in the zone that enables $V_f$ to be reached, it is filled more quickly. This means $V_f$ can be reached more rapidly, but it also means a greater quantity of altered glass.

In a situation where the media are extremely diffusive and reactive, the model highlights a special mathematical feature that distinguishes diffusion in a three-dimensional media from diffusion in one or even two dimensions. After a sufficiently long period of time, the flux of matter removed by diffusion can no longer decrease, and the glass alteration rate tends asymptotically towards a minimal value. This is simply due to the fact that the matter diffuses in a corona whose volume grows at the rate of the diffusion length cubed. The model then allows a comparison between the value of this asymptotic rate and that of the final alteration rate of the glass. With these parameters retained, it can be seen that the asymptotic rate is several orders of magnitude lower than the final rate $V_f$. In the case of a homogeneous medium, this rate is given by Eq. 1 (Equation for a homogeneous environment) when $t$ tends towards infinity.

## Package Lifetime: Time necessary to alter the entire package

Figure 7 gives the results concerning the package lifetime. For these parameter values, it can be seen that even when the metallic envelopes have a reactivity at $10^5$, the package lifetime is almost always greater than or equal to $2 \times 10^5$ years, and that it varies very little depending on the diffusion coefficient. This means that the transition from $V_0$ to $V_f$ occurs sufficiently rapidly to have little impact on the package lifetime. In these conditions, the final rate regime therefore sets the lifetime.

## Discussion

With the MOS, the goal was to calculate the impact that the reactive diffusion of silicon within the different layers of material surrounding glass could have on the alteration rate. This alteration rate potentially evolves between the initial hydrolysis rate (a situation where silicon diffuses extremely rapidly in the media) and the final alteration rate (a situation where silicon only diffuses very slowly in the media until reaching the $V_f$ value).

The first application exercises have led to the following general conclusions:

* The ability of the materials in the immediate vicinity of the package to consume the silicon coming from glass alteration has a high impact on the glass alteration rate. This consumption, which essentially occurs through the precipitation of secondary silicated phases, could in turn be limited by the clogging that inevitably accompanies the transformation of dense primary materials into newly-formed hydrated or oxidized minerals.
<figure><img src="image_10.png" /><figcaption>**Fig. 6** | Time necessary, as a function of the diffusion coefficient for silicon in the metallic containment envelopes *D* (at 20 °C), for the glass alteration rate to reach the final alteration rate (*V*<sub>i</sub>) for four different reactivities within these envelopes.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 7** | Number of years for the package lifetime as a function of the diffusion coefficient for silicon in the metallic envelopes *D* (at 20 °C). Five different silicon reactivities for the metallic envelopes (4 different curves. The curves for $\mathfrak{R} = 0$, $\mathfrak{R} = 10^3$, and $\mathfrak{R} = 10^4$ are very nearly overlapping).</figcaption></figure>

* The initial rate $V_0$ enables the alteration to be somewhat over-estimated, but it is not a realistic parameter for package alterations in their repository environment. Even a very reactive environment does not enable the glass block to be altered at $V_0$ for significant durations. This is particularly true when the package temperature is high because the activation energy associated with the surface reaction is higher than that associated with the transport mechanisms. Just the simple thickness of a container in stainless steel, or the long-term of its corrosion products, is enough to account for a rate drop of at least one order of magnitude.

Two major groups of perspectives are open at this stage in the model design. The first concerns all the studies associated with reinforcing the robustness of the MOS model. They consist of confirming the relevance of the model hypotheses listed in Table 1. Some of these hypotheses (T1, T3, R1, A1, A3, A5) are related to a simplification approach for reactive chemistry-transport codes. They could be explored using a glass alteration model coupled to a reactive transport code (e.g., GRAAL model/HYTEC code). As well as increasing MOS robustness, these studies will also provide arguments that can enable the identification and, consequently, the decrease

of unnecessarily high conservative margins associated with certain hypotheses.

The second group of perspectives relates to potential uses for the model: sensitivity studies, mastery, and propagation of the uncertainties. To assist with this, the MOS was coded in C++ so that it could be accessible for digital and statistical processing tools (URANIE platform, developed by the CEA²⁵). It will then become possible to determine precisely what constraints concern the intrinsic and relative values of the parameters, in order to meet a satisfactory package performance level (source term, fraction of the package altered at the end, lifetime...).

## Methods

### Resolution method for a heterogeneous environment

Equation 1 (Equation for a homogeneous environment) is a function of time that describes reactive transport in a semi-infinite media. There is no analytical resolution in the general case of several successive layers with different properties (different reactivity, porosity, and diffusion coefficients). Meshing must be used.

The complexity of the calculation given in Eq. 1 mainly comes from the fact that transport and reactivity are simultaneous. The silicon is transported and consumed at the same time, whether it is saturating the solution in the pores or saturating the reactivity of the surrounding materials. Thus the idea behind the simplified resolution of the MOS for an environment that is heterogeneous in reactivity and in transport properties was to decouple chemistry and transport. This idea is supported by the following points:

- The reactivity of silicon with the environment is fast, compared to its diffusive transport. This hypothesis, though over-estimating somewhat, is consubstantial with the simplifying hypothesis R1 (Table 1). If the reactional kinetics are fast, decoupling transport can be considered.

- Decoupling is justified by the large number of mesh nodes used, similar to what occurs in a reactive transport calculation code where the solution composition is homogeneous within a calculation mesh node.

In practice, the idea is:

1. to calculate the time necessary to fill and saturate the reactivity of a mesh node. This includes the quantity of dissolved silicon plus the quantity of silicon retained in newly formed solids.

2. to take into account the brake on diffusive transport created by a mesh node once its reactivity is exhausted, before filling the next mesh node.

This idea is illustrated in Fig. 8. In this diagram, the first mesh node fills at a maximal rate imposed by the maximal initial dissolution rate of the glass
<figure><img src="image_12.png" /><figcaption>**Table 3** | Physic al parameters of the package</figcaption></figure>

<table><thead><tr><th>Physical parameters of the package</th><th>Notation</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Block radius</td><td>R<sub>0</sub></td><td>m</td><td>0.21</td></tr><tr><td>Block height</td><td>H<sub>0</sub></td><td>m</td><td>1.1</td></tr><tr><td>Real geometric surface of the block</td><td>S<sub>0</sub></td><td>m<sup>2</sup></td><td>1.73</td></tr><tr><td>Internal porosity of the block</td><td>ω<sub>0</sub></td><td>SU</td><td>0.06</td></tr><tr><td>Effective fracturing ratio at V<sub>0</sub></td><td>T<sub>x0</sub></td><td>SU</td><td>5</td></tr><tr><td>Total fracturing ratio</td><td>T<sub>xf</sub></td><td>SU</td><td>40</td></tr><tr><td>Glass density</td><td>d<sub>v</sub></td><td>g m<sup>-3</sup></td><td>2.75 × 10<sup>6</sup></td></tr><tr><td>Mass fraction of Si in the glass</td><td>x<sub>Si</sub></td><td>g<sub>Si</sub> · g<sup>-1</sup><sub>glass</sub></td><td>0.21</td></tr><tr><td>Molar mass of Si</td><td>M<sub>Si</sub></td><td>g · mol<sup>-1</sup></td><td>28</td></tr></tbody></table>

The values given here are those used for the calculations in Chapter 3.

package (V₀). Later, the silicon removal rate becomes controlled by the diffusion of silicon through the mesh nodes whose reactivity has already been saturated. As the reactivity has ceased, it is then possible to use the simple diffusive transport equation through a limited spherical wall, in the permanent regime given in *Mathematics of diffusion*, page 89$^{26}$. This equation leads to Eq. 2 for the alteration rate controlled by diffusion through a spherical wall with internal and external radii [R₁, R₂] at concentrations within the limits C₁ and C₂.

$$V(t) = 4\pi D(C_2 - C_1) \left(\frac{R_1 R_2}{R_2 - R_1}\right) \frac{M_{\text{Si}}}{x_{\text{Si}}(1 - T_{x\text{Si}}^{\text{ret}})T_{x0}S_0} \quad (2)$$

This equation no longer includes the time, given that it concerns a permanent regime, nor the media reactivity $\mathcal{R}$, since it is exhausted, nor the porosity $\omega$. It is applied iteratively between the internal radius R₁ and an external radius R₂ that increases with the number of mesh nodes whose reactivity has been saturated. The concentrations at the limits C₁ and C₂ are respectively the concentration of silicon in contact with the package (C₁ = C<sub>sat</sub>) and the concentration in the distant medium (C₂ = C<sub>∞</sub>).

The approximation for a transitory regime using the equation for a permanent regime is classical in diffusive transport. In Cartesian 1D, it consists in approximating the concentration profile with a linear profile. For the diffusive flux, the approximation gives a systematic error that is a constant close to 1. In this calculation, however, a second approximation is added compared to the analytical equation. It concerns the choice of mesh node size. Overly large meshing would artificially maintain alteration throughout the time needed to saturate the entire volume of the mesh node before taking into account its diffusive braking. Inversely, overly-small meshing could lead to under-estimation of the diffusive flux, as the permanent regime equation corresponds to the weakest concentration gradient vertical to the package. An appropriate choice of mesh sizes is, therefore, important.

There is a special mesh size able to remove the same silicon flux as that coming from the package when glass is altering at the initial rate. This size is noted $\Delta L_0$, and depends on the radius. $\Delta L_0$ is the size of the reference mesh with inner radius R₁ and outer radius R₂ that verifies Eq. 3, where $\Delta L_0 = R_2 - R_1$:

$$V_0 = 4\pi D(C_2 - C_1) \left(\frac{R_1 R_2}{R_2 - R_1}\right) \frac{M_{\text{Si}}}{x_{\text{Si}}(1 - T_{x\text{Si}}^{\text{ret}})T_{x0}S_0} \quad (3)$$

Considering that for the current mesh, C₁ = 0 and C₂ = C<sub>sat</sub>, $\Delta L_0$ is given by Eq. 4:

$$\begin{align}
\Delta L_0 (R_1) &= 4\pi D C_{\text{sat}} \left( \frac{R_1 (R_1 + \Delta L_0)}{V_0} \right) \frac{M_{\text{Si}}}{x_{\text{Si}}(1 - T_{x\text{Si}}^{\text{ret}})T_{x0}S_0} \nonumber \\
&= \frac{\frac{4\pi D C_{\text{sat}} M_{\text{Si}}}{V_0 x_{x0} S_0 x_{Si}} (1 - T_{x\text{Si}}^{\text{ret}})^* R_1^2}{1 - \frac{4\pi D C_{\text{sat}} M_{\text{Si}}}{V_0 x_{x0} S_0 x_{Si}} (1 - T_{x\text{Si}}^{\text{ret}})^* R_1} \tag{4} \\
&\approx \frac{D C_{\text{sat}} M_{\text{Si}}}{V_0 x_{Si} (1 - T_{x\text{Si}}^{\text{ret}}) T_{x0}} \left( \frac{R_1}{R_0} \right)^2 \nonumber
\end{align}$$

$R_0$ is the radius of the glass block (Equation for a homogeneous environment).

The size $\Delta L_0(R_1)$ given by the approximate expression of Eq. 4 becomes smaller as the diffusion coefficient of the porous medium lowers, the concentration gradient becomes weaker, and the glass alteration speeds up.

Having a first meshing smaller than $\Delta L_0$ presents no interest: the medium is not thick enough to slow the glass alteration, which is controlled by the initial rate. Two successive meshings with sizes of $\Delta L_0$ enable the alteration rate to be slowed by a factor of 2 compared to the initial rate. With ten thousand meshes of this size, it is, therefore, possible to describe a rate drop of up to four orders of magnitude, i.e., the order of magnitude of the difference between the initial rate and the smallest final rate that may be proposed as input parameter in the MOS. The meshing adapts itself, as it presents the sizes of elements depending on the local physical parameters. The less diffusive media close to the package have finer meshing. An analytical demonstration of the resolution method is given in Supplementary Information Note 2.

It is necessary to calculate the brake on transport created by the different successive layers within which the reactivity has already stopped. In *Mathematics of diffusion*$^{26}$ we find the solution for a Cartesian 1D medium (Eq. 5), which gives the equivalent diffusion coefficient $D^{\text{eq}}$, for a layer of thickness $L^{\text{eq}}$ corresponding to the summation of the thicknesses of the layers having their individual diffusion coefficient.

$$\frac{L^{\text{eq}}}{D^{\text{eq}}} = \frac{L_1}{D_1} + \frac{L_2}{D_2} + \dots + \frac{L_n}{D_n} \quad (5)$$

Applied in spherical 1D and expressed under the iterative form implemented in the calculations, this equation becomes Eq. 6. It gives the equivalent diffusion coefficient for the diffusive medium with the same
<table><caption>**Table 4 | Parameters related to alteration in the vapor phase**</caption><thead><tr><th>Parameters for the arrival of water and of the vapor phase</th><th>Notation</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Water arrival time</td><td>t<sub>0</sub></td><td>an</td><td>...</td></tr><tr><td>Time 1</td><td>t<sub>1</sub></td><td>an</td><td>...</td></tr><tr><td>Fraction of package height submerged at time t<sub>1</sub></td><td>h<sub>1</sub></td><td>SU</td><td>...</td></tr><tr><td>Time 2</td><td>t<sub>2</sub></td><td>an</td><td>...</td></tr><tr><td>Fraction of package height submerged at time t<sub>2</sub></td><td>h<sub>2</sub></td><td>SU</td><td>...</td></tr><tr><td>Time at the end of saturation</td><td>t<sub>3</sub></td><td>an</td><td>...</td></tr><tr><td>Fraction of package height submerged at time s t<sub>3</sub></td><td>h<sub>3</sub></td><td>SU</td><td>1</td></tr></tbody></table>

thickness.

$$D_1^{eq} = D_1; D_n^{eq} = \frac{R_0^{-1} - R_n^{-1}}{\frac{R_0^{-1} - R_{n-1}^{-1}}{D_{n-1}^{SI}} + \frac{R_{n-1}^{-1} - R_n^{-1}}{D_n}} \quad (6)$$

This way, it is possible to solve the considered equations for a heterogeneous environment composed of several layers with different properties.

## Input parameters

This section presents all the input parameters. The package parameters are listed in Table 3.

The geometric surface of the cylindrical block is calculated using the radius and the height. This surface is multiplied by the maximum surface flux achievable in this geometry. This maximum surface flux is the one discharged in three spatial directions by a sphere with a radius equivalent to that of the cylindrical package. Multiplying the surface flux from a sphere by the total external surface of the cylinder is the overestimating choice made in the MOS. Everything happens as if we were altering three spherical packages. This assumption overestimates the reactivity of the environment. Such a significant margin could be reduced through modeling dedicated to the study of the G1 hypothesis.

It should be noted that the “fracturing ratio” is a parameter that does not only come from fracturing in the package. It also depends on the flow of water in the cracks and the initial alteration rate of the glass, and thus on the temperature. In fact, what is called the “fracturing ratio” is defined as the ratio between the quantity of altered glass measured and that which would be obtained if the block was monolithic. It corresponds to a magnitude measured during the alteration of scale one 400 kg glass block and reflects the brake on the transport of the matter that is within a fractured glass block. It is a function of the first-order glass alteration and of the ability of the glass block cracks to remove the matter coming from glass dissolution. When the solution is said to be saturated with respect to the passivating gel and the solution composition is no longer the first-order limitation on the alteration rate, all the surfaces contribute to the alteration in the same way. The fracturing ratio, which is also the ratio between the surface concerned and the geometric surface of the package, is then at its maximum level, equal to the block’s total fracturing ratio $T_{xf}$(SU), a factor of approximately 40. In contrast, for a block of glass submerged in a regularly renewed solution of pure water, the external surface alters at the maximal alteration rate for the glass, to which are added the contributions just from the cracks close to the external surface (in which the solution is sufficiently renewed). The internal cracks quickly constitute a very confined medium where the solution is loaded with elements from the glass alteration. Because of this, the alteration rate drops by several orders of magnitude, and the contribution from these internal surfaces to the total quantity of altered glass is low. The contribution from the surfaces that apparently alter at the initial rate (external surfaces) corresponds to approximately 4 times the geometric surface, giving an effective fracturing ratio of 5 in $V_0$. This figure is the result of the quantity of altered glass brought to the external geometric surface during the alteration

<table><caption>**Table 5 | Glass alteration parameters**</caption><thead><tr><th>Glass alteration parameters</th><th>Notation</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Temperature</td><td>T</td><td>°C</td><td>50</td></tr><tr><td>pH</td><td>pH</td><td>SU</td><td>8</td></tr><tr><td>Concentration at solution saturation</td><td>C<sub>sat</sub></td><td>mol<sub>Si</sub> · m<sup>-3</sup></td><td>1</td></tr><tr><td>Silicon retention ratio at V<sub>0</sub></td><td>T<sub>Si/Si</sub><sup>ref</sup></td><td>SU</td><td>0.2</td></tr><tr><td>Silicon retention ratio</td><td>T<sub>Si</sub><sup>ref</sup></td><td>SU</td><td>0.3</td></tr><tr><td>Dissolution constant for the initial rate</td><td>k<sub>disso</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>1.17 × 10<sup>8</sup></td></tr><tr><td>pH dependence</td><td>n<sub>disso</sub></td><td>SU</td><td>-0.4</td></tr><tr><td>Activation energy of the initial rate</td><td>Ea<sub>disso</sub></td><td>kJ · mol<sup>-1</sup></td><td>76</td></tr><tr><td>Initial rate (T, pH)</td><td>V<sub>0</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>Calculated = 0.474</td></tr><tr><td>Correction factor for the solution composition</td><td>R<sub>compsol</sub></td><td>SU</td><td>5</td></tr><tr><td>Dissolution constant for the final rate</td><td>k<sub>f</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>2.21 × 10<sup>4</sup></td></tr><tr><td>Activation energy of the final rate</td><td>Ea<sub>f</sub></td><td>kJ · mol<sup>-1</sup></td><td>53</td></tr><tr><td>Final rate (7)</td><td>V<sub>f</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>Calculated = 5.93 × 10<sup>-5</sup></td></tr><tr><td>Dissolution constant for the hydration rate</td><td>k<sub>h</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>1.75 × 10<sup>3</sup></td></tr><tr><td>Activation energy of the hydration rate</td><td>Ea<sub>h</sub></td><td>kJ · mol<sup>-1</sup></td><td>40</td></tr><tr><td>Glass alteration rate in the vapor phase</td><td>V<sub>hydr</sub></td><td>g · m<sup>-2</sup> · j<sup>-1</sup></td><td>Calculated = 5.94 × 10<sup>-4</sup></td></tr></tbody></table>

The values given here are those used for the calculation described in Chapter 3.

of a scale one block of glass in constantly-renewed pure water at 100 °C. This effective fracturing ratio at $V_0$ $T_{x0}$(SU) is the parameter applicable for an initial rate $V_0$ in glass alteration. The competition between reactive transport within the block and glass dissolution kinetics is the source of this fracturing ratio value

The parameters for the arrival of water and of the vapor phase are listed in Table 4. They are defined from the water arrival time, noted $t_0$, which concerns liquid water as well as water vapor. After $t_0$, a sectioned linear law is used to describe the evolution through time of the submerged height of the package, $H(t) = 2R_0h(t)$ (Fig. 9) where $h(t)$ is the fraction of submerged height (in the Supplementary Information Note 3, column W).

This gives:

$$h(t) = \frac{H(t)}{2R_0} = 0 \quad \begin{array}{l} si \ t \le t_0; \ h_1 \frac{t-t_0}{t_1-t_0} si \ t_0 \le t \le t_1; \ h_1 \\ + (h_2-h_1) \frac{t-t_1}{t_2-t_1} si \ t_1 \le t \le t_2; \ h_2 \\ + (h_3-h_2) \frac{t-t_2}{t_3-t_2} si \ t_2 \le t \le t_3; \ h_3 si \ t \ge t_3 \end{array} \quad (7)$$

The law $h(t)$ is therefore defined based on the parameters of time ($t_0$, $t_1$, $t_2$, $t_3$) and of the fraction of submerged height ($h_1$, $h_2$, $h_3$, corresponding respectively to the fractions submerged at the times

$t_1$, $t_2$, $t_3$). The package is completely submerged at time $t_3$, i.e., $h_3 = 1$.

During the saturation period between the times $t_0$ and $t_3$, the package is assumed to alter at the rate $V_{hydr}$ on the surfaces of the package above water, and at the rate $V(t)$ between $V_0$ and $V_f$ on the submerged surfaces. In other words, the total altered fraction of the package ($t$), sum of the altered glass
fractions in liquid water and under the vapor phase, is given by Eqs. 8 and 9:

$$ 
\begin{align*}
FTVA(t) &= FVA_{\text{vapor}}(t) + FVA_{\text{water}}(t) \\
FVA_{\text{vapor}}(t) &= \int_{t_0}^{t} V_{\text{hydr}} S_0 T_{xf} \frac{\Omega_{\text{Non-subm}}(h(t'))}{\pi H_0 R_0^2} dt' \\
FVA_{\text{water}}(t) &= \int_{t_0}^{t} V(t') S_0 T_{\times 0} \frac{\Sigma_{\text{subm}}(h(t'))}{2\pi(H_0 R_0 + R_0^2)} dt'
\end{align*}
$$

with:

$$ 
\begin{align}
\Omega_{\text{Non-subm}}(h) &= H_0 R_0^2 \left( \pi + 2(1-2h)(h-h^2)^{\frac{1}{2}} - \text{Arc cos}(1-2h) \right) \\
\Sigma_{\text{subm}}(h) &= 2(H_0 R_0 + R_0^2) \text{Arc cos}(1-2h) - 4R_0^2(1-2h)(h-h^2)^{\frac{1}{2}}
\end{align}
$$

$\Omega_{\text{Non-subm}}$ and $\Sigma_{\text{subm}}$ correspond respectively to the non-submerged internal volume (in Supplementary Information Note 3, column Y) and the submerged external surface of the package (Supplementary Information Note 3, column X), expressed in function of the fraction of submerged height $h$.

The glass alteration parameters are listed in Table 5.

The rates $V_0$, $V_f$ and $V_{\text{hydr}}$ are calculated depending on the temperature and the pH according to the following laws (Eqs. 10–12):

$$ V_0(T, \text{pH}) = k_{\text{disso}} \cdot 10^{-n_{\text{disso}}} \cdot \text{pH} \cdot \exp\left(-\frac{Ea_{\text{disso}}}{R \cdot T}\right) \cdot R_{\text{compsol}} \quad (10) $$

$$ V_f(T) = k_f \cdot \exp\left(-\frac{Ea_f}{R \cdot T}\right) \quad (11) $$

$$ V_{\text{hydr}}(T) = k_{\text{h}} \cdot \exp\left(-\frac{Ea_{\text{h}}}{R \cdot T}\right) \quad (12) $$

<figure><img src="image_13.png" /><figcaption>Fig. 9 | Diagram showing progressive saturation in liquid water.</figcaption></figure>

$R_{\text{compsol}}$ in Eq. 10 is the correction factor for the solution composition, equal to 1 in initially pure water and 5 in Bure porewater.

The *residual rate for glass alteration in a closed system* is called the *final rate*. This term was chosen to distinguish it from the *residual rate in the environment*, a term that has gained wider acceptance in recent publications.

The chemical and geometric parameters for the environment are listed in Table 6.

The name, internal and external radii for each of the zones, their porosity, and the silicon diffusion coefficient in the pores must be entered. Care must be taken not to define overly-thin size zones if their diffusion coefficient is very high: the size of each zone must be at least greater than $\Delta L_0$ (Cf. Equation 4, Resolution method for a heterogeneous environment) to ensure an accurate if somewhat over-estimated calculation (Cf. in the Supplementary Information Note 1 and Fig. 10).

The diffusion coefficient for silicon at temperature is calculated with the value given at 20 °C, using a polynomial law reproducing the Stokes-Einstein law between 0 and 90 °C (Eq. 13).

$$ D(T) = D(20 \,^{\circ}\text{C})(2.00 \times 10^{-4}(T - 20)^2 + 2.79 \times 10^{-2}(T - 20) + 1) \quad (13) $$

The reactivity coefficient $\mathfrak{R}$ expresses the way in which the reactivity of the medium maintains the difference in the chemical potential of silicon, accelerates its transport, and slows the progress of the dissolved silicon front in the media. $\mathfrak{R}$ is a dimensionless number. It is a partition coefficient between the solid phase and the liquid phase (Hypothesis R1). This formalism is useful because it leads to a simple analytical solution for the reactive transport equation (Eq. 1, Equation for a homogeneous environment). The physico-chemical reality of the interaction between silicon and the environmental materials is probably different depending on the material in the interaction, and potentially more complex than a partition coefficient.

For example, as concerns an iron silicate precipitation, the following approach is suggested: $\mathfrak{R}$ is defined as the ratio between the maximal concentration in the solid after the chemical reaction and the maximal concentration in solution that can be found in contact with the glass (Eq. 14). With this definition, hypothesis R1 must enable a slightly over-estimated modeling of the silicon and iron interaction. A digital application in the case of metallic iron allows calculation of the maximal possible value for $\mathfrak{R}$: the metallic iron has an iron concentration $C_{\text{Fe}}(\text{mol m}^{-3})$ of solid equal to its density divided by its molar mass. Thus, considering (i) $x_{\text{Si/Fe}} = 1$, (ii) a 10% porosity, and (iii) a silicon concentration at saturation close to 1 mol m$^{-3}$, (ivii) the fact that, due to their respective densities, iron transformed into iron silicate occupies at least twice the volume of metallic iron, we arrive at $\mathfrak{R}<10^6$ for iron. This means a maximum acceleration

<table><caption>Table 6 | Chemical and geometric parameters for the environment</caption><thead><tr><th>N°</th><th>Designation</th><th>Internal R</th><th>External R</th><th>Porosity</th><th>Diff Coef (20 °C)</th><th>Diff Coef (T)</th><th>Reactivity</th></tr></thead><tbody><tr><td></td><td>Unit</td><td>m</td><td>m</td><td>SU</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td><td>SU</td></tr><tr><td></td><td>Notation</td><td>R<sub>int</sub></td><td>R<sub>ext</sub></td><td>ω</td><td>D<sub>0</sub></td><td>D</td><td>γR</td></tr><tr><td>1</td><td>Container</td><td>0.21</td><td>0.28</td><td>0.1</td><td>10<sup>-12</sup> at 5 × 10<sup>-9</sup></td><td>Calculated 2 × 10<sup>-12</sup> at 10<sup>-8</sup></td><td>0–10<sup>5</sup></td></tr><tr><td>2</td><td>Infill zone</td><td></td><td></td><td></td><td></td><td>Calculated</td><td>0–10<sup>5</sup></td></tr><tr><td>3</td><td>Overpack</td><td></td><td></td><td></td><td></td><td>Calculated</td><td>0–10<sup>5</sup></td></tr><tr><td>4</td><td>Infill zone</td><td></td><td></td><td></td><td></td><td>Calculated</td><td>0–10<sup>5</sup></td></tr><tr><td>5</td><td>Tube guide</td><td></td><td></td><td></td><td></td><td>Calculated</td><td>0–10<sup>5</sup></td></tr><tr><td>6</td><td>Mortar (MREA)</td><td>0.28</td><td>0.32</td><td>0.1</td><td>10<sup>-11</sup></td><td>Calculated 2 × 10<sup>-11</sup></td><td>0</td></tr><tr><td>7</td><td>Oxidized clay</td><td>0.32</td><td>0.52</td><td>0.1</td><td>10<sup>-11</sup></td><td>Calculated 2 × 10<sup>-11</sup></td><td>0</td></tr><tr><td>8</td><td>Clay</td><td>0.52</td><td>1</td><td>0.1</td><td>10<sup>-11</sup></td><td>Calculated 2 × 10<sup>-11</sup></td><td>0</td></tr><tr><td>9</td><td>...</td><td>1</td><td>2</td><td>0.1</td><td>10<sup>-11</sup></td><td>Calculated 2 × 10<sup>-11</sup></td><td>0</td></tr><tr><td>10</td><td>...</td><td>2</td><td>3</td><td>0.1</td><td>10<sup>-11</sup></td><td>Calculated 2 × 10<sup>-11</sup></td><td>0</td></tr></tbody></table>

The values given here are those used for the calculation described in Chapter 3, and the three metallic layers are associated with a single layer.
<table><thead><tr><th colspan="3">Glass Package parameters</th></tr></thead><tbody><tr><td>Radius of the glass block</td><td>Ro</td><td>m</td><td>0.21</td></tr><tr><td>Height of the block</td><td>Ho</td><td>m</td><td>1.1</td></tr><tr><td>Actual geometric surface of the block</td><td>So</td><td>m2</td><td>1.73</td></tr><tr><td>Internal porosity of the block</td><td>w0</td><td>SU</td><td>0.06</td></tr><tr><td>Effective Fracturing Ratio</td><td>Txo</td><td>SU</td><td>5</td></tr><tr><td>Total Fracturing Ratio</td><td>Txf</td><td>SU</td><td>40</td></tr><tr><td>Glass density</td><td>dv</td><td>g/m3</td><td>2750000</td></tr><tr><td>Mass fraction Si</td><td>xSi</td><td>gSi/gverre</td><td>0.21</td></tr><tr><td>Molar mass Si</td><td>MSi</td><td>g/mol</td><td>28</td></tr></tbody></table><table><thead><tr><th colspan="10">Calculation DV</th></tr><tr><th colspan="2">% VA recherché</th><th colspan="2">10 % de VA</th><th colspan="2">Recherche du temps associé à un %VA</th></tr><tr><th colspan="2">Result →</th><th colspan="2">7.9 % AG</th><th colspan="2">21 thousands of years</th></tr><tr><th colspan="2">Water/steam distribution</th><th colspan="2">6.0 % in water</th><th colspan="2">1.9 % in vapor phase</th></tr><tr><th colspan="5">Paramètres de l'environnement</th><th colspan="2">Diff(20 °C Diff (T))</th><th colspan="2">Reactivity</th></tr><tr><th>N°</th><th>Nom</th><th>Internal Ra</th><th>External Ra</th><th>Porosity</th><th>m/m</th><th>m/m²</th><th>SU</th><th>SU</th></tr><tr><th></th><th></th><th>Ri</th><th>Re</th><th>SU</th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Stainless steel containt</td><td>0.21</td><td>0.23</td><td>w</td><td>0.1</td><td>Do 20°C</td><td>μ/p*(g/T)</td><td>R</td></tr><tr><td>2</td><td>Non-alloy steel overp</td><td>0.23</td><td>0.31</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>3</td><td>MREA</td><td>0.31</td><td>0.4</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>4</td><td>Modified Clay</td><td>0.4</td><td>0.5</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>5</td><td>Clay</td><td>0.5</td><td>0.6</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>6</td><td>Clay</td><td>0.6</td><td>0.7</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>7</td><td>Clay</td><td>0.7</td><td>0.8</td><td>0.1</td><td>1.0E-11</td><td>1.3E-11</td><td>1E+00</td><td></td></tr><tr><td>8</td><td>Clay</td><td>0.8</td><td>0.9</td><td>0.1</td><td>1.0E-12</td><td>1.3E-12</td><td>1E+00</td><td></td></tr><tr><td>9</td><td>Clay</td><td>0.9</td><td>4</td><td>0.1</td><td>1.0E-12</td><td>1.3E-12</td><td>1E+00</td><td></td></tr><tr><td>10</td><td>Clay</td><td>4</td><td>5</td><td>0.1</td><td>1.0E-12</td><td>1.3E-12</td><td>1E+00</td><td></td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>Fig. 10 | “MOS” tab on the calculation sheet. This tab shows the input parameters to be entered and the calculation results: evolution through time of the alteration rates, the quantities of glass altered, and the position of the diffusion front. The digital values displayed are examples, only intended to illustrate the calculation.</figcaption></figure>

factor of 10³ for the alteration rate when it is controlled by reactive diffusion into the environment. The same factor applies throughout the initial rate phase. For magnetite transforming into iron silicate, and with the same reasoning, $\Re$ is $<6 \times 10^5$.

$$ \mathfrak{R} = \frac{(1 - \omega) C_{\text{Fe}} x_{\text{Si}}}{\omega C_{\text{sat}}} \quad (14) $$

The same partition calculation between silicon in the solid phase and in the liquid phase can be applied for reactivity with the clay. For example, for clay that has a silicon consumption capacity per g of clay $CC = 0.1 \text{ mmol g}^{-1}$, density $\rho_A = 2 \text{ kg l}^{-1}$, and porosity $\omega = 0.1$, the digital application (Eq. 15) with $C_{\text{sat}} = 1 \text{ mmol l}^{-1}$ gives $\mathfrak{R} = 1800$.

$$ \mathfrak{R} = \frac{(1 - \omega)CC\rho_A}{\omega C_{\text{sat}}} \quad (15) $$

In the absence of chemical interaction, $\mathfrak{R}$ has a null value.

## Implementing the calculation

Despite the number of parameters, equations, and lines in the calculation sheet, implementing the calculation is simple as the resolution method enables each cell to be calculable by an equation. The calculation principle is explained here:

* Each calculation line represents a mesh node.

* The first columns are used to read the characteristics of the mesh in question: number, internal radius, and then the number, the external radius, the porosity, and the reactivity of the zone the mesh node belongs to. Logic tests using embedded “if” functions are used to identify these different magnitudes.

* The reference size $\Delta L_0$ is calculated (Eq. 4, Resolution method for a heterogeneous environment). This site is specific to each mesh node, depending on its position and the physico-chemical properties of the zone. It gives the thickness of the node in question and is used to calculate the internal radius of the following mesh node.

* Knowing the current diffusion coefficient in the mesh node and its size, it is possible to deduce the value of the equivalent diffusion coefficient through the medium consisting of the current node and of all the

previous ones (Eq. 6, Resolution method for a heterogeneous environment).

* Knowing the equivalent diffusion coefficient for all the previous nodes, the flow of silicon entering the mesh node can be found (Eq. 2, Resolution method for a heterogeneous environment).

* The volume of the spherical corona corresponding to the current mesh node is calculated. From this, the quantity of silicon retained in the corona can be deduced, both in the solution and the solids. Given the spherical geometry, this calculation slightly overestimates the quantities of solution and of reactive minerals effectively contained in the package periphery. This is hypothesis G1.

* From the flow of silicon entering the mesh node and the quantity contained in the corona, the time necessary to fill the current node can be calculated. This indexes the evolution of time in the calculation sheet.

* From the flow of silicon entering the mesh node, the glass alteration rate is deduced taking into account the silicon content in the glass and the silicon retention ratio in the alteration products. This rate is normalized at the geometric surface of the package. It is the alteration rate in water under the effect of the reactive diffusion of silicon. This is compared to the initial rate ($V_0$), to the final rate ($V_f$), and to the rate in vapor phase ($V_{\text{hydr}}$). These rates are each normalized at the geometric surface of the package in order to compare them, as, in fact, they apply to different surfaces: effective surface in initial rate ($S_{\text{acc}}=S_0 T_{x0}$), total surface in final rate, and vapor phase ($S_{\text{tot}}=S_0 T_{xf}$).

* The quantities of altered glass are calculated, and accumulate in time. They are deduced from the rates and from the fractions of the package surface to which they apply: the external surface for a submerged package for the rate under water, or the internal above-water surface for the vapor phase. The sum of the quantities altered under water and under the vapor phase gives a total quantity of altered glass. Knowing whether the matter altered during the vapor phase is available to be transported into the environment is a question that is not dealt with here.

The Supplementary Information Note 3 gives a complete presentation of the equations implemented in the Calculation tab of the MOS sheet, line by line and column by column. The designations and abbreviations used for
<figure><img src="image_15.png" /><figcaption>**Fig. 11 | “View water arrival” tab on the calculation sheet.** In this example, water arrives after 1000 years in the form of vapor. Liquid water completely and suddenly replaces the vapor phase after 5000 years. The digital values displayed are examples, only intended to illustrate the calculation.</figcaption></figure>

the variables in the Supplementary Information Note 3 and in the calculation sheet tabs are identical. The MOS model is resolved in a calculation sheet. The model is also resolved in a dedicated program written in C language, in order to simplify studies of uncertainty propagation.

Figure 10 shows the “MOS” tab on the calculation sheet. The upper half presents the input parameters, while the lower half gives graphs of the results. The graph on the left shows the rates under water through time ($V(t)$), the hydration rate ($V_{hydr}$), and the final rate ($V_f$). All are normalized at the geometric surface of the package in order to be able to compare them. The left graph also shows the advance of the diffusion front. The graph on the right gives the quantities of glass altered under water and in the vapor phase, as well as the sum of the two contributions. Logarithmic scales are used on all the axes, given the amplitude of the variations. The graphs and, in particular, the choice of their scales are updated automatically to enable the preservation of a correct view of the results as soon as the value of a parameter is changed. The initial time of the graphical representation is the moment water arrives at the package (in liquid or in vapor phase).

In the top center of the sheet (pinkish blocks) is a summary of the calculation in figures. It is possible to select a quantity of altered glass in the cell «FVA Search». FVA stands for Fraction de Verre Altéré, the altered glass fraction. The spreadsheet searches among the thousands of calculation lines for the closest value below or equal to this. It gives the total fraction value of altered glass (FTVA(%)) and the time at which it was reached. Unlike graph representations, this time takes into account the time $t_0$ passed before the water's arrival. It also gives the FVA under water (FVA\_water) and in the vapor phase (FVA\_vapor) at that time.

Figure 11 shows the tab “View water arrival” on the calculation sheet. It enables a display of the function describing the duration of the vapor phase depending on the values chosen for the times $t_1, t_2, t_3$ and the intermediary fractions of submerged height $h_1$ and $h_2$, knowing that the final fraction for the submerged height is $h_3$ equals 1.

## Conclusions

The MOS model originated in the desire to have a simple tool that explicitly integrated the contribution from the reactive diffusive environment to the alteration of a nuclear waste package in deep geological repository conditions.

The model enables a quantification of the environment's contribution to glass alteration. The structure of the model has been presented. At this stage, its predictive use should remain limited until the key initial steps of validation through reactive transport code coupled with glass alteration model, have been accomplished. Nevertheless, the model already provides a framework for considering the influence of each parameter over the relevant scales of time and space on the predominance of alteration mechanisms. It, therefore, helps anticipate the consequences of disposal concept choices. Finally, the model provides a way to quantify the extent of margins taken by simpler safety models. The MOS remains a simple and envelope model. It is not intended to replace best estimate models (reactive transport code + glass alteration model). It improves the connection between safety models and scientific models.

Although built for a specific purpose, the proposed model holds general significance for assessing the alteration of minerals and materials in their environment. Its uniqueness lies in the reversal of the classical paradigm, wherein the focus is on evaluating the long-term evolution kinetics of a mineral by considering its intrinsic characteristics rather than accounting for the reactivity and transport properties of its environment. The equation, derived from a simple mass balance, is formalized and can be easily adapted to other minerals. It remains straightforward in the case of a single layer with constant transport and reactivity properties.

The original resolution method for reactive transport in the case of multiple layers has been described in detail. It could be useful for other corrosion issues when it is a question of whether or not a diffusive medium can remove the matter in solution by a surface reaction.

## Data availability

The authors declare that the data supporting the findings of this study are available within the article and its supplementary information files.

## Code availability

The code generated during the current study is not publicly available due to confidentiality agreements with third parties. However, it is available from the corresponding author on reasonable request.

Received: 29 September 2023; Accepted: 21 July 2024;

Published online: 02 September 2024

## References

1. ANDRA. Dossier d'Autorisation de création de l'installation nucléaire de base (INB) CIGEO - Pièce 0 Présentation non technique. 36 (ANDRA CG-TE-PRE-AMOA-PUO-0000-21-0026/A, 2022).

2. ANDRA. Dossier d'options de sûreté - Partie après fermeture (DOS-AF). Report No. CG-TE-D-NTE-AMOA-SR2-0000-15-0062, 1-467 (ANDRA, 2016).

3. Ribet, I. Dossier Opérationnel verre. Modèles opérationnels et calculs d'intégration. Report No. RT DTCD/2004/04, (2004).

4. Ribet, I. et al. in *Advances for future nuclear fuel cycles Nimes*. 1-8 (2004).

5. Minet, Y. & Godon, N. in *American Ceramic Society-104th meeting*. 275-282 (2003).

6. Frugier, P. et al. SON68 Nuclear glass dissolution kinetics: current state of knowledge and basis of the new GRAAL model. *J. Nucl. Mater.* **380**, 8-21 (2008).

7. Minet, Y., Bonin, B., Gin, S. & Frugier, P. Analytic implementation of the GRAAL model: application to a R7T7-type glass package in a geological disposal environment. *J. Nucl. Mater.* **404**, 178-202 (2010).

8. Frugier, P., Minet, Y., Rajmohan, N., Godon, N. & Gin, S. Modeling glass corrosion with GRAAL. *npj Mater. Degrad.* **2**, 35 (2018).

9. Frugier, P., Chave, T., Gin, S. & Lartigue, J. E. Application of the GRAAL model to leaching experiments with SON68 nuclear glass in initially pure water. *J. Nucl. Mater.* **392**, 552-567 (2009).
10. Rajmohan, N., Frugier, P. & Gin, S. Composition effects on synthetic glass alteration mechanisms: Part 1. Experiments. *Chem. Geol.* **279**, 106–119 (2010).

11. Verney-Carron, A., Gin, S., Frugier, P. & Libourel, G. Long-term modeling of alteration-transport coupling: application to a fractured Roman glass. *Geochim. Cosmochim. Acta* **74**, 2291–2315 (2010).

12. Debure, M., Frugier, P., De Windt, L. & Gin, S. Borosilicate glass alteration driven by magnesium carbonates. *J. Nucl. Mater.* **420**, 347–361 (2012).

13. Jollivet, P., et al. Effect of clayey groundwater on the dissolution rate of the simulated nuclear waste glass SON68. *J. Nucl. Mater.* **420**, 508–518 (2012).

14. Debure, M., De Windt, L., Frugier, P. & Gin, S. HLW glass dissolution in the presence of magnesium carbonate: Diffusion cell experiment and coupled modeling of diffusion and geochemical interactions. *J. Nucl. Mater.* **443**, 507–521 (2013).

15. Godon, N., Gin, S., Rebiscoul, D. & Frugier, P. in *Water Rock Interaction WRI 14*. (2013).

16. Frugier, P., Fournier, M. & Gin, S. Modeling resumption of glass alteration due to zeolites precipitation. *Proc. Earth Planet. Sci.* **17**, 340–343 (2017).

17. IRSN. Projet de stockage Cigéo—Examen du Dossier d'Options de Sûreté. 243 (IRSN, 2017).

18. Gin, S., Ribet, I. & Couillard, M. Role and properties of the gel formed during nuclear glass alteration: importance of gel formation conditions. *J. Nucl. Mater.* **298**, 1–10 (2001).

19. Collin, M. et al. Structure of International Simple Glass and properties of passivating layer formed in circumneutral pH conditions. *npj Mater. Degrad.* **2**, 4 (2018).

20. Gin, S. et al. Origin and consequences of silicate glass passivation by surface layers. *Nat. Commun.* **6**, 6360 (2015).

21. Jollivet, P., et al. Investigation of gel porosity clogging during glass leaching. *J. Non-Crystalline Solids* **354**, 4952–4958 (2008).

22. Van der Lee, J., De Windt, L., Lagneau, V. & Goblet, P. Module-oriented modeling of reactive transport with HYTEC. *Comput. Geosci.* **29**, 265–275 (2003).

23. Gin, S. & Frugier, P. in *MRS 2002*. 175–182 (2003).

24. Neeway, J., Abdelouas, A., Grambow, B. & Schumacher, S. Dissolution mechanism of the SON68 reference nuclear waste glass: New data in dynamic system in silica saturation conditions. *J. Nucl. Mater.* **415**, 31–37 (2011).

25. Blanchard, J. P., Damblin, G., Martinez, J. M., Arnaud, G. & Gaudier, F. The Uranie platform: an open-source software for optimisation, meta-modelling and uncertainty analysis. *EPJ Nucl. Sci. Technol.* **5**, 4 (2019).

26. Crank, J. *The mathematics of diffusion*. 414 (1975).

27. ANDRA. Rapport Andra CG-TE-D-NTE-AMOA-SR1-0000-15-0060 – « Dossier d'options de sûreté - Partie exploitation ». (2016).

## Acknowledgements

The authors wish to thank the organizers of the Sernhac Seminar and its participants, a working group that lasted several days and from which this model was born: Bernard BONIN, Maxime FOURNIER, Pierre FRUGIER, Nicole GODON, Yves MINET, Sylvain PEUGET, Isabelle RIBET, Magaly TRIBET. This work was made possible by funding from CEA, ORANO, and EDF, and we thank them for their support.

## Author contributions

Pierre FRUGIER: conceptualization, equation resolutions, implementation, report writing; Nicole GODON: coordination, exploitation of the model and sensitivity studies, writing and proofreading; Yves MINET: checking equations and their resolution, writing and proofreading.

## Competing interests

The authors declare no competing interest.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41529-024-00496-0.

**Correspondence** and requests for materials should be addressed to Pierre Frugier.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024