Ion-exchange mechanisms and interfacial reaction kinetics during aqueous corrosion of sodium silicate glasses

The ion-exchange and associated interfacial reaction mechanisms of silicate glasses are critical in elucidating their aqueous corrosion behaviors, surface modification and property changes, hence have potential impact on both science and technology. This work reports findings of the atomic and nanoscale details of the glass–water interfacial reactions revealed by applying reactive force field (ReaxFF) based molecular dynamics (MD) simulations, from which the key mechanisms of the ion exchange, as well as the kinetics of associated interfacial reactions, are elucidated. It was found that the Na+ and H+ ion exchange can happen between two oxygen ions on a single silicon oxygen tetrahedron or adjacent tetrahedra. In addition, the clustered reaction of two non-bridging oxygens mediated by an adjacent water molecule was also identified. The latter reaction might be the main mechanism of water transport after initial surface reactions that consume the non-bridging oxygen species on the surface. Water molecules thus can play two roles: as an intermediate during the proton transfer processes and as a terminator of the clustered reactions. Statistical analyses were performed to obtain reaction kinetics and the results show that silanol formation is a more favored process than the silanol re-formation within the first 3 ns of interfacial reactions. The results obtained thus shed lights on the complex ion-exchange mechanisms during glass hydration and enable more detailed understanding of the corrosion and glass–water interactions of silicate glasses.


