NaCl aggregation in water at elevated temperatures and pressures: Comparison of classical force fields

At high temperature and pressure, water becomes supercritical, adopting features of both liquid and vapor.¹ The continuous dependence of physical properties on temperature and pressure in this region of the phase diagram enables the tuning of water–solute interactions. Industrial application of supercritical fluids takes advantage of this tunability to optimize an environmentally benign and reusable solvent for separation and oxidation.² For example, the recovery of expensive metals from waste streams of electronic and electrical equipment using supercritical water has shown promise.³ The properties of supercritical water are also essential for the transport of ions and the formation of ores in geothermal fluids.⁴ 

Our interest in subcritical to supercritical fluids stems from emerging applications in water desalination.^5,6 One option for supercritical desalination involves pressurizing and heating salty feed water to conditions under which it forms a supercritical water phase and a liquid phase consisting of concentrated brine. In the supercritical water phase, the solubility of sodium chloride (NaCl) can be low enough to meet drinking water standards.⁵ The phases are then physically separated to yield product steam/water and a salty brine for disposal. Another option for supercritical desalination involves compressing salt-laden feed water above the solution critical pressure with subsequent heating above the critical temperature. De-salting this single phase fluid near, or above, the critical temperature exploits changes in the solubility power of water as the temperature and density vary. Under ambient conditions, water acts as a strong dielectric to screen ionic charges. However, the dielectric constant decreases significantly near and above the critical temperature, which can be exploited to control solubility. For example, supercritical conditions have previously been used to precipitate NaCl by lowering its aqueous solubility.⁷ An additional benefit of this second option for supercritical desalination is that it yields solid salt(s) and purified water. As such, it could eliminate the need to dispose of concentrated brine, which is a processing expense encountered with most traditional desalination techniques. By manipulating the temperature and density of the supercritical fluid, it may also be possible to selectively precipitate salts from multi-component brines. This could enable the recovery of valuable co-products from natural brines (e.g., rare earth elements or lithium) leaving low-to-moderate density water containing very low levels of residual salt. However, to demonstrate the viability of generating valuable co-products by selective precipitation, understanding and controlling solute aggregation and nucleation is essential. Thus, our present focus is on developing that understanding in single phase subcritical and supercritical salt solutions, where it proves beneficial to first understand the aggregation behavior of the primary solutes before considering how the addition of other solutes alters aggregation, and ultimately, solubility.^8,9 With accurate force fields, molecular simulation could provide the required insights into atomic scale aggregation under subcritical to supercritical conditions.

While there are a number of source waters that can be used for desalination,^4,10 NaCl is a common contaminating solute. It is also a major constituent of water encountered in industrial applications. In addition to these practical considerations, the relative wealth of experimental data for aqueous NaCl and the development of force fields for both water and NaCl make this system an excellent candidate for simulation studies aimed at comparing salt/water force fields under subcritical to supercritical conditions. Recently, the importance of polarizability and hydrogen bonding in water has motivated development of force fields that include Gaussian charge distributions, hydrogen bonding interactions, fluctuating charges, and Drude oscillators.^11–17 Including polarizability into water force fields can result in better predictions along the near critical region of the vapor–liquid coexistence curve,^11,14,15,17 vapor phase properties,^11,13,17 dielectric constants,^13,17 and physically realistic behavior at elevated temperatures¹² while requiring relatively modest increases in computation cost.¹³ However, despite the obvious importance of polarizability and other quantum phenomena, classical Lennard-Jones (LJ) and Coulombic force fields are still often employed to study such systems since they are relatively computationally inexpensive. Moreover, the varying functional forms of polarizable models are not often included out-of-the-box in open source molecular dynamics (MD) simulation codes.¹⁸ There have been a number of simulation studies performed on NaCl in subcritical to supercritical water using classical force fields, including studies of clustering,^19–24 diffusion,^20,24–27 viscosity,^28,29 conductance,^24,30–33 dielectric permittivity,^34–39 potentials of mean force,^40–42 density,^29,43 and hydrogen bonding and solution structure.^{20,25,44–55} Although there are several comparisons of force field combinations under ambient conditions^{25,34,35,37,56–58} and one at temperatures up to 473 K,²⁹ there are no prior studies comparing the performance of classical force fields to each other near supercritical conditions.

Here, we compare classical force field combinations under near-critical and supercritical conditions using molecular dynamics (MD) simulations. Specifically, we compare five classical point charge NaCl force fields (Joung–Cheatham,⁵⁹ Dang 1995,⁶⁰ Kirkwood–Buff,⁶¹ Madrid 2019,⁶² and a scaled charge version of Joung–Cheatham) in combination with two rigid non-polarizable water force fields (TIP4P/2005⁶³ and SPC/E⁶⁴). Simulations at temperatures spanning subcritical to supercritical states along the 200 bar isobar are used to evaluate the force field combinations with NaCl weight fractions ranging from 0% to 20%. We present force field comparisons for the solution density, equilibrium ion clustering, hydrogen bonding structure, and static dielectric constant.

Several performance-based considerations were involved in selecting force fields for comparison. For rigid point-charge water force fields, TIP4P/2005⁶³ and SPC/E⁶⁴ were selected. Since TIP4P/2005 was parameterized to reproduce the phase behavior of water, one might expect it to be a fairly accurate water potential as supercritical conditions are approached. The SPC/E potential also does a relatively good job at reproducing the vapor coexistence curve. In a comparative study by Vega and Abascal, TIP4P/2005 and SPC/E both scored well relative to other rigid non-polarizable force fields.⁶⁵ In addition, a comprehensive review of available water force fields from ambient to supercritical conditions⁶⁶ indicated that TIP4P/2005 and SPC/E were among the best rigid non-polarizable force fields for water. It should be noted that while both water force fields accurately predict liquid phase densities, they are not as accurate in reproducing vapor phase densities (and pressures). Furthermore, both force fields underestimate the critical temperatures, pressures, and densities, although TIP4P/2005 is somewhat more accurate than SPC/E (see Table SI). For instance, experimentally T_c is 647 K, while TIP4P/2005 and SPC/E yield 640 K and 639 K, respectively.^66,67 Note that the simulated critical points can vary depending on the method used, e.g., others have reported a value of T_c as low as 630 K for SPC/E.⁶⁸ 

