# Advanced Monte Carlo method for simulating glass alteration: Application to aluminoborosilicate glasses

Jean-Marc Delaye, Sumit Tiwari, Emeric Brun, Paul Fossati, Stephane Gin

► To cite this version:

Jean-Marc Delaye, Sumit Tiwari, Emeric Brun, Paul Fossati, Stephane Gin. Advanced Monte Carlo method for simulating glass alteration: Application to aluminoborosilicate glasses. Journal of the American Ceramic Society, 2024, pp.1-19. 10.1111/jace.20167. cea-04738813

HAL Id: cea-04738813

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

Submitted on 15 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_1.png" /></figure>

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

# Advanced Monte Carlo method for simulating glass alteration: Application to aluminoborosilicate glasses

Jean-Marc Delaye<fnref n="1" /> ID
Stéphane Gin<fnref n="1" /> ID

Sumit Tiwari<fnref n="1" />

Emeric Brun<fnref n="2" />

Paul C. M. Fossati<fnref n="3" /> ID

## Correspondence

Jean-Marc Delaye, CEA, DES, ISEC, DPME, SEME, University of Montpellier, Marcoule, F-30207 Bagnols-sur-Cèze, France.

Email: jean-marc.delaye@cea.fr

## Funding information

CEA

## Abstract

Understanding the alteration mechanisms of nuclear glasses is paramount for safety assessment of deep geological disposal of nuclear waste. This endeavor poses significant challenges due to the experimental limitations in observing atomic-scale mechanisms, which are essential for developing macroscopic alteration models. Consequently, a novel Monte Carlo approach has been developed to simulate the alteration process of simplified nuclear glasses at the atomic level. This approach employs a dual lattice to represent both the solid material and the solution, thereby facilitating the incorporation of diffusion and hydrolysis-recondensation mechanisms. Preliminary findings underscore the critical influence of alteration conditions on glass behavior and elucidate the interplay between rate of hydrolysis at the gel-solution interface and the morphology of the outer layer formed within the gel. Intense hydrolysis promotes the formation of a porous, low-density layer, whereas slower hydrolysis fosters the development of a dense, cross-linked layer.

## KEYWORDS

glass alteration, Monte Carlo, passivating layer

# 1 | INTRODUCTION

In France, after the reprocessing of used fuel, nonrecyclable nuclear wastes are confined in borosilicate glasses, which are themselves intended to be disposed of in a deep geological repository. After several thousand years, groundwater will come into contact with the glass packages, which will then be altered, releasing the soluble radionuclides into the environment.<fnref n="1" />⁻⁴ Hence, it is crucial to identify the mechanisms controlling glass alteration. This enables the extrapolation of dissolution rates over

extended periods. Models are derived from this knowledge to calculate the source term and assess the safety of the disposal.⁵

Many studies have been conducted on aqueous corrosion of nuclear glasses over the previous decades. The references<sup>2,6–8</sup> are several recent reviews providing an overview of the current state of knowledge on the subject. It appears that not all aspects have been conclusively elucidated. For instance, the mechanisms governing glass passivation have not been definitively identified. It is likely that different mechanisms may be limited depending

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *Journal of the American Ceramic Society* published by Wiley Periodicals LLC on behalf of American Ceramic Society.

<footnote n="1">¹CEA, DES, ISEC, DPME, SEME, University of Montpellier, Bagnols-sur-Cèze, France</footnote>

<footnote n="2">²CEA, DES, ISEC, DMRC, University of Montpellier, Bagnols-sur-Cèze, France</footnote>

<footnote n="3">³CEA, Service de recherche en Corrosion et Comportement des Matériaux, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>
on the glass composition, alteration conditions (temperature, pH, and solution composition), and the progress of the reaction.<sup>9–14</sup> Many experimental approaches have been implemented to study the alteration behavior comprehensively, to analyze specific issues such as the morphology of alteration layers (electron microscopy),<sup>15–18</sup> their structures (Raman and Infrared spectroscopies, X-ray diffraction, and Nuclear magnetic resonance),<sup>19–23</sup> their porosity (small angle X-ray scattering, spectroscopic ellipsometry, and X-ray reflectometry),<sup>24,25</sup> concentration profiles between the solution and the pristine glass (time-of-flight secondary ions mass spectrometry),<sup>26,27</sup> and the amounts of aqueous species contained in the gels (Infrared spectroscopy and sum frequency generation).<sup>28</sup> However, experimental approaches are limited by the difficulty of investigating amorphous materials and dissociating elementary processes. Accessing information at a small scale focusing on an individual process remains a challenge.

Therefore, in parallel, modeling work has also been
performed to explore glass alteration mechanisms at the
atomic scale. Ab initio methods have been employed to
describe hydrolysis mechanisms and estimate correspond-
ing energies,<sup>29–33</sup> whereas classical molecular dynamics
or ReaxFF methods have been used to simulate model
gels,<sup>34–39</sup> describe water diffusion within them,<sup>40–43</sup> and
analyze their maturation over time. However, because of
the computational time, these methods cannot be directly
applied to simulate the behavior of glass over long periods,
and simulating complex glasses such as nuclear glasses
remains out of reach. Stochastic lattice models relying on
Monte Carlo method have been developed to bridge the
gap between highly precise yet temporally and spatially
constrained atomistic simulations and macroscopic reality.

The first Monte Carlo code (developed by Aertsens<sup>44</sup>) successfully replicated the effect of Na<sub>2</sub>O on the alteration behavior of a binary glass (SiO<sub>2</sub>-Na<sub>2</sub>O). Building on this initial attempt, a new approach was devised.<sup>45-47</sup> The alteration of simplified glasses (SiO<sub>2</sub>-B<sub>2</sub>O<sub>3</sub>-Na<sub>2</sub>O) was simulated by creating an interface between a glass, represented by a cubic lattice, in contact with a solution, localized on the same lattice. The mechanisms underlying the transition between initial and residual alteration rates have thus begun to be understood.

Then a more comprehensive algorithm was developed
at Commissariat à l'énergie atomique et aux énergies
alternatives (CEA, France) and later adapted at Pacific
Northwest National Laboratory (PNNL, USA) to reproduce
the behavior upon alteration of glasses containing sev-
eral oxides among the following: SiO<sub>2</sub>-Al<sub>2</sub>O<sub>3</sub>-B<sub>2</sub>O<sub>3</sub>-ZrO<sub>2</sub>-
Na<sub>2</sub>O-CaO.<sup>45,48</sup> These approaches notably deciphered the
role of Zr on the kinetics of formation of a dense film
between the pristine glass and the solution. An excel-
lent agreement was achieved between experimentation

and simulation regarding the morphology of gels.⁴⁹ More
recently, the same approach was utilized and digitized
into an independent code by PNNL to complement the
studies.⁵⁰⁻⁵² The effect of increasing Al₂O₃ content on the
alteration behavior was detailed as well as the impact of
the presence of boroxol rings. For low concentrations of
Al₂O₃, it is seen that an increasing Al₂O₃ content slows
down the release of B atoms, but ultimately, the amount
of B in solution increases because the protective film takes
longer to form. For high concentrations of Al₂O₃, the
decrease in the rate of B release is accompanied by a
decrease in its final concentration in solution as the glass
starts to dissolve into clusters and the protective film forms
more rapidly. Regarding boroxol rings, the clustering of B
leads to an increase in the glass dissolution rate.

In parallel, various models of silicon diffusion in the
pore water of the gel were compared to understand the
impact of diffusion of aqueous species in the gel on glass
alteration: homogenization of the solution, application of
a constant concentration gradient, or resolution of Fick's
law. Slower diffusion of silicon in the pore water of the
gel delays the formation of the external protective layer
enriched in silicon.⁵²

It is important to specify that all Monte Carlo
approaches developed so far for modeling glass alter-
ation rely on a unique and ordered lattice to model both
the glass and the solution. Moreover, only mechanisms
occurring at the glass-solution interface can be imple-
mented, namely, the possibility for glass atoms to be
released into solution and for the least soluble dissolved
atoms (especially Si and Al) to be redeposited at the
glass-solution interface. To our knowledge, the sole
endeavor to model the glass using a disordered lattice was
undertaken by Kerisit and coworkers.⁵³ This effort aimed
to simulate water diffusion into the vitreous network and
its hydrolysis. However, the challenge of modeling recon-
densation mechanisms on a disordered lattice precluded
the formation of a protective gel. To address this issue, an
approach using “patchy particles” was proposed.⁵⁴ The
principle involves representing the atomic system with
spherical particles interacting with short-range potentials
to accelerate relaxations. An application example was pro-
posed for hydrated silica. This approach could enable the
consideration of relaxation in Monte Carlo calculations
performed on genuinely disordered lattices.

More recently, a collaboration between PNNL and CEA
has addressed the question of the effect of glass irradiation
on the alteration rate.<sup>55,56</sup> Experimentally, it was observed
that pre-irradiation by heavy ions, which deposit elastic
energy within the glass, leads to a significant increase in
both the initial dissolution rate and the residual rate. All
attempts to reproduce these results using the Monte Carlo
approach, by adjusting the probabilities of mechanisms or
introducing pseudo-disorder in the form of Gaussian probability distributions, have failed. It has been concluded that the absence of accounting for water diffusion within the glass, independent of hydrolysis-recondensation mechanisms, which are currently impossible to implement in the current state of the algorithm, could be the cause of this failure because an acceleration of diffusion in an irradiated glass, which could not be simulated, was highly suspected.⁵⁵ This is why a new version of the Monte Carlo method was developed to introduce the possibility for aqueous species to diffuse within the solid without requiring, as it was previously the case, the release of species into solution beforehand. We have also transitioned to modeling glass hydrolysis bond by bond instead of atom by atom. As a result, the contact surface between the glass and the solution has increased, which we anticipate will facilitate the representation of radiation effects.

