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

Abstract

Phase behavior of carbon dioxide/water binary mixtures plays an important role in various CO₂-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 CO₂/H₂O 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.

1 Introduction

The interaction of CO₂ with H₂O is frequently seen in several subterranean processes (such as CO₂-based enhanced recovery and CO₂ storage). Phase behavior of CO₂/H₂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₂/H₂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₂/H₂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₂/H₂O binary mixtures (Aasen et al. 2017; Michelsen and Mollerup 2007). Numerous articles have addressed phase-composition modeling of CO₂/H₂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₂/H₂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 (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₂/H₂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₂/H₂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₂-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₂/H₂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₂/H₂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₂/H₂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 phase-density 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. (2012b, 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₂/H₂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. (2012b, 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. (2012b, 2013) to see if the use of this model can further improve phase-density predictions for CO₂/H₂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₂/H₂O mixtures shows that Parachor model is generally appropriate for the IFT estimation for VLE of CO₂/H₂O systems, but less suitable for the IFT estimation for LLE of CO₂/H₂O systems. This is primarily because an LLE of a CO₂/H₂O mixture has a smaller density difference than a VLE. Several empirical IFT correlations for CO₂/H₂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₂/CH₄/H₂O ternary systems based on mutual solubility, and this model performs well for CO₂/H₂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₂/H₂O mixtures tend to primarily focus on phase-composition modeling and pay less attention to phase-density calculations (especially for the CO₂-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₂/H₂O mixtures but also gives smooth and 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₂/H₂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. (2012b, 2013) volume translation model in phase-density calculations to check if applying this model can further improve the density-prediction accuracies. A new empirical IFT correlation for CO₂/H₂O mixtures is then proposed based on the reliable thermodynamic model that incorporates the Huron-Vidal mixing rule and the Abudour et al. (2013) volume translation model.

2 Methodology

2.1 PR EOS and α functions

In this study, PR EOS (Peng and Robinson 1976) is implemented because of its simplicity and more accurate liquid-density 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₂/H₂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 phase-density calculations. “Appendix A” shows the expressions of Twu α function and Gasem α function.

2.2 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₂/H₂O mixtures. Hence, in the first part of this study, we collect more phase equilibria data for CO₂/H₂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 liquid-phase densities for CO₂/H₂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”.

2.3 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. (2012b, 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. (2012b, 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”.

2.4 IFT correlations for CO₂/H₂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₂/H₂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₂. 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”.

2.5 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₂/H₂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 component-dependent 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:

$$\sigma = \left[ {\sum\limits_{i = 1}^{N} {\alpha_{i} } P_{i} (x_{i} \rho_{{\text{L}}}^{{\text{M}}} - y_{i} \rho_{{\text{V}}}^{{\text{M}}} )} \right]^{n}$$

(1)

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_{{{\text{H}}_{2} {\text{O}}}} = 52$, $P_{{{\text{CO}}_{2} }} = 78$) (Liu et al. 2016); $\rho_{{\text{L}}}^{{\text{M}}}$ is the molar density of liquid phase in mol/cm³; $\rho_{{\text{V}}}^{{\text{M}}}$ is the molar density of vapor phase in mol/cm³. 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₂-rich phase and aqueous phase:

$$n = C_{1} \ln K_{{{\text{CO}}_{2} }} + C_{2} \ln K_{{{\text{H}}_{2} {\text{O}}}} + C_{3}$$

(2)

where C₁, C₂, and C₃ are empirical coefficients; $K_{{{\text{CO}}_{2} }}$ and $K_{{{\text{H}}_{2} {\text{O}}}}$ are the equilibrium ratios (as known as K-values) of CO₂ and H₂O:

$$K_{i} = y_{i} /x_{i}$$

(3)

Since using one coefficient set for both α_i and n on the whole CO₂-rich-phse density range cannot converge after reaching the maximum iterations, we use two coefficient sets based on CO₂-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.

Table 1 Values of the correlation coefficients and α_i in Scenario #1

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:

$$\alpha_{i} = C_{1} \ln K_{{{\text{CO}}_{2} }} + C_{2} \ln K_{{{\text{H}}_{2} {\text{O}}}} + C_{3}$$

(4)

Specifically, when the CO₂-rich-phase density is greater than 0.2 g/cm³, $\alpha_{{{\text{H}}_{2} {\text{O}}}}$ can be simplified as:

$$\alpha_{{_{{{\text{H}}_{2} {\text{O}}}} }} = C_{1} \ln K_{{{\text{CO}}_{2} }} + C_{3}$$

(5)

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₂-rich-phse density range in Scenario #2 cannot converge after reaching the maximum iterations, we use two coefficients based on CO₂-rich-phase densities (the same as Scenario #1).

Table 2 Values of the correlation coefficients and α_i in Scenario #2

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 $\rho_{{{\text{CO}}_{2} {\text{-rich}}}} < 0.2$ g/cm³ and $\rho_{{{\text{H}}_{{2}} {\text{O }\text{-rich}}}} \ge 0.2$ g/cm³, respectively. In addition, based on the study by Chen and Yang (2019), introducing the reduced pressure of CO₂ can improve IFT predictions. Thus, we introduce the reduced pressure of CO₂ 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:

$$\sigma = \left[ {\mathop \sum \limits_{i = 1}^{N} \alpha_{i} P_{i} \left( {x_{i} \rho_{{\text{L}}}^{{\text{M}}} - y_{i} \rho_{{\text{V}}}^{{\text{M}}} } \right)} \right]^{4}$$

(6)

where the α_i term in the new correlation can be expressed as:

$$\alpha = C_{1} + (C_{2} p_{{\text{r}}} + C_{3} )\ln K_{{{\text{CO}}_{2} }} + (C_{4} p_{{\text{r}}} + C_{5} )\ln K_{{{\text{H}}_{2} {\text{O}}}}$$

(7)

where p_r is the reduced pressure of CO₂.

Table 3 lists the values of these coefficients determined by fitting the proposed correlation to the IFT training dataset.

Table 3 Coefficients in the α_i term for H₂O and CO₂

3 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).

3.1 Performance comparison of thermodynamic models in phase equilibrium calculations

Table 4 summarizes the measured phase equilibrium data of CO₂/H₂O mixtures over 278–378.15 K and 0.92–709.3 bar reported in the literature. Note that these experimental data were not included in the study by Aasen et al. (2017). Comparison between the measured and calculated phase-composition results is evaluated by the average absolute percentage deviation (AAD) defined as:

$${\text{AAD}} = \frac{1}{N}\sum\limits_{i}^{N} {\left| {\frac{{x_{{{\text{CAL}}}} - x_{{{\text{EXP}}}} }}{{x_{{{\text{EXP}}}} }}} \right|}_{i} \times 100\%$$

(8)

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₂ or H₂O in the aqueous phase (or the CO₂-rich phase), respectively.

Table 4 Phase equilibrium data of CO₂/H₂O mixtures employed in this study

Table 5 details the settings of the four thermodynamic models examined in this work. Table 6 summarizes the performance of different thermodynamic models in phase-composition predictions.

Table 5 Settings of four thermodynamic models examined in this work

Table 6 AAD of calculated mole fraction of CO₂ in the aqueous phase ($x_{{{\text{CO}}_{2} }}$

As shown in Table 6, although AAD for $x_{{{\text{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_{{{\text{H}}_{{2}} {\text{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 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.

Fig. 1

Measured and calculated pressure-composition data for CO₂/H₂O mixtures at T = 323.15 K (a) and T = 348.15 K (b). Solid circles are the experimental data from the studies by Briones et al. (1987) and Gillepsie and Wilson (1982)

3.2 Evaluation of thermodynamic models in density calculations

Table 7 summarizes the experimental aqueous-phase and CO₂-rich-phase densities of CO₂/H₂O mixtures over 278–478.35 K and 2.5–1291.1 bar documented in the literature. The pressure–temperature coverage of the phase density data collected from the literature are shown in “Appendix F”.

Table 7 Aqueous-phase ($\rho_{{{\text{H}}_{{2}} {\text{O}}}}$) and CO₂-rich-phase ($\rho_{{{\text{CO}}_{2} {\text{-rich}}}}$) density data of CO₂/H₂O mixtures employed in this study

Since Case 3 (Twu + HV) outperforms other thermodynamic models in phase-composition predictions for CO₂/H₂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₂-rich-phase density calculations.

Table 8 AAD of the calculated aqueous-phase density ($\rho_{{{\text{H}}_{{2}} {\text{O}}}}$) and CO₂-rich-phase density ($\rho_{{{\text{CO}}_{2} {\text{-rich}}}}$) by different thermodynamic models

As shown in Table 8, incorporation of VT into the thermodynamic framework can generally improve the phase-density prediction accuracy. Case 3–1 (Twu + HV + Abudour VT) provides the most accurate estimates of both aqueous-phase and CO₂-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.

Fig. 2

Predictions of aqueous-phase and CO₂-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 phase-density data from the study by Efika et al. (2016)

It can be seen from Fig. 2 that, regarding aqueous-phase 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 aqueous-phase density predictions; however, it fails to accurately predict CO₂-rich-phase densities. As a lighter phase, CO₂-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 accuracy (i.e., AAD of 2.62%). In contrast, applying constant VT in CO₂-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₂-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.

Fig. 3

Bar chart plots comparing the AAD in aqueous-phase (black) and CO₂-rich-phase (gray) density predictions by different models over 382.14–478.35 K and 35.3–1291.9 bar. Calculation results by the CPA EOS method are from the study by Tabasinejad et al. (2010)

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₂/H₂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₂/H₂O mixtures.

3.3 Evaluation of the newly proposed IFT correlation

Table 9 summarizes the experimental IFT data of CO₂/H₂O mixtures over 278.15–477.59 K and 1–1200.96 bar documented in the literature. Ideally, phase densities should be directly measured; however, only Chiquet et al. (2007), Kvamme et al. (2007), Bikkina et al. (2011), Bachu and Bennion (2009), and Shariat et al. (2012) applied measured phase densities in IFT calculations. 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).

Table 9 Measured IFT data for CO₂/H₂O mixtures used in this study

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.

Mean absolute errors (MAE), AAD, and coefficient of determination (R²) are used as performance measures. The expressions of MAE and R² are as follows:

$${\text{MAE}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left| {\sigma_{{{\text{EXP}},i}} - \sigma_{{{\text{CAL}},i}} } \right|$$

(9)

$$R^{2} = 1 - \frac{{\mathop \sum \nolimits_{i = 1}^{N} \left( {\sigma_{{{\text{EXP}},i}} - \sigma_{{{\text{CAL}},i}} } \right)^{2} }}{{\mathop \sum \nolimits_{i = 1}^{N} \left( {\sigma_{{{\text{EXP}},i}} - \overline{\sigma }_{{{\text{EXP}}}} } \right)^{2} }}$$

(10)

where σ_EXP is the measured IFT data in mN/m; σ_CAL is the calculated IFT in mN/m by different correlations; $\overline{\sigma }_{{{\text{EXP}}}}$ is the average of the measured IFTs in mN/m.

3.3.1 Performance of different IFT correlations

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.

Table 10 Technical Characteristics of different IFT models examined in this study

Table 11 Summary of the performance of different correlations in IFT estimations

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₂/H₂O mixtures over a wide range of temperatures and pressures.

Fig. 4

IFT predictions at different temperature conditions by different models. At T = 297.9 K (a), VLE is transformed to LLE at p = 64 bar. Model 1 (Parachor model) shows a more deteriorating performance when the vapor CO₂-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)

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 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₂ remains subcritical) are due to the transformation of VLE to LLE.

Fig. 5

taken from previous studies by Heuer (1957), Chun and Wilkinson (1995), Park et al. (2005), Akutsu et al. (2007), Bachu and Bennion (2009), Bikkina et al. (2011), Shariat et al. (2012), and Liu et al. (2016)

Plots of predicted IFTs vs. pressure by the newly proposed IFT correlation Model 5 at the temperature ranges of 278–368 K (a) and 378–478 K (b). The curves are plotted with an interval of 10 K. Experimental data are

3.3.2 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).

Fig. 6

Frequency distribution of the differences between the measured IFT data (i.e., the whole dataset including 778 data points) and calculated ones by Model 2 (refitted Chen and Yang (2019)’s correlation with two coefficient sets). Blue columns are instances, and the red curve is probability density function which follows Gaussian distribution with μ = 0.0941 and σ² = 3.3457

Fig. 7

Frequency distribution of the difference between the measured IFT data (i.e., the whole dataset including 778 data points) and calculated ones by Model 3 (this study). Blue columns are instances, and the red curve is probability density function which follows Gaussian distribution with μ = -0.2051 and σ² = 3.2781

P-value is used to check if one model is better than another one. Typically, the significance threshold α is 0.05; when $P > \alpha$, two models have the same performance. In contrast, when $P \le \alpha$, 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 < \alpha$ ($\alpha = 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.

4 Conclusions

The objective of this study is to screen and develop reliable models for describing the VLE, LLE, phase density, and IFT of CO₂/H₂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₂/H₂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₂/H₂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₂-rich-phase density. Abudour et al. (2013) VT method can slightly improve CO₂-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₂/H₂O mixtures. The new IFT correlation is only slightly more accurate than the refitted Chen and Yang (2019)’s correlation with two coefficient sets. But the new correlation yields smooth IFT predictions, avoiding the issue of discontinuous IFT predictions yielded by Chen and Yang (2019)’s correlation.

Appendices

Appendix A: PR EOS and α functions

The PR EOS (Peng and Robinson 1976) can be expressed as:

$$p = \frac{RT}{{v - b}} - \frac{a}{{v\left( {v + b} \right) + b\left( {v - b} \right)}}$$

(11)

where p is the pressure in bar; v stands for the molar volume in cm³/mol; T is the temperature in K; a is the attraction parameter with unit of bar cm⁶/mol², and b is the repulsion parameter with unit of cm³/mol; a and b can be determined by Eqs. (12) and (13):

$$a = 0.457535\frac{{R^{2} T_{{\text{c}}}^{2} }}{{p_{{\text{c}}} }}\alpha$$

(12)

$$b = 0.077796\frac{{RT_{{\text{c}}} }}{{p_{{\text{c}}} }}$$

(13)

where R is the universal gas constant in J/(mol K); T_c is the critical temperature in K; p_c is the critical pressure in bar; and α is the so-called alpha function.

The expression of Twu α function can be written as (Twu et al. 1991):

$$\alpha \left( {T_{{\text{r}}} } \right) = T_{{\text{r}}}^{{N\left( {M - 1} \right)}} {\exp}\left[ {L\left( {1 - T_{{\text{r}}}^{MN} } \right)} \right]$$

(14)

where T_r is the reduced temperature; L, M and N are compound-specific parameters. The values of these parameters regressed by Martinez et al. (2018) are used in this study.

Gasem α function can be expressed by (Gasem et al. 2001):

$$\alpha \left( T \right) = {\exp}\left( {\left( {A + BT_{{\text{r}}} } \right)\left( {1 - T_{{\text{r}}}^{{C + D\omega + E\omega^{2} }} } \right)} \right)$$

(15)

where the values of correlation parameters A through E are 2.0, 0.836, 0.134, 0.508 and −0.0467, respectively.

Appendix B: Summary of van der Waals one-fluid mixing rule and its BIPs

The van der Waals one-fluid mixing rule can be expressed as (van der Waals 1873):

$$a_{{\text{m}}} = \mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} z_{i} z_{j} \sqrt {a_{i} a_{j} } \left( {1 - kc_{ij} } \right)$$

(16)

$$b_{{\text{m}}} = \mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} z_{i} z_{j} \frac{{\left( {b_{i} + b_{j} } \right)}}{2}\left( {1 + kd_{ij} } \right)$$

(17)

where z_i is the molar fraction of the ith component in the mixture; a_i and b_i can be calculated by Eqs. (12) and (13); kc_ij and kd_ij are the BIPs that need to be fitted. In this study, the linear temperature-dependent BIP correlations from the study by Abudour et al. (2012a, b) are applied for CO₂/H₂O mixtures. Table

Table 12 BIPs correlations in the van der Waals mixing rule as obtained by Abudour et al. (2012a, b)

12 lists the BIP correlations obtained by Abudour et al. (2012a, b).

When vdW mixing rule is used in PR EOS, the fugacity coefficient can be written as:

$${\ln}\varphi_{i} = \frac{{bb_{i} }}{{b_{{\text{m}}} }}\left( {Z - 1} \right) - {\ln}\left( {Z - B} \right) - \frac{A}{2\sqrt 2 B}\left( {\frac{{2aa_{i} }}{{a_{{\text{m}}} }} - \frac{{bb_{i} }}{{b_{{\text{m}}} }}} \right){\ln}\left( {\frac{{Z + \left( {1 + \sqrt 2 } \right)B}}{{Z - \left( {1 + \sqrt 2 } \right)B}}} \right)$$

(18)

where

$$bb_{i} = 2\mathop \sum \limits_{j = 1}^{n} z_{j} \frac{{b_{i} + b_{j} }}{2}\left( {1 + kd_{ij} } \right) - b_{{\text{m}}}$$

(19)

$$aa_{i} = \mathop \sum \limits_{j = 1}^{n} z_{j} \sqrt {a_{i} a_{j} } \left( {1 - kc_{ij} } \right)$$

(20)

where Z is the compressibility factor. For PR EOS, Z can be calculated by Eq. (21).

$$Z^{3} - \left( {1 - B} \right)Z^{2} + \left( {A - 3B^{2} - 2B} \right)Z - \left( {AB - B^{2} - B^{3} } \right) = 0$$

(21)

where

$$A = \frac{{a_{{\text{m}}} p}}{{R^{2} T^{2} }}$$

(22)

$$B = \frac{{b_{{\text{m}}} p}}{RT} .$$

(23)

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):

$$b_{{\text{m}}} = \mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} z_{i} z_{j} \frac{{\left( {b_{i} + b_{j} } \right)}}{2}$$

(24)

$$a_{{\text{m}}} = b_{{\text{m}}} \left[ {\mathop \sum \limits_{i = 1}^{n} z_{i} \frac{{a_{i} }}{{b_{i} }} - \frac{{G_{{\infty}}^{{\text{E}}} }}{\varLambda}} \right]$$

(25)

where $G_{{\infty}}^{{\text{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):

$$G_{{\infty}}^{{\text{E}}} = RT\mathop \sum \limits_{i = 1}^{n} z_{i} \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{ji} b_{j} z_{j} {\exp}\left( { - \alpha_{ji} \tau_{ji} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{ki} \tau_{ki} } \right)}}$$

(26)

where

$$\tau_{ji} = \frac{{\Delta g_{ji} }}{RT}$$

(27)

$$g_{ii} = -\varLambda \frac{{a_{i} }}{{b_{i} }}$$

(28)

$$g_{ij} = - 2\frac{{\sqrt {b_{i} b_{j} } }}{{b_{i} + b_{j} }}\sqrt {g_{ii} g_{jj} } \left( {1 - k_{ij} } \right)$$

(29)

The generalized BIP correlations for τ_ij obtained by Aasen et al. (2017) are given below:

$$\frac{{g_{12} }}{{RT_{0} }} = 5.831 - 2.559\left( {\frac{T}{{T_{0} }}} \right)$$

(30)

$$\frac{{g_{21} }}{{RT_{0} }} = - 3.311 + 0.03770\left( {\frac{T}{{T_{0} }}} \right)$$

(31)

where T₀ = 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):

$${\ln}\varphi_{i} = \frac{{b_{i} }}{{b_{{\text{m}}} }}\left( {Z - 1} \right) - {\ln}\left( {Z - B} \right) - \frac{1}{2\sqrt 2 }\left( {\frac{{a_{i} }}{{b_{i} RT}} + \frac{{{\ln}\gamma_{i} }}{{\varLambda}}} \right){\ln}\left( {\frac{{Z + \left( {1 + \sqrt 2 } \right)B}}{{Z - \left( {1 + \sqrt 2 } \right)B}}} \right)$$

(32)

where lnγ_i is the activity coefficient of component i and can be expressed as (Zhao and Lvov 2016):

$$\begin{aligned} {\ln}\gamma_{i} = \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{ji} z_{j} b_{j} {\exp}\left( { - \alpha_{ji} \tau_{ji} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} z_{k} b_{k} {\exp}\left( { - \alpha_{ki} \tau_{ki} } \right)}} \\ + \mathop \sum \limits_{j = 1}^{n} \left[ {\frac{{b_{i} z_{j} {\exp}\left( { - \alpha_{ij} \tau_{ij} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} z_{k} b_{k} {\exp}\left( { - \alpha_{kj} \tau_{kj} } \right)}} \cdot \left( {\tau_{ij} - \frac{{\mathop \sum \nolimits_{l = 1}^{n} \tau_{lj} z_{l} b_{l} {\exp}\left( { - \alpha_{lj} \tau_{lj} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} z_{k} b_{k} {\exp}\left( { - \alpha_{kj} \tau_{kj} } \right)}}} \right)} \right] \\ \end{aligned}$$

(33)

The derivation of the expression of the activity coefficient in Huron-Vidal mixing rule is detailed in “Appendix H”.

Appendix D: Summary of volume translation models applied in this study

The constant volume translation can be expressed as (Peneloux et al. 1982; Jhaveri and Youngren 1988):

$$v_{{{\text{corr}}}} = v_{{{\text{EOS}}}} - \mathop \sum \limits_{i = 1}^{n} z_{i} c_{i}$$

(34)

where v_corr is corrected molar volume in cm³/mol; v_EOS stands for PR-EOS-calculated molar volume in cm³/mol; c_i is the component-dependent volume shift parameter which can be determined by Eq. (35) (Young et al. 2017).

$$c_{i} = s_{i} \times b_{i}$$

(35)

The values of s_i used by Liu et al. (2016) are applied in this study ($s_{{{\text{H}}_{2} {\text{O}}}} = 0.23170$ and $s_{{{\text{CO}}_{2} }} = - 0.15400$).

Abudour volume translation model can be expressed as (Abudour et al. 2012a, b):

$$v_{{{\text{corr}}}} = v_{{{\text{EOS}}}} + c - \delta_{{\text{c}}} \left( {\frac{0.35}{{0.35 + d}}} \right)$$

(36)

where $\delta_{{\text{c}}}$ is volume correction at the critical temperature in cm³/mol; d is the dimensionless distance function given by (Mathias et al. 1989; Abudour et al. 2012b):

$$d = \frac{1}{{RT_{{\text{c}}} }}\left( {\frac{{\partial p^{{{\text{PR}}}} }}{\partial \rho }} \right)_{T}$$

(37)

where ρ is the molar density in mol/cm³. The volume translation function proposed by Abudour et al. (2013) was extended to mixtures by the following equations (Abudour et al. 2013):

$$v_{{{\text{corr}}}} = v_{{{\text{EOS}}}} + c_{{\text{m}}} - \delta_{{{\text{cm}}}} \left( {\frac{0.35}{{0.35 + d_{{\text{m}}} }}} \right)$$

(38)

where (Abudour et al. 2013):

$$c_{{\text{m}}} = \left( {\frac{{RT_{{{\text{cm}}}} }}{{p_{{{\text{cm}}}} }}} \right)\left( {c_{{1{\text{m}}}} - \left( {0.004 + c_{{1{\text{m}}}} } \right){\text{e}}^{{ - 2d_{{\text{m}}} }} } \right)$$

(39)

$$c_{{1{\text{m}}}} = \mathop \sum \limits_{i = 1}^{n} z_{i} c_{{1{i} }}$$

(40)

$$d_{{\text{m}}} = \frac{1}{{RT_{{{\text{cm}}}} }}\left( {\frac{{\partial p^{PR} }}{\partial \rho }} \right)_{T} - \left( {\frac{1}{{RT_{{{\text{cm}}}} \rho^{2} }}} \right)\frac{{a_{v1}^{2} }}{{a_{11} }}$$

(41)

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_{{\text{c}}}$) (Abudour et al. 2012b):

$$c_{1} = 0.4266Z_{{\text{c}}} - 0.1101$$

(42)

The term d_m can be derived using the original PR EOS (Matheis et al. 2016):

$$d_{{\text{m}}} = \frac{{v^{2} }}{{RT_{{{\text{cm}}}} }}\left[ {\frac{RT}{{\left( {v - b} \right)^{2} }} - \frac{{2a\left( {v + b} \right)}}{{\left( { - b^{2} + 2bv + v^{2} } \right)^{2} }}} \right]$$

(43)

The volume correction of the given mixture at the critical point, δ_cm, can be determined by (Abudour et al. 2013):

$$\delta_{{{\text{cm}}}} = 0.3074\frac{{RT_{{{\text{cm}}}} }}{{p_{{{\text{cm}}}} }} - \mathop \sum \limits_{i = 1}^{n} \theta_{i} v_{{{\text{c}}i}}$$

(44)

where v_ci is the critical volume of component i; θ_i is the surface fraction of component i defined by (Abudour et al. 2013):

$$\theta_{i} = \frac{{z_{i} v_{{{\text{c}}i}}^{2/3} }}{{\mathop \sum \nolimits_{i = 1}^{n} z_{i} v_{{{\text{c}}i}}^{2/3} }}$$

(45)

The critical temperature of the mixture can be calculated via the following mixing rule (Abudour et al. 2013).

$$T_{{{\text{cm}}}} = \mathop \sum \limits_{i = 1}^{n} \theta_{i} T_{{{\text{c}}i}}$$

(46)

The critical pressure of the mixture can be determined by the correlation proposed by Aalto et al. (1996):

$$p_{{{\text{cm}}}} = \frac{{\left( {0.2905 - 0.085\omega_{{\text{m}}} } \right)RT_{{{\text{cm}}}} }}{{\mathop \sum \nolimits_{i = 1}^{n} \theta_{i} v_{{{\text{c}}i}} }}$$

(47)

where ω_m is the acentric factor of the mixture (Abudour et al. 2013):

$$\omega_{{\text{m}}} = \mathop \sum \limits_{i = 1}^{n} z_{i} \omega_{i}$$

(48)

where ω_i is the acentric factor of component i.

Appendix E: Summary of existing IFT correlations for CO₂/H₂O mixtures

Parachor model (Sugden 1930) can be expressed as below (Schechter and Guo 1998):

$$\sigma = \left[ {\mathop \sum \limits_{i = 1}^{n} P_{i} \left( {x_{i} \rho_{{\text{L}}}^{{\text{M}}} - y_{i} \rho_{{\text{V}}}^{{\text{M}}} } \right)} \right]^{4}$$

(49)

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_{{{\text{H}}_{2} {\text{O}}}} = 52$, $P_{{{\text{CO}}_{2} }} = 78$) (Liu et al. 2016); $\rho_{{\text{L}}}^{{\text{M}}}$ is the molar density of liquid phase in mol/cm³; $\rho_{{\text{V}}}^{{\text{M}}}$ is the molar density of vapor phase in mol/cm³.

Hebach et al. (2002)’s correlation can be expressed as:

$$\begin{aligned} \sigma = & k_{0} \left( {1 - {\exp}\left( {k_{1} \sqrt {dd} } \right)} \right) + k_{2} \cdot dd + k_{3} \cdot dd^{2} \\ + k_{4} \cdot dd^{3} + k_{5} {\exp}\left( {k_{6} \left( {dd - 0.9958} \right)} \right) \\ \end{aligned}$$

(50)

where (Hebach et al. 2002):

$$dd = \left( {\rho_{{{\text{H}}_{2} {\text{O}}}} - \rho_{{{\text{corr}}}} } \right)^{2}$$

(51)

$$\rho_{{{\text{corr}}}} = \left\{ {\begin{array}{ll} {\frac{{\rho_{{{\text{CO}}_{2} }} + b_{0} \left( {304 - T} \right)\left( {10 \times p} \right)^{{b_{1} }} }}{1000}\qquad 0.025\,{\text{g}}/{\text{cm}}^{3} < \rho_{{{\text{CO}}_{2} }} < 0.25\,{\text{g}}/{\text{cm}}^{3} } \\ {\qquad \qquad \qquad \rho_{{{\text{CO}}_{2} }}\qquad \qquad {\text{in other cases}}} \\ \end{array} } \right.$$

(52)

where $\rho_{{{\text{CO}}_{2} }}$ is the CO₂-rich-phase density in g/cm³; $\rho_{{{\text{H}}_{2} {\text{O}}}}$ is the aqueous-phase density in g/cm³; k₀ to k₆ and b₀ to b₁ are empirical coefficients. The units of T, p, and dd are K, bar, and g²/cm⁶, respectively. Table

Table 13 Coefficients in Hebach et al. (2002)’s correlation

13 lists the values of original and refitted coefficients.

Chen and Yang (2019)’s correlation is given as:

$$\sigma = C_{1} + \left( {C_{2} p_{{\text{r}}} + C_{3} } \right){\ln}K_{{{\text{CO}}_{2} }} + \left( {C_{4} p_{{\text{r}}} + C_{5} } \right){\ln}K_{{{\text{H}}_{2} {\text{O}}}}$$

(53)

where σ is IFT in mN/m; p_r is the reduced pressure of CO₂; C₁ to C₅ are empirical coefficients. To make fair comparisons, we refit these coefficients based on the IFT database employed in this study. Table

Table 14 Refitted coefficients in Chen and Yang’s correlation

14 summarizes the values of these refitted coefficients.

Appendix F: Pressure–temperature coverage of phase-density and IFT data collected from the literature.

Figure 

Fig. 8

Pressure–temperature coverage of phase-density and IFT data collected from the literature. The solid curves stand for pure-CO₂ (left) and pure-H₂O (right) saturation curves, respectively

8 depicts the pressure–temperature coverage of phase-density and IFT data collected from the literature over 278–478.35 K and 2.5–1291.1 bar, and 278.15–477.59 K and 1–1200.96 bar, respectively.

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 

Fig. 9

Identification of the outliers at T = 278.15–288.15 K (a) and T = 293.15–298.15 K (b). Outliers are from the studies by Chun and Wilkinson (1995) and Park et al. (2005)

9 shows the identification of the outliers from the collected data over 40–60 bar and 278.15–298.15 K, while Fig. 

Fig. 10

Identification of the outliers at moderate temperature (307.15–314.15 K) conditions. Outliers are from the study by Bachu and Bennion (2009)

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 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:

$${\ln}\gamma_{i} = \frac{1}{RT}\frac{{\partial G_{{\infty}}^{{\text{E}}} }}{{\partial z_{i} }}$$

(54)

where the excess Gibbs free energy can be expressed as (Huron and Vidal 1979):

$$G_{{\infty}}^{{\text{E}}} = RT\mathop \sum \limits_{i = 1}^{n} z_{i} \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{ji} b_{j} z_{j} {\exp}\left( { - \alpha_{ji} \tau_{ji} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{ki} \tau_{ki} } \right)}}$$

(55)

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:

$$\begin{aligned} \frac{{G_{{\infty}}^{{\text{E}}} }}{{RT}} &= z_{1} \cdot \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}} \\ &\quad+ z_{2} \cdot \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j2} b_{j} z_{j} {\exp}\left( { - \alpha_{j2} \tau_{j2} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)}} \\ \end{aligned}$$

