A New Force Field for OH– for Computing Thermodynamic and Transport Properties of H2 and O2 in Aqueous NaOH and KOH Solutions

The thermophysical properties of aqueous electrolyte solutions are of interest for applications such as water electrolyzers and fuel cells. Molecular dynamics (MD) and continuous fractional component Monte Carlo (CFCMC) simulations are used to calculate densities, transport properties (i.e., self-diffusivities and dynamic viscosities), and solubilities of H2 and O2 in aqueous sodium and potassium hydroxide (NaOH and KOH) solutions for a wide electrolyte concentration range (0–8 mol/kg). Simulations are carried out for a temperature and pressure range of 298–353 K and 1–100 bar, respectively. The TIP4P/2005 water model is used in combination with a newly parametrized OH– force field for NaOH and KOH. The computed dynamic viscosities at 298 K for NaOH and KOH solutions are within 5% from the reported experimental data up to an electrolyte concentration of 6 mol/kg. For most of the thermodynamic conditions (especially at high concentrations, pressures, and temperatures) experimental data are largely lacking. We present an extensive collection of new data and engineering equations for H2 and O2 self-diffusivities and solubilities in NaOH and KOH solutions, which can be used for process design and optimization of efficient alkaline electrolyzers and fuel cells.


INTRODUCTION
Modeling aqueous alkaline solutions is of interest for a broad array of manufacturing and separation processes. 1,2 Aqueous alkaline solutions containing potassium hydroxide (KOH) and sodium hydroxide (NaOH) are often used for electrolysis and in fuel cells due to their high ionic conductivities and low cost. 3−7 NaOH and KOH have solubilities exceeding 18 mol/ kg in water at 293 K (and above) 8,9 and significantly influence the thermophysical properties of the solution. 10 The interplay between different thermophysical properties (e.g., densities, viscosities, and ionic conductivities) of aqueous NaOH and KOH solutions influences the product gas purity, and the energy (and Faradaic) efficiency of alkaline-water electrolyzers. 11−13 Knowledge of the thermodynamic and transport properties of hydrogen (H 2 ) and oxygen (O 2 ) gas in aqueous NaOH and KOH solutions is therefore highly relevant for optimization and process design of electrolyzers. 12,14 Modeling electrolyte systems is a challenging endeavor because of the strong long-range ionic interactions, which make the solutions highly nonideal. 1,15−17 Electrolyte solutions are commonly modeled using semiempirical equations of state and molecular based simulations. 1,15−26 Semiempirical equations provide a rapid and convenient method for the prediction of thermophysical properties. 15 The quality of these equations depends on the availability of accurate experimental and simulation data. 18−22 For aqueous alkaline solutions, experimental data for self-diffusivities and solubilities of H 2 and O 2 at high concentrations (above 4 mol/kg), temperatures (323− 373 K), and pressures (above 50 bar) is lacking, especially in the case of aqueous NaOH solutions. 14,27 These temperatures (ca. 353 K) and concentrations (4−12 mol/kg electrolyte solution) are especially relevant for alkaline electrolyzers. 3,7,28−30 Molecular simulations (i.e., molecular dynamics (MD) and Monte Carlo (MC)) can be used as a complementary approach to experiments 31 to provide insight at conditions for which experimental data are limited and difficult to obtain due to high temperatures, pressures, and the corrosiveness of the solution (in case of strong alkaline solutions).
Molecular simulations of electrolyte systems can be studied using either ab initio simulations or force-field-based methods. 1,32,33 Ab initio simulations have the potential to more accurately describe the structure and solvation of the ions, 32,34 but these simulations are computationally expensive and are limited to systems comprising hundreds of atoms for time scales of the order of pico-seconds. To precisely calculate transport properties of fluids, long simulations of several nanoseconds are essential. 24,35 To account for ion−ion and ion−water interactions at high electrolyte concentrations, water molecules need to be modeled explicitly, 36,37 which makes the computations more costly. To overcome both the time and system-size restrictions of ab initio calculations, forcefield-based methods are usually preferred for large-scale production of thermophysical data.
Force fields for aqueous electrolytes can be polarizable or nonpolarizable. 38 38,47 and used in combination with the TIP4P/2005 water model. 38 These force fields yield reasonable predictions for densities, dynamic viscosities, and self-diffusivities of aqueous electrolytes with a scaled charge of 0.85 for concentrations up to 4 mol/kg salt. 38 However, the dynamic viscosities computed using the Madrid-2019 force field deviate from experiments at higher molalities. To address this, Vega and co-workers have developed a new force field called the Madrid-Transport with a scaled charge of 0.75. 47 This force field can accurately predict dynamic viscosities of aqueous NaCl and KCl solutions up to their solubility limit. 47 Despite the importance of alkaline systems, there is no Madrid-force field for OH − to accurately predict densities and dynamic viscosities of aqueous NaOH and KOH systems. Existing OH − force fields are often used to simulate the solvation energy 48−50 and structure, 51−56 and cannot be used directly in combination with the TIP4P/2005 water model and the Madrid-force fields 38,47 for Na + and K + ions as they do not use scaled charges of -0.85 or -0.75.
Here, we propose several nonpolarizable two-site OH − force fields with scaled charges of −0.85 and −0.75, respectively. One of the newly proposed OH − force fields with a scaled charge of −0.75 yields accurate predictions for both densities and dynamic viscosities of aqueous NaOH and KOH solutions for concentrations ranging from 0 to 8 mol/kg, at temperatures ranging from 298 to 353 K. We use this force field to compute the self-diffusivities of H 2 and O 2 in aqueous NaOH and KOH solutions using MD. Solubilities of these gases as functions of concentrations, and temperatures, and pressures are computed using Continuous Fractional Component Monte Carlo (CFCMC) simulations. 57−59 Our data, obtained from molecular simulations, are compared to available experimental data on H 2 and O 2 in KOH solutions. Our simulations can adequately describe the trends observed in experiments for variations in both concentration and temperature. The selfdiffusivities and solubilities of H 2 and O 2 in NaOH and KOH solutions are then fitted to semiempirical engineering equations. These engineering equations can be used for process modeling, and for optimizing electrolyzers and fuel cells. 12 This paper is organized as follows. In section 2, details on the force fields are provided, and the molecular simulation (MD and MC) techniques are explained. In section 3, force field optimization of OH − is discussed, and the results for viscosities, H 2 and O 2 self-diffusivities, and solubilities at temperatures ranging from 298 to 353 K are provided. Our conclusions are summarized in section 4.  Figure S1). This force field is used for computing self-diffusivities of H 2 in NaOH and KOH solutions. The Marx model yields significantly more accurate H 2 solubilities than the Vrabec model in pure TIP4P/2005 water (see Figure S1) and is used for computing H 2 solubilities in NaOH and KOH solutions. For the K + and Na + ions, the Madrid-Transport (+0.75) 47 and Madrid-2019 (+0.85) 38 force fields are used (parameters listed in Table 1). For OH − , several force fields are proposed in this work. The details for OH − force field are discussed in section 3.1. All force fields considered in this work are rigid. All interaction parameters for the TIP4P/2005 water, H 2 , and O 2 models are provided in the Supporting Information (Tables S1−S3). The Lennard-Jones (LJ) and Coulombic interactions are considered for modeling the intermolecular interactions. The Lorentz−Berthelot mixing rules 65 . Analytic tail corrections for energies and pressures are applied to the LJ part of the potential. The cutoff radius for both LJ and Coulombic potentials is set to 10 Å. The particle−particle particle-mesh (PPPM) 66,70 method is used for long-range electrostatic interactions with a relative error of 10 −5 .