The Joung–Cheatham (JC) force field for NaCl (which was originally parameterized for use with the SPC/E water potential) was chosen due to its popularity and parameterization.⁵⁹ The JC force field was designed to reproduce the hydration free energies, lattice energies, and lattice constants across all of the alkali halide pair combinations. The JC force field was used in combination with both the SPC/E and TIP4P/2005 water potentials. The combination of TIP4P/2005 with JC has been used previously and found to have performance comparable to that of JC with SPC/E.^56,69 We also examined the 1995 model of Dang (D95), which was designed as a nonpolarizable force field for use with SPC/E.⁶⁰ This force field is similar to the more common Smith–Dang model⁷⁰ but does not include any polarization terms. Although this is an older model, it is widely used in the literature and, thus, warrants inclusion in any comparison of classical NaCl force fields. In addition to SPC/E, we also examined this force field in TIP4P/2005 water. We also chose to examine the Kirkwood-Buff force field (KBFF) of Weerasinghe and Smith, which was designed to reproduce the experimental Kirkwood–Buff integrals and activity coefficients.⁶¹ Although KBFF predicts the experimental solubility as a function of NaCl concentration very well at ambient conditions, its accuracy at higher temperatures is uncertain.⁷¹ The Madrid 2019 (M19) NaCl force field is a scaled charge (±0.85 vs ±1.00) force field that was initially parameterized for use with the TIP4P/2005 water force field.⁶² M19 is a newer force field parameterized to reproduce solution densities while avoiding the propensity for excessive clustering exhibited by many older NaCl force fields. Instead of a mixing rule like Lorentz–Berthelot (LB) to define unlike atom/ion force field interaction terms, the M19 force field has fitted potential cross terms (ϵ_ij and σ_ij). The final force field used in this study was a scaled charge version of the Joung–Cheatham (SJC) NaCl force field with TIP4P/2005 water. The SJC force field serves as a control for comparison to the M19 force field. The intent is to separate the influence of scaling the NaCl charges to ±0.85 from the influence of the LJ cross terms. This combination of force fields represents a broad sample of commonly used NaCl and water models that were parameterized to reproduce different properties.

Classical MD simulations were performed on sodium chloride (NaCl) in water at concentrations of 2 wt. %, 10 wt. %, 15 wt. %, and 20 wt. % NaCl. Simulations contained 1000 water molecules and either 6, 34, 54, or 77 molecules of NaCl, present as Na⁺ and Cl⁻ ions. These combinations were chosen to span low-to-moderate NaCl concentrations where one would not anticipate the formation of large ionic clusters or nucleation of solid NaCl in liquid phases. In addition, simulations were performed at select state points with 5000 water molecules and equivalent NaCl concentrations to verify that a system size of 1000 water molecules is sufficient (see supplementary material, Fig. S1).

We evaluated eight classical force field combinations that rely on Lennard-Jones potentials

and Coulomb potentials

These potentials are a function of the distance r_ij separating two atoms/ions i and j, the depth of the interaction well ϵ_ij, the distance at which the potential is zero σ_ij, the electronic charges q_i and q_j, and the permittivity of the vacuum ϵ₀. Specific force field parameters are provided in Table I. All of the force fields, except the Madrid 2019 NaCl force field, use the Lorentz–Berthelot (LB) mixing rules for unlike atoms i and j,^72,73

whereas the ij interaction for the Madrid 2019 salt force field involves fitted parameters.

Simulations were performed along the 200 bar isobar for all of the force field combinations. A pressure of 200 bar is well within the range that can safely be accessed experimentally. For each concentration, simulations were conducted at 450 K–750 K in 50 K increments. In addition, simulations were performed along isobars spanning the range of 100 bar–600 bar in 100 bar increments to investigate the influence of pressure on NaCl aggregation for the JC NaCl force field with SPC/E and TIP4P/2005. Simulations of pure SPC/E and TIP4P/2005 water were also performed along isobars at 100 bar–400 bar as well to provide reference points for how the salts altered water properties.

The GROMACS 2019.2^18,74–78 simulation package was used to perform all MD simulations. Atomic positions and velocities were propagated using leap-frog integration of Newton’s equations of motion with a time step of 2 fs.⁷⁹ Configurations from the simulation trajectory were saved every 0.5 ps. The Verlet cutoff scheme was used to construct the neighbor lists.⁸⁰ Electrostatics were handled using Particle Mesh Ewald (PME) summations.^81,82 Both LJ and PME interactions used static cutoff distances of 14 Å for simulations with initial simulation box dimensions less than or equal to 40 Å and a distance that is 40% of the initial simulation box for simulations with dimensions greater than 40 Å.^81–83 These cutoffs are necessary because some of the state points end up being low density phases. Tables of the exact cutoff distances used for each replica and state point are provided in the supplementary material. Temperature coupling was achieved with velocity re-scaling and a time constant τ_T of 1.0 ps.⁸⁴ Isotropic pressure coupling was performed for the isobaric–isothermal ensemble (NpT) simulations using the Berendsen⁸⁵ barostat for an initial 1 ns equilibration. A second 1 ns equilibration with the Parrinello–Rahman⁸⁶ barostat was followed by 19 ns production simulations with a time constant τ_p of 4.0 ps and a compressibility of 4.9 × 10⁻⁵ bar⁻¹. Each state point had four independent simulations performed: the first starting with randomly placed free ions in solution and a high starting density, the second with a low starting density, the third with randomly placed ion pairs in solution, and the last with all of the ions in a single large cluster. The Packmol program was used to generate the initial configurations.⁸⁷ 

