Supercritical Water: A Simulation Study to Unravel the Heterogeneity of Its Molecular Structures

Abstract

(1) Background: In the quest to accurately model the radiolysis of water in its supercritical state, a detailed understanding of water’s molecular structure, particularly how water molecules are arranged in this unique state, is essential. (2) Methods: We conducted molecular dynamics simulations using the SPC/E water model to investigate the molecular structures of supercritical water (SCW) over a wide temperature range, extending up to 800 °C. (3) Results: Our results show that at a constant pressure of 25 MPa, the average intermolecular distance around a reference water molecule remains remarkably stable at ~2.9 Å. This uniformity persists across a substantial temperature range, demonstrating the unique heterogeneous nature of SCW under these extreme conditions. Notably, the simulations also reveal intricate patterns within SCW, indicating the simultaneous presence of regions with high and low density. As temperatures increase, we observe a rise in the formation of molecular clusters, which are accompanied by a reduction in their average size. (4) Conclusions: These findings highlight the necessity of incorporating the molecular complexity of SCW into traditional track-structure chemistry models to improve predictions of SCW behavior under ionizing radiation. The study establishes a foundational reference for further exploration of the properties of supercritical water, particularly for its application in advanced nuclear technologies, including the next generation of water-cooled reactors and their small modular reactor variants that utilize SCW as a coolant.

1. Introduction

Alongside hydroelectric power, nuclear energy holds significant potential for reducing global greenhouse gas emissions compared to coal- or gas-powered plants, making it an essential alternative to fossil fuels in the fight against climate change. The importance of this perspective was underscored during the IAEA International Conference on Climate Change and the Role of Nuclear Power, which took place in Vienna, Austria, during 7–11 October 2019 [1]. In recent years, significant advancements have been made in develo** fourth-generation, advanced nuclear reactors, notably supercritical water-cooled reactors (SCWRs) that are central to this study. These reactors are among six innovative designs being pursued for commercial applications by the “Generation IV International Forum” (GIF) research and development collaboration [2,3,4,5,6,7,8,9,10]. The SCWR, a class of high-temperature, high-pressure, water-cooled reactors, operates above the thermodynamic critical point of water—temperatures exceeding 373.95 °C and pressures above 22.1 megapascals (MPa) for light water (H₂O) [11]. With core inlet and outlet temperatures set at 300 °C and 600 °C, respectively, and a nominal pressure of 25 MPa, the SCWR is an extremely energy-efficient system, achieving thermodynamic cycle efficiencies of over ~45%, significantly surpassing the ~28–32% efficiencies of current conventional pressurized water reactors [2,8,9]. This high efficiency offers substantial economic benefits through lower-cost electricity generation. SCWRs are anticipated to be commercially deployed between 2030 and 2040 [2].

Over the past fifteen years, several countries have recognized the potential of SCWR technology and actively engaged in its research. This interest has driven the development of “small modular reactor” (SMR) variants of SCWRs that utilize SCW cooling [12,13], catering to the need for smaller, more flexible reactor designs that large-scale SCWRs cannot meet. These reactors aim to significantly reduce capital costs through modularization and offer potential for enhanced safety features that could minimize the consequences of severe accidents. Noteworthy among such initiatives is the ECC-SMART project (“European-Canada-China-Small Modular Supercritical Water Reactor Technology”), launched in September 2020 [14]. This project marks a major international collaboration, with twenty partners dedicated to advancing these energy systems and integrating them as central elements in their future energy strategies.

In this context, there is an urgent need to identify the optimal chemistry and materials under SCWR conditions [8,15]. This research is especially pivotal for the development of small modular SCWRs (SCW-SMRs), which are being considered for deployment in small remote communities and develo** countries [7].

A major challenge in understanding the water chemistry of SCWRs and their SMR variants arises from limited knowledge about in-core radiolysis and the specific radiolytic species these systems produce. SCWRs operate under significantly higher temperatures and pressures—approximately 300–600 °C and 25 MPa—compared to current generation water-cooled reactors, which typically operate at temperatures of 250 to 330 °C and pressures around 10 MPa [8,15,16]. This distinct operating environment of SCWRs results in water chemistry that is markedly different from that of conventional water-cooled reactors. Bridging this knowledge gap is crucial, as it directly impacts critical factors such as material corrosion and degradation, along with the transport of corrosion products [8]. Thus, addressing this gap is essential for ensuring the long-term viability and safety of SCWRs and SCW-SMRs.