METHODOLOGY
The OCTP tool 71 is used in LAMMPS to calculate the transport properties. The simulations are initially equilibrated in the NPT and NVT ensembles for a period of ca. 2 ns. Production runs (in NVT) of 10−50 ns are used to calculate dynamic viscosities and self-diffusivities. To obtain an ensemble mean and a standard deviation, each calculation is repeated 5 times with a different random seed for the initial velocity. The Nose−Hoover thermostat and barostat 66,72,73 are used, with a coupling constant of 100 and 1000 fs, respectively. The modifications of the Nose−Hoover for rigid bodies, proposed by Kamberaj,74 are used in LAMMPS. The densities and transport properties are calculated in a simulation box containing 700 H 2 O molecules. The corresponding numbers of NaOH and KOH molecules, in combination with the respective molarities are provided in Tables S4 and S5. All initial configurations are created using the PACKMOL software. 75 Two gas molecules (H 2 or O 2 ) (corresponding to infinite dilution) are used to calculate self-diffusivities of the gases in the aqueous NaOH and KOH solutions. All selfdiffusivities, computed from mean-square displacements, 71 are corrected for finite-size effects using the Yeh−Hummer equation: 76 where D i MD and D i denote the self-diffusivities calculated by MD and corrected for finite-size effects for species i, respectively, k B is the Boltzmann constant, T is the temperature (in K), ξ is a dimensionless number equal to 2.837298 for a cubic simulation box, and L is the length of the simulation box. 76,78 The dynamic viscosities (η), obtained from the MD simulations, do not have finite-size effects. 78 57,84,85 All molecules are considered as rigid, and only intermolecular LJ and Coulombic interactions are considered. A cutoff radius of 10 Å is used for both the LJ and Coulombic interactions. The Ewald summation with a relative precision of 10 −6 is used for the electrostatics. Analytic tail corrections for energies and pressures are applied to the LJ part of the potential. 66 Periodic boundary conditions are imposed in all directions. All MC simulations contained 300 water molecules. For all the compositions considered, the corresponding numbers of NaOH and KOH molecules and the respective molarities are provided in Table S4 and S5. Simulations are carried out at temperatures of 298, 323, 333, 343, and 353 K, at H 2 and O 2 pressures of 1, 50, and 100 bar.
To calculate the excess chemical potentials and solubilities of H 2 and O 2 , "fractional" molecules are introduced. In contrast to "whole" or normal molecules, the interactions of "fractional" molecules with other molecules are scaled with a continuous order parameter λ (in the range of [0, 1]): 57,86 λ = 0 indicates no interactions between the fractional and whole molecules (ideal gas), while λ = 1 indicates full interactions, corresponding to a "whole" or normal/unscaled molecule. For more details regarding the scaling of the interactions of fractional molecules the reader is referred to refs 87−89. A single fractional molecule of H 2 (and O 2 ) is used to calculate the excess chemical potentials of the respective molecules in the solution. All other molecules in the simulation are whole molecules. The Wang−Landau algorithm 90,91 is used to construct a biasing weight function for λ (W(λ)). The biasing weight function helps in overcoming possible energy barriers in λ-space, to ensure a flat observed probability distribution. 83 100 bins are used to obtain a histogram of λ values, thereby computing the probability of occurrence for each λ value. The Boltzmann average of any parameter (A) can be computed using 83 The infinite dilution excess chemical potential (μ ex,∞ ) can be related to the Boltzmann sampled probability distribution of λ (p(λ)) using 57,83 where p(λ = 1) and p(λ = 0) are the Boltzmann sampled probability distribution of λ at 1 and 0, respectively. The molarity based Henry coefficient (H) is defined as 83 in which f i is the fugacity of a solute in the gas phase, m i is the molarity of the gas in the solution (mol/L), and m 0 is set to 1 mol/L. The infinite dilution excess chemical potential of H 2 and O 2 can be related to the molarity based Henry coefficient using 83 where R is the universal gas constant. For all simulations, 4 × 10 5 equilibration cycles are carried out followed by 4 × 10 5 production cycles. A cycle refers to N number of trial moves, with N corresponding to the total number of molecules, with a minimum of 20. Trial moves are selected with the following probabilities: 1% volume changes, 35% translations, 29% rotations, 25% λ changes, and 10% reinsertions of the fractional molecules at random locations inside the simulation box. The maximum displacements for volume changes, molecule translations, rotations, and λ changes are adjusted to obtain ca. 50% acceptance of trial moves. For each condition (concentration, temperature, pressure), 20 independent simulations are performed. The final Boltzmann probability distributions of λ are averaged in blocks of 4 to obtain 5 independent averaged distributions. For all averaged distribu-tions, the excess chemical potentials and Henry coefficients are calculated to obtain a mean value and the standard deviation. All the raw data for the MD and MC simulations are shown in Tables S6−S11.  Figure 1 shows the variation of densities as functions of electrolyte concentrations for both NaOH and KOH. By adjusting the value of σ OO , it is possible to obtain an excellent agreement for all the different models. All the densities obtained deviate less than 2% from experimental fits found in literature. A larger negative charge on O (q O ) results in a larger optimum σ OO parameter, to counteract the strong attractive Coulombic interactions. The experimental fits of Olsson 94 (for densities and viscosities of aqueous NaOH), Gilliam 93 (densities of aqueous KOH), and Guo 95 (viscosities of aqueous KOH) are used and shown as lines in Figure 1.

