key: cord-0518369-54omvg2p
authors: Tran, Leann; Skvara, Jiri; Smith, William R.
title: Atomistic Simulation Framework for Molten Salt Vapor-Liquid Equilibrium Prediction and its Application to NaCl
date: 2022-03-07
journal: nan
DOI: nan
sha: f69b2f3996fdac2ce0f43349c8639f9ef7ffa642
doc_id: 518369
cord_uid: 54omvg2p

Knowledge of the vapor-liquid equilibrium (VLE) properties of molten salts is important in the design of thermal energy storage systems for solar power and nuclear energy production applications. The high temperatures involved make their experimental determination problematic, and the development of both macroscopic thermodynamic correlations and predictive molecular-based methodologies are complicated by the requirement to appropriately incorporate the chemically reacting vapor-phase species. We derive a general thermodynamic-based atomistic simulation framework for molten salt VLE prediction and show its application to NaCl. Its input quantities are temperature-dependent ideal-gas free energy data for the vapor phase reactions, and density and residual chemical potential data for the liquid. If these are not available experimentally, the former may be predicted using standard electronic structure software, and the latter by means of classical atomistic simulation methodology. The framework predicts the temperature dependence of vapor pressure, coexisting phase densities, vapor phase composition, and vaporization enthalpy. It also predicts the concentrations of vapor phase species present in minor amounts (such as the free ions), quantities that are extremely difficult to measure experimentally. We furthermore use the VLE results to obtain approximations to the complete VLE binodal dome and the critical properties. We verify the framework for molten NaCl, for which experimentally based density and chemical potential data are available in the literature. We then apply it to the analysis of NaCl simulation data for two commonly used atomistic force fields. The framework can be readily extended to molten salt mixtures and to ionic liquids.

The vapor-liquid equilibrium (VLE) properties of molten salts are important in nuclear energy production [1] [2] [3] [4] and solar power systems [5] [6] [7] . Due to the high temperatures involved, the experimental determination of these properties is problematic, and most studies for even the simplest case of pure molten alkali halides date prior to 1980 (e.g., Janz et al. [8] ).

Atomistic molecular simulation is a promising alternative approach for molten salt property prediction, but the vast majority of such studies to date have focussed on density and transport properties (for recent reviews, see Wang et al. [9] and Sharma et al. [10] ), and its use for VLE prediction has been limited. In a 1992 study, Panagiotopoulos [11] applied the Gibbs Ensemble Monte Carlo algorithm [12, 13] to the restricted primitive model (RPM) force field (FF) with its diameter fitted to the NaCl liquid density. Guissani and Guillot [14] in 1994 performed molecular dynamics (MD) simulations in the N V E ensemble (N is the number of particles, V the volume, and E the internal energy) using the more realistic Born-Huggins-Mayer-Fumi-Tosi (BHMFT) FF [15, 16] with parameters previously reported by Lewis and Singer [17] . They calculated the VLE binodal curve from an empirical equation of state fitted to simulation data for sets of P v isotherms exhibiting van-der-Waals loops (P is the pressure and v is the molar volume). The same FF and variants thereof were studied in 2007 by Rodrigues and Fernandes [18] (RF), who calculated chemical potentials by the integration of N V T MD simulation isotherms and determined the VLE binodal curve from the equality of the simulated chemical potentials and pressures of the coexisting phases. Abramo et al. [19] recently used a correlation approach to fit BHMFT temperature-dependent ionic size parameters to MD N V T isothermal compressibility simulations, and used these to calculate the resulting VLE envelope using the same methodology as that of RF. Most recently, Kussainova et al. [20] (KMYYP) performed MD simulations of NaCl vapor pressures and coexisting phase densities using two modern electrolyte force fields: that of Joung and Cheatham [21] , and the Madrid FF of Benavides et al. [22] . They used a treatment that combined N V T liquid chemical potential simulations in conjunction with simulations for a vapor phase model of an ideal gas of ion pairs.