(56)

Taking the partial derivative of the first term in the right hand side of Eq. (56) yields:

$$\begin{aligned} \frac{\partial }{{\partial z_{1} }} & \left( {z_{1} \cdot \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}}} \right) \\ &\quad= \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}} + z_{1} \\ &\qquad\cdot \left( {\frac{{\tau_{11} b_{1} {\exp}\left( { - \alpha_{11} \tau_{11} } \right) \cdot \left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)} \right)}}{{\left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)} \right)^{2} }} - \frac{{b_{1} {\exp}\left( { - \alpha_{11} \tau_{11} } \right) \cdot \left( {\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)} \right)}}{{\left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)} \right)^{2} }}} \right) \\ \end{aligned}$$

(57)

Taking the partial derivative of the second term in the right hand side of Eq. (56) yields:

$$\begin{aligned} \frac{\partial }{{\partial z_{1} }} & \left( {z_{2} \cdot \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j2} b_{j} z_{j} {\exp}\left( { - \alpha_{j2} \tau_{j2} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)}}} \right) \\ &\quad= z_{2} \cdot \left( {\frac{{\tau_{12} b_{1} {\exp}\left( { - \alpha_{12} \tau_{12} } \right) \cdot \left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)} \right)}}{{\left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)} \right)^{2} }} - \frac{{b_{1} {\exp}\left( { - \alpha_{12} \tau_{12} } \right) \cdot \left( {\mathop \sum \nolimits_{j = 1}^{n} \tau_{j2} b_{j} z_{j} {\exp}\left( { - \alpha_{j2} \tau_{j2} } \right)} \right)}}{{\left( {\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)} \right)^{2} }}} \right) \\ \end{aligned}$$