Reported properties for NaCl clustering, water hydrogen bonding, and density were calculated over 19 ns from the NpT production simulations. A custom analysis code for clustering and hydrogen bonding was developed in-house (see details in Sec. II B). We chose to perform the cluster and hydrogen bond analysis on the production run from the NpT ensemble simulations to ensure that no artifacts from constraining the density were present in our results. However, it turns out that the results are invariant to the choice of NpT vs NVT ensemble (see Figs. S2 and S3 in the supplementary material). The final configurations from each of the 200 bar NpT simulation replicas were then scaled to the average volume of the respective simulations and used as starting configurations for 50 ns canonical ensemble (NVT) simulations. Dielectric properties were calculated from the 50 ns simulations in the NVT ensemble.

Following Stillinger,⁸⁸ we define the NaCl cluster as a group of ions that are linked by a sequence of nearest neighbor associations. A nearest neighbor association is taken to be two ions within a given cutoff distance r_c of one another. For each temperature, pressure, NaCl wt. %, and force field combination, r_c was set as the distance to the first minimum in the radial distribution function (RDF) rounded to two significant figures.

Case-specific selection of r_c is necessary for each force field combination and state point due to variations in the RDFs and broadening of the first peak with increasing temperature (see Fig. 1). Changing the water force field from SPC/E to TIP4P/2005 does not significantly alter r_c for the JC, D95, and KBFF combinations (Fig. 1). However, changing the salt force field does have an impact. The complete table of specific r_c values used for the NaCl cluster cutoff distance is reported in the supplementary material, Tables SIV and SV. The locations of the first RDF maxima exhibit less variance than the locations of the first minima (supplementary material, Tables SVI and SVII). Different force field combinations at the same state points in Fig. 1 have RDF maxima distances that differ at most by 0.12 Å and minima distances that differ at most by 0.3 Å.

The cluster size s is defined as the sum of Na⁺ ions and Cl⁻ ions in a given cluster. Cluster size histograms are calculated for each NpT ensemble simulation, and we report the average values from four independent simulations. The fraction of ions f_s in a given cluster of size s is calculated as

where n_s,i is the number of clusters of size s identified in the frame i, N_frame is the number of trajectory frames analyzed, and N_ions is the total number of ions in the simulation.

A hydrogen bond between two water molecules is defined as satisfying the following three geometric conditions: The distance between the two oxygen atoms d_O⋯O must be less than or equal to 3.3 Å, and the distance between an oxygen atom and a hydrogen atom d_O⋯H must be less than or equal to 2.5 Å. The angle between the first molecule’s covalently bonded OH pair and the proposed hydrogen bond to the second molecule’s oxygen, θ_OH⋯O, must also be greater than or equal to 95.7°.⁸⁹ A schematic of the hydrogen bond definition is shown in Fig. 2.

Static dielectric constants ε are computed for the water molecules using the GROMACS gmx dipoles utility. ε is calculated for each of the 50 ns NVT ensemble simulations using the dipole fluctuation definition obtained from linear response theory,^83,90

where ϵ₀ is the permittivity of the vacuum and k_B is the Boltzmann factor. $M2−M2$ is the mean squared fluctuations in the total system dipole moment vector M contained within a volume V at a temperature T. The total system dipole is taken to be the sum of the product between the atomic charges q_i and their respective positions r_i over all N water atoms. Long simulations are required since $M$ and, thus, also ε are slow to converge. A typical example of the convergence of the static dielectric constant is shown in supplementary material, Fig. S6.

This method of calculating ε only takes into account dielectric contributions from the water molecules. This is a common approximation when calculating the static dielectric response and has been applied previously to electrolyte solutions.^34,91 When ions are present in solution, the static dielectric constant from Eq. (5) is not equivalent to the complex frequency-dependent permittivity measured experimentally.⁹¹ This formula does not account for dielectric and conductivity contributions from ion pairs and clusters, which may contain dipoles and/or other electrostatic moments, which complicate the calculation of conductivity.^92,93 The static dielectric constant obtained from Eq. (5) does, however, still provide some insight into how the introduction of ions alters the orientation of water molecules and their response under an applied field.

It should be noted that some of the state points at 200 bar exhibited volume changes that indicate instability. Specifically, simulations of SPC/E water with 20 wt. % JC and KBFF at 700 K exhibited larger volume fluctuations that are likely indicative of proximity to the vapor–liquid phase envelope (see the supplementary material, Fig. S7). Simulations of TIP4P/2005 water with 20 wt. % JC and KBFF at 700 K show sustained volume jumps that suggest that the state point is very close to the liquid–vapor coexistence curve and vacillating between high and low density states. For clarity, these data points are omitted in the presentation of our results.

In order to probe the quantitative accuracy of the various force field combinations, experimental data are required. Due to the limited availability of quantitative experimental data on structural properties for aqueous NaCl at elevated temperatures, we have compared the densities to the equation of state (EOS) developed by Driesner et al.^94,95 Driesner’s EOS provides an accurate P–V–T relationship, as a function of NaCl concentration (x). In our comparisons, it is used as a surrogate for the experimental P–V–T–x data. We note that a direct and comprehensive comparison to experimental pressure-composition and/or temperature-composition phase diagrams would be useful; however, the required simulations are very time-consuming and beyond the scope of the present study.

Figure 3(a) shows the predicted density for each of the force field combinations along the isobar at 200 bar. A plot of the density differences between the Driesner EOS and the simulated densities is shown in Fig. 3(b) for ease of comparison. All of the force field combinations succeed in qualitatively reproducing the trend in the liquid phase densities. However, the simulated liquid densities are consistently lower than the EOS, with the difference between the simulated and EOS densities increasing with increasing temperature up to the phase transition. This may be partially due to the behavior of the water force fields, which exhibit densities that are lower than the IAPWS 95 EOS⁶⁷ (see the supplementary material, Fig. S8). Both water force fields are known to underpredict liquid densities and overpredict vapor densities (and, thus, underpredict vapor pressures).^{65,66,68,96,97}