A main challenge for VLE molten salt modeling is the appropriate treatment of the vapor phase (see, for example, Zhao et al. [6] ), which requires knowledge of the species present. For NaCl and other alkali halides, the vapor phase is experimentally known to be dominated by the neutral ionpair monomer and its dimer [23] [24] [25] [26] [27] [28] . Since molten salt vapor pressures are very low, the vapor phase may be treated as an ideal gas, and its composition is determined by the relevant chemical reaction equilibria, which must be taken into account in the VLE calculations. These were not considered in any of the aforementioned simulation studies.

The purpose of this paper is to provide a predictive thermodynamic-atomistic simulation framework for molten salt VLE properties, whose novelty is the integration of ideal-gas chemical reaction free energy quantities with atomistic liquid-state chemical potential and density atomistic simulation data. This data is required as functions of temperature at the standard state pressure P 0 = 1 bar.

For many molten salts (e.g., alkali halides), the ideal-gas data are readily available in the NIST-JANAF [29] or other thermochemical compilations. If it is unavailable for a salt of interest, it can be calculated by means of electronic structure methodology such as that incorporated in software such as Gaussian [30] .

The framework also predicts the VLE vaporization enthalpy as a function of temperature, T , in addition to the VLE binodal curve and critical properties. The enthalpy has not been considered in prior simulation studies, and its experimental values are sparsely available in the literature (for example, it is not included in the comprehensive review of molten salt properties by Janz et al. [8] ).

Finally, the framework also predicts the concentrations of vapor phase species that are present in minor amounts in a computationally efficient manner.

The paper is organised as follows. The next section describes the thermodynamic framework. In the subsequent Results and Discussion section, we first validate the framework using experimental NaCl density data [31] and chemical potential data from the NIST-JANAF compilation [29] to calculate the vapor pressure as a function of temperature. We then show the framework's application to the NaCl FFs previously considered by Kussainova et al. [20] , using their liquid density and chemical potential simulation data. We calculate the NaCl vapor pressure, the vaporization enthalpy and the VLE phase diagram and critical properties. This is followed by a section giving an uncertainty analysis of our calculation procedures, and then a final Conclusions section.

For simplicity and generality, we consider a 1-1 electrolyte AB comprised of an A + cation and a B − anion.

A. Vapor Pressure, P * Treating the liquid as an incompressible fluid between P 0 = 1 bar and P * (the experimental isothermal compressibility of NaCl(liq) is only 0.343 per GPa [32] ), its mean ionic chemical potential at (T, P ) can be written as

where µ 0 i (T, P 0 ) is the molar standard chemical potential of species i at T and P 0 , R is the universal gas constant, IG denotes the ideal gas standard state, and µ res;N P T AB,liq is the liquid mean ionic residual chemical potential with respect to the ideal-gas model in the N P T ensemble.

The vapor phase may consist of multiple species, as identified from experimental knowledge or from electronic structure calculations. Based on experimental information for NaCl and other alkali halides [23] [24] [25] [26] [27] [28] , we begin by assuming only the AB monomer and the (AB) 2 dimer species in the vapor phase. We will show how to verify (or negate) this assumption by considering additional species after the initial P * (T ) calculation.

At the low pressures typical of molten salt VLE, the vapor phase species chemical potentials can be modelled by their ideal-gas form:

where y i is the mole fraction of species i. The VLE condition for the vapor pressure P * (T ) arising from the equality of the AB liquid and vapor phase chemical potentials is then ∆G 0 r1 (T, P 0 ) + µ res,N P T AB,liq (T, P 0 ) RT + P 0 ρ AB,liq (T, P 0 )RT P * P 0 − 1 = ln

where

The vapor-phase monomer mole fraction y(T, P * ) (henceforth dropping its subscript for notational convenience) is constrained by the equilibrium condition for its dimerization reaction:

where

P * (T ) is determined from the numerical solution of Eqs. (3) and (5).