INTRODUCTION
Glass as an ancient yet modern material has transformed our lives in many ways, with applications ranging from window panes, cookwares, and containers in our daily lives to high technology applications such as display, optical communication, biomedical devices, and nuclear waste disposal 1,2 . In many of these applications, glass will be in contact with water, either the liquid or vapor phase, in their lifetime and the property changes due to these interactions are thus of critical scientific and technological importance. In practice, various levels of interactions exist ranging from top surface modification, such as during glass processing and packaging to long term corrosion, such as in pharmaceutical and nuclear waste disposal. Recent studies have shown that water incorporation can effectively improve the strength and lower the brittleness of glass, increasing the composition space of glass materials [3][4][5] . However, understanding glass-water interactions is challenging due to the complex interfacial reactions that happen in the nanometer scales and varying time scales ranging from nano-seconds to years, which makes direct observations very challenging 1,6 . Advanced characterization methods, such as NMR, SIMS, XPS, FTIR, and Raman spectroscopy, coupled with isotope substitution and other techniques, have provided insights to the water-glass reactions [5][6][7][8][9][10][11] . Recent developments of first principles and classical atomistic simulations 12,13 , especially those with reactive potentials 14,15 , have enabled detailed investigations of the surface and interfacial reactions with atomic level resolution [16][17][18] .
Glass dissolution is a complex process that depends on both intrinsic material properties, such as glass composition, structure, surface condition, thermal history, and extrinsic parameters such as temperature, solution pH value, and solution composition 6,7 . The commonly accepted stages of corrosion of multicomponent glasses, based on intensive studies of corrosion of boroaluminosilicate glasses for nuclear waste disposal 19 , include the initial fast ion-exchange reactions, followed by hydrolysis reactions of bridging oxygen bonds between network formers, then the formation of porous silica rich amorphous gel layer by reorganization of the remaining network or through the precipitation of the dissolved species 1 . When the condition is correct, this silica rich gel layer can serve as a passivation layer that separates the bulk glass from the solution, thus provide a low residual dissolution rate 6 . The actual dissolution behavior can deviate from these steps depending on the glass compositions and solution conditions. Hydration, hydrolysis, and ion-exchange are the three fundamental processes during glass dissolution 19 . The three processes can all present in the corrosion of complex glasses and the rate of one affects the other two. In modified silicate glasses, due to lack of channels or openings of molecular water diffusion, hydrolysis and re-condensation are the main mechanisms of water transport. The ion-exchange process leads to the selective leaching of silicate glasses and accounts for the initial fast dissolution of glass (as discussed earlier) 6 . The rate and extent of ion-exchange depend on the solution chemistry (such as the pH value) and glass structure (network structures and type of modifier cations). In addition, temperature and pressure-induced modifications on local structure can also affect the reaction processes 3,20 . Despite the importance of these reactions, detailed atomistic mechanisms are still lacking due to the complex reactions in the amorphous matrix. The goal of this work is to understand the fundamental reaction mechanism such as ion-exchange and hydration using sodium silicate glass as a model system and archetype for more complex multicomponent silicate glasses. To accomplish this, the latest reactive potentials 14 were employed to perform atomistic simulations of sodium silicate glass-water interfacial reactions using parallel high performance computing. Based on the dynamic simulation data, both the mechanisms and kinetics of these reactions are evaluated.
In addition to the experimental techniques, various computational simulation methods covering wide time sacles, such as ab initio molecular dynamics (AIMD), classical molecular dynamics (MD), and Monte Carlo (MC) simulations, have been applied to investigate the glass structures and properties as well as glass-water interfacial reactions 18,[21][22][23][24][25][26] . Among these methods, the MD method is widely used in glass simulations since it can provide atomistic level details of medium-range-order structures with relatively lower computational cost, whereas AIMD is accurate and capable of describe chemical reactions but very expensive to investigate glassy materials, where a large system size is usually required due to the complex structure features. The early MD simulations of the sodium silicate glass systems dated back to Soules's report about the structures of sodium silicate glasses in 1979 27 . Detailed structures and comparisons with experimental results have been reported by subsequent studies including Garofalini and Levine 28 , Du and Cormack 29 , and Vessal and coworkers 30 . The various aspects of the glass structures including pair distribution function, bond angle distribution function, ring size distribution, Q n speciation, etc. have been investigated. Although the structures of the bulk silicate glass have been well investigated through MD simulations, the glass-water interfacial reaction is more challenging due to lack of available empirical potential parameters that are capable to describe the glass-water chemical reations. Du and Cormack 31 proposed a set of potentials to investigate the hydroxylation process of the silica glass surface. However, this set of potentials is not able to simulate water dissociation reactions due to the different charges of oxygen ions used in glass and the hydroxide group. The CLAYFF force field developed by Cygan, Liang, and Kalinichev 32 was originally designed for hydroxide, oxyhydroxide, various clay phases, and more complex systems including the surfaces and interfaces associated with the clays and cementrelated phases. However, water in this set of potentials is also not dissociable hence not able to fully describe glass-water interfacial reaction. Mahadevan and Garofalini 33 proposed another set of potentials to simulate the water dissociation, and recently this set of potentials has been extended to systems including H, O, Si, and Na 15 . This extension enables the simulation of sodium silicate glass-water interfacial reactions. Another potential that can simulate the water dissociation and silicate glass structures is the reactive force field (ReaxFF) proposed by van Duin and coworkers 34 . This potential is a bond order-based empirical potential fitted to the data form first-principles calculations, which can simulate the chemical reactions and bond formation/breakage. Recently, refinement of the ReaxFF parameters for the Na-Si-O-H systems by incorperating extensive first principles data of bulk sodium silicate crytal phases, sodium ion reaction diffusion energetics, and glass-water intercctions have been reported 14,35 . This refined potential enables us to simulate the interfacial reaction between sodium silicate glass and water. Comparing with accurate but computationally much more expensive AIMD simulations, MD simulations with the ReaxFF can generate glass models with similar structural features for sodium silicate glasses 14 , while being able to access larger systems and longer simulaiton times. With the recent Na-Si-O-H potentials, several following up studies focused on glass-water simulations have been carried out 16,36 and these studies showed that it is able to provide detailed understandings about the glass-water interfacial reactions.
Although most of the practical glass compositions are multicomponent, sodium silicate glass serves as an archetype of these glasses and a detailed understanding of sodium silicate glass-water can provide insights of glass-water interfacial reactions of more complex glass systems. In our earlier studies, we investigated the initial stage of water-sodium silicate glass interactions using the ReaxFF based MD simulations 16 . It was found that the reactions can be divided into three stages, depending on the location relative to the top glass surface. Different from silica surfaces, it was found that reaction with nonbridging oxygen (NBO − ) bonded to a silicon to form silanol (Si-OH) groups is the main reactions. Deeper into the surface, proton hopping between a Si-OH and Si-NBO − groups serves as the main means of water transport in sodium silicate glasses. In between the two regions, both reactions exist 16 . This has helped to answer questions, such as what is the main reactions in sodium silicate glasses when in contact with water and what is the mechanism of water transport in these glasses. However, deeper understanding of glass corrosion requires the fundamental mechanism of the ion-exchange process. There are several mechanisms proposed on how water transport into silicate glasses. For silica glass, Tomozawa and coworkers proposed that molecular water diffusion is the dominant mechanism of water transport in silica glass 37,38 . On the other hand, a study from Sterpenich and Libourel proposed that the ion-exchange process between hydrogen and network modifier ions such as alkali ions plays the major role of the hydration process of the glass 39 . Although the exact mechanism remains controversial, advanced simulation work like in this work can shed light on the fundamental steps of glass corrosion.
The purpose of this work is to use the reactive potential (ReaxFF) based MD simulations to study the sodium silicate glass-water interfacial reactions and find the detailed mechanisms during the ion-exchange process, which exhibits the silanol formation and re-formation processes. Furthermore, the statistics of the interfacial reactions are studied by calculating the concentrations of involved species and frequency of each ionexchange process as a function of simulation time.