Directly measuring the chemistry within reactor cores beyond the critical point of water is exceedingly challenging, if not impossible. This difficulty primarily arises from the SCWR coolant being exposed to intense neutron and γ-radiation fields during its circulation through the reactor core. These radiation fields play a major role in contributing to the formation of various troublesome reactive oxidizing species, including hydroxyl radicals (^•OH), hydrogen peroxide (H₂O₂), oxygen (O₂, as a decomposition product of H₂O₂), and the superoxide anion/hydroperoxyl radical (O₂^•−/HO₂^•, depending on pH level). These species critically affect the chemical environment, operational efficiency, and aging of the reactor (see, e.g., [8,15,16,17,18,19,20,21,22,23]), potentially accelerating corrosion processes in in-core materials, particularly fuel cladding. This can lead to fuel failures and the release of fuel fragments and fission products into the coolant, which further affects the transport of radioactive materials out of the core into downstream pi** components, increasing radiation exposure to reactor maintenance personnel. Furthermore, understanding the behavior of this water is crucial for predicting and mitigating corrosion processes like embrittlement, pitting corrosion, and stress corrosion cracking [8,24,25]. Addressing these issues is vital for maintaining the integrity of reactor components.

There is still limited experimental information on radiolytic yields and reaction rates in the SCW regime or at temperatures ranging from those of current subcritical water reactors to the critical temperature (see, e.g., [8] and references cited therein). Consequently, theoretical modeling and computer simulations have become essential tools for predicting the effects of water radiolysis under these extreme conditions, and their impact on materials [15,16]. In fact, the current understanding of the potential effects of water radiolysis in SCWRs or their SMR versions relies almost entirely on such models [8].

Water at high temperatures exhibits structural differences from water at ambient conditions. For instance, our recent molecular dynamics simulations demonstrated that liquid water undergoes a significant structural change around 150 °C [26]. Interestingly, this structural transformation could explain some apparent anomalies observed experimentally in the radiolysis of water at elevated temperatures. These include the thermalization distance of secondary subexcitation electrons produced in large quantities in the radiolysis of water and the rate constant for the bimolecular reaction of two hydrated electrons, both of which experience a sudden and sharp drop near 150 °C [26].

In the quest to accurately model the radiolysis of water in its supercritical state, a detailed understanding of water’s molecular structure, particularly how water molecules are arranged in this unique state, is crucial. This knowledge is especially relevant for develo** Monte Carlo track-structure chemistry simulations, which rely on a precise representation of the molecular nature of the medium—all current models of radiation tracks employ a continuum approach to water structure, thereby neglecting the actual molecular structure [27]. An accurate description of this structure is key to advancing radiolysis studies, as it enables more reliable predictions of the yields of radiolytic species in nuclear reactor coolants [26,28,29].

Unlike ordinary water, the hydrogen bond (H-bond) network in SCWRs becomes unstable due to the high operating temperatures and pressures. This instability results in the frequent breaking and re-forming of hydrogen-bonded structures, leading to substantial density fluctuations within the water. As a result, a heterogeneous medium forms, characterized by areas of high density comprising locally hydrogen-bonded molecules arranged in diverse tetrahedral configurations, suggestive of a ‘liquid-like’ state, and simultaneously by areas of low density where molecules exist in a ‘gas-like’ state, not bonded to each other. A number of studies have helped to shed light on key parameters (see, e.g., [30,31,32,33,34,35,36,37,38,39,40,41]), yet significant gaps in knowledge persist, especially regarding the spatial arrangement of water molecules in such extreme conditions. In this context, molecular dynamics (MD) simulations have proven to be particularly valuable. They skillfully generate realistic configurations of water molecules, accurately reflecting the conditions dictated by the high temperatures and pressures characteristic of SCW. In the field of radiation science, this detailed molecular picture is essential for a deeper comprehension of how radiation energy is deposited within SCWRs [42], an aspect that remains incompletely understood.

This study builds upon prior research, where we probed the heterogeneous properties of SCW at diverse densities at 400 °C [28,29]. We now aim to employ MD simulations to conduct a detailed analysis of the structural characteristics of SCW over a temperature range of 400 to 800 °C while maintaining a constant pressure of 31 MPa at 400 °C and 25 MPa between 500 and 800 °C. These particular conditions closely correspond with the operational parameters established for SCWRs and their SCW-SMR variants [7,12,14].

2. Results