Pitzer [28] has argued that in addition to the ring form of the dimer in the vapor phase, inclusion of its linear isomer (NaCl) 2 is needed at higher temperatures. Rather than including both species individually in the VLE calculation, they can be incorporated in the analysis without changing the form of Eq. (5) by the "lumping technique" described in Section 9.7 of Smith and Missen [33] , which expresses the single (lumped) standard chemical potential of the two (AB) 2 isomers as µ (AB) 2 ,ring+linear (T, P ) = µ 0 (AB) 2 ,ring+linear (T, P 0 ) + RT ln

The molar fraction z of the ring dimer within the two-member dimer group is given by [33] 

For future reference, the lumped dimer ideal-gas species standard enthalpy arising from Eq. (7) is 

As previously mentioned, we can calculate the compositions of additional vapor-phase species by first calculating P * assuming that the monomer and dimers are the only species present, and subsequently performing a vapor-phase reaction equilibrium calculation at the resulting (T, P * ) including the additional species (e.g., A(g), A 2 (g), B(g), B 2 (g), A + (g), B − (g)) (this requires knowledge of their standard chemical potential data). This strategy can both serve to verify the neglect of the minor species in the original P * calculation, and provide an approximation to their concentrations.

The calculation may be rapidly performed using any readily available Gibbs energy minimization algorithm [33, 34] . (Iif the monomer and dimer compositions change significantly from the original values in the post hoc calculation, then the original P * calculation must be modified to incorporate the additional species at the outset.)

Finally, we note in passing that Eq. (3) can be written using the identity

to give the alternative form

Standard chemical potentials may be expressed in many forms. Molar standard Gibbs formation energies, ∆G f i (T, P 0 ) are commonly used, but in this paper we use the following expression:

where gef i (T, P 0 ) is the Gibbs Energy Function, ∆H f i (T, P 0 ) is the species standard molar enthalpy of formation and

∆H 0 r1 (298.15, P 0 ) and ∆H 0 r2 (298.15, P 0 ) data are given in the Supporting Information (SI). Ideal-gas species molar enthalpy values are expressed as

and liquid species enthalpy values h i (T, P ) are expressed as

where we have neglected the pressure dependence.

Using the Gibbs-Helmholtz equation for the chemical potentials, the differential of Eq. (3) on the VLE saturation curve is

where y is evaluated at (T, P * ) and

Similarly, the differential of the reaction equilibrium condition of Eq. (5) on the VLE saturation

Eliminating dy from Eqs. (16) and (18) gives the temperature derivative of the vapor pressure:

where

Eq. (20) can be written in the Clapeyron form

or alternatively as

where Z is the compressibility factor. ∆H vap AB (T ) is the enthalpy change per mole of per mole of vapor formed (from (2 − y) mol liquid), given by

A. Validation for NaCl using experimental data

In this section, we validate our thermodynamic framework using experimental data for the required simulation quantities. We henceforth refer to this as the "NJ" (for "NIST-JANAF") methodology.

µ res,N P T NaCl,liq (T, P 0 ) values were calculated from Eq. (1) and species Gibbs Energy Function data in the NIST-JANAF compilation [29] from

We smoothed the data for use in Eqs. (3)-(6) by fitting each of ∆G 0 r1 (T, P 0 )/T , ∆G 0 r2 (T, P 0 )/T and µ res,N P T (T, P 0 )/T to the functional form

We used the interval (1200 K, 2500 K) for the NJ ∆G 0 (T, P 0 ) and µ res, N P T (T, P 0 ) data, in conjunction with the density data of Kirschenbaum et al. [31] (see Table III ), which we fitted to the functional form

over the interval (1149 K, 3200 K). The regression coefficients and the standard regression errors are given in the SI. Table I and for the other cases, it is provided in the SI. At about 2300 K and higher temperatures where the vapor pressure exceeds 10 bar (and no experimental data exist), the vapor phase ideal-gas approximation introduces errors that will increase with temperature. These could be reduced by incorporating second virial coefficients, but this is beyond the scope of the current study. experimental results as follows. blue: Kvande [27] , orange: Horiba [35] , green: Fiock [36] , red: Stull [37] , purple:

von Wartenberg et al. [38] and Ruff [39] , brown: Barton [23] . Table I . NaCl VLE properties calculated by the framework of Section II A using different liquid data, with an ideal-gas vapor phase of monomer plus ring and linear dimers. In all cases, ∆G 0 r1 (T, P 0 ) and ∆G 0 r2 (T, P 0 )

were obtained from the NIST-JANAF compilation [29] , augmented with the Pitzer [28] linear dimer data.

NJ: Kirshenbaum density data [31] and NIST-JANAF µ res N P T NaCl (T, P 0 ) data; JC: Joung-Cheatham FF [21] density and µ res N P T NaCl (T, P 0 ) data; MFC-MSC: MSC density data and the hybrid method for calculating µ res N P T NaCl (T, P 0 ) for the Madrid FF [22] (denoted by Kussainova et al. [20] as the "Madrid rerun"). Table I with the vapor phase containing the indicated additional vapor-phase species. The additional ideal-gas thermochemical data were obtained from the NIST-JANAF compilation [29] . Table II shows vapor-phase ideal-gas reaction equilibrium calculations that include additional species in the vapor phase at the example temperatures 1200 K and 2100 K and the corresponding NJ P * (T ) values in Table I . In addition to providing useful additional information, these indicate that the concentrations of the additional species are very small, and that they thus have an insignificant effect on the originally calculated VLE results. The atomic Na and Cl species have the next largest concentrations, and the ion concentrations are vanishingly small at 1200 K and larger but still negligible at 2100 K. with the experimental data, and will be treated henceforth as the benchmark for comparison with subsequently shown results obtained from simulation data. The green dashed curve, which includes only the vapor-phase monomer, somewhat under-predicts the experimental data and the NJ curve at all temperatures. The black dashed curve, including the monomer and the ring dimer, is in good agreement with the experimental data and the NJ results at temperatures below about 1400 K, but deviates increasingly from them at higher temperatures, consistent with the observation of Pitzer [28] .

In a recent study, Kussainova et al. [20] (KMYYP) performed MD density and µ res N V T simulations for NaCl, using the nonpolarizable FF of Joung and Cheatham [21] (JC) and variants of the Madrid scaled-charged FF of Benavides et al. [22] (MSC). They used this data to apply a thermodynamic framework different from that of Section II to calculate P * (T ).

Finding the MSC liquid chemical potential to significantly overestimate the experimental NIST-JANAF [29] values, they followed an earlier suggestion of Vega [40] and calculated its chemical potential by insertion of a fully charged ion pair into configurations of an MSC molecular dynamics simulation trajectory, which they called the "Madrid rerun" approach. We refer to this in the following as the "MFC-MSC" methodology.

In the next two subsections we briefly review the KMYYP density and chemical potential results and their agreement with experiment, and in the subsequent two subsections we apply the framework of Section II A to calculate P * (T ) and the vaporization enthalpy ∆H vap (T ) (not calculated in the KMYYP study) using their simulation data in conjunction with ideal-gas thermochemical data from the NIST-JANAF compilation [29] .

We used the regression form of Eq. (28) in Section III A to smooth the KMYYP JC and MSC density simulation data [20] over the temperature range for which the data was available (the regression coefficients and standard errors of the fits are given in the SI). The smoothed curves and the underlying data are shown in Fig. 2 in comparison with the experimental data of Kirshenbaum et al. [31] fitted to the same functional form. Both force fields give somewhat low densities in comparison with experiment, and the Joung-Cheatham [21] (JC) simulation results give the best agreement. 