RESULTS AND DISCUSSION
The dynamic viscosities of aqueous NaOH and KOH solutions calculated using FF1−FF4 are shown in Figure 2. It can be observed that the choice of the total charge (q OH ), and the resulting σ OO has a significant influence on the viscosities, especially at higher concentrations in which the influence of ion−ion interactions become more important. The influence of ion size on the viscosities and densities is shown in Figure S2. In case of FF2 (with q OH = −0.85), the dynamic viscosity is overestimated by more than a factor 3 compared to the experimental fit for the highest concentration of NaOH. For aqueous KOH, the FF2 model overestimates the dynamic viscosity by around 40% at the highest concentration of KOH. The FF1, FF3, and FF4 models with q OH = −0.75 show a much better agreement with the experimental fit. The findings of the Madrid-Transport model for aqueous NaCl and KCl solutions 47 also show that a scaled charge of 0.75 leads to better predictions of transport properties (especially at concentrations above 4 mol/kg salt) compared to a scaled charge of 0.85. Overall, the FF1 model shows the best agreement with the experimental viscosities and densities. For this reason, only the results of the FF1 model will be used and discussed further in this work.
where g w is the anion/cation−O W RDF, r is the radial distance, r min is the position of the first minimum in the RDF, and ρ w is the number density of water in the solution. Our results show a first peak at approximately 2. 13 53 The reported hydration numbers (in the first shell) are in the ranges of 4−8 and 6−8 for Na + and K + , respectively. 96 For OH − , the results show a first peak at approximately 2.75 Å for OH − −O W , with hydration numbers of 4.8 and 5.9 for KOH and NaOH, respectively, at a molality of 5 mol/kg. Other molecular simulations in literature report a first peak ranging from 2.3 to    Table  3. Even though for Na + and K + the self-diffusivities at infinite dilution are close, this is not the case for OH − (underestimated by a factor ca. 5). For reasonable values of σ OO , ϵ OO , and q O , we could not obtain OH − self-diffusivities close to the values reported by experiments 97 without causing significant deviations from experimental densities and viscosities. This result is expected as classical OH − models cannot capture the details of the solvation of OH − in water and the proton transfer mechanism, which lead to anomalously high OH − mobilities as discussed by Tuckerman et al. 34,98 As such, our model, similarly to other classical force fields, is not suitable for predicting OH − diffusivities of NaOH and KOH. Since electrical conductivities vastly depend on the mobility of the OH − ions in the solution, the new OH − model presented here is unable to accurately predict electrical conductivities of aqueous NaOH and KOH solutions. Although our classical force field cannot capture the proton transfer mechanism, it can correctly predict the dynamic viscosities of the electrolyte solutions. As the aim of this study is to study the transport properties and solubilities of H 2 and O 2 gas in aqueous NaOH and KOH electrolytes, correct predictions of densities and viscosities are sufficient. Developing an OH − force field by taking into account the proton transfer mechanism and accurate OH − mobilities is beyond the scope of this work as quantum mechanical based force fields will be required.