RESULTS
The ion-exchange process in silicate glass dissolution The ion-exchange process during silicate glass dissolution is mainly proton (H + )-alkali ion, e.g., sodium (Na + ) exchange. The process can be complex due to the different local environment around the exchanged atoms and short time scales that make direct observation challenging. As ion exchange is a dynamic process, we chose to observe trajectories from MD simulations by monitoring two important distances: (a) distance between the hydrogen and oxygen atoms in a silanol group, and (b) distance between the non-bridging oxygen (NBO − ) atom (associated with the same silicon) and the nearby Na + in the initial configuration. The changes of these distances provide a detailed understanding of ion exchange at the sodium silicate glass-water interfaces. Thereby, we devoted to track a typical reaction site over the time scale of nanoseconds. The initial and final configurations of the monitored units are shown in Fig. 1a and b, respectively.
Evolution of the distances between the selected oxygen and H/ Na atoms have been monitored and plotted as a function of simulation time (as shown in Fig. 2a). In the beginning of the process, the H + belongs to a water molecule in the bulk water region hence the initial distance of O-H is more than 15 Å. This water molecule moves inside the water box until it finally enters the glass region at around 0.15 ns. It interacts with the nearby SiO 4 tetrahedral unit and possible proton transfer happens as indicated by the fluctuation of the O-H distance (Fig. 2a). This journey continues for another 2.1 ns, until the distance decreases to 4-5 Å as the water molecule diffuses close to the reaction site where NBO − and Na + are located. The O-Na distance plot shows a quite similar trend: initially, the Na + is away from the monitored oxygen atom; after 0.1 ns, it comes close to the NBO − and vibrates around it with a distance of 2.3 Å, which is a typical Na + -NBO − distance in silicate glasses. This vibration process continues for about 2.2 ns, until the next stage of the reaction happens. After a residence time of around 2.3 ns, significant changes occur in both O-Na (increase) and O-H (decrease) distances, indicating an ionexchange process happened between the monitored Na + -H + ion pair. After the reaction, a silanol group (Si-OH) is formed and the group is stable, whereas the sodium atom leaves the reaction region, after staying close by at around 3 Å for a certain period of time (around 0.4 ns). Observation of this process is provided in the visualization of the trajectory shown in the Fig. 2b, where colored Na + and H + ions showing their locations as a function of time (red for 0 ns and blue for 3.0 ns). The top panel of Fig. 2b shows that the H + ion diffuses into the glass and finally stays there (as a silanol group). The Na + ion moves around in the glass matrix in the beginning, and after the ion exchange it diffuses toward the surface and glass-water interface (illustrated as the dark blue arrow in the bottom panel of Fig. 2b).