At very low NaCl concentrations, the salt force fields do not exhibit much influence. At 2 wt. % NaCl, there is a clear split between the water force fields, which is best shown in Fig. 3(b). The force field combinations with SPC/E water exhibit lower liquid densities than TIP4P/2005 for the same temperatures and pressures, which is what one would expect from the behavior of the pure water force fields.^65,66 Upon reaching 10 wt. % NaCl, the influence of the salt force fields on the densities is more pronounced.

The T+JC combination shows the smallest deviation from the Driesner EOS. Interestingly, the combination of TIP4P/2005 with KBFF performs almost as well as T+JC and better than SPC/E with KBFF, even though the KBFF NaCl force field was parameterized with SPC/E.⁶¹ The T+SJC and S+D95 force field combinations are less accurate in matching the EOS density at higher concentrations. The T+SJC and T+M19 scaled charge force fields exhibit a somewhat different dependence on temperature. The degree of aggregation is similar for these two scaled charge force fields (see Fig. 7), which would imply that the differences in density between these scaled charge models are unlikely to be attributed to the number and/or type of species in solution. This difference is most likely due to differences between the LJ parameters and optimized cross terms of T+M19 and the original JC force field LJ parameters (developed with charges of ±1), which were kept for T+SJC.

A phase change from liquid to vapor (or possibly a low density supercritical phase) is characterized by an abrupt change in density, which can be seen at temperatures around 650 K along the 200 bar isobar. In Fig. 3(a), such phase changes are seen for simulations with 10 wt. % NaCl between 650 K and 700 K for all of the force field combinations. At higher NaCl concentrations, the phase behavior of the force field combinations differs. The scaled charge force field combinations (T+M19 and T+SJC) exhibit a phase change near 700 K at 15 wt. % and between 700 K and 750 K at 20 wt. %, whereas the combinations with full charge salt models exhibit a phase change at somewhat lower temperatures. Moreover, as concentration increases, we find that all models show a small increase in the approximate temperature where vaporization occurs, indicating that the liquid phase becomes more stable at high temperatures with increasing NaCl concentration at 200 bar. Based on the Driesner equations of state, the critical temperatures increase from pure water at 647 K to 673 K, 738 K, 790 K, and 851 K for 2 wt. %, 10 wt. %, 15 wt. %, and 20 wt. % NaCl, respectively. Simultaneously, the critical pressures increase from pure water at 221 bar to 277 bar, 463 bar, 640 bar, and 845 bar.⁹⁸ While the qualitative behavior we see is consistent with a higher critical temperature, it is unlikely that any of these models quantitatively reproduce the P–V–T–x phase diagram, especially at higher concentrations of NaCl.

Isobars for the simulation density at 10 wt. % JC NaCl, dissolved in SPC/E and TIP4P/2005 water, are displayed in Fig. 3(c). Changes in pressure over the range of 100 bar–600 bar have a minimal impact on the liquid phase densities due to the incompressibility of the solution. As the temperature increases, the density decreases for both the TIP4P/2005 and SPC/E force fields with JC NaCl, in agreement with the Driesner EOS. The system transitions to a low density phase between 600 K and 650 K at 100 bar and between 650 K and 700 K at 200 bar. The higher pressures do not exhibit the sharp transition in density seen at 100 bar and 200 bar, indicating that the liquid phase is stable to higher temperatures at higher pressures. Both force field combinations exhibit similar qualitative trends and produce lower densities than the Driesner EOS. The differences in the densities are displayed in Fig. 3(d). In both cases, the difference decreases with increasing pressure, although this trend is more pronounced for S+JC than T+JC.

The clustering behavior of ions in solution is known to be sensitive to the force field parameters used in the simulation.²⁶ Figure 4 quantifies NaCl aggregation along (a) the 200 bar isobar and (b) the 700 K isotherm for 10 wt. % NaCl modeled by the S+JC force field. As the temperature increases along the isobar, the simulated densities decrease (Sec. III A). Significant density decreases are accompanied by a distinct increase in ion clustering. At 200 bar, the fraction of free ions (f₁, s = 1) decreases by 30% over the range of 600 K–650 K [Fig. 4(a)]. However, f₁ decreases much more dramatically when the temperature increases from 600 K to 700 K (90%) or 750 K (95%), where a change to a low density phase occurs [vapor or low density supercritical phase, see Fig. 3(a)]. The decrease in free ions coincides with the formation of larger cluster sizes, especially for the low density states at 700 K and 750 K. The histograms demonstrate that these larger clusters are not limited to ion pairs (where s = 2). A single large amorphous cluster that may be approaching the critical nucleus size exists in equilibrium with a few small clusters (free ions and/or ion pairs), which appears to be typical for high temperature, low density phases. Cluster size distributions of this nature are the result of the high salt concentration and finite size effects imposed by the simulation.⁹⁹ This suggests that 10 wt. % NaCl may be above the saturation concentration at these conditions. This is to be expected since salt solubility in low density steam is very low.^9,94,100,101 Sherman and Collings also observed the formation of larger NaCl clusters in higher temperature simulation of Smith–Dang NaCl salt with SPC/E.¹⁹ As the pressure decreases along the 700 K isotherm [Fig. 4(b)], the fraction of free ions decreases by 6%, 19%, 40%, and 92% relative to 600 bar. At low densities (200 bar, 700 K), there is also a preference for charge neutral clusters with sizes such as s = 2, 4, and 6. The preference for ion pairs over free ions is also present at pressures of 300 bar–400 bar. This could be the result of reduced electrostatic screening in the low density phases.

Figure 5 shows cluster size histograms for all of the force field combinations with 10 wt. % NaCl at 200 bar and (a) 650 K (liquid phase) and (b) 700 K (low density phase). The degree of ion clustering increases with decreasing density for all of the force field combinations. The extent of clustering is strongly influenced by the salt force field employed, while the water force field plays a lesser role. The scaled charge salt force fields (M19 and SJC) show less clustering than the full charge force fields (JC, D95, and KBFF).