The µ total NaCl,liq (T, P 0 ) liquid data shown in Fig. 2 of Kussainova et al. [20] and calculated from Eqs. (3) and (4) of their paper use NIST-JANAF standard Gibbs formation free energies for Na + and Cl − in conjunction with their density and µ res,N V T NaCl,liq simulation data. We obtained µ res,N P T NaCl,liq (T, P 0 ) simulation values for each FF from Eq. (26) using their total µ total NaCl,liq (T, P 0 ) values. We smoothed the resulting µ res,N P T NaCl,liq (T, P 0 ) values (see the SI for details) using the functional form of Eq. (27) . The results are shown in Fig. 3 for the JC and for both Madrid force fields (the original scaled charge version and the hybrid MFC-MSC "Madrid rerun" calculation), where they are compared with the corresponding experimental quantity obtained from the NIST-JANAF compilation [29] . 6) incorporating the monomer and both dimers in the vapor phase using the KMYYP density and chemical potential data of Figs. 2 and 3 in conjunction with the NIST-JANAF [29] and Pitzer [28] ideal-gas data for ∆G 0 r1 (T, P 0 ) and ∆G 0 r2 (T, P 0 ). They are compared with the KMYYP results (open symbols joined by dotted curves) shown in Fig. 4 of their paper [20] and with the NJ solid blue curve, replicated from Fig. 1 for comparison. Our JC results are about an order of magnitude . higher than and almost parallel to the NJ curve (the latter aspect will be seen in the next section to have important implications for the vaporization enthalpy results). Our MFC-MSC results, indicated by the red dash-dot curve, are lower than the NJ curve at low temperatures and increase linearly with temperature until they intersect the NJ curve at about 2000 K. Although these trends with respect to the NJ curve are similar to the trends followed by their respective µ res N P T results in Fig.   3 , differences between the curves are magnified in Fig. 4 due to the logarithmic term involving P * in Eq. (3). The best overall agreement with the NJ curve is achieved by the MFC-MSC approach.

Interestingly, no solutions to Eqs. (3)-(6) were found to exist in the case of the original scaled-charge Madrid FF (MSC). This is a consequence of the large positive deviation of its µ res N P T from the NJ results in Fig. 4 .

Our results in Fig. 4 using the framework of Section II are seen to differ markedly from the KMYYP [20] results, due to different approaches are used to calculate the ideal-gas ∆G 0 r1 quantity, and an incorrect factor of 2 in Eq. (6) of their paper. We obtained ∆H vap NaCl (T ) from (25) using the liquid residual enthalpy from differentiation of the regression equation for µ res N P T NaCl,liq in Eq. (27):

Values of ∆H vap NaCl (T ) of. (25) are given in Table I and shown in Fig. 5 . The experimentally-based NJ curve shows a weak temperature dependence, and the JC curve is almost flat and lies very close to it. This good agreement is due in large part to the parallelism of the corresponding P * (T ) curves in Fig. 4 The MSC-MFC curve is seen to lie well above the other results.

Experimental data for ∆H vap NaCl (T ) is sparse, and unlike the analysis of Section II B, have been derived indirectly from regressions of ln[P * (T )] curves under various assumptions. Thus, assuming only monomer species in the vapor phase, Fiock and Rodebush [36] calculated a constant value of 180 kJ·mol −1 over the temperature range (1180 K-1429 K); Kelley [41] re-analyzed their data and obtained ∆H vap NaCl = 4.184(52.800 − 0.0069T ) kJ·mol −1 , which gives decreasing values of 186-180 kJ·mol −1 over the same temperature range. Barton and Bloom [23] give a value of 180.9 ± 1.7 kJ·mol −1 at 1339 K.

Except near the critical point, the coexisting vapor densities are very small with respect to the liquid densities, and special care must be taken in modeling them simultaneously in order to avoid the prediction of negative vapor densities. We thus used the method of Rodrigues and Fernandes [18] , which models the liquid-vapor coexistence curve using the following scaling law equations:

where T c and ρ c are respectively the critical temperature and density, (E 0 , E 1 , F 0 , F 1 ) are parameters, and τ = 1 − T /T c , leading to the following expressions for the densities:

Our analysis proceeds by first fitting the density data to Eq. (30) to determine (T c , E 0 , E 1 ), and then using the resulting T c value to determine (ρ c , F 0 , F 1 ) by fitting the data to Eq. (31). The first regression can be treated as a separable least-squares problem with the nonlinear parameter T c and the linear parameters (E 0 , E 1 ), and the second regression is linear in the parameters (F 0 , F 1 ).

Following the determination of T c , we fitted the P*(T) data to the Antoine equation to determined the critical pressure, P c .

The Antoine parameter values are given in the SI High T expt data NJ, low T expt data JC, low T sim. data MFC-MSC, low T sim. data  Figure 6 . NaCl VLE phase diagram. Curves are fits to the indicated data points using the scaling law methodology of Rodrigues and Fernandes [18] and the stars indicate the critical points of the last four rows of Table IV. is given in Table III . The resulting critical properties are given ini the first row of Table IV .

We first re-analyzed the Kirshenbaum et al. [31] VLE data set using the scaling-law-based methodology of Rodrigues and Fernandes [18] . The results in the second row of Table IV (regression coeffi- cients (E 0 , E 1 , F 0 , F 1 ) corresponding to each of the last 4 rows are given in the SI) indicate that the uncertainty intervals overlap those of Kirshenbaum et al. [31] in the first row, and show a considerably smaller ρ c uncertainty range. P c in the second row is blank because Kirshenbaumet al.did not provide P * (T ) in their paper.

To investigate the effect of using only low-temperature VLE data where the vapor phase exhibits ideal-gas behavior, we repeated the Rodrigues and Fernandes [18] analysis using our experimentallybased NJ VLE data at P * < 10 bar from that the critical property uncertainty intervals overlap the (relatively large) uncertainty intervals of Kirshenbaum et al. [31] in the first row. However, they lie slightly lower than those of the second row.

Finally, the last two rows of Table IV , based on results from the framework of Section II A using the simulation data of Kussainova et al. [20] , both give T c values well below those of the experimentallybased results of the first two rows. The JC value of ρ c is very high and the MFC-MSC value is very low.

The VLE coexistence curves are shown in Fig. 6 . The Joung-Cheatham force field gives better overall results than the hybrid MFC-MSC methodology.

D. Uncertainty analysis

The uncertainties ∆ ln(P * /P 0 ) in ln(P * ) calculated from Eqs. (3)-(6) arise from uncertainties in the input quantities {x i } = {∆G r1 (T, P 0 )/RT, ∆G r2 (T, P 0 )/RT, µ res,N P T NaCl,liq (T, P 0 )/RT, ρ liq (T, P 0 )}. The sensitivity coefficients ∂ ln(P * /P 0 )/∂x i are derived in the SI. We use the NIST-JANAF thermochemical data for ∆G 0 r1 (T, P 0 ) and ∆G 0 r2 (T, P 0 ) and consider these quantities to be exact. The density sensitivity coefficient is very small and its contribution to the uncertainty may be neglected, and the sensitivity coefficient with respect to µ res, N P T (T, P 0 ) gives an uncertainty in ln P * of ∆ ln[P * ] = ∂ ln P * ∂(µ res,N P T /RT ) ∆(µ res, N P T (T, P 0 )/RT ) = (2 − y)∆(µ res, N P T (T, P 0 )/RT )

We use chemical potential simulation standard deviations reported by Kussainova et al. [20] , and our calculated values of y at each temperature. We find that the JC and MFC-MSC P * values vary slightly with T and are approximately 6% and 18%, respectively. They are indicated by the error bars in Fig. 4 .

For the uncertainties of ∆H vap AB (T, P * ) of Eq. (25), we used (see the SI for the derivation)

where V is the variance-covariance matrix of the µ res, N V T (T, P 0 )/T regression of Eq. (27),

and

We found that the JC and MFC-MSC ∆H vap AB (T, P * ) values have uncertainties (one standard deviation) of approximately 0.9 kJ·mol −1 and 4.5 kJ·mol −1 , respectively. These are indicated as error bars in Fig. 5 