Ion exchange during silanol formation
The ion-exchange process in the glass-water interfacial reaction is a complex process and involves multiple species (reactants) and dissociation/formation of water molecules. The ion-exchange process can also be observed in key surface reaction mechanisms, such as silanol formation and re-formation processes. To understand further details of the ion-exchange process during the silanol formation process, the trajectories between 2.0 and 2.8 ns  Representative Na + and H + ion-exchange reactions at the sodium silicate glass-water interface. a The configuration before the ionexchange process, showing an NBO − and Na + nearby; and b the configuration after the exchange and formation of silanol group. The O, Si, H, and Na are colored by red, orange, cyan, and purple, respectively. The larger spheres represent for the monitored atoms, and distances between O and monitored Na/H atoms are also shown.
(see Fig. 2) were investigated and results plotted in Fig. 3. Between 2.0 and 2.3 ns, the local environment of the reaction site consists of a Si-NBO − , an adjacent Na + ion associated with the NBO − , a nearby hydroxide group and a nearby silanol group (Fig. 3a). Both the hydroxide group and silanol group result from the reaction that the water molecule, where the monitored proton initially stays, undergoes with a nearby Si-NBO − during its journey from bulk water to the reaction site, as indicated by the shorter O-H distance in this period in Fig. 2a. This local structure exists for about 0.7 ns with minor changes due to the vibration, until the proton from the adjacent silanol group transfers back to the hydroxide group at around 2.3 ns, forming a water molecule (Fig. 3b). This change breaks the temporary stability of the local structure, which triggers the next step of the ion-exchange process. The water molecule then dissociates and one proton of it reacts with Si-NBO − to form a silanol (Si-OH) group, leaving the hydroxide group (OH − ) behind while the Na + -NBO − distance increases in the process (Fig. 3c). Eventually, a stable silanol group is formed and the ion-exchange process is completed while the Na + (and OH − ) leave the focused region for further possible reactions (Fig. 3d).
Ion exchange in the self-exchange silanol re-formation process Two different reactions at the sodium silicate glass-water interfaces have been revealed in our previous work 16 : (1) silanol formation controlled by the water diffusion and (2) silanol reformation controlled by the ion-exchange process, with the later reaction being more complicated. To better understand the details of this reaction, configurations at different reaction stages are monitored. Firstly, we observe proton transfer between NBO − atoms in the same silicon tetrahedron unit during the silanol reformation process (Fig. 4). This process starts from the configuration containing a silanol group (Si-O bot H A ) and a hydroxide group (Fig. 4a). The H A then dissociates from the silanol group and forms a water molecule together with an adjacent hydroxide group (Fig. 4b). After a certain length of simulation, the water molecule dissociates and the resulting proton (H A ) forms bond with the other NBO − atom of this [SiO 4 ] tetrahedron unit (O top ) to form a new silanol group (Si-O top H A as shown in Fig. 4c). Similar process can also happen between a water molecule and silicon units with two NBO − atoms, and the two hydrogens can switch their positions with each other (either in the silanol or hydroxide group). This process is called the self-exchange silanol reformation process: the proton transfer takes place within a single silicon tetrahedron. Thereby, we can conclude that Q 2 silicon (with two NBO − atoms) is very reactive and a possible transit point of water moving into glass inside.
In addition to the proton transfer, the self-exchange silanol reformation also involves adjacent Na + ions and subsequent ionexchange processes with the protons. To further investigate these processes, the Na + ions around the two non-bridging oxygen atoms (named as O top and O bot ) have been identified and important distances during these ion-exchange processes were measured, as shown in Fig. 5. Four Na + ions, two for each NBO − , are involved in the processes: Na1 and Na2 are around O bot , and Na3 and Na4 are around O top , respectively ( Fig. 5a and b). In addition, the protons from water molecule are denoted as H A and H B , and the whole reaction can be divided into four stages, as shown in Fig. 5c and d (colored by yellow, green, cyan, and gray, respectively).
In the first stage, the hydrogen atoms are in a water molecule, which travels among the water box for about 0.1 ns and then enters the glass matrix (yellow region of Fig. 5c and d). At the same time, the Na + becomes around the target NBOs as the distances between Na + and NBOs become around 2.1-2.3 Å (the bottom panels).
The second stage starts at around 0.3 ns, and the configuration, which is right before the reaction happened is shown in Fig. 5a. As the water molecule reaches the reaction site, it first reacts with the O bot atom. Then the H A together with O bot form a silanol group confirmed by the short H A -O bot distance, which is around 1.0 Å (green region, Fig. 5c, top). At the end of this region, the H A and H B frequently switch the position with each other undergoing a proton exchange process, suggesting the system becomes unstable and further process will take place. On the other hand, the Na3-O bot and Na4-O bot distances also increase starting from 0.3 ns until the end of stage II (Fig. 5c, bottom), indicating that these Na + ions and the protons from the water molecule undergo ion-exchange processes.    Fig. 5d). The Na3-O top becomes longer than the previous stage during the same time range, indicating the Na3 is undergoing a complete ion-exchange process. However, the Na4-O top distance shows a quite different trend, that it is apparently shorter during 1.28-1.70 ns and later gradually increases. This is due to the proton exchange between H A and H B during 1.28-1. The fourth stage starts at around 2.3 ns and the configuration right before the stage is similar to Fig. 5b. The two protons left the reaction site as suggested by the increase of H-O distances (gray region, Fig. 5c and d, top). The distances between Na + and O show more diversity: Na1, Na3, and Na4 move back to the associated NBOs after traveling around for about 0.2 ns, while the Na2 initially stays close to the O bot but later moves away from the reaction site (gray region, Fig. 5c and d, bottom). The exact processes in this stage are quite complicated since more water molecules and other surrounding units are involved as the reaction continues.
In the final configuration at 3 ns (Fig. 6), the H A atom is found to be a part of nearby silanol group (Si S -OH A ), while the H B remains in the hydroxide groups dissociated from the original water molecule. Another water molecule (H C -O N -H D ) from the environment reacts with O bot forming a new silanol group, while the remaining hydroxide group (O N H D − ) attracts Na1. The Na3 and Na4 remain around the associated oxygen atoms due to limited local void space in the bulk side of glass, while Na2 moves away from the reaction site as more available space exist in the surface side of glass. Na2 is eventually captured by a nearby NBO − , waiting for another possible ion-exchange process. It is worth noting that the sodium ions stay near the reaction sites as the ionexchange process is going on. This is due to the residual time after initial reaction, as the time after initial dissociation and before both protons (A and B) leave the reaction site is less than 1 ns. According to an earlier work by Hahn and van Duin 36 using the same ReaxFF potential of surfaces of sodium silicate glasses, sodium ions tend to reside in the local area after initial dissociation from the reaction site. The maximum residence time of the Na + ions can be up to 1 ns before its distance from the reaction site exceeds 5 Å after initial dissociation. This can be explained by the fact that sodium ion diffusion in silicate glasses is usually considered through a process similar to the vacancy mediated mechanism in crystals 40,41 . In addition, the void space surrounding the reaction sites is quite limited, which results in a much slower movement of Na + ions than protons which are much smaller in size. Therefore, the Na + ions are still traveling locally while H + comes into the reaction sites. Although Na1, Na3, Na4 do not move much, the Na2 is about 4 Å away from the O bot . This can be an indication of the Na + /H + ion exchange, although clearer observations will need longer simulations.
Ion exchange in the adjacent-exchange silanol re-formation process Another type of the silanol re-formation process is proton transfer between NBOs in different silicon tetrahedron units. This process starts with a configuration containing a silanol group and a Si-NBO − unit (Si-O L H A and Si-O R − , Fig. 7a). After vibration of the silanol group for certain length of simulation, H A becomes close to O R thus shared by the two NBOs, forming a clustered linkage Si-O L -H A -O R -Si (Fig. 7b). The hydrogen atom is finally captured by the O R atom and forms a new silanol group, while the previous silanol group becomes a Si-NBO − unit (Fig. 7c). This process is called the adjacent-exchange silanol re-formation process: the proton transfer takes place between two nearby silicon tetrahedra.
The ion-exchange process also takes place in the adjacentexchange silanol re-formation process and this is shown in detail in Fig. 8. The process involves the proton (H A ) and two Na + ions  (Na1 and Na2), where Na1 is around the left NBO − (O L ) and Na2 is around the right NBO − (O R ) before this process, respectively. The H A atom initially belongs to a water molecule, which travels within the water box for the first 0.2 ns, as indicated by the large H-O distances in the yellow region of Fig. 8c. This water molecule then enters the glass region and interacts with nearby units during the path in the next 2.1 ns (green region, Fig. 8c), until it finally comes to the reaction site as illustrated in Fig. 8a. The adjacent-exchange silanol re-formation then begins, as the H A -O L and H A -O R distances both significantly decrease and become close to each other for a long time (cyan region, Fig. 8c). The correlation between these two distances confirms that H A in this region is moving back and forth between O L and O R during this stage, as illustrated in Fig. 7.
When the proton transfer process starts, the subsequent ionexchange processes occur. After H A dissociates from the water molecule and moves into the center site of O L and O R , the Na1-O L and Na2-O R distances both increase (cyan regions, Fig. 8d and e). The Na2-O R distance exhibits a fluctuating behavior due the reversibility of silanol re-formation, e.g., it increases when H A is associate with O R and decreases when H A goes back to the O L . In addition, the hydroxide group dissociated from the initial water molecule, O W -H B in Fig. 8a and b, accompanies with the Na2 ion starting from 2.7 ns, which slightly increases the overall distance of Na2-O R . On the other hand, the fluctuation of Na1-O L distance is smaller due to the change of environment, e.g., another nearby hydroxide group (the O N H C group in Fig. 8b) comes close to Na1. The attractive interaction with O N H C group slightly increases the Na1-O L distance at the beginning of the approaching process, and it eventually attracts the Na1 ion and left the reaction region together after 2.7 ns. It is worth noting that the sudden increase of the Na1-O L distance and the decrease of Na2-O W occur at the same time (at around 2.7 ns) even though Na1 and Na2 are sufficiently separated.
The process can be explained by following steps: at first, the O N H C group approaches to and finally associates with Na1 (2.30-2.42 ns, in Fig. 8d). Then Na1 and O N H C group leave the reaction region from the surface side, as there is more void space in the surface side than in the bulk glass side, which allows those species to move away. The initial speed of the Na1 removal is slow during 2.42-2.70 ns, and this process results in a less negative charge of O L (from −0.967 at 2.42 ns to −0.814 at 2.70 ns). Evolution of the O L charge in this time range is plotted in Fig. S1 in Supplementary information. This change makes the attraction between O L and H A weaker, resulting in a close H A -O R distance thus an increase of the Na2-O R distance. On the other hand, the attraction between Na2 and O W H B group from the initial water decreases the Na2-O W distance. The sudden increase of the H A -O L distance at 2.7 ns significantly enhances this process, which leads to the Na2-O W distance finally reaches to around 2.1 Å as shown in Fig. 8e.