This article aims to describe this new Monte Carlo approach and the experiments used for parameter adjustments. The initial results aimed at reproducing alteration experiments on SiO₂–Al₂O₃–B₂O₃–Na₂O type glasses. This system incorporates the main elements of nuclear glass and already exhibits complex behaviors depending on the Al₂O₃ or B₂O₃ contents.⁵⁷⁻⁶⁰ We will explore how the alteration conditions impact the nature of the interface formed between the gel and the solution. The subsequent discussion will revisit the links between the rate of hydrolysis at the gel-solution interface and the ability to form a dense, cross-linked external layer.

## 2 | METHODOLOGY

### 2.1 | Experiments

The Monte Carlo calculations described in this article were adjusted based on alteration experiments conducted on three glasses with SiO₂–B₂O₃–Al₂O₃–Na₂O composition. The results of these experiments have been detailed elsewhere.⁵⁹ Here, we summarize only the “long-term” experiments, as they served as references for determining the probabilities of the Monte Carlo mechanisms.

The three sodium aluminoborosilicate glasses were prepared from raw materials (SiO₂, H₃BO₃, Na₂CO₃, and Al₂O₃) in a Pt-Rh crucible. The powders were melted at 1450°C for 3 h and 20 min then poured onto a metal plate. After grinding, the glass was reheated to 1450°C and transferred to a graphite crucible for refining at 550°C for 1 h, followed by cooling to room temperature.

The obtained glasses were analyzed by SEM-EDX for Si, Al, and Na elements, and by ICP-OES for all elements. The nominal compositions are given in Table 1. Chemical analyses confirmed the nominal compositions within ±5%.

<table><caption>TABLE 1 Nominal compositions (%mol) of the glasses used for the alteration experiments.</caption><thead><tr><th></th><th>SiO<sub>2</sub></th><th>B<sub>2</sub>O<sub>3</sub></th><th>Na<sub>2</sub>O</th><th>Al<sub>2</sub>O<sub>3</sub></th></tr></thead><tbody><tr><td>SBNA3.5</td><td>60.5</td><td>20.1</td><td>15.9</td><td>3.5</td></tr><tr><td>SBNA4</td><td>64.9</td><td>17.3</td><td>13.7</td><td>4.1</td></tr><tr><td>SBNA6</td><td>66.8</td><td>15.9</td><td>11.3</td><td>6</td></tr></tbody></table>

*Note:* S, B, N, and A stand for Si, B, Na, and Al, respectively. The following numbers represent the mol% of Al₂O₃ in the corresponding glass.

Rectangular parallelepiped coupons (20 × 20 × 1 mm³) were then prepared. All six faces were polished using various polishing grades up to a finishing with a diamond paste (1 µm) to ensure a roughness of less than 5 nm. Powder sizes of 120–125 µm and 40–100 µm were prepared using a standard procedure: grinding, sieving, and ultrasonic cleaning in ethanol and acetone to remove fine particles.

Alteration experiments were then conducted for 6 months. Leachings were carried out in PFA reactors. Glass powders were immersed in 50 mL of water (18.2 MΩ), corresponding to a surface-to-volume ratio (S/V) of 50 cm⁻¹. The pH was maintained at 9.0 ± 0.1 using a LiOH solution (10⁻² M). The temperature was set at 90°C. The high S/V ratio allowed for a rapid saturation of the silica solution. To limit evaporation, the reactors were placed in an over-container made of PFA containing a few mL of deionized water. Solution samples of 0.5 mL were regularly taken using a syringe. They were filtered (0.2 µm) and diluted in 2 mL of ultrapure HNO₃ (0.5 N). Samples were taken every 7 days for the first month and every 14 days for the following 5 months. The samples were analyzed by ICP-OES to measure the concentrations of Si, B, Na, and Al (uncertainties of ±5%).

Normalized mass losses (NL(i)) and equivalent alteration thicknesses (ETh(i)) were calculated using the following classical formulas:

$$NL(i)_t = NL(i)_{t-1} + \frac{[C(i)_t - C(i)_{t-1}] \cdot V_t}{S \cdot x_i} \quad (1)$$

with C(i)ₜ being the concentration of element i in solution at time t, Vₜ the volume of the solution at time t, S the surface area of the glass determined by BET, and xᵢ the fraction of element i in the glass. Normalized mass losses are expressed in g/m² as follows:

$$ETh(i) = \frac{NL(i)}{\rho} \quad (2)$$

with ρ being the density of the glass.

The uncertainties on the equivalent thicknesses are around 10%, taking into account the uncertainty in the measurements of the solution concentrations and that of the surface-to-volume ratio.
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Illustration of the initial network utilized to initiate a Monte Carlo calculation. The solid network is depicted on the right-hand side, consisting of Si (yellow), Al (green), and B (blue) atoms. An instance of the insertion of two [3]B atoms as a pair is demonstrated. On the left-hand side, the water molecules (red) are situated on the liquid network, positioned at the midpoint between two nodes on the “solid network.”</figcaption></figure>

**2.2 | Description of the Monte Carlo method**

The novelty of this Monte Carlo algorithm lies in its utilization of two distinct and interlinked networks, one representing the solid and the other the solution.

Regarding the network representing the solid, a cubic lattice, similar to those utilized in previous Monte Carlo methods described in the literature, is employed. Each node corresponds to a vertex of a cube, with connections to six neighboring nodes. To accurately model the tetrahedral coordination of elements found in glass structures (such as Si, Al, and tetrahedrally coordinated B = [4]B), certain bonds between neighboring nodes are selectively removed.

When a glass composition of the SiO₂–Al₂O₃–B₂O₃–Na₂O type is chosen, the requisite quantities of Si, Al, and B atoms are randomly distributed across the nodes of the lattice. Depending on the concentration of Na₂O, B atoms may exhibit either tri- or tetrahedral coordination. In-line with experimental findings, it is assumed that Na atoms preferentially compensate for Al atoms (as %Na₂O > %Al₂O₃, all the Al atoms are tetra-coordinated in the glasses used here), with any remaining Na atoms stabilizing tetrahedrally coordinated B atoms. Any uncompensated B atoms adopt tri-coordination. To incorporate

tri-coordinated B atoms ([3]B), the following procedure is employed: two [3]B atoms are placed adjacently on neighboring sites, and the bond between these nodes is then removed. This procedure mirrors which used in previous Monte Carlo algorithms. The Na atoms are not physically present on the solid network, but the Na₂O concentration influences the quantities of four-coordinated Al and B atoms, as well as three-coordinated B atoms.

As a result, the network encompasses both tetrahedrally coordinated species (Si, Al, [4]B) and tri-coordinated species ([3]B). An illustrative example of the ordered network representing the solid portion of the system is depicted in Figure 1. Throughout the manuscript, this structured and cubic lattice will be referred to as the “solid network.”

A significant departure from previous Monte Carlo algorithms is the introduction of a secondary subnetwork, henceforth referred to as the “liquid network,” to incorporate the solution and simulate its interactions with the glass. The nodes of the “liquid network” are positioned at the midpoint of bonds connecting two nodes on the solid network, effectively intertwining the liquid network with the solid network. The water molecules, denoted by red squares in Figure 1, are initially placed on the nodes of the liquid network within the volume corresponding to the empty portion of the solid network. This process creates an
interface between the solid and the solution, resulting in their initial separation.

As the solid and liquid networks are interconnected, the diffusion of water within the solid network, independent of hydrolysis or recondensation mechanisms, is permitted. As said before, in the previous Monte Carlo methods, water diffusion into the glass was solely achieved through bond hydrolysis and the release of elements into the solution.

A Monte Carlo simulation involves progressing through various stages representing elementary mechanisms occurring at the solid-solution interface. Each event is assigned a specific probability to regulate its frequency. Below, we outline the eight mechanisms considered.

1. Diffusion of a water molecule within the liquid network: This process capitalizes on the interplay between the two intertwined networks. A probability, termed wjump, governs the frequency of a water molecule's movement from one site to a neighboring site within the liquid network. Notably, a water molecule can only jump to an adjacent site that is unoccupied. Initially, the solution and the glass are solely in contact through their interface. However, due to diffusive jumps, water molecules can rapidly occupy interstitial sites within the solid network. Consequently, when a water molecule has jumped into the glass, it can be positioned at the midpoint of a bond connecting two atoms on the solid network, enabling the hydrolysis of this bond. Furthermore, when a water molecule from the primary solution traverses into the solid network, a new water molecule is introduced at the initial site. This mechanism treats the primary solution as an infinite reservoir of water molecules.

However, it quickly became apparent that the characteristic time for the jump probability of diffusing water molecules within the glass (associated to 1/wjump) could be significantly lower than the characteristic time for hydrolysis-recondensation mechanisms (associated respectively to 1/wbreak and 1/wrefor). In consequence, the use of a simple probability wjump applied to all possible diffusive jumps was insufficient to reproduce the amount of released boron in solution. Therefore, the code was modified to replace wjump with another parameter called nc_vois, which regulates the number of times the water molecule diffusion loop is called during a Monte Carlo step. This allowed for an increase in the number of diffusive events and consequently increased the amount of boron released in solution, in accordance with experimental observations. Therefore, at each Monte Carlo step, the subroutine for identifying all possible water molecule jumps is executed nc_vois times. All feasible jumps are then applied. For instance, if the final site of a possible event is occupied