Figure 1 illustrates the variation in water density for SCW at 400 °C across a pressure range of 24 to 35 MPa. Notably, at approximately 31 MPa, our molecular dynamics simulations align remarkably well with experimental data [43], exhibiting a minimal relative error of 1.3%. This high degree of accuracy at 31 MPa renders it an exceptionally suitable pressure for our computational analysis. In contrast, deviations from experimental densities become evident at other pressures, particularly at 25 MPa, where the error significantly increases, thereby reducing its applicability for our studies (at 400 °C). As a result, we have determined that a pressure setting of 31 MPa offers a more dependable and appropriate choice for conducting precise calculations.

To explore the molecular heterogeneity of SCW at 400 °C and the optimal pressure of 31 MPa, we analyzed the radial distribution function, g(r), with respect to the center-of-mass (COM) of water molecules. This function is defined as follows [28,29]:

$$g\left(r\right)\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\frac{〈n\left(r\hspace{0.17em} \pm \hspace{0.33em} \frac{\delta r}{2}\right)〉}{4 \pi r^2 \rho \delta r\hspace{0.17em}},$$

(1)

where $n\left(r\hspace{0.33em}\pm \hspace{0.33em}\frac{\delta r}{2}\right)$ represents the number of molecules within a spherical shell spanning from radius r to $\left(r\hspace{0.33em} \pm \hspace{0.33em} \frac{\delta r}{2}\right)$, essentially within a volume of $4\hspace{0.17em} \pi \hspace{0.33em}{r}^2 \delta r$. The notation $〈\dots 〉$ signifies that this average is computed over a duration of 10 ns in the MD simulation. This function plays a pivotal role in quantifying the likelihood of locating water molecules at a precise distance r from a given reference water molecule.

Figure 2 reveals a prominent peak in g(r) at a distance of 2.9 Å under the defined thermodynamic conditions of 400 °C and 31 MPa. This pronounced peak closely resembles those observed in water at temperatures below the critical point, specifically at 200, 300, and 350 °C, as documented in [26]. The origin of this peak can be traced to the first hydration shell encircling a water molecule, attributable primarily to the presence of strong hydrogen bonding within the water molecular structure. As the distance r increases beyond this point, g(r) shows a noticeable dip, implying the formation of void spaces. The secondary peak at ~5.5 Å indicates the presence of a neighboring molecular cluster, likely comprising 4.5 water molecules [29]. As the distance further increases, the interactions among water molecules reach a state of equilibrium, indicative of a balancing of molecular forces. This pattern, marked by alternative peaks and wells, highlights the molecular heterogeneity of SCW at 400 °C and a pressure of 31 MPa [44].

To enhance our understanding of molecular structures, we conducted a study at significantly higher temperatures in the supercritical state, specifically at 500–800 °C, and at a pressure of 25 MPa (as mentioned in the Introduction, this pressure corresponds with the operational environment chosen for SCWRs and their SCW-SMR variants). This investigation allowed us to identify and analyze the distinct molecular configurations present at these temperatures. In Figure 3, we compare the experimental and simulated densities at 25 MPa pressure. Small discrepancies are noted, which decrease as the temperature increases. Further studies to investigate local phenomena can thus be undertaken.

Examination of these structures using the COM radial distribution function uncovers a uniformity in the interactions among water molecules beyond the 2.9 Å peak position, as shown in Figure 4. Moreover, in contrast to the observations made in Figure 2, under conditions of lower pressure and temperatures exceeding 400 °C, a certain instability becomes apparent in these interactions within the 4 Å to 8 Å range, denoted by the absence of a well. This lack of stability can be ascribed to the weakening of the structural integrity of the water network, which leads to the breakdown of hydrogen bonds [45,46]. As these weaker intermolecular forces break down, the molecular structure becomes less stable, leading to increased instability.

To unveil this behavior, an investigation into the variation of oxygen–oxygen coordination numbers in response to temperature changes provides a way to uncover the development of water molecule clusters at 500–800 °C. A key observation is the progressive decline in the O–O coordination number as temperature increases from 500 °C to 800 °C, spanning a significant range from ~1.4 at the lowest temperature to ~0.5 at the highest (Figure 5). Moreover, our findings reveal a corresponding increase in the number of clusters with increasing temperature, if we consider the cluster based on a neighbor cutoff distance of 3.5 Å, which is the maximum distance between molecules to be considered as part of the same cluster (Figure 6). This trend indicates a direct correlation between rising temperature and the reduction of coordination numbers. At lower temperatures, water molecules have less kinetic energy and cluster more closely together, forming more stable structures through hydrogen bonding. However, as temperature increases, the increased thermal motion leads to the disintegration of molecular clusters by breaking their hydrogen bonds. This disruption causes the clusters to fragment into a greater number of smaller clusters or isolated molecules, as illustrated in the snapshots shown in Figure 6. In other words, this observed decrease in cluster size, along with the reduced number of molecules linked to the coordination number at higher temperatures, could be explained by a decrease in hydrogen bonds, which is directly correlated with the O–O coordination number in the water system.