Continuous reactions that led to multiple silanol formation
In the mechanisms reported earlier, it leads to the formation of one silanol group. A more complicated reaction involving proton transfer processes and step reactions that led to two silanol formation have also been observed. This is schematically shown in Fig. 9f. This process starts from the configuration with two silicon tetrahedron units and a nearby water molecule (Fig. 9a). The initial distances between hydrogen and both oxygen atoms are more than 5 Å. At around 0.01 ns, the water molecule comes close to the top silicon tetrahedron, one of the hydrogen atoms then interacts with the top oxygen atom, forming a Si-O-H-O-H cluster (Fig. 9b). This hydrogen atom later dissociates from the water molecule and, together with the top silicon unit, forms a silanol group (Fig. 9c). The newly formed silanol group oscillates between two silicon units and can be any position between configurations c and d. In addition, this process (a to c) can be reversed depending on surrounding species, e.g., water molecule and hydroxide group. For example, configuration b is also observed at around 0.12 ns, indicating that the silanol group can interact with the nearby hydroxide group. This reversible process continues for about 0.35 ns, until the hydroxide group captures another hydrogen from outside and forms a new water molecule (Fig. 9d). The brand-new water molecule has less attraction comparing with the bottom Si-NBO − unit, thus the top silanol tends to get closer to the bottom Si-NBO − unit, as shown in configuration Fig. 9d. After about 0.02 ns, the final proton exchange process immediately occurs: the bottom NBO − finally captures the hydrogen atom from the top silanol group, forming a bottom silanol group and leaving the top oxygen atom as a NBO − (Fig. 9e). The top Si-NBO − then reacts with the water molecule to form another silanol group, with a hydroxide group dissociated from the water molecule (Fig. 9f). After this process, the bottom silanol group becomes stable as the distance between the monitored hydrogen and bottom oxygen is almost constant.
The whole process is controlled by the surrounding species of the silicon units: the hydroxide group introduces the reversibility of the silanol formation, while the water molecule terminates the reversible process by converting the available NBO − to another silanol group. The ion exchanges of this mixed process take place during a very short period, but it is even more complicated as both the silanol formation and the adjacent-exchange silanol re-formation processes are involved.

Statistics of the interfacial reactions
Different processes contribute to the glass-water interfacial reactions as discussed in the above section. The overall reaction speed of the interfacial reactions, therefore, can be varied due to the different dominant process in various reaction stages. Both the silanol formation and silanol re-formation processes contribute to the progress of sodium silicate glass-water interfacial reactions, and three different reaction stages are observed in the previous study 16 . In the earlier reaction state (mainly takes place at the near-surface region of the glass alteration layer), the majority part of the reactions is dominated by the water diffusion together with subsequent silanol formation process which can be described as: This process is a fast reaction as the concentration changes of both reactants and products are significant within the first 0.1 ns, as shown in Fig. 10. At the same time, the Na + ion, which initially associated with the Si-NBO − unit, also participates in this process as it undergoes an ion-exchange process with the proton from the water. When the reactions continue in deeper glass matrix, in addition to the silanol formation process (Eq. 1), some of the protons from the newly formed silanol groups react with inner NBOs to proceed the silanol re-formation process, which can be expressed as: Si cs --OH þ Si cb ---NBO À ! Si cs ---NBO À þ Si cb ---OH (2) where the Si cs and Si cb represents for close-to-surface and closeto-bulk silicon units, respectively. This process can be conducted through either self-exchange or adjacent-exchange proton transfer process, as illustrated in Figs. 4 and 7, and it also involves ion process. The stoichiometry of this process is 1:1:1:1, indicating the concentrations of both the reactants and products keep as constants during this process. Both the silanol formation (Eq. 1) and re-formation (Eq. 2) processes contribute to the second stage of the interfacial reactions at the middle of the glass alteration layer. The third stage of the reaction (mainly takes place at the near-bulk region of the glass alteration layer) is dominated by the silanol re-formation process, since the water diffusion and subsequent silanol formation process at this region are rare. Since the silanol re-formation process is reversible until terminated by the other hydrous species, e.g., water molecule, the reaction speed is relatively slower at this stage. As the reaction is a dynamic process, the boundaries of these stages in the alteration layer of the glass-water interface are not fixed; in fact, all the three boundaries tend to propagate deeper into glass with different speed as reaction time increases.