due to a previous jump, the event is not applied. Moreover, considering the diffusion coefficients vary by a factor of 30 between SBNA3.5 and SBNA6 (if considering equivalent thickness after 16 days), adjustments in the equivalence between real-time and one Monte Carlo step are necessary to accurately represent this range (further details are provided in the text).

At every Monte Carlo step, the water molecules are categorized into two groups: the main solution, comprising water molecules interconnected (i.e., located on neighboring sites) and in contact with the initial solution. The second group consists of isolated water molecules (located on the liquid network), which are disconnected from the main solution, indicating there is no continuous path connecting these water molecules to the main solution. Further details on the identification of water molecules belonging to the main solution are provided in the Supporting Information. By default, the remainders are classified as isolated water molecules.

2. Hydrolysis of bonds within the glass network is a critical aspect of the simulation. The probability of bond opening is denoted as wbreak. Although theoretically, different probabilities could be utilized depending on the chemical nature of the X-Y bonds (where X and Y represent Si or Al), in practice, only one value is presently defined. Within a given glass composition, the probability of breaking Si-O-Si, Si-O-Al, and Al-O-Al bonds is uniform. This probability, however, is modified when the glass composition changes. When a water molecule is situated at the midpoint of an X-Y bond containing at least one B atom, hydrolysis occurs immediately, that is, with a probability of 1.0. This approach reflects our assumption that B atoms undergo hydrolysis more readily compared to Si or Al atoms.

The opening of an X-Y bond corresponds to the chemical reaction (3). In the Monte Carlo network, no oxygen (O) atoms are present, which is why we refer to X-Y bonds. However, in reality, these bonds correspond to X-O-Y bonds:

$$X-O-Y + H_2O \rightarrow X-O-H + Y-O-H \quad (3)$$

As a consequence of the hydrolysis mechanisms, the glassy structure undergoes progressive depolymerization, leading to the formation of non-bridging oxygens resulting from silanol or aluminol groups.

3. Bond reformation is governed by a probability parameter called wrefor. However, this mechanism is not applied to the bonds initially removed on the cubic network to introduce tri-coordinated B atoms. Additionally, a bond cannot be reformed around a B atom. When an X-Y bond (where X and Y represent Si or Al)
is reformed, it corresponds to the following chemical reaction:

$$X-O-H + Y-O-H \rightarrow X-O-Y + H_2O \quad (4)$$

After a bond reformation mechanism, a new water molecule is introduced into the liquid network and can once again diffuse into the interstitial sites within the solid network. As *wrefor* does not significantly affect the results in the first order, its value is kept constant at 1/1000 for all calculations. This indicates that a chemical bond has a one in 1000 chance of reforming.

4. Release of a glass element into the main solution occurs as the glassy network undergoes progressive depolymerization. It is possible that an atom (Si, Al, or B) becomes completely surrounded by broken bonds. Such atoms can be released into the solution. w

For Si and Al atoms, two scenarios are considered. If these atoms completely surrounded by broken bonds are in contact with the main solution, they are released into the solution. Otherwise, if they are not in contact with the main solution, they remain in their positions. As for B atoms completely surrounded by broken bonds, they are immediately released into the solution regardless of their locations, whether in contact with the main solution or not. This choice thus considered that diffusion of hydrolyzed species is much faster than the other mechanisms. To date, the possibility for hydrolyzed species to diffuse slowly within the liquid network before reaching the main solution is not implemented.

5. Redeposition of dissolved elements at the glass-water interface. When Si and Al atoms are released into the solution, their concentration in the solution increases. This mechanism involves considering the possibility for a Si or Al atom in solution to be redeposited at the glass-water interface. The glass-water interface is defined as the collection of sites within the solid network in contact with the main solution. The probability for this redeposition event is denoted as *wred*. Let us consider the case of Si atoms. If $C_{Si}$ is the concentration of Si in the main solution (i.e., the percentage of sites in the main solution occupied by Si atoms), the quantity of Si atoms redeposited during a Monte Carlo step is given by $wred \times C_{Si}$. To estimate the percentage of sites occupied by Si atoms in the main solution, we considered a virtual volume equivalent to 33,333 layers of 30 × 30 sites above the glass (it represents the solution), resulting in a surface-to-volume ratio of 1000 cm⁻¹. As probabilities will be calibrated against experimental data later on, the surface-to-volume ratio serves as an intermediate value in our calculations, where adjusting it essentially modifies the redeposition probability. The redeposited Si atoms are randomly placed on the available sites at

the glass-water interface. The same procedure applies to the Al atoms. This step allows for the gradual formation of an alteration layer enriched with Si and Al between the pristine glass and the main solution. In this study, *wred* is different from *wrefor* because the events considered by *wred* involve the formation of bonds between a Si or Al atom dissolved in the solution (surrounded by OH groups) and an Si or Al atom in the solid (surrounded by closed bonds), whereas the events considered by *wrefor* involves the reformation of bonds between atoms within the solid network (both surrounded by closed bonds). There is no expectation for the energies represented by *wred* and *wrefor* to be identical.

6. Vacancy migration. Jump of vacancies within the solid network occurs when B, Si, and Al atoms surrounded by broken bonds are released into the solution, leaving empty sites within the solid network, known as vacancies. The possibility for vacancies to diffuse within the solid network was implemented with a probability referred to as *wvacan*. Different probabilities can be considered depending on the change in the number of neighboring vacancy pairs after the jump. Specifically, a higher probability may be applied to promote the clustering of vacancies, that is, when the number of neighboring vacancy pairs increases after a jump. This facilitates the gradual formation of larger pores. It is important to note that when a vacancy jumps, the broken bonds are reorganized around the exchange sites to maintain the same number of broken bonds. A vacancy jump is never associated with hydrolysis, reformation, or redeposition reactions. The possibility for vacancy site migration within the solid network will be fully explored in later studies on gel ripening. Currently, this capability is deemed noncritical, and the parameter *wvacan* is maintained constant at 1/250. This implies that a vacancy migration event has a probability of 1/250 of occurring.

7. Jump of water molecules within the solid network. The phenomenon of water molecule migration within the solid network occurs when a water molecule within the liquid network moves into an adjacent vacancy in the solid network. This process occurs when a water molecule is positioned adjacent to a vacancy site within the solid structure. In such cases, the vacancy is promptly occupied by the water molecule, thereby causing the water molecule to be extracted from the liquid network. Consequently, the site vacated by the water molecule becomes accessible for subsequent diffusion mechanisms. Water molecules on the solid network have the potential to migrate through exchanges with atoms located on the solid network. This option of the code will be employed later to simulate gel ripening,
but it is not utilized in the current study, which con- centrates on the events at the glass-solution interface.

8. Dissolution of solid clusters formed in the main solu- tion. The dissolution of solid clusters formed within the main solution involves the elimination of clus- ters composed of atoms from the solid network that have become isolated within the main solution. Upon detection of these clusters, they are promptly dissolved, liberating the Si, Al, and B atoms into the solution. Con- sequently, the concentrations of Si, Al, and B in the main solution are elevated accordingly.

To outline, the implementation of the new Monte Carlo algorithm involves several steps as follows. Initially, the structure is prepared comprising the solid portion, the solution portion, and the interface between them. The Monte Carlo steps are then applied, each encompassing the following sub-steps:

- Identifying water molecules belonging to the main solu- tion and isolated water molecules (first pass).

- Jump of water molecules within the liquid network. The loop is repeated *nc_vois* times.

- Assessing potential breakage events for chemical bonds within the solid network, executed based on the proba- bility *wbreak*.

- Determining events for bond reformation, applied according to the probability *wrefor*.

- Re-evaluating water molecules belonging to the main solution and isolated water molecules (second pass).

- Release of Si, Al, and B atoms into solution. Completely hydrolyzed Si and Al atoms are released if in contact with the main solution, whereas completely hydrolyzed B atoms are released regardless.

- Reassessing water molecules belonging to the main solution and isolated water molecules (third pass).

- Removing isolated solid clusters within the main solu- tion.

- Identifying events for Si and Al redeposition at the glass- water interface, executed based on the probability *wred*.

- Reassessing water molecules belonging to the main solution and isolated water molecules (fourth pass).

- When a vacancy on the solid network is adjacent to a water molecule in the neighboring site within the liquid network, filling the vacancy with water. If mul- tiple neighboring water molecules exist, one is ran- domly selected. And subsequently removing this water molecule from the liquid network.

- Identifying events for vacancies (filled or empty) to jump within the solid network, with execution based on the probability *wvacan*. The promotion of vacancy clustering can be achieved using various probabilities.

Table 2 provides a summary of the different probabilities and their corresponding functions.

<table><caption>TABLE 2 Probabilities and parameter of the Monte Carlo code.</caption><thead><tr><td>Probability/ Parameter</td><td>Role</td></tr></thead><tbody><tr><td>nc_vois</td><td>Number of times the loop controlling the jumps of water molecules on the liquid network is called</td></tr><tr><td>wbreak</td><td>Opening of chemical bonds</td></tr><tr><td>wrefor</td><td>Reformation of chemical bonds</td></tr><tr><td>red</td><td>Redeposition of Si and Al atoms at the glass-water interface</td></tr><tr><td>wvacan</td><td>Jump of vacancies within the solid network</td></tr><tr><td>srempt</td><td>Parameter introduced to limit the quantity of water molecules within the solid</td></tr></tbody></table>