The increase in the number of clusters or the reduction in cluster size not only corroborates the formation of more atomic clusters at higher temperatures but also reflects a significant alteration in the molecular structure of water. Additionally, investigations into hydrogen bonds have revealed a decrease in their prevalence with increasing temperature, indicating significant alterations in the molecular dynamics of water under varying thermal conditions.

3. Materials and Methods

In our simulations, water molecules are modeled as rigid entities, with their intramolecular degrees of freedom held constant. Interaction among these molecules is treated through the extended simple point charge (SPC/E) pair potential [47]. This rigid water model is renowned for its relative simplicity and computational efficiency, making it a prevalent choice for MD simulations of water-based systems. It differs from the original simple point charge (SPC) model [48] by integrating a corrective term that reconciles the variances in polarization self-energy observed between a molecule in its liquid state and one in a vacuum.

In the SPC/E representation, a water molecule is depicted through three interaction sites: one centered on the oxygen atom and the others on the two hydrogen atomic nuclei. The O–H bond length is consistently set at 1 Å, and the H–O–H angle is maintained at the precise tetrahedral angle of 109.47°. Charge distribution is modeled with the oxygen carrying a charge of q_O = −0.8476|e| and each hydrogen having q_H = −q_O/2 = +0.4238|e|, where |e| represents the magnitude of the electron’s charge. This arrangement ensures charge neutrality and accurately accounts for the Coulomb interactions within the system:

$$U_{\mathrm{Coulomb}}\hspace{0.33em}=\hspace{0.17em}\hspace{0.33em}\frac{1}{4 \pi \hspace{0.17em}{\epsilon }_0}\hspace{0.33em}{\displaystyle \sum _{i,j}\hspace{0.17em}\frac{q_i\hspace{0.17em}{q}_j}{r_{ij}}},$$

(2)

where the summation is over all distinct pairs i and j, ε₀ is the vacuum permittivity, q_i and q_j are the partial charges on atoms i and j, and r_ij is the non-bonding distance between atoms i and j on two different molecules.

The van der Waals forces are described using a “12-6” Lennard-Jones (LJ) potential, which only concerns oxygen atoms. This potential quantifies the interaction energy between pairs of oxygen atoms, each belonging to distinct water molecules, and depends on their separation distance, denoted as r_OO:

$$U_{\mathrm{LJ}}\hspace{0.33em}= \hspace{0.33em}4\hspace{0.17em}\hspace{0.17em}{\epsilon }_{\mathrm{LJ}}\hspace{0.17em}\hspace{0.17em}\left[{\left(\frac{\sigma }{r_{OO}}\right)}^{12}\hspace{0.17em}-\hspace{0.17em}\hspace{0.33em}{\left(\frac{\sigma }{r_{OO}}\right)}^6\right],$$

(3)

where ε_LJ = 0.6502 kJ/mol is the LJ potential well depth and σ = 3.166 Å is twice the van der Waals radius of oxygen.

As mentioned above, the SPC/E model further incorporates an average self-polarization energy correction term, detailed as follows [47]:

$$U_{\mathrm{pol}}\hspace{0.33em}=\hspace{0.33em}\frac{{\left(\mu \hspace{0.33em}-\hspace{0.33em}{\mu }_0\right)}^2}{2\hspace{0.17em}\alpha },$$

(4)

where μ is the effective dipole moment of polarized water in the liquid state, which is 2.351 Debye (D), μ₀ = 1.85 D denotes the dipole moment of the water molecule in its isolated, gas phase form, and α = 1.45 ų is the isotropic polarizability of water.

The total intermolecular interaction potential, U, for the SPC/E water model is therefore represented as the sum:

$$U\hspace{0.33em} =\hspace{0.33em}\hspace{0.17em}{U}_{\mathrm{Coulomb}}\hspace{0.33em}+\hspace{0.33em}\hspace{0.17em}{U}_{\mathrm{LJ}}\hspace{0.33em}+\hspace{0.33em} U_{\mathrm{pol}}$$