Fig. 9
Continuous reactions that led to the formation of two silanol groups in adjacent tetrahedra involving proton transfer.
Configurations a-f are the snapshots at different reaction stages (with the time labeled). The top side of each configuration is closer to the glass-water interfacial boundary.
The reaction rate can be evaluated by the concentration change of either reactants or products, and it can help determine when and which of the stages starts. Based on the definition of the reaction rate, the rate value R can be obtained as decreased concentration of a reactant or increased concentration of a product in a unit of time. Although multiple processes exist in each stage of the glass-water reaction, only the silanol formation process (expressed by Eq. 1) changes the reactant/product concentration. Both self-exchange and adjacent-exchange silanol re-formation processes control the overall reaction progress since these ion-exchange processes go back and forth due to their reversibility. However, they have minor effects on the concentration changes of the reactants/products, only when new NBOs are formed in the water-rich region that the subsequent cluster reactions take place. Therefore, the reaction rate of the overall interfacial reaction can be approximately predicted by calculating the concentration changes of reactants/products of the silanol formation process. It thus can be calculated as shown in the following equation: Two stages with different rates of the water concentration change can be observed according to Fig. 10: In the first stage, the concentrations of both H 2 O and NBO − undergo a significant change (the first 0.25 ns of the reaction). This indicates that the dominant processes at early reaction stage is mainly molecular water diffusion (or direct contact between the water box and bulk glass at the very beginning) and the subsequent silanol formation process, and this process has a high reaction speed. The silanol reformation processes in this region are negligible comparing with the silanol formation process. In the second stage, the much slower reaction rate is attributed to the fact that molecular water diffusion becomes less due to less void space than in the previous stage. At the same time, the self-exchange/adjacent-exchange silanol reformation processes increase after initial silanol formation in the first stage. The mixture of both silanol formation and selfexchange/adjacent-exchange silanol re-formation processes (like the continuous reaction shown in Fig. 9) leads to the reduced overall reaction speed. This is because the proton involved in the silanol re-formation process moves back and forth, which reduces the speed of the proton propagation into glass. As a result, the possibility of new NBO − formation in the water-rich region is reduced, which leads to a lower frequency of the subsequent silanol formation process. One thing to mention is that the transition between the two stages is a gradual process rather than an immediate one, as the alteration layer is continuously developed while the reaction is going on. It will also get affected by the local structure features, such as available void space and the local concentrations of the involved species. It can be expected that there would be another stage that with much slower changing rate of the concentrations than the previous two stages if the simulation lasts much longer until reaching a constant reaction rate. In this case, the major process in the third stage would become the silanol re-formation process. The silanol formation process in this stage would be rare since the hydration process of the established alteration layer is almost completed in previous stages. It would only happen when new NBOs are formed in the water-rich region through certain proton transport processes (as shown in Figs. 4 and 7). This process would allow more protons or water molecules to come into the alteration layer and compensate the newly formed close-to-surface NBO − . This stage, therefore, would result in the further growth of the glass alteration layer into the glass matrix. With longer time of simulations, this stage may be accessible. It might be useful to extract the reaction constant at different temperature from which the reaction energy barrier can be calculated for these reactions. However, due to coupling of multiple reaction mechanisms, it is difficult to obtain reliable reaction rates and the barriers.
To understand the dominant reaction in each stage during the interfacial reaction, the frequencies of the silanol formation process and two types of silanol re-formation processes were analyzed by counting the number densities of the hydrogens which transit between involved species (H 2 O and silanol group for the silanol formation process, or Si cs -OH and Si cb -OH for the silanol reformation processes) in two consecutive frames (time interval of 10 ps). The numbers of reactions happening at certain time with a time scale of 10 ps can thus be obtained, as shown in Fig. 11. The silanol formation has the highest frequency among the three processes. This is followed by the adjacent-exchange silanol reformation process, while the self-exchange silanol re-formation process has the lowest frequency of all due to the limited number of the Si-(NBO − ) 2 units in the glass structure. In addition, at the very early stage of the reaction (<0.025 ns), the number of the silanol formation process is significantly higher than the value in later reaction. This is due to the initial silanol formation process happens directly between the water molecules and the surface NBOs, which involves minimal amount of water diffusion. Water diffusion then plays a more important role in further reactions and it make these reactions more difficult to take place. Moreover, higher temperature increases the frequencies of all three reactions, especially for the silanol formation process and the adjacent-exchange silanol reformation process. The accelerated reaction leads to a deeper penetration of the proton, which make it possible for the proton to encounter with more Si-(NBO − ) 2 units at deeper glass region. This leads to a slightly increased number of the self-exchange silanol reformation process. The statistics of the accumulated numbers of each process at different reaction temperatures are listed in Table 1. These results suggest that the silanol formation process is the favored reaction comparing with the silanol re-formation processes.
The more complicated continuous reactions (Fig. 9) consist of multiple processes by combining both silanol formation and  re-formation processes. It can have various pathway to achieve the final reaction, many of them can be even more complicated than the one presented in Fig. 9. Therefore, it is difficult to quantify the exact number of the reactions happened.