Densities and Viscosities.
It is important to show that the NaOH and KOH models (FF1 OH − model, and the Madrid-Transport models of Na + , and K + ) 47 can accurately predict the temperature-dependence of densities and viscosities. Figure 4 shows the densities and viscosities at different temperatures for both NaOH and KOH solutions. The agreement between MD simulations and experimental fits is excellent for aqueous KOH. For aqueous NaOH solutions, the results of densities are overestimated by ca. 2% and for dynamic viscosities by ca. 20% at the highest concentration (molality 8 mol/kg). Despite this, the trends of densities and viscosities for variations of electrolyte concentration and temperature are well-predicted by the MD simulations using the new force fields. Densities and viscosities show a much weaker dependence on pressure (in the range of 1 to 100 bar) compared to temperature (in the range of 298 to 353 K) due to the incompressibility of the liquid phase. The variations of densities and viscosities as a function of pressure are shown in Figure S4. 3.3. Self-Diffusivities of H 2 and O 2 in Aqueous NaOH and KOH. The finite size-corrected self-diffusivities (using eq 1) of H 2 and O 2 in aqueous NaOH and KOH solutions calculated using MD simulations at various temperatures are shown in Figure 5 Table 4. The experimental diffusivities of O 2 in NaOH solution at 296 K (black) provided by Zhang et al. 14 are plotted as points. The FF1 OH − model, in combination with the TIP4P/2005 water model, 42 the Bohn O 2 model, 62 the Vrabec H 2 model, 63 and the Madrid-Transport Na + and K + models, 47 is used for the MD simulations.    Figure S5    The values for f 0 (10 −1 (L/mol)), f 1 (10 −4 (L/mol) K −1 ), f 2 (10 −3 (mol/L)),   Figure  S5.