(58)

As such, Eq. (54) can be expressed as:

$$\begin{aligned} \frac{1}{RT}\frac{{\partial G_{{\infty}}^{E} }}{{\partial z_{1} }} = \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}} \\ + \frac{{z_{1} b_{1} {\exp}\left( { - \alpha_{11} \tau_{11} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}} \\ \cdot \left( {\tau_{11} - \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j1} b_{j} z_{j} {\exp}\left( { - \alpha_{j1} \tau_{j1} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k1} \tau_{k1} } \right)}}} \right) \\ + \frac{{z_{2} b_{1} {\exp}\left( { - \alpha_{12} \tau_{12} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)}} \\ \cdot \left( {\tau_{12} - \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{j2} b_{j} z_{j} {\exp}\left( { - \alpha_{j2} \tau_{j2} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{k2} \tau_{k2} } \right)}}} \right) \\ \end{aligned}$$

(59)

Using letter i to replace number 1 leads to a general expression of activity coefficient (Huron and Vidal 1979):

$$\begin{aligned} {\ln}\gamma_{i} &= \frac{1}{RT}\frac{{\partial G_{{\infty}}^{{\text{E}}} }}{{\partial z_{i} }} \\ &= \frac{{\mathop \sum \nolimits_{j = 1}^{n} \tau_{ji} b_{j} z_{j} {\exp}\left( { - \alpha_{ji} \tau_{ji} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{ki} \tau_{ki} } \right)}} \\ &\quad+ \mathop \sum \limits_{j = 1}^{n} \left[ {\frac{{b_{i} z_{j} {\exp}\left( { - \alpha_{ij} \tau_{ij} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{kj} \tau_{kj} } \right)}}\left( {\tau_{ij} - \frac{{\mathop \sum \nolimits_{l = 1}^{n} \tau_{lj} b_{l} z_{l} {\exp}\left( { - \alpha_{lj} \tau_{lj} } \right)}}{{\mathop \sum \nolimits_{k = 1}^{n} b_{k} z_{k} {\exp}\left( { - \alpha_{kj} \tau_{kj} } \right)}}} \right)} \right] \\ \end{aligned}$$

(60)