DISCUSSIONS
The interfacial reactions between sodium silicate glass and water at different temperatures have been simulated using the Na-Si-O-H ReaxFF parameters 14 , and detailed analyses on the ion-exchange processes during the interfacial reaction have been performed. By monitoring the distances between oxygen and the exchanged ions, the time when the exchange process happens can be determined. Configurations during the process have been captured and the related distances have been verified, as shown in Figs. 2 and 3. The ion-exchange process is carefully monitored by taking snapshots of the configurations at different stages, consequently, the sodium leaching process is found, as shown by the arrow in the bottom panel of Fig. 2b. Two different types of the proton transfer process have been noticed: (A) the proton transfer between two NBOs in the same silicon tetrahedron unit, and (B) the proton transfer process between two NBOs belonging to different silicon tetrahedron units. Both processes play an important role in the interfacial reactions between sodium silicate glass and water, and the mixed processes of these two as well as the silanol formation process results in the overall proton propagation into the glass and subsequent glass corrosion. In addition, the ion-exchange processes were observed in both types of silanol re-formation processes. Since Na + ion diffusion in sodium silicate glasses is through a mechanism similar to the vacancy mediated process observed in crystals 40 , and there is a barrier of 0.6-0.7 eV for sodium ion diffusion in these glasses, this makes the local movement a more favored event as compared to the longerdistance diffusion for the Na + ions, especially in glasses with high sodium content like the composition in this study. As a result, Na + ion diffusion is much slower than the proton diffusion in the bulk glass region, which leads to the reversible proton transfer processes (both self-exchange and adjacent-exchange silanol re-formation, Figs. 4 and 7) frequently happened processes. On the other hand, the Na + ions at the surface side have more space to move, which allows them to move far away from the reaction site. This step is always together with a nearby hydroxide group, and it can stabilize the silanol group thus lead the whole reaction go further.
The three types of proton transfer processes, as well as the subsequent ion-exchange processes, are simplified into a twodimensional illustration as shown in Fig. 12.
Moreover, the clustered reaction involving both silanol formation and adjacent-exchange silanol re-formation processes, as shown in Fig. 9, is found as one of the most important processes for the propagation of the hydrogen-related species (silanol groups and water molecules). Since the silanol re-formation process is reversible when it occurs alone, additional hydrogen atoms (or water molecule) will significantly affect the lifetime of the silanol group. Therefore, the local concentration of the water molecules/hydroxide groups will change the whole reaction process: the water molecules can either play a role of the intermediate during the singular reaction or be the terminator of the cluster reaction.
As the practical dissolution reactions go further, the release of Na + and OH − will increase, which will enhance the pH value of the solution. Therefore, the reaction mechanisms in the system from this work can be slightly different from those in the glass-acidic solution systems due to the additional H + /H 3 O + concentration. The majority difference affected would be the silanol group, as the silanol formation process (Eq. 1) would release hydroxide groups and changes the local pH value continuously as the reactions develop. The H + /H 3 O + will consume the hydroxide group thus accelerate the silanol formation process and thus the overall reaction.
The overall interfacial reaction might be divided into two stages with different reaction rates, according to the changes of the H 2 O and NBO − concentrations, while an additional stage is expected in longer time simulations which is much closer to the real situation in experiments. The different dominant processes in each stage lead to the differences in reaction rate. The statistics analyses on the frequencies of the three ion-exchange processes show that the silanol formation process has the highest frequency over the whole simulation, followed by the adjacent-exchange silanol re-formation process. The self-exchange silanol re-formation process has the smallest frequency due to the limited existence of the Si-(NBO − ) 2 units in the glass. Higher reaction temperature will increase the frequencies of all three ion-exchange processes, which leads to an accelerated overall reaction.
In conclusion, the sodium silicate glass-water interfacial reactions have been studied at different temperatures with recently refined Na-Si-O-H ReaxFF potentials, and mechanistic understanding of the ion-exchange process of the glass-water interactions have been obtained by detailed analyses of the trajectories of the MD simulations up to 3 nano-seconds. The ion-exchange process was found to be mediated by water molecules and two types of proton transfer processes were observed during silanol re-formation, which serves as a means of water transport in the glasses. The Na + and H + ion exchange is found to be a fundamental step in both types of silanol re-formation processes, and diffusion of Na + ion from the reaction site (e.g., toward the surface) can stabilize the silanol group (as a reaction product) and increase the possibility of further reactions deeper in the bulk glass. The clustered reaction process in which a proton is stabilized in the middle of two NBOs for a certain period of time has been observed as an important step of the proton propagation after the initial silanol formation process. Two roles of water molecules played during the ion-exchange reactions have been identified: it can serve as an intermediate product during the silanol re-formation, or as a terminator of the clustered reaction process. Importantly, it was found that the local glass structure will significantly affect the direction of the reversible reaction, which leads to different speed of the glass corrosion. The statistical analyses on the frequencies of the ion-exchange processes revealed that the silanol formation process happens most frequently and this is followed by adjacent-exchange silanol re-formation while self- Fig. 12 Schematics of the three types of ion-exchange processes during the sodium silicate glass-water interfacial reactions. a Silanol formation, b self-exchange silanol re-formation, c adjacentexchange silanol re-formation. The small green arrows indicate the target directions of the proton movement. The purple, cyan, red, orange circles represent for the sodium, hydrogen, oxygen, and silicon atoms, respectively. exchange silanol formation is the least frequent, during the early stage of the interfacial reactions. It is expected as the reaction progresses, the silanol re-formation reactions will become dominating. The adjacent site silanol re-formation happens more frequently than the self-exchange silanol reformation, both of which will play more important role with longer reaction time and as the reactants move further into the bulk glass region. These reactive dynamic simulations hence provide direct atomic level evidence of ion exchange processes happened during aqueous corrosion of silicate glasses. The results also shed lights on the fundamental steps of the ion-exchange processes of sodium silicate glass-water interfacial reactions that will lead to deeper mechanistic understanding of corrosion of silicate glass materials.