The fraction of free ions (f₁, s = 1) in solution is a general indicator of the relative degree of clustering between the different force field combinations. Force field comparisons at the same state point are consistently made relative to T+M19, the force field combination that has shown the least clustering and, thus, the highest f₁. Relative to the T+M19 force field combination, f₁ at 650 K decreases by 8% (T+SJC), 35% (S+JC), 38% (S+KBFF), 47% (T+JC), 53% (T+KBFF), 65% (S+D95), and 72% (T+D95) for the respective force field combinations shown in Fig. 5(a). With an increase in temperature to 700 K, the densities decrease [see Figs. 3(a)] and f₁ for the full charge force fields decreases by 69%–96% relative to their values at 650 K. The scaled charge force fields show less sensitivity to the decreased density, and f₁ only decreases by 69% and 70% for M19 and SJC, respectively.

All of the force fields in the low density phase exhibit a preference for neutral clusters. Ion pairs (s = 2) are 1.5–2.2 times more prevalent than free ions at 700 K. The D95 force field combinations are much more prone to forming large clusters, which is best exhibited by the cluster size distribution shifting to primarily the largest clusters possible in the system at 700 K (Fig. 5(b), inset). On the other hand, the scaled charge force field combinations, T+M19 and T+SJC, yield almost exclusively small clusters sizes (s ≤ 10). Snapshots of the simulations for T+M19 and T+JC are given in Fig. 6 to show a representative sample of the ensemble averaged distributions at 700 K and 200 bar.

In Fig. 7, we show the clustering trends as a function of the salt concentration in the liquid phases at 450 K and 650 K. There is a separation in the clustering behavior of the force fields. The salt force field primarily determines the extent of aggregation, but the extent of aggregation varies slightly depending on the water force field. Salt force field combinations with SPC/E water exhibit less NaCl clustering. Increasing the temperature from 450 K to 650 K increases the propensity for ions to cluster.

All of the force field combinations show a decrease in f₁ as the wt. % of NaCl increases. As the salt content is increased from 2 wt. % to 20 wt. % at 450 K, the force field combinations JC and KBFF with both water models, as well as T+M19 and T+SJC, show similar trends with f₁ decreasing by 31%–45% in response to the increase in concentration. On the other hand, the D95 force field combinations exhibit a greater degree of ion clustering, with a lower f₁ than any other force field combination, and a concentration induced decrease of 67% and 72% for S+D95 and T+D95, respectively. Increasing the concentration from 2 wt. % to 20 wt. % at 650 K results in decreasing f₁ by 52%–74% across all force field combinations.

Ion pairs (s = 2) are the first ion clusters to form in solution. At 450 K and 2 wt. % salt, there are almost no mid-sized or larger clusters formed for any of the force field combinations (see the two bottom panels on the left of Fig. 7). As the salt concentration increases to 10 wt. %, the fraction of ions in the form of ion pairs (f₂) increases from 0.04–0.10 to 0.15–0.24 for T+JC, T+KBFF, T+M19, T+SJC, S+JC, and S+KBFF. The increase is significantly less pronounced in the D95 salt combinations (from 0.19 to 0.27 and 0.14 to 0.26, for T+D95 and S+D95, respectively) due to their propensity to form larger clusters (s ≥ 3).

At 650 K, differences between the force field combinations are evident as the concentration increases from 2 wt. % to 10 wt. %. The fraction of ions bound as ion pairs increases by 24% and 21% for the scaled charge salt force fields M19 and SJC but decreases by 29%–55% for the full charge NaCl force field combinations. This is likely a reflection of NaCl solubility differences from the salt force fields. The scaled charge force fields primarily form ion pairs at concentrations where the other force field combinations predict the formation of mid-sized and/or large clusters. The D95 salt force field predicts even more mid-sized and large clusters at the same concentrations. As the concentration is increased from 10 wt. % to 20 wt. % at 650 K, all of the force field combinations exhibit decreases in the fraction of ions found in ion pairs.

Under conditions where the salt concentration is below the saturation limit and/or the system size is small enough to suppress the formation of a critical nucleus, mid-sized clusters (3 ≤ s ≤ 10) should be present and could even dominate the cluster size distribution. At 450 K, all of the force field combinations predict an almost negligible fraction of ions in the mid-sized cluster range at 2 wt. %. This changes dramatically as the concentration is raised to 10 wt. %, with the fraction of ions in mid-sized clusters increasing by at least an order of magnitude to 0.03–0.28. As the salt concentration is increased from 10 wt. % to 20 wt. %, the force fields then exhibit comparatively moderate increases of 1.6–4.4 times their values at 10 wt. %.

At 650 K, three different trends emerge. The D95 force field combinations exhibit an initial increase in the fraction of ions in mid-sized clusters as the concentration increases from 2 wt. % to 10 wt. %. This is followed by a decrease in mid-sized clusters as the concentration further increases from 10 wt. % to 20 wt. %. This is due to the formation of a long tail in the cluster size distribution toward larger cluster sizes. Persistent large clusters are not formed in these simulations. The full charge force fields with JC and KBFF as well as both scaled charge force fields exhibit monotonic increases in the fraction of ions in mid-sized clusters as the salt content increases from 2 wt. % to 20 wt. %. The JC and KBFF combinations show a larger portion of ions in mid-sized clusters than the scaled charge force fields, with some larger clusters at higher concentrations. The scaled charge force field combinations exhibit increases in mid-sized clusters with increasing concentration and almost no large (s > 10) clusters.

In summary, the propensity for NaCl to aggregate is dominated by the choice of ion force field. This becomes particularly apparent at 650 K where the system may be close to the critical point. Force field combinations with SPC/E water tend to result in slightly less NaCl clustering than with TIP4P/2005 water. The force field combinations can be ordered by their propensity to aggregate at 650 K as T+D95 ≥ S+D95 > T+KBFF ≥ T+JC ≥ S+KBFF ≈ S+JC > T+SJC ≥ T+M19.