(5)

The SPC/E force field, despite its simplicity, effectively replicates many of the thermodynamic and physical properties of liquid water, including diffusion coefficients as well as the dielectric constant over a wide range of temperatures and pressures [49,50]. It also reproduces, with good agreement, the liquid–vapor coexistence curve for water and predicts critical parameters, which compare favorably to the experimental values [51]. Furthermore, it has also demonstrated proficiency in accurately characterizing the structures and properties of water at high temperature and pressure conditions up to the supercritical state [26,28,29,39,40,41,42].

The simulations were conducted in accordance with the established procedure previously employed by Metatla et al. [28,29]. In summary, each simulation involved a cubic system measuring 81 Å on each side, comprising N = 7000 rigid water molecules, which corresponds to a total of 21,000 atoms. This system was extended in all three directions through periodic boundary conditions. The selected system size is deemed suitable for adequately representing the thermodynamic limit, thereby ensuring a proper representation of the configuration space. For handling long-range Coulomb interactions, the Ewald summation method [52] was employed and a cutoff radius of 9.5 Å was applied to the LJ interactions.

To preserve the rigid tetrahedral geometry of water molecules, the SHAKE iterative algorithm [53] was employed, applying a relative geometric tolerance of 10⁻⁴. Newton’s equations of motion were integrated using the velocity Verlet algorithm [52,54], and the integration time step was fixed at 1 femtosecond (fs). The simulations extended over 10⁶ steps, equivalent to a duration of 1 nanosecond (ns). All MD calculations were conducted using the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) simulation package [55].

The temperature ranged from 400 to 800 °C, with the pressure held constant at 31 MPa at 400 °C and maintained at 25 MPa for temperatures between 500 and 800 °C. The initial configurations were devised using the BIOVIA Materials Studio Amorphous Cell building tool [56], set to simulate normal ambient water conditions of 25 °C and a density of 1 g/cm³. These configurations were then subjected to a standard energy minimization process for equilibration. This involved employing a mix of steepest descent and conjugate gradient algorithms, carried out across 5000 steps. Subsequent to this minimization, the systems were progressively heated to their respective target temperatures, increasing in increments of 25 °C, within the canonical NVT ensemble, where the number of molecules, the volume of the simulation cell, and the temperature remain constant. Following the heating phase, the systems were equilibrated for 100 ps at an external pressure of 25 MPa. The Nosé-Hoover thermostat algorithm [57,58] was used to maintain the systems at the desired temperature, and the Andersen barostat algorithm [59] was employed for pressure regulation. Finally, the MD simulations were performed in the isobaric–isothermal NPT ensemble, maintaining a constant number of molecules at constant pressure and temperature, until the system achieved an equilibrated density, a process completed over the course of 1 ns.

The molecular configurations corresponding to the range of temperatures studied, and at a constant pressure of 25 MPa, have been determined. The clusters were identified based on a selection criterion that considers the separation distance between adjacent atoms [60]. We chose a threshold of 2.9 Å for the separation distance to identify these clusters. This particular value was selected because it matches the highest molecule count in our center-of-mass (COM) radial distribution function (RDF) diagrams. This pattern is consistently observed across all the temperatures and pressures we studied.

4. Conclusions

The study investigated the heterogeneity of water in its supercritical state through SPC/E molecular dynamics simulations. The results indicated that at 400 °C, there is a good correlation between density values derived from theoretical and experimental methods at a pressure of 31 MPa. This finding underscores the distinct heterogeneity of water, as demonstrated by its molecular structure and the center-of-mass radial distribution function. The latter exhibits distinct molecular peaks and wells, indicating the presence of voids. A noteworthy observation is the consistency of the peak positions at 2.9 Å across all studied temperatures and pressures, mirroring patterns observed in the subcritical state. Moreover, it was observed that an increase in temperature correlates with a rise in the number of clusters formed, along with a reduction in hydrogen bonding.

This study indicates that traditional SCW radiolysis models, which consider water as a continuum, may require updates to incorporate these findings. Implementing such revisions would enhance the accuracy of predictions regarding SCW behavior under radiation. This could be important in applications like nuclear power generation and waste processing, where water is frequently utilized as a coolant and can reach supercritical states. This is especially the case for advancing our knowledge and capabilities in the context of SCWRs and their SCW-SMR variants. By better understanding the molecular interactions and transformations under these conditions, researchers and engineers can more accurately predict and mitigate the effects of radiolysis, leading to safer and more efficient reactor designs.