Simulation details and glass structure generation
The glass composition investigated in this paper is 22.8 Na 2 O-77.2 SiO 2 with a density of 2.365 g/cm 3 from experiment 16 . The whole simulation of this work is divided into three parts: the bulk glass and surface generation, the bulk water generation, and the glass-water interfacial reactions. This procedure followed the previously adopted protocol 16 , and details of the model creation can be found below. The total interaction energy in ReaxFF potential consists of eight terms, including the Coulomb energy, the van der Waals energy, bond energy, valence angle energy, torsion angle energy, over-coordination penalty energy, under-coordination penalty energy, and lone-pair energy 14 . These energy terms enable the ReaxFF potential to monitor bond formation and breakage during a chemical reaction, thus can be applied to investigate the glass-water interfacial reaction. The potential parameters were obtained by fitting to the first principle data of bulk crystalline and glass sodium silicate structures, sodium ion diffusion energy barriers, sodium ion solvation and Na-OH dissociation in water 14 . Nosé-Hoover thermostat and barostat 42,43 were used in the MD simulations, and the damping parameters were chosen as 100 and 1000 times of the timestep used (1.0 fs for glass formation and 0.1 for glass-water reaction), respectively. The simulations were carried by the LAMMPS simulation package 44 .
The bulk glass with 8001 atoms were obtained using simulated melt and quench process with partial charge pairwise potential 29 . The initial glass structure was relaxed following the procedures described in earlier work 16 with recently refined ReaxFF parameters 14 . After the glass formation, two surfaces were cleaved by elongating the simulation box in one direction, and the surfaces were subsequently annealed at 300 K for 100 ps with the ReaxFF potential 14 . The glass formation and surface relaxation processes were under an isobaric-isothermal ensemble (NPT). The timestep used in the processes of glass formation and surface creation was 1.0 fs.

Glass-water interface models and interfacial reaction simulations
A water box with 3000 H 2 O molecules was created using the ReaxFF potential with same size of glass surface in both x and y directions, and the total volume was fixed to maintain the density to be 1.0 g/cm 3 . The canonical ensemble (NVT) was used during the relaxation at 300 K. The interfacial model was then created by putting this water box on top of the relaxed glass surface with a small gap (3-5 Å) between glass and water in both surface regions. This process was carefully performed with a smaller timestep of 0.1 fs to ensure accuracy of integration for hydrogenrelated species and avoid unfavorable reactions at the interface.
After the interfacial model was created, the system was relaxed at the reaction temperatures under NPT ensemble for 10 ps. The reaction temperatures were set to 300 K, 350 K, 400 K, and 450 K based on the previous simulation results using similar potential parameters 16 . The usage of higher than experimental temperatures was to achieve faster reactions, so reaction kinetics can be obtained while maintaining the basic physical process of the system. It was shown in our earlier study that too high a temperature (>500 K) will lead to formation of unrealistic byproducts such as molecular oxygen, etc 16 . This initial MD run minimizes the artifact of the initial gaps between water box and glass surface by immediately relaxing the simulation box with reasonable density. The reaction simulation was then further performed under NVT ensemble at the reaction temperatures up to 3.0 ns. The NVT reaction process was used to simplify the analyses that the simulation size is constant; therefore, the positions of atoms can be comparable during the whole simulation. The volume of the sampling block, the x and y directions of simulation box size and the z direction of selected bin size (1 Å), was used in the concentration calculations. Fig. S2 (Supplementary information) shows the pair distribution functions of the sodium silicate glass before and after reactions after 3 ns simulations at 450 K. The overall distribution functions of the glass remain essentially unchanged although modifications of the glass structure happened in the period due to glass-water interactions and reactions. Figure 13 shows the sodium silicate glass-water interfacial model after 3 ns reactions at 450 K. Water has reacted with the top glass surface and pockets of water transport deeper into the bulk glass.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors on reasonable request. Fig. 13 Snapshot of the sodium silicate glass-water interface after reactions for 3 ns at 450 K. The image shows water penetration after the reactions and pockets of water transport deeper into the glass. The balls in red, orange, purple, and cyan represent for the oxygen, silicon, sodium, and hydrogen atoms, respectively. The surface of area containing hydrogen atoms is also presented as the cyan cloud.