The structure of water and the influence of ions on that structure have long been an area of study and contention.¹⁰² In a seminal work, Hofmeister demonstrated that different ions have varying abilities to assist in the salting out of egg white proteins.^103,104 This inspired the development of the Hofmeister series of ions, ordered in terms of their ability to enhance or diminish the solubility of proteins.¹⁰⁴ This effect on protein solubility has been attributed to the notion that ions have long range effects on the water structure.¹⁰² There is some debate on how far the influence of ions on the water structure extends into the solution and how different ions may exert that influence.

Figure 8 highlights how increasing the temperature from 450 K to 750 K decreases the average number of hydrogen bonds per water molecule $nHB$⁠. This trend is in qualitative agreement with what has been observed for pure water, where increasing thermal energy disrupts the water structure.^1,105–107 It is also in agreement with previous work that examined hydrogen bonding in high temperature aqueous NaCl solutions.^51–54 Force field combinations with SPC/E water have a slightly lower $nHB$ value than combinations with TIP4P/2005 water. This discrepancy increases with increasing temperature and is particularly apparent in Fig. 9(a) at 2 wt. % NaCl and 650 K.

Previous research has found a correlation between decreasing density and decreasing $nHB$ for supercritical water.^53,107–110 In a comparison of water force fields, Vega and Abascal found that both SPC/E and TIP4P/2005 exhibit lower densities than those in the experiment.⁶⁵ Furthermore, SPC/E exhibits notably lower densities than TIP4P/2005 at elevated temperatures (see the coexistence curves in Refs. 68, 96, and 97). This suggests that at 2 wt. % NaCl, the force fields with SPC/E water are likely approaching the vapor–liquid phase envelope at lower temperatures than seen with TIP4P/2005 and, thus, show larger changes in density. This is consistent with the larger differences in density for SPC/E force field combinations when compared to the Driesner EOS (see Fig. 3(b) at 2 wt. % NaCl). If $nHB$ is correlated with the density, then the growing difference in $nHB$ with increasing temperature, e.g., as observed when comparing SPC/E and TIP4P/2005 force field combinations in Fig. 9(a) at 2 wt. % NaCl, could be attributed to the SPC/E force field combinations being closer to the phase envelope.

Increasing NaCl concentration further also decreases $nHB$ at temperatures of 450 K–600 K (Fig. 8). At 650 K, the $nHB$ for all of the force fields exhibit a change in trend with increasing NaCl concentration. The initial increase in $nHB$ as the concentration is increased from 0 wt. % to 10 wt. % NaCl could be correlated with proximity to the two phase region, an effect that lessens as more salt is added to the system and the critical temperature and pressure increase (see Sec. III A). The initial increase in $nHB$ with increasing salt concentration may, therefore, be related to the state point at 650 K and 200 bar being further away from the phase envelope upon addition of salt. At 700 K–750 K, the density is low for the JC, D95, and KBFF NaCl force field combinations that exhibit an increase in $nHB$ with increasing salt concentration. The scaled charge force fields (M19 and SJC) show a dramatic increase in $nHB$ with increasing salt concentration at 700 K. This can be largely attributed to the transition from being at low density for 2 wt. %–10 wt. % NaCl to higher density at 15 wt. %–20 wt. % (see Fig. 3).

Figure 9(a) demonstrates that the average number of hydrogen bonds per water molecule $nHB$ at 2 wt. % NaCl is dominated by the water force fields. Most noticeably at 650 K, force field combinations with SPC/E show a lower $nHB$ than those with TIP4P/2005. This could be due to SPC/E having a lower critical temperature and density than TIP4P/2005, as was discussed earlier. At 20 wt. % NaCl, the choice of NaCl force field has a larger impact. This is most notable at 700 K, where the scaled charge force fields exhibit significantly more hydrogen bonding, likely due to the higher densities for these model combinations (see Fig. 3).

Figures 9(a) and 9(b) can be used to probe whether or not differences in $nHB$ are density driven. It appears that the water force field dominates the trend at low concentrations (2 wt. % NaCl) and the trend in solution densities matches that of $nHB$⁠. At 20 wt. % NaCl, the trends in $nHB$ seen in Fig. 9(a) are no longer as strongly correlated with the solution density trends shown in Fig. 9(b), especially at temperatures of 650 K and below. A reason for the differences in the hydrogen bonding trends may be that they are a byproduct of the differences in phase equilibria that arise from differences in the water and NaCl force fields.

The extent to which ions disrupt the water structure has also been an area of extensive study. A number of inquiries suggest that the disruption to the hydrogen bonding network in water is limited to the first hydration shell surrounding the ions,^50,111–116 although this effect may be concentration dependent.⁵⁵ One proposed explanation is that the ions act as hydrogen bond acceptors/donors.¹¹⁴ By replacing one of the hydrogen bonds, the ions then disrupt the water–water hydrogen bonding network within the first solvation shell without affecting the extended hydrogen bond network.

Localized disruption to water hydrogen bonding was calculated for 2 wt. % and 20 wt. % NaCl at 450 K and 650 K (supplementary material, Fig. S10). First solvation shell water molecules were defined as those with an oxygen atom that was within a distance r_io of an ion. This distance was chosen based on the radial distance to the first minimum in the ion–oxygen radial distribution function. The $nHB$ for waters both within and outside the first solvation shell were averaged separately. The $nHB$ for the first solvation shell and bulk water matched to within the standard deviations of the independent simulations (supplementary material, Fig. S10). If ions occupy a significant fraction of the hydrogen bond acceptor/donor positions, one would expect the average number of hydrogen bonds between water molecules to be measurably lower in the first solvation shell. Since there is no statistically significant difference between $nHB$ in the first solvation shell and $nHB$ in the remaining water, it implies that the ionic influence on hydrogen bonding in the first solvation shell does not dominate. This may be due to only a small number of water molecules forming hydrogen bonds with an ion in the first solvation shell. Jedlovszky et al. found that only the water molecules closest to the solute were able to assume orientations that replaced a hydrogen bond.⁴⁴ Water molecules that were further from the central ion, but still within the first solvation shell, were free to assume independent orientations.

