Toward accurate density and interfacial tension modeling for carbon dioxide/water mixtures

Phase behavior of carbon dioxide/water binary mixtures plays an important role in various CO2-based industry processes. This work aims to screen a thermodynamic model out of a number of promising candidate models to capture the vapor–liquid equilibria, liquid–liquid equilibria, and phase densities of CO2/H2O mixtures. A comprehensive analysis reveals that Peng–Robinson equation of state (PR EOS) (Peng and Robinson 1976), Twu α function (Twu et al. 1991), Huron–Vidal mixing rule (Huron and Vidal 1979), and Abudour et al. (2013) volume translation model (Abudour et al. 2013) is the best model among the ones examined; it yields average absolute percentage errors of 5.49% and 2.90% in reproducing the experimental phase composition data and density data collected in the literature. After achieving the reliable modeling of phase compositions and densities, a new IFT correlation based on the aforementioned PR EOS model is proposed through a nonlinear regression of the measured IFT data collected from the literature over 278.15–477.59 K and 1.00–1200.96 bar. Although the newly proposed IFT correlation only slightly improves the prediction accuracy yielded by the refitted Chen and Yang (2019)’s correlation (Chen and Yang 2019), the proposed correlation avoids the inconsistent predictions present in Chen and Yang (2019)’s correlation and yields smooth IFT predictions.