Solubilities of H 2 and O 2 in
Aqueous NaOH and KOH. In Figure 6, the H 2 and O 2 solubilities obtained using CFCMC calculations are shown as functions of NaOH and KOH concentrations. In this figure, only the results at 298 and 333 K are shown as solubilities (especially at higher electrolyte concentrations) vary only weakly in the temperature range of 298−353 K. The solubilities of H 2 and O 2 at 298, 323, 333, 343, and 353 K are shown in Figure S7.
As a comparison the experimental data provided by Walker et al. 99 on the solubilities of H 2 and O 2 in aqueous KOH are fitted and plotted in Figure 6a,b. This experimental data are also in agreement with the experiments of Davis et al. 102 for O 2 solubilities (at 298 and 333 K) and with the Sechenov model. 101 The Sechenov model 100 where C G and C G,0 are the solubility of the gas in the electrolyte and pure water at 1 bar, respectively. f 0 and f 1 are fitting constants. The temperature dependence of the parameter C G,0 can be fitted as where f 2 −f 4 are additional fitting parameters. The optimized fitting parameters for MC simulations in this work and the experimental data of Walker et al. for H 2 and O 2 solubilities are shown in Table 5. Equation 8 provides an excellent fit for both the simulation results found in this work and the experimental data present in the literature, as shown in Figure S7.

CONCLUSIONS
The self-diffusivities and solubilities of H 2 and O 2 in aqueous NaOH and KOH solutions are modeled using MD and CFCMC simulations. A new two-site nonpolarizable OH − force field (FF1 model) is proposed with a scaled charge of −0.75, which matches with the TIP4P/2005 water and the Madrid-Transport models for Na + and K + . Although our classical force field cannot capture the proton transfer mechanism, which influences the OH − diffusivities, it can predict the densities, dynamic viscosities, and the salting out of H 2 and O 2 in aqueous NaOH and KOH solutions. Excellent agreement is observed between simulation and experimental data for both densities and dynamic viscosities of NaOH and KOH for a concentration range of 0−6 mol/kg and a temperature range of 298−353 K. This model is used to generate self-diffusivity and solubility data for H 2