However, the first solvation shell and bulk water molecules both show the same decrement in $nHB$⁠, and therefore, a more likely explanation is that the hydrogen bonding disruption at these concentrations extends beyond the first solvation shell. This is in qualitative agreement with several experimental studies. Okhulkov and Gorbaty performed x-ray diffraction experiments that produced pair correlation functions.¹¹⁷ The addition of NaCl to liquid and supercritical water altered these correlation functions, which was interpreted as being indicative of hydrogen bond disruption. Femtosecond elastic second harmonic experiments performed by Chen et al. at 297 K showed water orientation order perturbations at concentrations as low as 6 × 10⁻⁶ wt. % NaCl.¹¹⁸ These perturbations increased and then started to level off at 3 × 10⁻⁴ wt. % NaCl. This leveling off was attributed to overlapping electrostatic spheres of influence from the ions that could lead to spatial correlations in the water structure. Kondoh et al. used dielectric spectroscopy to demonstrate that dissolved ions weaken the hydrogen bonding of water.¹¹⁹ They too hypothesized that the disruption of the water structure extends beyond the first solvation shell of an ion.

Dielectric screening provided by the solvent is important in determining ion aggregation/clustering and solubility. Thus, differences in the predicted static dielectric constant of the water were examined. Although we have not included contributions from the ions, this approximation still provides useful insight into how the introduction of ions alters the dielectric response of the water.^34,35,39 In Fig. 10, we show the static dielectric constants (ε) for the water molecules in simulations at 200 bar and 450 K–750 K as a function of NaCl concentration. Force field combinations with SPC/E water have a higher dielectric constant than those with TIP4P/2005. This is in agreement with prior simulations of the pure solvent at ambient conditions, where the respective static dielectric constants were found to be 68 and 55.3.⁶⁶ Both water models underpredict the dielectric constant under ambient conditions, which is 78.5, by a substantial amount, although SPC/E performs somewhat better. As the salt concentration increases, the static dielectric constant decreases for all of the force field combinations. This is known as the dielectric decrement and has been studied extensively in both simulations and experiments under ambient conditions.^{34–39,120–125} Although some previous work has examined the dielectric response of aqueous NaCl at elevated temperatures, the highest temperatures studied by the experiment¹²⁴ and simulation¹²⁶ are not very close to supercritical conditions (333 K and 473 K, respectively).

Force field combinations involving SPC/E showed slightly less ion clustering than their TIP4P/2005 counterparts (see Fig. 7). This is consistent with a higher dielectric solvent (SPC/E) providing better electrostatic screening between charged species. However, there is a subtle correlation between larger NaCl aggregates and the static dielectric constant of the water. For example, as seen in Fig. 7 at 450 K, the D95 force field combinations lead to a higher degree of ion aggregation. The ε values in Fig. 10 for 450 K show that the D95 force field with both water models has a higher ε, while the JC, KBFF, M19, and SJC force fields have similar and lower ε. Force fields that form more (and larger) ion clusters show higher dielectric constants than those where ions are present predominantly as free ions, which is consistent with the recent results of Saric et al.³⁵ Electrostatic screening provided by oppositely charged species within the same cluster may result in reduced disruption to the dielectric response of the water molecules. At higher temperatures, the dielectric constants calculated using the various force fields start to differ, for example, ε for the TIP4P/2005 water with 20 wt. % NaCl; ε values ranked by the force field are D95 > KBFF > JC. This again aligns with the propensity to form clusters for the full change force fields. The M19 and SJC force fields exhibit the opposite trend, although they have a very similar propensity to cluster as well as dielectric constants.

Interestingly, at 650 K, the dielectric constants appear to increase with increasing salt concentration. Recall that for pure water, 650 K is approaching the critical isotherm for both water force fields. Experimentally, addition of salt results in an increase in both the critical temperature and pressure.⁹⁴ The observed increase in the dielectric constant at 650 K may partly result from the decrease in simulation volume (increase in density) associated with moving away from the phase envelope as the salt concentration increases. At 700 K–750 K, almost all of the force field combinations are in a vapor or low density supercritical fluid phase and, thus, exhibit low dielectric constants. The exception is the scaled charge models with TIP4P/2005, which remain in a liquid phase at 700 K and 20 wt. % (see Fig. 3) and, therefore, have somewhat larger dielectric constants.

To provide a foundation for applying molecular simulation to supercritical desalination involving salt precipitation from a single phase fluid, we have performed a series of NpT and NVT ensemble MD simulations with five different salt force fields (JC, D95, KBFF, M19, and SJC) combined with TIP4P/2005 and/or SPC/E water. Simulations were performed at 200 bar for temperatures ranging from 450 K to 750 K. The results suggest that increasing temperature at constant pressure results in increased NaCl clustering for all force field combinations studied. Simulations were also performed along isobars ranging from 100 bar to 600 bar for JC NaCl salt with SPC/E and TIP4P/2005 water. Increasing the pressure at constant temperature results in decreased NaCl clustering.

All of the force field combinations reproduce the qualitative trends in density as represented by the Driesner equation of state for aqueous NaCl solutions. Combinations of salt and water force fields that are not specifically parameterized for the intended use should be treated with caution. A recent work by Dopke et al.⁶⁹ has shown that ion force field parameters optimized for use with a specific water model are not necessarily accurate when used with other water force fields and some errors can be introduced through inaccuracies in the combining rules. However, the results for T+JC and T+KBFF demonstrate that it is possible to replicate the density quite well, although other thermodynamics properties may not be accurately reproduced (such as the static dielectric constant). Comparisons to other experimental data would be needed to assert that T+JC is always the better force field combination (of the ones examined here). Despite the variation in clustering behavior predicted by the different force field combinations, there was no apparent correlation to the density predictions over the range of salt concentrations examined.