To compare the results of Monte Carlo calculations with experiments, it is necessary to convert a Monte Carlo step into an effective duration. The method by which this conversion is carried out is explained in Section 3.

One additional aspect needs clarification to complement the previous explanation. A parameter, *srempt*, is introduced to regulate the amount of water within the solid. If we did not regulate this, due to the diffusion events controlled by *nc_vois*, all sites of the liquid network within the glass would eventually be occupied. This parameter determines the maximum percentage of sites on the liquid network that can be occupied by water molecules, considering only those not belonging to the main solution. This means that only the water molecules not linked to the main solution are considered. Specifically, the probability for a water molecule to jump within the liquid network is equal to the factor (*srempt*—*srempli*)/*srempt*. Here, *srempli* represents the current percentage of occupied sites by water molecules. Consequently, the jump probability reaches zero when the maximum capacity for water molecule occupancy is reached. In such cases, no additional water molecules from the main solution can diffuse on the liquid network. However, water molecules continue to diffuse deeper within the liquid network, causing *srempli* to decrease (because *srempli* represents the concentration of water molecules on the liquid network, and the volume occupied by the water molecules increases in this scenario), thereby allowing the penetration of main solution-sourced water molecules to resume. In practice, a stationary state is achieved after several thousand steps. Without imposing a limit on the quantity of water molecules within the liquid network, it would become rapidly saturated, which is physically unrealistic. Here, the parameter *srempt* is consistently maintained at 10% throughout all calculations.

During a Monte Carlo calculation, the water molecules initially present in the main solution gradually infiltrate the liquid network. Additionally, atoms from the glass can undergo hydrolysis and subsequently be released into the
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Equivalent thickness of dissolved glass calculated from the amounts of released Si into solution for glasses SBNA3.5, SBNA4, and SBNA6 during the 6-month alteration experiments.</figcaption></figure>

solution. This process enables the tracking of the quantities of elements (Si, Al, or B) released into the solution over time. These quantities will be compared with experimental data, as illustrated below.

3 | RESULTS

**3.1** **Experiments**

The results of the glass alteration experiments for SBNA3.5, SBNA4, and SBNA6 over a duration of 6 months are presented in a former article.⁵⁹ It is noteworthy that the pH remains stable throughout the experiments, and the three glasses studied form a gel but no crystallized secondary phases. Figure 2 illustrates the calculated equivalent thickness based on the release of Si into solution over time. It is observed that a saturation plateau is quickly reached, with its altitude decreasing as the Al content increases. This outcome has been attributed to an increase in the hydrolysis energy of bonds around Si when they are in proximity to Al.⁶¹ Concurrently, a decrease in the energies of bond reformation for these same bonds has been highlighted. The inset in Figure 2 depicts an enlargement of the curves between 0 and 500 h. The plateau is reached more rapidly in glasses containing higher Al contents (SBNA4 and SBNA6).

In these experiments, B is not retained in the gels and
can be considered as a tracer of alteration. This means that
the equivalent thickness of B corresponds to the thickness
of dissolved glass and, assuming isovolumetric alteration,
it also corresponds to the gel thickness. The equivalent
thicknesses calculated from the concentrations of B in
solution are plotted in Figure 3 as a function of the square
root of time. It is indeed expected that under conditions

<figure><img src="image_5.png" /><figcaption>FIGURE 3 Equivalent alteration thickness measured from the amounts of released B into solution for glasses SBNA3.5, SBNA4, and SBNA6 during the 6-month alteration experiments.</figcaption></figure>

of Si saturation, gel growth leads to a release of B propor-
tional to the square root of time. Indeed, for the SBNA4
and SBNA6 glasses, the releases initially increase linearly
with the square root of time before bending, indicating
a behavior reflecting an evolution of transport properties
within the gel. The same conclusion likely applies to glass
SBNA3.5, despite the absence of one intermediate point at
1.8 days due to experimental issues. However, the curva-
ture of the curve occurs more quickly for glass SBNA3.5
compared to the other two glasses.

Deviation from a monotonic diffusive regime reflects the
formation of an increasingly passivating gel between the
glass and the solution. The formation of this passivating
gel is faster for glass SBNA3.5 than for glasses SBNA4 and
SBNA6. One possible explanation is the increasing rigidity
of the gel as the Al content increases, limiting the dynamics
of bonds (i.e., the frequency of bond openings and closures)
and thus the possibilities of local rearrangements. Indeed,
as mentioned earlier, increasing the Al₂O₃ content leads to
higher hydrolysis energies of bonds around Si and lower
energies of bond reformation. This can lead, on the one
hand, to more difficult bond openings in the gel and, on
the other hand, to faster bond reformation. Hence, there
is greater difficulty in constructing a passivating gel that
requires greater degrees of freedom for bond opening at
one location and closure at another after bond rotation.

**3.2** | Monte Carlo simulations

3.2.1 | Parameter tuning based on experiments

Initially, a series of calculations were conducted by vary-
ing the parameters *wbreak*, *wred*, and *nc_vois* over wide
ranges of values to find the set of parameters that best reproduces the experimental results for glass SBNA6. The simulation box size was as follows: 1755 000 sites (a grid of 30 × 30 × 1950) on the solid network, comprising 1059 980 Si, 190 418 Al, and 504 602 B. The calculations were systematically run for 10⁶ steps. The Z-axis is perpendicular to the solid-solution interface. The total size of the box, including both the solid and the liquid, was 30 × 30 × 2000, meaning that the last 50 layers were reserved for the initial solution. The initial position of the glass-water interface is at Z = 1950. Periodic boundary conditions are applied along the X and Y axes.

The best fit was achieved with the following values: wbreak = 0.003, wred = 0.42 nc_vois = 16. A value of wred = 0.42 means that a redeposition event of Si or Al was applied with the following p_redep probability:

$$p_{redep} = wred \times C_{X} \qquad (5)$$

where C_X represents the solution concentration of Si or Al (these concentrations are between 0 and 1), depending on the atom being considered.

A value of nc_vois = 16 means that the loop for water molecule diffusion mechanisms was applied 16 times, whereas the loops for other mechanisms were applied only once. This significant acceleration of diffusion compared to other mechanisms was necessary to best replicate the experimental values, particularly the release of B into solution.

Figure 4 compares the quantities of Si and B released into the solution as a function of the square root of time between simulation and experiment. The simulation fairly accurately reproduces the experiment in terms of the amount of Si released over time, whereas reaching saturation value occurs more quickly with the Monte Carlo calculation. The initial slope of B release is well reproduced by the simulation, but not the nonlinearities that subsequently appear in the experimental curve, that could be due to measurement uncertainties. The entirety of the experimental curves is not simulated, as Monte Carlo calculations are limited to 10⁶ steps due to computational time constraints. At steady state, the final concentrations of Si and Al in the solution, expressed as percentages of sites, are 0.11% and 0.02%, respectively.

The equivalence between a Monte Carlo simulation step and real time was calculated as follows: First, we identified the number of Monte Carlo steps required to reach Si saturation, which was determined to be 120 000 steps. Then, we equated this value with the experimentally determined time required to reach Si saturation. An experimental saturation time of 191.8 h was considered here. This equivalence resulted in a duration of one Monte Carlo step being 5.75 s.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Equivalent thicknesses calculated from the amounts of Si and B released into solution for glass SBNA6 with parameter settings wbreak = 0.003, wred = 0.42, and nc_vois = 16. Monte Carlo simulation: blue. Experiment: orange.</figcaption></figure>

Figure 5 shows the profile of the quantity of Si+Al and that of the degree of cross-linking (expressed as the percentage of closed bonds on the solid network, that is, bonds connected two atoms among Si, Al, and B) as a function of depth within the alteration gel. The quantity of Si+Al and the degree of cross-linking were calculated layer by layer on the solid network (with one layer corresponding to a plane of 30 × 30 sites containing 1800 closed bonds at maximum). The depth was converted into Angstroms considering that the interval between two successive layers of atoms is 3.5 Å. The curves corresponding to the case that best reproduce the experiment are shown in orange (other parameter settings are also displayed in the figure, which will be discussed later). The number of timesteps for the different calculations is not strictly equal due to some calculations being interrupted prematurely by unexpected stops of the computing cluster. A steady state is reached with respect to water diffusion inside the gel, as the degree of cross-linking remains constant throughout the glass. Toward the exterior of the gel, at the interface with the solution, there is a layer enriched in Si+Al, which corresponds to a peak in the degree of cross-linking. A more compact zone enriched in Si+Al has thus formed toward
<figure><img src="image_7.png" /></figure>

FIGURE 5 Monte Carlo simulations for glass SBNA6. Top: profile of the quantity of Si+Al at the interface between the gel and the solution. Bottom: profile of the degree of cross-linking (i.e., percentage of closed bonds in each layer parallel to the glass-solution interface) at the interface between the gel and the solution. Multiple cases corresponding to different parameter settings are shown (see text for explanations). For all the calculations, *nc_vois* = 16. The numbers of Monte Carlo timesteps are as follows: 778 500 (*wbreak* = 0.0015), 1000 000 (*wbreak* = 0.003), and 635 100 (*wbreak* = 0.004).

the exterior of the gel through the interplay of hydrolysis and redeposition mechanisms. The morphology of this zone is shown in Figure 6 (along with those obtained with other parameter settings). In this figure, water molecules belonging to the main solution are depicted in grey, and water molecules trapped in the vacancies on the solid network are depicted in red.

A series of calculations were then conducted to test the role of the parameters *wbreak* and *wred* on both the SBNA6 and SBNA4 glasses, using the same parameter sets to highlight the influence of glass composition. These results are detailed in Section 3.2.2. From all these calculations, it was possible to identify the *wbreak*, *wred*, and *nc_vois* parameters leading to the best reproduction of the experiments for glass SBNA4. The optimal parameter set is as follows: *wbreak* = 0.0035; *wred* = 0.3; and *nc_vois* = 16.

Figure 7 depicts the equivalent thicknesses measured from the amounts of Si and B released into solution as a function of the square root of time. The Monte Carlo

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

FIGURE 6 Monte Carlo simulations for glass SBNA6. Morphologies of alteration gels formed outside the gel for different parameter settings. Si: yellow, water molecules in the vacancies of the solid network: red, water molecules from the main solution: gray. The numbers of Monte Carlo timesteps are as follows: 778 500 (*wbreak* = 0.0015), 1000 000 (*wbreak* = 0.003), and 635 100 (*wbreak* = 0.004).

time step for this simulation corresponded to a duration of 0.775 s. Once again, we observe the slightly too rapid attainment of the steady state for the quantity of Si in solution. However, the slope of B release is very close to the experimental curve. The entirety of the experimental curves was not simulated due to computational time constraints. Here, the Monte Carlo calculation was run for 10⁶ steps. At steady state, the final concentrations of Si and Al in the solution, expressed as percentages of sites, are 0.20% and 0.025%, respectively.

Figure 8 illustrates the profiles of the quantities of Si+Al and the degree of cross-linking at the interface between the gel and the solution as a function of depth. Again, a peak in the Si+Al profile and in the degree of cross-linking is observed toward the exterior of the gel for the optimal
<figure><img src="image_9.png" /><figcaption>FIGURE 7 Equivalent thicknesses calculated from the amounts of Si and B released into solution for glass SBNA4 with the parameter setting *wbreak* = 0.0035; *wred* = 0.3; and *nc_vois* = 16. Monte Carlo simulation: blue. Experiment: orange.</figcaption></figure>

parameter setting (shown in gray in Figure 8). It indicates the presence of a denser layer as depicted in Figure 9 (other simulations with different parameter settings are shown in Figures 8 and 9, which will be discussed later).

To reproduce the experimental results corresponding to glass SBNA4, it was necessary to increase the parameter regulating the hydrolysis mechanisms (*wbreak*) and reduce the parameter regulating the redepositions (*wred*) in order to increase the amount of Si released into solution. As the Monte Carlo time step corresponds to a shorter duration compared to glass SBNA6, the curve of B release into solution also remains close to the experimental values.

### 3.2.2 | Parametric study

To better understand the role of various parameters, we varied *wbreak* and *wred* within intervals around the reference values obtained for the SBNA6 glass. Thus, for SBNA6 glass, *wbreak* was varied between 0.0015 and 0.004 while keeping *wred* at 0.42. The results are shown in Figure 10. As expected, an increase in the amount of released Si into the solution is observed as the value of *wbreak* increases, indicating accelerated hydrolysis. At steady state, the final

<figure><img src="image_10.png" /><figcaption>FIGURE 8 Monte Carlo simulations for glass SBNA4. Top: profile of the quantity of Si+Al at the interface between the gel and the solution. Bottom: profile of the degree of cross-linking (i.e., percentage of closed bonds in each layer parallel to the glass-solution interface) at the interface between the gel and the solution. Multiple cases corresponding to various parameter settings are shown (see text for explanations). The number of Monte Carlo timesteps is set to 1000 000 for each calculation.</figcaption></figure>

concentrations of Si (Al) in the solution, expressed as percentages of sites, are respectively equal to 0.032% (0.005%), 0.11% (0.02%), and 0.17% (0.031%) for *wbreak* = 0.0015, *wbreak* = 0.003, and *wbreak* = 0.004.

The study of Si+Al profiles and the degree of cross-linking (Figure 5) reveals the formation of an outer layer enriched in Si and Al at the gel-solution interface (with a peak at the glass-solution interface), but the height of this peak decreases as the value of *wbreak* increases. This observation is further discussed.

However, despite these differences in behavior at the interface between the gel and the solution, the impact on the quantity of released B into the solution remains marginal, and the curves are close to each other.

Interestingly, the case of SBNA4 glass shows a change in alteration behavior when the calculation parameters change. The parameter *wbreak* was varied between 0.0015 and 0.004 while keeping the value of *wred* at 0.42. For the highest values of *wbreak*, peaks of Si release are observed, corresponding to solid alteration by clusters at certain moments (Figure 11). The peaks of Si release are significant
<figure><img src="image_11.png" /><figcaption>FIGURE 9 Monte Carlo simulations for glass SBNA4. Morphologies of alteration gels formed outside the gel for different parameter settings. Si: Yellow, water molecules in the vacancies of the solid network: Red, water molecules from the main solution: Gray. The number of Monte Carlo timesteps is set to 1000 000 for each calculation.</figcaption></figure>

because, after detailed analysis of the results, it appears
that a poorly adherent layer forms on the surface due to
the high roughness of the interface between the solid and
the solution. This layer is poorly adherent in the sense that
it is connected to the underlying solid by a small num-
ber of bonds. Consequently, this layer easily transforms
into isolated clusters from the solid, as it only requires a
few bonds to break. Therefore, it is regularly dissolved.
In Figure 11, at steady state, the final concentrations of
Si (Al) in the solution, expressed as percentages of sites,
are respectively equal to 0.036% (0.004%), 0.06% (0.007%),
0.086% (0.010%), 0.117% (0.014%), 0.148% (0.019%), and
0.208% (0.026%) for *wbreak* = 0.0015, *wbreak* = 0.002,
*wbreak* = 0.0025, *wbreak* = 0.003, *wbreak* = 0.0035, and
*wbreak* = 0.004.

Increasing *wbreak* while keeping *wred* constant leads to
a decrease in the overconcentration of Si+Al appearing at
the interface between the gel and the solution (Figure 12).
It even disappears completely for the highest value of
*wbreak* = 0.004. Figure 13 shows the different morpholo-
gies of the gel for the various values of *wbreak*. The outer

layer enriched in Si+Al becomes less distinct and increas-
ingly porous as *wbreak* increases until it gives way to a
highly porous zone when *wbreak* = 0.004.

Similar behavior was observed for SBNA3.5 glass (see
Figure S1) with the onset of cluster hydrolysis occurring
earlier than for SBNA4 glass. The amount of B₂O₃ and
Na₂O is higher in SBNA3.5 glass, which therefore con-
tains less SiO₂ and Al₂O₃. It becomes easier in this case
to hydrolyze the material through clustering because the
number of bonds to break to isolate a volume of solid is
lower.

Starting from the optimal parameterization for SBNA4
glass, a parametric study was conducted by varying the
parameter *wred*. When *wred* decreases, meaning the rate
of redeposition on the gel surface is reduced, there is a sig-
nificant increase in the quantities of Si and B released into
the solution (results are shown in Figure S2). This result
was expected as hydrolysis is less and less compensated by
redeposition. Furthermore, as *wred* decreases, the Si+Al-
enriched outer layer becomes less distinct (Figure 8). A
cross-section of the gels is presented in Figure 9.
<figure><img src="image_12.png" /><figcaption>FIGURE 10 Equivalent thicknesses calculated from the quantities of Si and B released into solution for SBNA6 glass as a function of the time step number for different values of wbreak. In all calculations, wred = 0.42 and nc_vois = 16.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 11 Equivalent thicknesses calculated from the quantities of Si and B released into solution for SBNA4 glass as a function of the time step number for different values of wbreak. In all calculations, wred = 0.42 and nc_vois = 16.</figcaption></figure>

# 4 | DISCUSSION

## 4.1 | Effect of composition and rate of hydrolysis at the gel-solution interface

One result that previous simulation methods failed to achieve is the modification of the glass alteration mode when the parameters of Monte Carlo calculations change, particularly when the hydrolysis or redeposition rates of species are modified. For SBNA4 and SBNA3.5 glasses, when the hydrolysis rate increases (*wbreak*), a transition is observed from an alteration mode that leads to the formation of an outer layer enriched in Si+Al to a mode of alteration more dominated by clusters, resulting in the appearance of a porous external zone without enrichment in Si+Al.

This transition occurs at a lower value of *wbreak*, indicating a lower hydrolysis rate, for SBNA3.5 glass compared to SBNA4 glass. It is possible that the external roughness of the glass, which increases when the content of soluble elements (here, B and Na atoms) increases, plays a role in the alteration behavior. SBNA3.5 glass contains more B and Na than SBNA4 glass. Consequently, due to the ability of water to penetrate inside the glass, an increasingly rough external layer appears as water molecules penetrate, subse-

quently facilitating the hydrolysis of less soluble elements (Si and Al) by clusters rather than atom by atom.

However, the composition of the glass is not the sole criterion at play, as the evolution of *wbreak* also strongly impacts the alteration behavior of the glass at constant composition. When hydrolysis is accelerated without modifying the redeposition rate of Si and Al elements (*wred*), the transition toward a more cluster-dominated alteration eventually occurs. This has been observed for SBNA3.5 and SBNA4 glasses as discussed earlier. It may be necessary to further increase *wbreak* to observe this phenomenon for SBNA6 glass. This transition is gradual. Initially, a decrease in the intensity of Si+Al enrichment in the outer part of the gel is observed, accompanied by a decrease in the peak of cross-linking degree. Then, it becomes impossible to redeposit sufficient Si and Al to allow the formation of this denser layer because hydrolysis has become too rapid, and the Si+Al-enriched outer layer disappears. When the hydrolysis rate remains low (indicated by low values of *wbreak*), there is more time for the formation of an external layer enriched with Si and Al, making the upper part of the glass surface more stable. This process has been previously highlighted by Kerisit et al.⁵¹

Thus, it appears that the difference between the outflow flux of elements (hydrolysis) and the inflow flux
<figure><img src="image_14.png" /><figcaption>FIGURE 12 Monte Carlo simulations for glass SBNA4. Top: profile of the quantity of Si+Al at the interface between the gel and the solution. Bottom: profile of the degree of cross-linking (i.e., percentage of closed bonds in each layer parallel to the glass-solution interface) at the interface between the gel and the solution. Multiple cases corresponding to various parameter settings are shown (see text for explanations). The calculations are performed with *wred* = 0.42 and *nc_vois* = 16. The number of Monte Carlo timesteps is set to 1000 000 for each calculation.</figcaption></figure>

(redeposition) also plays a significant role in determining
the glass alteration behavior. By increasing *wbreak* with-
out modifying *wred*, thereby widening this difference, a

transition toward a more cluster-releasing alteration mode
is observed. One possible explanation is as follows: As
the outflow flux of elements increases, the roughness of
the interface between the glass and the solution increases,
making it more challenging to form an outer layer enriched
in Si and Al because more voids need to be filled for this
layer to cover the entire interface. Redeposition, which
remains unchanged, is no longer sufficient. The glass, with
an increasingly rough outer layer, eventually undergoes
more cluster-releasing alteration.

Confirmation of this explanation is provided by the
behavior of SBNA4 glass when the parameters are set to
the following values: wbreak = 0.005, wred = 0.4, and
nc_vois = 8. It is observed (Figure 14) that in this case,
the experiment is relatively well reproduced concerning
the releases of Si and B, but the alteration occurs through
the release of numerous clusters. In this example, the dif-
ference between the outflow and inflow fluxes of elements
is quite close to the reference adjustment, so that the net
releases of Si and B into solution closely reproduce the
experiment. However, the alteration behavior is still modi-
fied. That is, an excessive outflow flux of elements prevents
the formation of a denser outer layer despite compensa-
tion by an increased inflow flux. Not only the difference
between the outflow and inflow fluxes of elements but also
the overall activity at the gel-solution interface, defined as
the sum of the outflow and inflow fluxes of elements, thus
impacts the glass alteration behavior. In Figure 14, at steady
state, the final concentrations of Si and Al in the solution,
expressed as percentages of sites, are 0.308% and 0.038%,
respectively.

This potential impact of interface morphology depen-
dent on the hydrolysis rate had already been discussed
in a previous article based on the former Monte Carlo
method.⁶² In this work, two borosilicate glasses, one with-
out Al₂O₃ and the other containing Al₂O₃, were simulated

<figure><img src="image_15.png" /><figcaption>FIGURE 13 Monte Carlo simulations for glass SBNA4. Morphologies of alteration gels formed outside the gel for different parameter settings. The calculations are performed with *wred* = 0.42 and *nc_vois* = 16. Si: yellow, water molecules in the vacancies of the solid network: red, water molecules from the main solution: gray. The number of Monte Carlo timesteps is set to 1000 000 for each calculation.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 14 Equivalent thicknesses calculated from the amounts of Si and B released into solution for glass SBNA4 with the parameter setting *wbreak* = 0.005; *wred* = 0.4; and *nc_vois* = 8. Monte Carlo simulation: blue. Experiment: orange.</figcaption></figure>

by Monte Carlo. It was found that the formation of the external protective film enriched in Si and Al took longer to form on the glass containing Al₂O₃ because the roughness of the interface between the gel and the solution was greater (due to the fact that in this work, Al was considered as less soluble than Si). We confirm in the present article the connection between external roughness and difficulty in forming a dense external layer.

## 4.2 | Links to the experiment

In this section we attempt to connect the previous observations from Monte Carlo calculations with experimental observations. The calculations thus reveal two types of alteration behavior depending on the glass composition and hydrolysis and redeposition parameters. Alteration behavior gradually changes along with the morphology of the interface between the gel and the solution.

Experimentally, recent works by Taron et al.,⁵⁸ Damodaran et al.,⁵⁹ Vienna et al.,⁶³ and Reiser et al.⁵⁷ have shown a transition in glass alteration behavior as the Al₂O₃ content increases. For zero or low Al₂O₃

content, rapid release of Si is observed followed by a sharp deceleration of alteration. In these cases, both Si and B concentrations quickly reach a plateau. This behavior has been explained, namely, according to TOF-SIMS characterizations data, by the precipitation of an outer precipitated layer made of amorphous silica, which plays a protective role by triggering a sudden decrease in the alteration rate.⁶²

Additionally, it has been observed that the lower the concentration of Al₂O₃, the easier the hydrolysis of Si, with the presence of Al in the immediate vicinity of Si atoms leading to an increase in the average hydrolysis energy.⁶¹

Monte Carlo calculations show that alteration is more likely to occur through the release of numerous clusters when hydrolysis or activity at the gel-solution interface is more significant, or when the content of B₂O₃+Na₂O is higher. As seen, a decrease in Al₂O₃ content leads to easier hydrolysis of Si. It also results in higher concentrations of B₂O₃+Na₂O if the SiO₂ concentration remains constant. Cluster release upon reaching the steady state, which is a specificity of Monte Carlo, nevertheless reflects the difficulty in forming a protective zone enriched in Si and Al because the interface between the gel and the solution becomes too rough. In this case, there remains the possibility of forming a protective zone on the gel by precipitating silica from the solution.

Unfortunately, the mechanisms currently implemented in this version of the Monte Carlo method do not allow for the simulation of silica precipitation in solution. However, the comparison of Monte Carlo results with experiments suggests the following sequence of mechanisms in the case of low Al₂O₃ content: release of a larger quantity of Si followed by silica precipitation at the gel-solution interface. This release is facilitated by the difficulty in forming an external zone enriched in Si and Al. A mechanism compatible with that observed in Monte Carlo simulations has been experimentally observed in a glass containing three oxides (without Al) altered for 33 days in a solution enriched with ²⁹Si and ¹⁸O. An enrichment of the external surface of the gel in ²⁹Si was observed, indicating a mechanism of rapid redeposition of the Si present in solution. This confirms that hydrolysis leading to significant release of Si atoms promotes their subsequent redeposition on the surface.⁶²

However, when the Al₂O₃ content increases, the average hydrolysis energy also increases, and the outflow flux at the gel-solution interface decreases.⁶¹ Experimentally, a change in alteration behavior is observed. There is no longer silica precipitation at the gel surface but rather reorganization of the gel through local mechanisms, resulting in a slower decrease in the alteration rate.
Monte Carlo calculations show a modification of the alteration behavior with the formation of an outer layer enriched in Si+Al when hydrolysis at the gel-solution interface is weaker (at constant redeposition rate). The formation of this Si+Al-enriched layer is observed close to the gel surface because the redeposition mechanism occurs there, and no internal gel reorganization mechanism is currently implemented in the method. However, if the rate of hydrolysis at the gel surface decreases, and as a result, the quantity of hydrolysis-recondensation mechanisms will decrease, it is possible that the formation of an outer layer enriched in Si+Al propagates toward the interior of the structure due to a rebalancing between the characteristic relaxation times (the time required to reform a bond) at the surface and in the bulk. This would lead to the appearance of a thicker protective layer, but conversely, it would be slower to form. The Monte Carlo results thus suggest the following mechanism in the case of a higher Al₂O₃ content: If hydrolysis remains sufficiently fast, a Si+Al-enriched layer may form relatively quickly toward the exterior of the gel. However, if hydrolysis slows down, the characteristic times between external and internal reorganizations may balance out, leading to the formation of a Si+Al-enriched zone over a greater thickness. This externally enriched layer of Si+Al, with increased cross-linking, would logically play a passivating role associated with the reduction of the alteration rate over time.

The existence of a dense outer layer in coherence with the Monte Carlo calculations has been observed in some experiments, notably during the alteration of a 4-oxide aluminoborosilicate glass for 21 years under static conditions⁶² or after the alteration of a 4-oxide borosilicate glass for 100 days.¹⁶

## 5 | CONCLUSIONS

A new Monte Carlo method has been developed to model the alteration of alumino-borosilicate glasses to overcome the limitations of previous approaches used so far. This new approach allows for the simulation of water diffusion within the glass network in addition to the classical mechanisms of hydrolysis and redeposition, which was not possible previously. This is achieved through the use of two shifted subnetworks, with one hosting the solid and the other the solution.

Incorporating this new diffusion mechanism has revealed new alteration behaviors dependent on the hydrolysis rate at the gel-solution interface. When the hydrolysis rate remains low, an external layer enriched in Si and Al gradually develops at the gel-solution interface during alteration. As it increases, alteration

occurs more through the release of clusters, preventing the formation of an outer layer enriched in Si and Al. Instead, an outer zone that is more porous and jagged is formed.

The glass composition plays a significant role in the alteration behavior. The transition from the first behavior to the second occurs at a lower hydrolysis rate threshold when the glass is more enriched in B₂O₃ and Na₂O.

More broadly, these findings highlight the importance of the morphology of the interface between the glass and the solution in the ability to form an external layer enriched in Si and Al. This morphology depends on both the hydrolysis and redeposition rates and the glass composition.

These innovative works will be continued by expanding the list of implemented mechanisms. In particular, it will be important to introduce the possibility for B atoms to diffuse within the gel before being released into solution.

## ACKNOWLEDGMENTS

This study was funded by CEA. The authors are grateful to the CEA LMAC team for the solution analyses.

## CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

## ORCID

Jean-Marc Delaye ID https://orcid.org/0000-0002-1311-642X

Paul C. M. Fossati ID https://orcid.org/0000-0001-8230-6422

Stéphane Gin ID https://orcid.org/0000-0002-1950-9195

## REFERENCES

1. Lee WE, Ojovan MI, Stennett MC, Hyatt NC. Immobilisation of radioactive waste in glasses, glass composite materials and ceramics. Adv Appl Ceram. 2006;105(1):3–12. https://doi.org/10.1179/174367606j81669

2. Gin S, Delaye J-M, Angeli F, Schuller S. Aqueous alteration of silicate glass: state of knowledge and perspectives. NPJ Mater Degrad. 2021;5:42(1). https://doi.org/10.1038/s41529-021-00190-5

3. Frugier P, Gin S, Minet Y, Chave T, Bonin B, Godon N, et al. SON68 nuclear glass dissolution kinetics: current state of knowledge and basis of the new GRAAL model. J Nucl Mater. 2008;380(1–3):8–21. https://doi.org/10.1016/j.jnucmat.2008.06.044

4. Weber WJ, Ewing RC, Angell CA, Arnold GW, Cormack AN, Delaye JM et al. Radiation effects in glasses used for immobilization of high-level waste and plutonium disposition. J Mater Res. 1997;12(8):1946–78. https://doi.org/10.1557/jmr.1997.0266

5. Poinssot C, Gin S. Long-term behavior science: the cornerstone approach for reliably assessing the long-term performance of nuclear waste. J Nucl Mater. 2012;420(1–3):182–92. https://doi.org/10.1016/j.jnucmat.2011.09.012

6. Frankel GS, Vienna JD, Lian J, Scully JR, Gin S, Ryan JV, et al. A comparative review of the aqueous corrosion of glasses, crys-
talline ceramics, and metals. NPJ Mater Degrad. 2018;2:15(1).
https://doi.org/10.1038/s41529-018-0037-2

7. Frankel GS, Vienna JD, Lian J, Guo X, Gin S, Kim SH, et al. Recent advances in corrosion science applicable to disposal of high-level nuclear waste. Chem Rev. 2021. https://doi.org/10.1021/acs.chemrev.0c00990

8. Vienna JD, Ryan JV, Gin S, Inagaki Y. Current understanding and remaining challenges in modeling long-term degradation of borosilicate nuclear waste glasses. Int J Appl Glas Sci. 2013;4(4):283-94. https://doi.org/10.1111/ijag.12050

9. Strachan D, Neeway JJ, Pederson L, Schreiber DK, Mitroshkov A, Zhu Z, et al. On the dissolution of a borosilicate glass with the use of isotopic tracing—insights into the mechanism for the long-term dissolution rate. Geochim Cosmochim Acta. 2022;318:213-29. https://doi.org/10.1016/j.gca.2021.12.004

10. Gin S, Guo X, Delaye JM, Angeli F, Damodaran K, Testud V, et al. Insights into the mechanisms controlling the residual corrosion rate of borosilicate glasses. NPJ Mater Degrad. 2020;4(1).
https://doi.org/10.1038/s41529-020-00145-2

11. Geisler T, Nagel T, Kilburn MR, Janssen A, Icenhower JP, Fonseca ROC, et al. The mechanism of borosilicate glass corrosion revisited. Geochim Cosmochim Acta. 2015;158:112-29.
https://doi.org/10.1016/j.gca.2015.02.039

12. Hellmann R, Cotte S, Cadel E, Malladi S, Karlsson LS, Lozano-Perez S et al. Nanometre-scale evidence for interfacial dissolution-reprecipitation control of silicate glass corrosion. Nat Mater. 2015;14(3):307-11. https://doi.org/10.1038/nmat4172

13. Tribet M, Mir AH, Gillet C, Jegou C, Mougnaud S, Hinks JA, et al. New insights about the importance of the alteration layer/glass interface. J Phys Chem C. 2020;124(18):10032-44.
https://doi.org/10.1021/acs.jpcc.0c02121

14. Bunker BC, Arnold GW, Day DE, Bray PJ. The effect of molecular-structure on borosilicate glass leaching. J Non Cryst Solids. 1986;87(1-2):226-53. https://doi.org/10.1016/S0022-3093(86)80080-1

15. Mir AH, Jan A, Delaye J-M, Donnelly S, Hinks J, Gin S. Effect of decades of corrosion on the microstructure of altered glasses and their radiation stability. NPJ Mater Degrad. 2020;4:11(1). https://doi.org/10.1038/s41529-020-0115-0

16. Cailleteau C, Angeli F, Devreux F, Gin S, Jestin J, Jollivet P, et al. Insight into silicate-glass corrosion mechanisms. Nat Mater. 2008;7(12):978-83. https://doi.org/10.1038/nmat2301

17. Gin S, Jollivet P, Fournier M, Berthon C, Wang Z, Mitroshkov A, et al. The fate of silicon during glass corrosion under alkaline conditions: a mechanistic and kinetic study with the International Simple Glass. Geochim Cosmochim Acta. 2015;151:68-85.
https://doi.org/10.1016/j.gca.2014.12.009

18. Hopf J, Eskelsen JR, Chiu M, Ievlev AV, Ovchinnikova OS, Leonard D, et al. Toward an understanding of surface layer formation, growth, and transformation at the glass-fluid interface. Geochim Cosmochim Acta. 2018;229:65-84. https://doi.org/10.1016/j.gca.2018.01.035

19. Angeli F, Villain O, Schuller S, Ispas S, Charpentier T. Insight into sodium silicate glass structural organization by multinuclear NMR combined with first-principles calculations. Geochim Cosmochim Acta. 2011;75(9):2453-69. https://doi.org/10.1016/j.gca.2011.02.003

20. Molieres E, Angeli F, Jollivet P, Gin S, Charpentier T, Majérus O, et al. Chemical durability of lanthanum-enriched borosilicate

glass. Int J Appl Glas Sci. 2013;4(4):383-94. https://doi.org/10.1111/ijag.12054

21. Kaya H, Gin S, Vogt BD, Kim SH. Impact of aqueous solution pH on network structure of corrosion-induced surface layers of boroaluminosilicate glass. J Am Ceram Soc. 2022;105(11):6581-92. https://doi.org/10.1111/jace.18606

22. Kaya H, Ngo D, Smith NJ, Gin S, Kim SH. Network structure in alteration layer of boroaluminosilicate glass formed by aqueous corrosion. J Non Cryst Solids. 2021;553. https://doi.org/10.1016/j.jnoncrysol.2020.120494

23. Geisler T, Dohmen L, Lenting C, Fritzsche MBK. Real-time in situ observations of reaction and transport phenomena during silicate glass corrosion by fluid-cell Raman spectroscopy. Nat Mater. 2019;18(4):342-48. https://doi.org/10.1038/s41563-019-0293-8

24. Rebiscoul D, der Lee A, Rieutord F, Né F, Spalla O, El-Mansouri A, et al. Morphological evolution of alteration layers formed during nuclear glass alteration:: new evidence of a gel as a diffusive barrier. J Nucl Mater. 2004;326(1):9-18. https://doi.org/10.1016/j.jnucmat.2003.10.015

25. Ngo D, Liu H, Sheth N, Lopez-Hallman R, Podraza NJ, Collin M, et al. Spectroscopic ellipsometry study of thickness and porosity of the alteration layer formed on international simple glass surface in aqueous corrosion conditions. NPJ Mater Degrad. 2018;2(1). https://doi.org/10.1038/s41529-018-0040-7

26. Jollivet P, Angeli F, Cailleteau C, Devreux F, Frugier P, Gin S. Investigation of gel porosity clogging during glass leaching. J Non Cryst Solids. 2008;354(45-46):4952-58. https://doi.org/10.1016/j.jnoncrysol.2008.07.023

27. Gin S, Collin M, Jollivet P, Fournier M, Minet Y, Dupuy L, et al. Dynamics of self-reorganization explains passivation of silicate glasses. Nat Commun. 2018;9. https://doi.org/10.1038/s41467-018-04511-2

28. Ngo D, Liu Hongshen, Chen Z, Kaya H, Zimudzi TJ, Gin S, et al. Hydrogen bonding interactions of H₂O and SiOH on a boroa-luminosilicate glass corroded in aqueous solution. NPJ Mater Degrad. 2020;4:1(1). https://doi.org/10.1038/s41529-019-0105-2

29. Zapol P, He H, Kwon KD, Criscenti LJ. First-principles study of hydrolysis reaction barriers in a sodium borosilicate glass. Int J Appl Glas Sci. 2013;4(4, SI):395-407. https://doi.org/10.1111/ijag.12052

30. Xiao YT, Lasaga AC. Ab initio quantum mechanical studies of the kinetics and mechanisms of quartz dissolution: OH-catalysis. Geochim Cosmochim Acta. 1996;60(13):2283-95. https://doi.org/10.1016/0016-7037(96)00101-9

31. Jabraoui H, Charpentier T, Gin S, Delaye J-M, Pollet R. Behaviors of sodium and calcium ions at the borosilicate glass-water interface: gaining new insights through an ab initio molecular dynamics study. J Chem Phys. 2022;156:134501(13). https://doi.org/10.1063/5.0087390

32. Jabraoui H, Charpentier T, Gin S, Delaye J-M, Pollet R. Atomic insights into the events governing the borosilicate glass-water interface. J Phys Chem C. 2021;125(14):7919-31. https://doi.org/10.1021/acs.jpcc.1c00388

33. Pelmenschikov A, Strandh H, Pettersson LGM, Leszczynski J. Lattice resistance to hydrolysis of Si-O-Si bonds of silicate minerals: ab initio calculations of a single water attack onto the (001) and (111) beta-cristobalite surfaces. J Phys Chem B. 2000;104(24):5779-83. https://doi.org/10.1021/jp994097r
34. Ohkubo T, Gin S, Collin M, Iwadate Y. Molecular dynamics simulation of water confinement in disordered aluminosilicate subnanopores. Sci Rep. 2018;8:3761. https://doi.org/10.1038/s41598-018-22015-3

35. Rimsza JM, Du J. Nanoporous silica gel structures and evolution from reactive force field-based molecular dynamics simulations. NPJ Mater Degrad. 2018;2:18(1). https://doi.org/10.1038/s41529-018-0039-0

36. Rimsza JM, Du J. Interfacial structure and evolution of the water-silica gel system by reactive force-field-based molecular dynamics simulations. J Phys Chem. 2017;121(21):11534-43. https://doi.org/10.1021/acs.jpcc.7b02734

37. Deng L, Miyatani K, Amma S, Suehara M, Ono M, Yamamoto Y, et al. Reaction mechanisms and interfacial behaviors of sodium silicate glass in an aqueous environment from reactive force field-based molecular dynamics simulations. J Phys Chem. 2019;123(35):21538-47. https://doi.org/10.1021/acs.jpcc.9b05030

38. Mahadevan TS, Du J. Atomic and micro-structure features of nanoporous aluminosilicate glasses from reactive molecular dynamics simulations. J Am Ceram Soc. 2021;104(1):229-42. https://doi.org/10.1111/jace.17465

39. Taron M, Delaye J-M, Gin S. A classical molecular dynamics simulation method for the formation of “dry” gels from boro-aluminosilicate glass structures. J Non Cryst Solids. 2021;553:120513. https://doi.org/10.1016/j.jnoncrysol.2020.120513

40. Collin M, Gin S, Dazas B, Mahadevan T, Du J, Bourg IC. Molecular dynamics simulations of water structure and diffusion in a 1 nm diameter silica nanopore as a function of surface charge and alkali metal counterion identity. J Phys Chem. 2018;122(31):17764-76. https://doi.org/10.1021/acs.jpcc.8b03902

41. Jabraoui H, Gin S, Charpentier T, Pollet R, Delaye J-M. Leaching and reactivity at the sodium aluminosilicate glass-water interface: insights from a ReaxFF molecular dynamics. J Phys Chem. 2021;125(49):27170-84. https://doi.org/10.1021/acs.jpcc.1c07266

42. Mahadevan TS, Du J. Surface reaction, water diffusion and mechanical properties of hydrated nanoporous calcium aluminosilicate gel structures. J Non Cryst Solids. 2023;621. https://doi.org/10.1016/j.jnoncrysol.2023.122604

43. Kagan M, Lockwood GK, Garofalini SH. Reactive simulations of the activation barrier to dissolution of amorphous silica in water. Phys Chem Chem Phys. 2014;16(20):9294-301. https://doi.org/10.1039/c4cp00030g

44. Aertsens M, VanIseghem P. Modelling glass dissolution with a Monte Carlo technique. In: Murphy WM, Knecht DA, editors. Sci. Basis Nucl. Waste Manag. XIX. Vol. 412. 1996. p. 271-78.

45. Ledieu A, Devreux F, Barboux P. Monte Carlo simulations of borosilicate glass corrosion: predictions for morphology and kinetics. J Non Cryst Solids. 2004;345:715-19. https://doi.org/10.1016/j.jnoncrysol.2004.08.152

46. Ledieu A, Devreux F, Barboux P. The role of aluminium in the durability of alumino-borosilicate glasses. Phys Chem Glas. 2005;46(1):12-20.

47. Ledieu A, Devreux F, Barboux P, Minet Y. Contribution of Monte Carlo modeling to understanding the alteration of nuclear glasses by water. Nucl Sci Eng. 2006;153(3):285-300. https://doi.org/10.13182/NSE06-A2614

48. Devreux F, Ledieu A, Barboux P, Minet Y. Leaching of borosilicate glasses. II. Model and Monte-Carlo simulations. J Non Cryst Solids. 2004;343(1-3):13-25. https://doi.org/10.1016/j.jnoncrysol.2004.06.007

49. Cailleteau C, Devreux F, Spalla O, Angeli F, Gin S. Why do certain glasses with a high dissolution rate undergo a low degree of corrosion? J Phys Chem C. 2011;115(13):5846-55. https://doi.org/10.1021/jp111458f

50. Kerisit S, Pierce EM. Monte Carlo simulations of the dissolution of borosilicate glasses in near-equilibrium conditions. J Non Cryst Solids. 2012;358(10):1324-32. https://doi.org/10.1016/j.jnoncrysol.2012.03.003

51. Kerisit S, Ryan JV, Pierce EM. Monte Carlo simulations of the corrosion of aluminoborosilicate glasses. J Non Cryst Solids. 2013;378:273-81. https://doi.org/10.1016/j.jnoncrysol.2013.07.014

52. Kerisit S, Pierce EM, Ryan JV. Monte Carlo simulations of coupled diffusion and surface reactions during the aqueous corrosion of borosilicate glasses. J Non Cryst Solids. 2015;408:142-49. https://doi.org/10.1016/j.jnoncrysol.2014.07.020

53. Kerisit S, Du J. Monte Carlo simulation of borosilicate glass dissolution using molecular dynamics-generated glass structures. J Non Cryst Solids. 2019;522:119601. https://doi.org/10.1016/j.jnoncrysol.2019.119601

54. Kerisit S, Mahadevan T, Du J. Patchy particle model of hydrated amorphous silica. J Non Cryst Solids. 2021;556:120555. https://doi.org/10.1016/j.jnoncrysol.2020.120555

55. Jan A, Delaye J-M, Gin S, Kerisit S. Monte Carlo simulation of the corrosion of irradiated simplified nuclear waste glasses. J Non Cryst Solids. 2019;519:6-13. https://doi.org/10.1016/j.jnoncrysol.2019.05.025

56. Jan A, Delaye J-M, Gin S, Kerisit S. Molecular dynamics simulation of ballistic effects in simplified nuclear waste glasses. J Non Cryst Solids. 2019;505:188-201. https://doi.org/10.1016/j.jnoncrysol.2018.11.021

57. Reiser JT, Parruzot B, Gin S, Bonnett JF, Smoljan CS, Seymour LM, et al. Effects of Al:Si and (Al plus Na):Si ratios on the static corrosion of sodium-boroaluminosilicate glasses. Int J Appl Glas Sci. 2022;13(1):94-111. https://doi.org/10.1111/ijag.16109

58. Taron M. Simulation à l'échelle nanoscopique du transport réactif: application à la dissolution des verres nucléaires. PhD Thesis. 2022.

59. Damodaran K, Gin S, Naranayasamy S, Delaye J-M. On the effect of Al on alumino-borosilicate glass chemical durability. NPJ Mater Degrad. 2023;7:46(1). https://doi.org/10.1038/s41529-023-00364-3

60. Neeway JJ, Reiser JT, Kerisit SN, Reyes RA, Daniel RC, Smith GL et al. Effect of composition on the corrosion behavior of 24 statistically-designed alkali-borosilicate waste glasses. J Nucl Mater. 2023;586:154674. https://doi.org/10.1016/j.jnucmat.2023.154674

61. Damodaran K, Delaye J-M, Kalinichev AG, Gin S. Deciphering the non-linear impact of Al on chemical durability of silicate glass. Acta Mater. 2022;225:117478. https://doi.org/10.1016/j.actamat.2021.117478

62. Gin S, Mir AH, Jan A, Reyes RA, Daniel RC, Smith GL, et al. A general mechanism for gel layer formation on borosilicate glass under aqueous corrosion. J Phys Chem C. 2020;124(9):5132-44. https://doi.org/10.1021/acs.jpcc.9b10491

63. Vienna JD, Crum VJ. Non-linear effects of alumina concentration on product consistency test response of waste glasses. J Nucl
Mater. 2018;511:396-405. https://doi.org/10.1016/j.jnucmat.2018.09.040

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Delaye J-M, Tiwari S, Brun E, Fossati P, Gin S. Advanced Monte Carlo method for simulating glass alteration: Application to aluminoborosilicate glasses. J Am Ceram Soc. 2024;1-19. https://doi.org/10.1111/jace.20167