Introduction
The interaction of CO 2 with H 2 O is frequently seen in several subterranean processes (such as CO 2 -based enhanced recovery and CO 2 storage). Phase behavior of CO 2 /H 2 O mixtures under subterranean conditions plays a great role in affecting the overall efficiency of these processes. Thus how to accurately model the phase behavior of CO 2 /H 2 O mixtures becomes drastically important. Overall, an appropriate combination of cubic equation of state (CEOS), mixing rule in CEOS, α function, volume translation, and interfacial tension (IFT) model should be determined to well capture the vapor-liquid equilibrium (VLE), liquid-liquid equilibrium (LLE), phase density, and IFT of CO 2 /H 2 O mixtures.
Due to their simplicity and good reliability, CEOSs such as SRK EOS (Soave 1972) and PR EOS (Peng and Robinson 1976) are the most widely used thermodynamic models for the phase behavior modeling of CO 2 /H 2 O binary mixtures (Aasen et al. 2017;Michelsen and Mollerup 2007). Numerous articles have addressed phase-composition modeling of CO 2 /H 2 O mixtures. Two types of methods, ϕ-ϕ (fugacity-fugacity) approach and γ-ϕ (activity-fugacity) approach (Trusler 2017;Zhao and Lvov 2016), are often applied in such modeling processes. Because γ-ϕ approach has a discontinuity issue in the phase diagram near the critical region (Zhao and Lvov 2016), this work focuses on ϕ-ϕ based methods. Valtz et al. (2004) found that the most accurate model is PR EOS (Peng and Robinson 1976), Mathias-Copeman α function (Mathias and Copeman 1983), and Wong-Sandler mixing rule (Wong and Sandler 1992) with average absolute percentage deviation (AAD) of 5.4% in reproducing the measured phase composition data for CO 2 /H 2 O mixtures. However, the temperature and pressure ranges used by Valtz et al. (2004) were narrow (278.2-318.2 K and 4.64-79.63 bar, respectively). In addition, the parameters in the Wong-Sandler mixing rule 1 3 (Wong and Sandler 1992) are given as discrete values at different isotherms. Zhao and Lvov (2016) applied PRSV EOS (Stryjek and Vera 1986) and the Wong-Sandler mixing rule to calculate phase compositions, obtaining an AAD of 7.12% in reproducing the measured phase composition data of CO 2 /H 2 O mixtures over a wide range of temperatures and pressures. Similar to Valtz et al. (2004)'s study, in the study by Zhao and Lvov (2016), the parameters in the Wong-Sandler mixing rule are provided as discrete values at different isotherms, instead of generalized correlations; their model is inconvenient to use since one has to make extrapolations based on the provided values when making predictions at conditions different from those given by Zhao and Lvov (2016). Abudour et al. (2012a) applied van der Waals (1873) one-fluid (vdW) mixing rule with several temperature-dependent BIP correlations in PR EOS in determining phase compositions of CO 2 /H 2 O mixtures. With the tuned BIPs, their model yielded good accuracy (i.e., AAD of 5.0%) in aqueous phase-composition predictions but lower accuracy (i.e., AAD of 13.0%) in CO 2 -rich phase-composition predictions, respectively.
A recent comprehensive study by Aasen et al. (2017) revealed that the most accurate thermodynamic model (among the ones examined by them) in phase-composition and phase-density predictions for CO 2 /H 2 O mixtures is PR EOS, Twu α function (Twu et al. 1991), Huron-Vidal mixing rule, and constant volume translation. This model only yields an AAD of 4.5% in phase-composition calculations and an AAD of 2.8% in phase-density calculations for CO 2 / H 2 O mixtures. Aasen et al. (2017), Valtz et al. (2004), and Zhao and Lvov (2016) also pointed out that more advanced models [e.g., the Cubic-Plus-Association (CPA) EOS (Kontogeorgis et al. 1996)] do not guarantee an improvement in the phase-composition predictions for CO 2 /H 2 O mixtures.
With regards to phase-density calculations, CEOS based methods tend to overestimate liquid-phase molar volumes. A detailed discussion of this issue can be found in the studies by Matheis et al. (2016) and Young et al. (2017). In order to address this problem, Martin (1979) introduced the volume translation concept in CEOS to improve liquid-phase volumetric predictions. Peneloux et al. (1982) developed volume translation schemes in SRK EOS for pure substances. Jhaveri and Youngren (1988) applied volume translation into PR EOS, leading to the improvement of liquid phasedensity predictions. A thorough comparison of different types of volume translation methods can be found in Young et al. (2017)'s work. According to the study by Young et al. (2017), the temperature-dependent volume translation method developed by Abudour et al. (2012bAbudour et al. ( , 2013 provides the most accurate estimates on liquid-phase densities without thermodynamic inconsistencies (e.g., crossover of pressure-volume isotherms). Aasen et al. (2017) applied constant volume translation to phase-density calculations of CO 2 /H 2 O mixtures and achieved a significant improvement in density-prediction accuracies. However, a more accurate volume translation function, the one proposed by Abudour et al. (2012bAbudour et al. ( , 2013, was not applied in Aasen et al. (2017)'s study; furthermore, it should be noted that Aasen et al. (2017) used GERG-2008(Kunz and Wagner 2012 and EOS-CG (Gernert and Span 2016) calculated densities as reference densities instead of experimental data. In this study, we apply the volume translation method by Abudour et al. (2012bAbudour et al. ( , 2013 to see if the use of this model can further improve phase-density predictions for CO 2 /H 2 O mixtures; these predictions are compared to the measured density data documented in the literature. Parachor model (Sugden 1930) is one of the most widely applied models in predicting mixtures' IFT (Schechter and Guo 1998). However, its accuracy heavily relies on the density difference between the two coexisting phases in a VLE or an LLE. Our experience in using Parachor model to calculate IFT of CO 2 /H 2 O mixtures shows that Parachor model is generally appropriate for the IFT estimation for VLE of CO 2 /H 2 O systems, but less suitable for the IFT estimation for LLE of CO 2 /H 2 O systems. This is primarily because an LLE of a CO 2 /H 2 O mixture has a smaller density difference than a VLE. Several empirical IFT correlations for CO 2 / H 2 O mixtures have been proposed in the literature. However, most of these correlations are only applicable to a limited temperature and pressure range (Zhang et al. 2016). Hebach et al. (2002) proposed a new correlation which correlated IFT with phase densities. Hebach et al. (2002)'s model is suitable over a wide range of temperature and pressure conditions, although the prediction accuracy decreases with an increase in temperature or pressure. Chen and Yang (2019) proposed a new empirical IFT correlation for CO 2 /CH 4 /H 2 O ternary systems based on mutual solubility, and this model performs well for CO 2 /H 2 O binary mixtures. However, our experience in applying Chen and Yang (2019)'s model shows that some breaking points can be observed in the predicted IFT curves under some conditions, hampering its ability in providing consistent and smooth IFT predictions. In addition, using two sets of BIPs (as applied in Chen and Yang (2019)'s study) in the aqueous phase and non-aqueous phase can lead to thermodynamic inconsistency issue near the critical region as demonstrated by Li and Li (2019).
The discussion above reveals that the previous studies of phase behavior modeling of the CO 2 /H 2 O mixtures tend to primarily focus on phase-composition modeling and pay less attention to phase-density calculations (especially for the CO 2 -rich phase). Whereas, phase-density is one important property in VLE and LLE since IFT calculations and flow simulations can heavily rely on such property. As for the IFT modeling, we are currently lacking a reliable IFT correlation that not only pays due tribute to the phase composition and density of CO 2 /H 2 O mixtures but also gives smooth and 1 3 consistent IFT predictions over a wider range of temperature/pressure conditions.
In this study, we first conduct a thorough literature review to select the most promising thermodynamic models that can well capture the VLE and LLE of CO 2 /H 2 O mixtures. Then, we conduct phase-composition calculations by using PR EOS, Twu α function, and Huron-Vidal mixing rule [as suggested by Aasen et al. (2017)], and validate the accuracy of this model by comparing the calculated VLE and LLE phase compositions to the measured ones. Then, we introduce Abudour et al. (2012bAbudour et al. ( , 2013

PR EOS and α functions
In this study, PR EOS (Peng and Robinson 1976) is implemented because of its simplicity and more accurate liquiddensity predictions compared with SRK EOS (Aasen et al. 2017). The expression of PR EOS is detailed in "Appendix A".
With regards to α functions in PR EOS, Twu α function and Gasem α function are used in this study. Compared with the other types of α functions, the Twu α function can describe the thermodynamic properties of polar compounds more accurately and perform better in the supercritical region (Young et al. 2016). In addition, according to the study by Aasen et al. (2017), Twu α function coupled with PR EOS and Huron-Vidal mixing rule yields the most accurate phase-composition estimations on CO 2 /H 2 O mixtures among the models evaluated by them. Therefore, we select Twu α function as one of the evaluated α functions in this study.
The Gasem α function provides more accurate representation of supercritical phase behavior (Gasem et al. 2001). Besides, based on the study by Abudour et al. (2013), Gasem α function coupled with PR EOS, vdW mixing rule and Abudour volume translation yields the most accurate liquid-phase-density predictions for the chemical compounds examined by their study. Therefore, we select the Twu α function and the Gasem α function in VLE/LLE and phasedensity calculations. "Appendix A" shows the expressions of Twu α function and Gasem α function.

Mixing rules
Mixing rules have a great impact on phase equilibrium calculations. Huron and Vidal (1979) proposed a new expression by considering the excess Gibbs energy for CEOS, which made more accurate the phase-composition predictions for mixtures containing polar substances. Furthermore, according to the comprehensive study by Aasen et al. (2017), the most accurate thermodynamic model among the ones examined by them is PR EOS coupled with Twu α function and Huron-Vidal mixing rule, which provides an AAD of 4.5% in reproducing the phase-composition data measured for CO 2 /H 2 O mixtures. Hence, in the first part of this study, we collect more phase equilibria data for CO 2 /H 2 O mixtures to verify the performance of the model suggested by Aasen et al. (2017). These additional experimental data are not included in the study by Aasen et al. (2017).
The vdW mixing rule is one of the most commonly used mixing rules in petroleum industry (Pedersen et al. 2014). Although vdW mixing rule is originally developed for nonpolar systems, the vdW mixing rule coupled with the tuned BIPs can be reliably used for describing the phase behavior of mixtures containing polar components (e.g., water). Besides, based on the study by Abudour et al. (2012a), Gasem α function with vdW mixing rule and their temperature-dependent volume translation function provided a promising means to well reproduce the measured liquidphase densities for CO 2 /H 2 O mixtures. Therefore, in this study, we also employ the model suggested by Abudour et al. (2012a) to test if it outperforms the model suggested by Aasen et al. (2017). The expressions of vdW mixing rule and Huron-Vidal mixing rule and their BIPs are detailed in "Appendices B and C".

Volume translation models
Volume translation is used to overcome the inherent deficiency of CEOS in liquid-phase-density predictions. In order to improve liquid-phase-density calculations, Peneloux et al. (1982) developed a constant volume translation model in SRK EOS, while Jhaveri and Youngren (1988) developed a constant volume translation model in PR EOS. Abudour et al. (2012bAbudour et al. ( , 2013 revised the temperature-dependent volume translation function to improve both saturated and single-phase liquid density calculations. Furthermore, unlike other temperature-dependent volume translation models, the volume translation model developed by Abudour et al. (2012bAbudour et al. ( , 2013 does not yield thermodynamic inconsistency issues unless at extremely high pressures. Therefore, we select the constant and Abudour et al. volume translation models in this study for phase-density predictions. The expressions of these two volume translation models are detailed in "Appendix D".

IFT correlations for CO 2 /H 2 O mixtures
In this study, we select Parachor model (Sugden 1930), Hebach et al. (2002)'s correlation, and Chen and Yang (2019)'s correlation to predict IFT of CO 2 /H 2 O mixtures. The Parachor model (Sugden 1930) is one of the most widely used methods in determining mixtures' IFT. It correlates IFT with phase compositions and molar densities of each phase. Parachor is a component-dependent constant. The expression of the Parachor model is shown in "Appendix E". Hebach et al. (2002)'s correlation correlates IFT with temperature, pressure, and phase densities. Phase compositions are not included in their correlation. To make fair comparison, we also refit coefficients in their correlation based on the IFT database employed in this study. Values of the original and refitted coefficients as well as the Hebach et al. (2002)'s correlation are shown in "Appendix E".
Chen and Yang (2019)'s correlation correlates IFT with phase equilibrium ratios (K-values) and the reduced pressure of CO 2 . Unlike the Parachor model and the Hebach et al. (2002)'s correlation, the density of the two equilibrating phases is not one input in the Chen and Yang (2019)'s correlation. Chen and Yang (2019)'s correlation, they proposed four groups of coefficient sets, i.e., one coefficient set (using one coefficient set on the whole pressure range) with or without the reduced pressure term, and two coefficient sets (dedicated to the pressure ranges of p ≤ 73.8 and p > 73.8 bar) with or without the reduced pressure term. Since using the reduced pressure term can improve prediction accuracy (Chen and Yang 2019), we introduce the reduced pressure term in this study. Similarly, we refit these coefficients based on the IFT databased employed in this study to make fair comparison. Values of the original and refitted coefficients as well as the Chen and Yang (2019)'s correlation are detailed in "Appendix E".

IFT correlation proposed in this study
Before we finalize our IFT correlation, we tried several scenarios to find the optimal one to correlate the IFT of CO 2 / H 2 O mixtures. Since the Parachor model is one of the most widely used models in mixtures' IFT predictions, we revise the original Parachor model by introducing a componentdependent correction term α i ; furthermore, we replace the constant exponential term in the original Parachor model by correlating it with several physical properties (e.g., equilibrium ratios). The new IFT correlation can be expressed as follows: where σ in the interfacial tension; α i is the introduced component-dependent correction term; x i and y i are the mole fractions of component i in liquid and vapor phases, respectively; P i is the Parachor value of component i ( P H 2 O = 52 , P CO 2 = 78 ) (Liu et al. 2016); M L is the molar density of liquid phase in mol/cm 3 ; M V is the molar density of vapor phase in mol/cm 3 . N is the number of component; n is the exponent.
First, the component-dependent correction term α i is set as a constant for each component, and the exponential term n can be expressed by the equilibrium ratios of CO 2 -rich phase and aqueous phase: where C 1 , C 2 , and C 3 are empirical coefficients; K CO 2 and K H 2 O are the equilibrium ratios (as known as K-values) of CO 2 and H 2 O: Since using one coefficient set for both α i and n on the whole CO 2 -rich-phse density range cannot converge after reaching the maximum iterations, we use two coefficient sets based on CO 2 -rich-phase densities. Table 1 listed the values of these coefficients and α i determined by fitting the proposed correlation (abbreviated as Scenario #1) to the IFT training dataset.
Since using constants to represent α i leads to a larger AAD compared with the refitted Chen and Yang (2019)'s correlation (i.e., 8.8746% vs. 7.8520%, respectively), we correlate equilibrium ratios to α i to see if it can improve the IFT predictions. The expression of n in this scenario (abbreviated as Scenario #2) is the same as that in Scenario #1. The expression for α i is given as: Specifically, when the CO 2 -rich-phase density is greater than 0.2 g/cm 3 , H 2 O can be simplified as:  Table 2 listed the values of these coefficients determined by fitting the proposed correlation to the IFT training dataset. Since using one coefficient set for α i and n on the whole CO 2 -rich-phse density range in Scenario #2 cannot converge after reaching the maximum iterations, we use two coefficients based on CO 2 -rich-phase densities (the same as Scenario #1).
We find that using correlations to represent α i can slightly improve the IFT predictions (i.e., AAD of 8.31% in Scenario #2 vs. AAD of 8.87% in Scenario #1). Besides, we find that the value of n is around 4 over a wide range of temperature/pressure conditions in all scenarios (i.e., its value only slightly changes with the change of equilibrium ratios); therefore, we set the value of n as 4 for simplicity.
However, our experience in applying Scenarios #1 and #2 shows that some breaking points can be observed in the correlated IFT curves due to the fact that two different sets of coefficients are adopted under the conditions of CO 2 -rich < 0.2 g/cm 3 and H 2 O -rich ≥ 0.2 g/cm 3 , respectively. In addition, based on the study by Chen and Yang (2019), introducing the reduced pressure of CO 2 can improve IFT predictions. Thus, we introduce the reduced pressure of CO 2 in the expressions of α i and use one coefficient set to see if these settings can further improve the prediction accuracies without yielding inconsistent IFT predictions. Based on the calculation results, the following IFT correlation yields the lowest AAD among the ones examined in this study: where the α i term in the new correlation can be expressed as: where p r is the reduced pressure of CO 2 . Table 3 lists the values of these coefficients determined by fitting the proposed correlation to the IFT training dataset.

Results and discussion
The values of critical pressure (p c ), critical temperature (T c ), acentric factor (ω), molecular weight (M), critical compressibility factor (Z c ) used in this study are retrieved from the NIST database (Lemmon et al. 2011).  Aasen et al. (2017). Comparison between the measured and calculated phase-composition results is evaluated by the average absolute percentage deviation (AAD) defined as:

Performance comparison of thermodynamic models in phase equilibrium calculations
where AAD is the average absolute percentage deviation; N is the number of data points; x CAL and x EXP are the calculated and measured mole fraction of CO 2 or H 2 O in the aqueous phase (or the CO 2 -rich phase), respectively. Table 5 details the settings of the four thermodynamic models examined in this work. Table 6 summarizes the performance of different thermodynamic models in phasecomposition predictions.
As shown in Table 6, although AAD for x CO 2 of Case 3 (Twu + HV) is slightly higher than that of Case 4 (Gasem + HV), i.e., 4.73% of Case 3 vs. 3.64% of Case 4, Case 3 (Twu + HV) significantly outperforms the other models in y H 2 O predictions, i.e., AAD of 10.37% of Case 3 vs. > 16% of other cases. Thus, given the overall performance, Case 3 (Twu + HV) is found to be the best model Table 2 Values of the correlation coefficients and α i in Scenario #2  in phase-composition predictions. Figure 1 compares the performance of different models at T = 323.15 K and T = 348.15 K. As can be seen from these two figures, the thermodynamic model Case 3 (Twu + HV) can well capture the trend exhibited by the measured solubility data over a wide pressure range. Since Case 3 (Twu + HV) outperforms other thermodynamic models in phase-composition predictions for CO 2 / H 2 O mixtures, we only focus on the performance of Case 3 coupled with volume translation in phase-density predictions. Table 8 summarizes the performance of different volume translation models in both aqueous-phase and CO 2 -rich-phase density calculations.  Spycher et al. (2003). We directly use these data mentioned in their paper for convenience i N is 2 for x CO 2 and 3 for y H 2 O , respectively j Only experimental data for CO 2 /pure water are selected in the study by Zhao et al. (2015) T, K p, bar    As shown in Table 8, incorporation of VT into the thermodynamic framework can generally improve the phasedensity prediction accuracy. Case 3-1 (Twu + HV + Abudour VT) provides the most accurate estimates of both aqueousphase and CO 2 -rich-phase density, yielding AAD of 2.90% in reproducing the measured phase-density data. Figure 2 further visualizes some of the calculation results by these three different models at different pressure/temperature conditions. It can be seen from Fig. 2 that, regarding aqueousphase density predictions, the performance of Case 3-2 (Twu + HV + Constant VT) improves dramatically as temperature rises. As shown in Fig. 2e, f, at high temperature conditions, Cases 3-2 yields the most accurate aqueousphase density predictions; however, it fails to accurately predict CO 2 -rich-phase densities. As a lighter phase, CO 2 -rich-phase density can be accurately predicted without the use of volume translation functions. Applying Abudour VT method is able to only slightly improve the prediction  Briones et al. (1987) and Gillepsie and Wilson (1982)  1 3 accuracy (i.e., AAD of 2.62%). In contrast, applying constant VT in CO 2 -rich-phase density predictions can lead to larger errors than the case without the use of VT. Figure 3 compares the performance of different models in terms of their accuracy in phase-density predictions over 382.14-478.35 K and 35.3-1291.9 bar. Note that the results of CPA EOS model from the work by Tabasinejad et al. (2010) focuses on the same pressure and temperature ranges. As can be seen from Fig. 3, although the CPA EOS model can accurately predict the aqueous-phase density, it tends to be less accurate in determining the CO 2 -rich-phase density. Overall, the thermodynamic model Cases 3-1 (Twu + HV + Abudour VT) give an accuracy comparable to the more complex CPA EOS model.

Evaluation of thermodynamic models in density calculations
In addition, according to the study by Aasen et al. (2017), CPA EOS model yields higher percentage errors (AAD) in reproducing phase-composition data for CO 2 /H 2 O mixtures compared with Case 3 (PR EOS + Twu + HV), i.e., 9.5% vs. 4.5% (Aasen et al. 2017). Therefore, overall, Case 3-1 (Twu + HV + Abudour VT) is a more accurate model in both phase-composition and phase-density predictions for CO 2 / H 2 O mixtures. In order to expand our IFT database, IFT data with precisely determined phase densities are also included in our IFT database. The collected IFT data are randomly placed into two bins: a training dataset (including 589 data points) and a validation dataset (including 189 data points). Results in Sects. 3.1 and 3.2 reveal that the thermodynamic model using PR EOS, Twu α function, Huron-Vidal mixing rule, and Abudour et al. (2013) VT yields the most accurate estimates on both phase compositions and densities. Therefore, the aforementioned thermodynamic model provides reliable phase-composition and phase-density predictions that can be fed into the proposed IFT correlation.

Evaluation of the newly proposed IFT correlation
Mean absolute errors (MAE), AAD, and coefficient of determination (R 2 ) are used as performance measures. The expressions of MAE and R 2 are as follows: where σ EXP is the measured IFT data in mN/m; σ CAL is the calculated IFT in mN/m by different correlations; EXP is the average of the measured IFTs in mN/m. Table 10 shows the details of the different IFT models examined in this study. Table 11 summarize the performance of different correlations in IFT estimations. As can be seen, the most accurate IFT model is Model 3 proposed in this study, although it only shows a marginal edge over Model 2. Figure 4 visually compares the measured IFTs vs. pressure and the calculated ones by different IFT models at selected temperatures. As shown in these plots, in general, Model 3 (this study) outperforms other empirical correlations over a wide range of temperatures and pressures. It can be also observed from these plots that breaking points appear in the predicted IFT curves at p = 73.8 bar by Model 2 (Refitted Chen and Yang (2019)'s correlation with two coefficient sets). Such discontinuous IFT prediction can be attributed to the fact that two different sets of coefficients are adopted under the conditions of p ≤ 73.8 and p > 73.8 bar, respectively, in Chen and Yang (2019)'s correlation. Although using one coefficient in Chen and Yang (2019)'s correlation (e.g., Model 5) can avoid such discontinuous IFT predictions, it yields larger percentage errors. Therefore, Model 3 (this study) is the best model in IFT predictions for CO 2 /H 2 O mixtures over a wide range of temperatures and pressures. Figure 5 illustrates how the IFTs predicted by Model 3 (this study) vary with pressure at different temperatures. It can be observed from Fig. 5 that the new IFT correlation Aqueous-phase/CO 2 -rich-phase density, g/cm 3 Aqueous-phase/CO 2 -rich-phase density, g/cm 3 Aqueous-phase/CO 2 -rich-phase density, g/cm 3 Aqueous-phase/CO 2 -rich-phase density, g/cm 3

Performance of different IFT correlations
Aqueous-phase density CO2-rich-phase density Aqueous-phase density CO2-rich-phase density Aqueous-phase density CO2-rich-phase density Aqueous-phase density CO2-rich-phase density Aqueous-phase density CO2-rich-phase density Fig. 2 Predictions of aqueous-phase and CO 2 -rich-phase density by Case 3-1 (Twu + HV + Abudour VT, dashed line), Case 3-2 (Twu + HV + constant VT, dotted line) and Case 3 (Base case, solid line) at different temperature conditions. The circles are the measured phasedensity data from the study by Efika et al. (2016) provides smooth and consistent IFT predictions at different pressures and temperatures. Overall, Model 3 proposed in this study yields accurate and consistent IFT predictions over the wide range of temperatures and pressures, although it yields relatively higher percentage errors at higher temperature conditions (e.g., T = 478 K) compared with that at lower temperature conditions (i.e., T < 378 K). It is interesting to observe from Fig. 5a that when the pressure is less than around 15 bar and the temperature is between 278.15 and 368.15 K, an increase in temperature leads to a decrease in the predicted IFT under the same pressure. In comparison, when the pressure is larger than around 15 bar, an increase in temperature leads to an increase in the predicted IFT. At higher temperatures of 378.15-478.15 K, an increase in temperature always results in a decline in the predicted IFT under the same pressure, as seen in Fig. 5b. Most of the measured IFTs documented in the literature follow this trend (Akutsu et al. 2007;Chalbaud et al. 2009;Chiquet et al. 2007;Chun and Wilkinson 1995;Da Rocha et al. 1999;Georgiadis et al. 2010;Hebach et al. 2002;Heuer 1957;Hough et al. 1959;Kvamme et al. 2007;Khosharay and Varaminian 2014;Liu et al. 2016;Park et al. 2005;Pereira et al. 2016;Shariat et al. 2012), except for the studies by Bachu and Bennion (2009) and Bikkina et al. (2011), i.e., an increase in temperature leads to an increase in IFT at a temperature range of 373.15-398.15 K in the study by Bachu and Bennion (2009), and an increase in temperature leads to an increase in IFT over 298.15-333.15 K in the study by Bikkina et al. (2011). Again, the sharp drops in the IFT curves at lower temperatures (where CO 2 remains subcritical) are due to the transformation of VLE to LLE.  b,c,e These data are already summarized by Park et al. (2005) and Shariat et al. (2012). We directly use these data mentioned in their papers for convenience d,f,g Some experimental data appear to be outliers and hence excluded for further analysis due to the significant deviation from other experimental data at similar temperature and pressure conditions (see Appendix G)

Statistical significance tests of IFT correlations
As shown in Table 11, the AADs yielded by Model 2 (refitted Chen and Yang (2019)'s correlation) and Model 3 (this study) are on the same scale. Therefore, it is necessary to conduct statistical significance tests to check if the marginal edge of Model 3 over Model 2 is statistically significant. Figure 1 shows the frequency distribution of the differences between the measured IFT data (i.e., whole dataset including 778 data points) and calculated ones by Model 2, while Fig. 7 shows the same information for Model 3. As can be seen from Figs. 6 and 7, the distribution of the deviations generated by the two models can be considered to follow Gaussian distributions. As such, paired one-tailed t-tests are applied as the statistical significance test method (Japkowicz and Shah 2011). P-value is used to check if one model is better than another one. Typically, the significance threshold α is 0.05; when P > , two models have the same performance. In contrast, when P ≤ , it is reasonable to say that one model is significantly better than another one (Japkowicz and Shah 2011).
P-value of Model 2 (refitted original Chen and Yang (2019)'s correlation with two coefficient sets) and Model 3 (this study) is 0.0069, P < ( = 0.05 ); therefore, it is reasonable to say that Model 3 statistically outperforms Model 2. In addition, the new model does not give discontinuous IFT predictions, while Chen and Yang (2019)'s IFT model bears such issue.

Conclusions
The objective of this study is to screen and develop reliable models for describing the VLE, LLE, phase density, and IFT of CO 2 /H 2 O mixtures. Based on the comparison between the experimental data and the calculated ones from different models, we can reach the following conclusions: 1. The most accurate method to represent CO 2 /H 2 O VLE and LLE is PR EOS, Twu α function, and Huron-Vidal mixing rule, which only yields AAD of 5.49% and 2.90% in reproducing measured CO 2 /H 2 O phase-composition data and phase-density data over a temperature range of 278-378.15 and 278-478.35 K and over a pressure range of 6.9-709.3 and 2.5-1291.1 bar, respectively. 2. Applying either constant or Abudour et al. (2013) VT method can significantly improve aqueous-phase density calculations. In addition, when the temperature is higher than 373 K, constant VT method can yield lower error in reproducing measured phase-density data than Abudour et al. (2013) VT method; 3. Constant VT method cannot improve the prediction accuracy of CO 2 -rich-phase density. Abudour et al. (2013) VT method can slightly improve CO 2 -rich-phase density predictions, but such improvement is more obvious at low to moderate temperature conditions. 4. The new IFT correlation based on the aforementioned PR EOS model outperforms other empirical correlations with an overall AAD of 7.77% in reproducing measured IFT data of CO 2 /H 2 O mixtures. The new IFT correla-  Model 1 (Parachor model) shows a more deteriorating performance when the vapor CO 2 -rich phase changes to a liquid phase. Experimental data are from the studies by Kvamme et al. (2007), Georgiadis et al. (2010), Liu et al. (2016), and Shariat et al. (2012) tion is only slightly more accurate than the refitted Chen and Yang (2019)  where

Appendix C: Summary of Huron-Vidal mixing rule and its BIPs
In the Huron-Vidal mixing rule, the following equations are applied to calculate a m and b m (Huron and Vidal 1979): where G E ∞ is the excess Gibbs energy at infinite pressure, and Λ is an EOS-dependent parameter. For PR EOS, Λ = 0.62323 (Huron and Vidal 1979).
The excess Gibbs energy corresponding to the non-random two-liquid (NRTL) (Zhao and Lvov 2016;Wong and Sandler 1992) model can be expressed by (Aasen et al. 2017;Huron and Vidal 1979): where The generalized BIP correlations for τ ij obtained by Aasen et al. (2017) are given below: where T 0 = 1000 K is the reference temperature.
When the Huron-Vidal mixing rule is used in PR EOS, the fugacity coefficient can be calculated by (Zhao and Lvov 2016): where lnγ i is the activity coefficient of component i and can be expressed as (Zhao and Lvov 2016): The derivation of the expression of the activity coefficient in Huron-Vidal mixing rule is detailed in "Appendix H".
where (Abudour et al. 2013): where T cm , p cm and δ cm are the critical temperature, critical pressure and volume correction of the mixture at the critical point, respectively. c 1i is the specie-specific parameter of component i and has a linear relationship with critical compressibility ( Z c ) (Abudour et al. 2012b): The term d m can be derived using the original PR EOS (Matheis et al. 2016): The volume correction of the given mixture at the critical point, δ cm , can be determined by (Abudour et al. 2013): where v ci is the critical volume of component i; θ i is the surface fraction of component i defined by (Abudour et al. 2013): The critical temperature of the mixture can be calculated via the following mixing rule (Abudour et al. 2013).
The critical pressure of the mixture can be determined by the correlation proposed by Aalto et al. (1996): where ω m is the acentric factor of the mixture (Abudour et al. 2013): where ω i is the acentric factor of component i.

Appendix E: Summary of existing IFT correlations for CO 2 /H 2 O mixtures
Parachor model (Sugden 1930) can be expressed as below (Schechter and Guo 1998): where x i and y i are the mole fractions of component i in liquid and vapor phases, respectively; P i is the Parachor value of component i ( P H 2 O = 52 , P CO 2 = 78 ) (Liu et al. 2016); M L is the molar density of liquid phase in mol/cm 3 ; M V is the molar density of vapor phase in mol/cm 3 . Hebach et al. (2002)'s correlation can be expressed as: where (Hebach et al. 2002):  where CO 2 is the CO 2 -rich-phase density in g/cm 3 ; H 2 O is the aqueous-phase density in g/cm 3 ; k 0 to k 6 and b 0 to b 1 are empirical coefficients. The units of T, p, and dd are K, bar, and g 2 /cm 6 , respectively. Table 13 lists the values of original and refitted coefficients. Chen and Yang (2019)'s correlation is given as: where σ is IFT in mN/m; p r is the reduced pressure of CO 2 ; C 1 to C 5 are empirical coefficients. To make fair comparisons, we refit these coefficients based on the IFT database employed in this study. Table 14 summarizes the values of these refitted coefficients.
Appendix F: Pressure-temperature coverage of phase-density and IFT data collected from the literature.

Appendix G: Experimental data selection in IFT database employed in this study
The collected IFT data are further screened to remove any obvious outliers. Figure 9 shows the identification of the outliers from the collected data over 40-60 bar and 278.15-298.15 K, while Fig. 10 shows the identification of outliers from the collected data over 100-270 bar and 307.15-314.15 K. As seen in Fig. 9a, the measured IFT data by Chun and Wilkinson (1995) and Park et al. (2005) fall into the range of 5-8 mN/m over 278.15-288.15 K and 40-60 bar, which are significantly lower than the measured IFT data (i.e., around 22-28 mN/m) obtained by other studies under similar conditions. Figure 9b indicates that the measured IFT data by Chun and Wilkinson (1995) and Park et al. (2005) fall into the range of 10-14 mN/m over 293.15-298.15 K and 50-70 bar, which are significantly lower than the measured (53) = C 1 + C 2 p r + C 3 lnK CO 2 + C 4 p r + C 5 lnK H 2 O IFT data (i.e., around 28-31 mN/m) obtained by other studies under similar conditions. Figure 10 indicates that the measured data by Bachu and Bennion (2009) fall into the range of 16-19 mN/m over 307.15-314.15 K and 120-270 bar, which are significantly lower than the measured IFT values (i.e., around 30 mN/m) obtained by other studies under similar conditions. No outlier exists at other temperature and pressure conditions. We have removed these outliers in the IFT regression analysis.

Appendix H: Derivation of activity coefficient in the fugacity expression when Huron-Vidal mixing rule is used.
Similar to the approach used by Wong and Sandler (1992), the activity coefficient of component i can be expressed by the following formula: where the excess Gibbs free energy can be expressed as (Huron and Vidal 1979): To make the derivation process more intuitive, we can set i = 1 (the first component) and n = 2 (two compounds in the system). Then the excess Gibbs free energy can be expressed as: Taking the partial derivative of the first term in the right hand side of Eq. (56) yields: Taking the partial derivative of the second term in the right hand side of Eq. (56) yields: As such, Eq. (54) can be expressed as: Using letter i to replace number 1 leads to a general expression of activity coefficient (Huron and Vidal 1979):