Contrasting the various force field combinations shows that the salt force field exerts a stronger influence than the water force field over the extent of NaCl clustering, although force field combinations involving SPC/E water exhibited slightly less clustering than those with TIP4P/2005 water. The propensity for ion aggregation across the force field combinations can be ranked at 650 K as T+D95 ≥ S+D95 > T+KBFF ≥ T+JC ≥ S+KBFF ≈ S+JC > T+SJC ≥ T+M19. Comparisons between the SJC and M19 force fields indicate that the scaling of the ionic charges from ±1.0 to ±0.85 drives a strong decrease in the degree of clustering exhibited by a salt force field. Differences in the Lennard-Jones parameters played a lesser role in clustering for the scaled charge force fields used in this study. Qualitatively, all of the force field combinations predict that while the number of free ions decreases with increasing temperature and concentration, the resulting aggregates present in solution are not just simple ion pairs (Na₁Cl₁) as supercritical conditions are approached. Larger clusters form.

Hydrogen bonding decreases with increasing temperature for all force field combinations. There is also a clear dependence on salt concentration. As NaCl concentration increases, hydrogen bonding generally decreases. Salt force fields combined with SPC/E water yield a lower average number of hydrogen bonds than the same salt force fields with TIP4P/2005 water. At low NaCl concentrations (2 wt. %), the number of hydrogen bonds is dominated by the water force field. At high NaCl concentrations (20 wt. %), the NaCl force field has a stronger influence on the average number of hydrogen bonds, indicating that the extent of hydrogen bonding does not depend only on the choice of water force field.

Calculated static dielectric constants for the water show that the water force field primarily determines the dielectric constant, especially at lower salt concentrations. This is not surprising since the method employed to calculate the dielectric constant only accounts for contributions from the water. All of the force field combinations qualitatively capture the trends in dielectric decrement reported under ambient conditions. There are, however, rather subtle differences between the salt force fields when combined with the same water force field. These differences would seem to indicate that ion clustering (i.e., the formation of larger ion clusters) can increase the solvent dielectric constant. This could be due to electrostatic screening within the cluster reducing the influence of cluster charges on the surrounding water.

Accurate simulations of salt precipitation processes under supercritical conditions require a salt force field that replicates the correct clustering and nucleation behavior. Dielectric response and conductance can serve as indirect measures of ionic aggregation and could be used to validate which force fields best replicate ion clustering. Given the degree of NaCl aggregation we have observed in these simulations, calculations of the dielectric response should take into account contributions from the ionic aggregates in solution as well as the water solvent. This requires a more rigorous treatment of the dielectric and conductivity properties under these conditions, as has been outlined by Rudas et al.⁹³ There is also currently a critical lack of experimental data for the dielectric response and conductance in salt solutions at elevated temperatures and pressures. Future experimental work should address this need. Recent studies of NaCl nucleation using rigid force field models have illuminated a two-step nucleation mechanism for NaCl beyond the spinodal, involving amorphous salt clusters.^127–129 Furthermore, simulations indicate that the nucleation rate is sensitive to the force field used and that polarizable force fields are apparently necessary to accurately predict nucleation rates.^130,131 Additional work is necessary to explore NaCl nucleation under supercritical conditions.

This study was focused on exploring the potential use of simulation in the design of single phase supercritical desalination with selective salt precipitation, with the goal of establishing confidence levels in the classical force fields. Subsequent efforts should be directed toward predicting solubility limits (beginning with just NaCl in water), which will require calculation of chemical potentials of salt and salt–water aggregates. To ultimately describe multi-component brines involving a wide suite of metal ions and associated ligands, additional classical force field models and mixing rules must also be developed and tested. A combination of simulation and experiment will ultimately be needed to determine whether or not the value of co-products plus the savings realized by avoiding brine disposal costs justifies the investment associated with the elevated temperature and pressure operations required for supercritical desalination. Finally, we note that molecular simulation may also prove useful in the design of supercritical desalination processes involving a concentrated liquid brine in equilibrium with a supercritical water phase. With accurate force fields, simulation could be used to locate the phase boundaries and critical points for such systems as a function of salt composition, temperature, and pressure. This may require the use of more accurate polarizable force field models. We also note that accurate chemical potential calculations require either specialized simulation software (e.g., osmotic ensemble Monte Carlo¹³²) or a very large number of simulations.⁵⁶ Moreover, accurately predicting phase boundaries near the critical point and salt solubility remains a challenge for any classical force field.

The supplementary material provides additional data and discussion to support arguments made in the main paper. This includes a comparison of the cluster analysis, hydrogen bonding, and densities for systems with 5000 vs 1000 water molecules (and equivalent concentrations of NaCl). We also compared the results of performing our cluster and hydrogen bond analysis in the NVT ensemble, rather than the NpT ensemble. The supplementary material further includes alternative representations of radial distribution functions (RDFs) and the specific cutoff distances used to define ion clusters. A discussion is provided on how state points were evaluated for stability when in proximity to the critical point and additional data for state point densities to contextualize the salt-containing solution densities with their pure solvent counterparts. Differences between the Driesner and Anderko–Pitzer equations of state for the water–NaCl system are explored. Attempts to identify the source of the differences in hydrogen bonding across force field combinations and potential correlations between hydrogen bonding and the number of charges within an ionic cluster are also discussed.

This work was supported by the Laboratory Directed Research and Development Program at Los Alamos National Laboratory under Project No. 20190057DR. This research also used resources provided by the Los Alamos National Laboratory Institutional Computing Program, supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. T.J.Y. acknowledges support from the Director’s Postdoctoral Fellow Program under Award No. 20190653PRD4. The authors also wish to thank the anonymous referees for their useful suggestions, which improved the manuscript.

The data that support the findings of this study are available from the corresponding author upon reasonable request.