An approximate uncertainty interval at the 95% confidence level for the T c nonlinear regression of Eq. (30) can be determined by calculating the upper (T upper ) and lower (T lower ) T values that satisfy [43] 

where S is the regression's residual sum of squares, S(T c ) is its minimum at the calculated T c value, n is the number of data points, m is the number of parameters, and F 0.05 (m, n − m) is the value of the F distribution at the α = 0.05 level. The calculated T upper and T lower values satisfying Eq. (39) were used to determine the uncertainty intervals given in Table IV .

Following the calculation of the T c uncertainty, we performed linear regressions of Eq. (31) at each of T lower , T upper , T c to determine the upper and lower 95% confidence limits for the parameter ρ c in each case. The ρ c uncertainty intervals in Table IV indicate upper and lower extremes of these limits.

We have presented a combined thermodynamic-atomistic simulation framework to calculate the vapor-liquid equilibrium (VLE) properties of a molten salt at temperatures for which the vapor pressure P * (T ) is low (less than 10 bar) and shown its results for NaCl. It is based on the assumption that the vapor phase exhibits ideal-gas behavior, which is an excellent approximation for molten salts over a salt-specific low-temperature range.

The framework provides methodology to calculate the temperature dependencies of vapor pressure, coexisting liquid and vapor-phase densities, vaporization enthalpy, and vapor-phase species concentrations (including those of very small concentration). We also show how to use the liquid and vapor coexistence densities to calculate approximations to the entire VLE coexistence curve and the critical properties. We provide uncertainty estimates for the vapor pressure, vaporization enthalpy, and critical properties.

The framework requires knowledge of the temperature dependence of ideal-gas standard Gibbs energy changes for the vapor-phase reactions (in the case of NaCl, these are the monomer ionization and dimerization reactions), and data for the liquid densities and residual chemical potentials. If the ideal-gas data is not available, it can be calculated using electronic structure software such as Gaussian [30] , and the liquid chemical potential and density data can be calculated by classical force-field-based simulation methodology.

We validated the framework for NaCl by using experimentally based data from standard sources.

We found that the assumption of only monomer and dimer species in the vapor phase is sufficient to provide results in excellent agreement with the experimental vapor phase measurements up to 2100 K. alculations are given for additional vapor-phase species present in very small concentrations; for example, at 1200 K and 2100 K, the Na + and Cl − ion mole fractions are respectively the order of 10 −9 and 10 −6 . We also use the calculated VLE T and coexisting phase densities to determine the VLE coexistence curve and the molten salt's critical properties, using the scaling law methodology of Rodrigues and Fernandes [18] . We found that the use of only low-temperature experimentally based VLE data where the vapor phase exhibits ideal-gas behavior gives slightly low values of the critical properties.

We then applied the framework to NaCl simulation data from the literature [20] for liquid chemical potential and density data determined from the Joung-Cheatham force field [21] and from the Madrid force field of Benavides et al. [22] . We found that these force fields are inadequate to provide accurate predictions of the NaCl VLE properties, although the former provides good estimates of the vaporization enthalpy. This means that molten salt force fields whose parameters are tuned to ambient conditions (such as those considered here) may not be suitable for accurate predictions of high-temperature VLE data.

The most important aspect of the framework is its thermodynamically rigorous incorporation of the vapor-phase species and their chemical reactions. For NaCl, the most important of these are the monomer ionization and dimerization reactions. The latter reaction was not considered in previous VLE simulation studies [14, 18, 20] .

Provided that the identities of the vapor-phase species are known or can be determined, the framework of this paper may be used to calculate the VLE properties of other molten salts, and of room-temperature ionic liquids.

The Supplementary Material contains uncertainty analysis derivations and data tables.

Gen3 CSP Summit

Fluid Phase Equilb

Chemical Reaction Equilibrium Analysis: Theory and Algorithms

and Compute Canada (www.computecanada.ca) are gratefully acknowledged. The authors thank Professor Athanassios Z. Panagiotopoulos for kindly providing numerical data for figures in reference [20] .