An Accurate Model to Calculate CO 2 Solubility in Pure Water and in Seawater at Hydrate–Liquid Water Two-Phase Equilibrium

: Understanding of CO 2 hydrate–liquid water two-phase equilibrium is very important for CO 2 storage in deep sea and in submarine sediments. This study proposed an accurate thermodynamic model to calculate CO 2 solubility in pure water and in seawater at hydrate–liquid water equilibrium (HL W E). The van der Waals–Platteeuw model coupling with angle-dependent ab initio intermolecular potentials was used to calculate the chemical potential of hydrate phase. Two methods were used to describe the aqueous phase. One is using the Pitzer model to calculate the activity of water and using the Poynting correction to calculate the fugacity of CO 2 dissolved in water. Another is using the Lennard–Jones-referenced Statistical Associating Fluid Theory (SAFT-LJ) equation of state (EOS) to calculate the activity of water and the fugacity of dissolved CO 2 . There are no parameters evaluated from experimental data of HL W E in this model. Comparison with experimental data indicates that this model can calculate CO 2 solubility in pure water and in seawater at HL W E with high accuracy. This model predicts that CO 2 solubility at HL W E increases with the increasing temperature, which agrees well with available experimental data. In regards to the pressure and salinity dependences of CO 2 solubility at HL W E, there are some discrepancies among experimental data. This model predicts that CO 2 solubility at HL W E decreases with the increasing pressure and salinity, which is consistent with most of experimental data sets. Compared to previous models, this model covers a wider range of pressure (up to 1000 bar) and is generally more accurate in CO 2 solubility in aqueous solutions and in composition of hydrate phase. A computer program for the calculation of CO 2 solubility in pure water and in seawater at hydrate–liquid water equilibrium can be obtained from the corresponding author via email.


Introduction
There has been a growing concern on the potential impact of rising greenhouse gas levels in the atmosphere. CO 2 is the greatest contributor to the greenhouse effect. Consequently, carbon capture and storage (CCS) becomes a promising option to reduce anthropogenic CO 2 emissions into the atmosphere. Geological storage and ocean storage are the most common storage methods for large CO 2 volumes. Proposed sites for CO 2 geological storage range from mined salt caverns, saline aquifers, depleted hydrocarbon reservoirs, unmineable coal seams, and so on [1]. Storage of CO 2 in depleted oil/gas reservoirs and unmineable coal beds is more economical than other CO 2 storage strategies due to additional oil and gas recovery, characterized geological conditions, existing installed equipment, and so on. Saline aquifers are favored due to their extended storage capacity. Storages of CO 2 in depleted oil/gas reservoirs and in saline aquifers are more technically feasible [2,3]. Tens of pilot and commercial projects for CO 2 storage in depleted oil/gas reservoirs or in saline aquifers have been launched to date [3,4]. Injecting CO 2 into the deep ocean to form CO 2 hydrate is attractive since the CO 2 storage capacity of deep oceans is huge [5]. However, this method is more controversial than other CO 2 storage methods due to the risk of leakage, local acidification of seawater, the correspondingly negative impact on the marine environment, and so on [4,6,7]. Recent studies relevant to ocean storage of CO 2 tried to develop numerical models for better evaluation of storage efficiency and risk assessment and investigate the determination and quantification of ocean site selection criteria [7,8]. In a hybrid method between ocean and geological storage, CO 2 would be injected into submarine sediments below the deep ocean floor to avoid the potential leakage of CO 2 and harm to marine ecosystems [9].
In addition to the deep ocean, CO 2 injected into shallow zones in deep-sea sediments or cold saline aquifers will dissolve in seawater or brine and then convert into CO 2 hydrate [10][11][12]. Alternatively, Ohgaki et al. [13] proposed the method to replace CH 4 in methane hydrate in the seafloor sediments with CO 2 , which provides twin advantages: sequestrating CO 2 and enhancing CH 4 hydrate recovery. A recent study found that the direct use of a CO 2 + N 2 gas mixture (flue gas), instead of pure CO 2 , replacing CH 4 can greatly enhance the overall CH 4 recovery rate and reduce the costs [14][15][16].
Gas hydrates, or clathrate hydrates, are non-stoichiometric crystalline solids that consist of a hydrogen-bonded network of host water molecules and enclathrated guests molecules such as CH 4 , C 2 H 6 , N 2 , and CO 2 [17]. They are thermodynamically stable at a low temperature and high pressure. Within the stable region in deep sea or subsurface aquifers, CO 2 hydrate coexists with seawater or brine while being absent of a CO 2 -rich phase. Therefore, understanding CO 2 hydrate-liquid water two-phase equilibrium is very important for CO 2 storage.
The stability boundary of CO 2 hydrate, corresponding to the pressure (P)-temperature (T) conditions for CO 2 hydrate-liquid water-CO 2 vapor (HL W V) three-phase equilibrium, has been studied extensively [18]. Based on the van der Waals-Platteeuw (vdW-P) model [19], many thermodynamic models were developed to calculate P-T conditions for HL w V three-phase equilibrium of gas hydrates (including CO 2 hydrate). In recent years, great improvements on the vdW-P model included the calculation of interaction potentials between gas molecules and water molecules using ab initio methods or quantum mechanics [17,18,[20][21][22][23][24][25][26][27][28], taking into account effects of lattice distortion and guest-guest interactions (which is essential for modeling multiple occupancy of small guests) by extension of the vdW-P model or using alternative methods [29][30][31][32][33][34][35].
In contrast to numerous studies for three-phase equilibrium of CO 2 hydrate, studies on CO 2 hydrate-liquid water two-phase equilibrium are relatively limited. Aya et al. [76], Yang et al. [77], Lee et al. [78], Servio and Englezos [79], Zhang [80], Someya et al. [81], Kim et al. [82], Zhang et al. [83], Sun et al. [84], and Geng et al. [85] measured CO 2 solubility in water at hydrate-liquid water equilibrium (HL W E). All measurements show that CO 2 solubility at HL W E increase with increasing temperature. In regards to the effect of pressure, available experimental data show weak pressure dependence of CO 2 solubility at HL W E although there are some discrepancies among them.
Yang et al. [77], Hashemi et al. [86], and Zhang et al. [83] proposed models to calculate CO 2 solubility in water at HL W E. All these models adopted the van der Waals-Platteeuw hydrate model [19] to describe the hydrate phase. Differences between them are how to calculate the Langmuir constants and how to describe the aqueous phase. The model of Yang et al. [77] used an empirical formula to calculate the Langmuir constants and used the lattice fluid equation of state to calculate thermodynamic properties of the aqueous phase. Yang's model agrees with the experimental data reported in the same work. However, data measured by Yang et al. [77] are larger than those measured by other studies. Hashemi et al. [86] adopted the Langmuir constants reported by Parrish and Prausnitz [87] and used the Trebble-Bishnoi EOS [39] to describe the aqueous phase. Zhang et al. [83] applied the lattice distortion theory developed by Holder and co-workers [88][89][90] to calculate the Langmuir constants and used Margules equations [91] to describe the aqueous phase. Calculations of this model agree with the experimental data measured by the same study. However, the lattice distortion theory overestimates the occupancy fraction of CO 2 in small cage. Recently, Sun et al. [92] combined the Chen-Guo hydrate model [93,94] and the Patel-Teja EOS [38] to calculate CO 2 solubilities in water and aqueous NaCl solution under HL W E. All these models predict that CO 2 solubility at HL W E increases with the increase of temperature, which is consistent with available experimental data. About pressure dependence of CO 2 solubility at HL W E, models of Yang et al. [77], Zhang et al. [83], and Sun et al. [92] suggested a negative trend (i.e., CO 2 solubility at HL W E decreases with the increasing pressure at constant temperature), whereas the model of Hashemi et al. [86] gives a positive trend. Bergeron et al. [95,96] made thermodynamic analyses on temperature and pressure dependence of CO 2 solubility at HL W E. Although SAFT EOS and CPA EOS have been widely used to model three-phase equilibrium of gas hydrates as mentioned above, they have not been used to model HL W E of gas hydrates to date. Sun and Duan [97] applied the method to calculate Langmuir constants from ab initio intermolecular potentials to model HL W E of CH 4 hydrate. However, this method has not been used to model HL W E of CO 2 hydrate.
The purpose of this study is to develop a new model to predict CO 2 solubility at HL W E in marine environments by taking into account the effect of temperature, pressure and salinity together. Van der Waals-Platteeuw basic hydrate model [19] is used to calculate the chemical potential of hydrate phase. The method to calculate Langmuir constants from angle-dependent ab initio intermolecular potentials, which can well predict cage occupancies of CO 2 molecule in small and large cages of sI hydrate [18], is followed by this study. The aqueous solution is described with two methods. One is using the Pitzer model [42] to calculate the activity of water and using the Poynting correction to calculate the fugacity of CO 2 dissolved in water. Another is using the Lennard-Jones-referenced Statistical Associating Fluid Theory (SAFT-LJ) EOS [55] to calculate the activity of water and the fugacity of dissolved CO 2 . The next section will introduce the principle of the thermodynamic model proposed by this study. Section 3 compares the experimental data with predictions of this model for CO 2 solubility in aqueous solution at HL W E. Finally, some conclusions are drawn.

Model Framework for HL W E
According to principles of thermodynamics, the fugacities or chemical potentials of species in the various phases must be equal at phase equilibrium. For CO 2 hydrate-liquid water two-phase equilibrium, the following equation must be satisfied: where µ is the chemical potential, ∆µ is the difference of chemical potential, the subscript W represents H 2 O, the superscript H represents hydrate phase, β represents the hypothetical empty hydrate lattice, and L represent liquid water (aqueous phase). ∆µ H W is calculated from the van der Waals-Platteeuw model [19] based on the statistical mechanics: where ν i is the number of i-type cages (also called "cavities") per water molecule, θ ij is the fractional occupancy of i-type cavities with j-type guest molecules, and C ij is the Langmuir constant of gas component j in i-type cavity. For sI hydrate, the number of small cavities, v 1 , equals to 1/23, the number of large cages, v 2 , equals to 3/23. f j in Equation (3) is the fugacity of gas component j in the hydrate phase, which equals to the fugacity of gas component in the aqueous phase at HL W E. Our previous study [18] presented a method to calculate the Langmuir constant from angle-dependent ab initio intermolecular potentials, which not only can represent P-T conditions of hydrate-liquid water-vapor three-phase equilibrium of pure CH 4 hydrate and CO 2 hydrate, but can also well-predict cage occupancies of CH 4 or CO 2 in small and large cages of sI hydrate. Thus, this method was followed by this study. It is the first time to apply Langmuir constants calculated from ab initio intermolecular potentials to model HL W E of CO 2 hydrate. In order to facilitate the application, the Langmuir constant of CO 2 in small and large cage of structure I (sI)-type hydrate, C Is,CO 2 and C Il,CO 2 , calculated from angle-dependent ab initio intermolecular potentials was fit using empirical equations as follows: C Is,CO 2 = exp(− 25.354129 C Il,CO 2 = exp(−23.632139 + 3697.1701/T).
Since ions (e.g., Na + , Cl − ) do not enter cages of the hydrate lattice, this study neglects the influence of ions on Langmuir constants as previous studies did. ∆µ L W in Equation (1) is calculated from the thermodynamic correction equation improved by Holder et al. [98]: where a W is the activity of water, ∆µ 0 W (T 0 , 0) is the reference chemical potential difference at the reference temperature, T 0 , usually taken to be 273.15 K, and zero pressure. ∆h L W is the difference of enthalpy between hydrate and liquid water. The variation of ∆h L W with temperature is described with the following equation: where ∆h 0 W (T 0 ) is the reference enthalpy difference at T 0 , ∆C L p is the difference of isobaric molar heat capacity. The variation of ∆C L p with temperature is linear: ∆µ 0 w and ∆h 0 W determined by our previous study [18] and ∆C 0 p and b determined by Parrish and Prausnitz [87] is adopted by this model.
∆V L W in Equation (6) is the difference of the molar volume between hydrate and liquid water. For sI hydrate, it equals to 4.6 cm 3 /mol after neglecting small variations of molar volumes of hydrate lattice and liquid water with pressure.
All solutes dissolved in water, such as CO 2 and ions, will change the activity of water. The fugacity of CO 2 dissolved is affected by ions. This study adopted two methods to calculate f CO 2 and a W . One is using the Pitzer model [42] to calculate the activity of water and using the Poynting correction to calculate the fugacity of CO 2 dissolved in water, as we did for HL W E of CH 4 hydrate [97]. Another is using the SAFT-LJ EOS [55] to calculate the activity of water and the fugacity of dissolved CO 2 .

Method 1: Poynting Correction + Pitzer Model
At hydrate-liquid water two-phase equilibrium, vapor CO 2 or liquid CO 2 phase is absent. The fugacity of CO 2 in hydrate phase is equal to that in aqueous phase. The latter can be calculated from the Poynting correction equation [99,100] as following: where V aq CO 2 means the partial molar volume of CO 2 in aqueous solution, P sat is defined as the pressure required to obtain a given solubility, x CO 2 . At P sat , CO 2 -rich vapor (or liquid) coexists with H 2 O-rich liquid (aqueous phase) so that P sat can be calculated from CO 2 solubility model developed by Duan et al. [101] for vapor-liquid equilibrium (VLE). f sat CO 2 in Equation (9) represents the fugacity of CO 2 at P sat , which is calculated from the EOS of Duan et al. [102]. V aq CO 2 can be evaluated from PVTX and solubility data of the CO 2 -H 2 O system based on thermodynamic methods [103]. The value of V aq CO 2 determined by this study is equal to 31.5 cm 3 /mol at temperatures below 300 K.
The Pitzer model [42] was used to calculate the activity of water in aqueous electrolytes solution containing dissolved CO 2 . It is well accepted that the Pitzer model can calculate thermodynamic properties of aqueous electrolyte solutions over a wide P-T-m range (from the saturation pressure to 1000 bar, from the Eutectic temperature to 523~573 K, from dilute solution to saturated solution) with high accuracy when sufficient data are available to constrain the regression of Pitzer ion-interaction parameters. Temperature-dependent ion-ion interaction parameters determined by Spencer et al. [104] are suitable for aqueous solutions containing Na + , K + , Mg 2+ , Ca 2+ , Cl -, and SO 4 2− at temperatures below 298 K. These parameters and CO 2 -ion interaction parameters determined by Duan et al. [101] were adopted by Sun and Duan [18] in the model for three-phase equilibrium of CO 2 hydrate. Although some approximations had to be made to evaluate CO 2 -ion interaction parameters (such as neglecting the ionization of the carbonic acid and neglecting interactions between CO 2 molecules), three-phase equilibria of CO 2 hydrate in the CO 2 -H 2 O-NaCl system and CO 2 -seawater system have been calculated with high accuracy. Thus, this study tried to check the performance of the Pitzer model with parameters determined by Spencer et al. [104] and Duan et al. [101] in modeling HL W E of CO 2 hydrate. The details about the Pitzer model is omitted here to keep the paper concise. The readers can find them in the paper of Duan and Sun [105].

Method 2: SAFT-LJ EOS
The Statistical Association Fluid Theory (SAFT) EOS is one type of advanced EOS based on statistical mechanics. It describes the intermolecular association (hydrogenbonding force) based on first-order thermodynamic perturbation theory of Wertheim [46][47][48] so that it is suitable for fluid system containing associating molecules (e.g., water, methanol). Depending on the type of reference fluid selected, there are various SAFT-type EOSs, such as SAFT-VR (Variable Range) [106], SAFT1 [107], PC (Perturbed Chain)-SAFT [108], and SAFT-LJ (Lennard-Jones) [109]. The SAFT-LJ EOS improved by Sun and Dubessy [54,55], which takes into account multi-polar interaction and ion-ion interaction besides association, allows for calculation of VLE and PVTx properties of the CO 2 -H 2 O-NaCl system at temperatures below 573 K. It was adopted by this study to calculate thermodynamic properties of aqueous solutions. As we know, it is the first time to apply the SAFT-type EOS to modeling HL W E. The improved SAFT-LJ EOS is defined in terms of the dimensionless residual Helmholtz energy a res as follows: where a represents the dimensionless molar Helmholtz energy, the superscript seg, chain, assoc, polar, and ion represent the contributions from the repulsive and the attractive forces between Lennard-Jones segments, the chain formation, the intermolecular association, the multi-polar interaction, and ion-ion interaction, respectively. a chain and a assoc in Equation (10) are calculated from Wertheim's first-order thermodynamic perturbation theory [46][47][48], a seg , a polar , and a ion is calculated based on the Kolafa-Nezbeda equation [110], the Gubbins-Twu equation [111], and the primitive model of the mean-spherical approximation [112]. Detailed equations (Equations (A1)-(A30)) for a seg , a chain , a assoc , a polar , and a ion are listed in Appendix A and can also be found in Sun and Dubessy [54,55]. Parameters for pure H 2 O, CO 2 and NaCl, binary parameters for interactions between H 2 O and NaCl, and binary parameters for interactions between CO 2 and NaCl determined by Sun and Dubessy [55] were followed by this study (these parameters are given in Tables A1 and A2 in Appendix A). Two binary parameters for interactions between H 2 O and CO 2, k 1,ij and k 2,ij (their definitions can be found in Equations (A12) and (A13) in Appendix A) were re-evaluated by this study in order to improve the accuracy of the SAFT-LJ EOS at temperatures below 300 K. New parameters are listed in Table 1. CO 2 solubility at H-L W equilibrium at a given P and T can be calculated by iteration of Equation (1). For Method 1, an initial guess for the value of x CO 2 is made firstly. Then, P sat is calculated from CO 2 solubility model of Duan et al. [101] for VLE, and f sat is calculated from EOS of Duan et al. [102]. The third step is to calculate f CO 2 (from Equation (9)) and a W (from the Pitzer model [42]) at P. Finally, f CO 2 and a W are substituted into Equations (2), (3), and (6) to calculate ∆µ H W and ∆µ L W . If ∆µ H W equals to ∆µ L W , the current value of x CO 2 is the answer. If ∆µ H W does not equal to ∆µ L W , we must change the value of x CO 2 and repeat calculation until the answer is found. For Method 2, the first step is to guess an initial value for x CO 2 , too. Then, both f CO 2 and a W at P are calculated from the improved SAFT-LJ EOS [55]. After that, they are substituted into Equations (2), (3), and (6)

Result and Discussion
Firstly, this study compared predictions of this model for H-L w equilibrium of the CO 2 -H 2 O system with available experimental data. Table 2 lists the absolute average deviations (AAD) of CO 2 solubilities calculated by this model from each solubility data sets. Figures 1-3 compare the calculated CO 2 solubility isobars, isopleths, and isotherms with experimental data, respectively. Both calculation results of Method 1 and those of Method 2 are given in Table 2 and Figures 1-3. It can be seen that the difference between Method 1 and Method 2 is small at pressures below 300 bar. The difference between two methods increases with pressure in that Method 2 gives larger pressure dependence of CO 2 solubility than Method 1. Figures 1-3 show that this model agrees well with most of experimental data sets. CO 2 solubilities calculated by this model are systematically lower than experimental data reported by Yang et al. [77] by 3-11.7%. A comparison of different experimental data sets shows that CO 2 solubilities measured by Aya et al. [76] are greater than those of Geng et al. [85] by 4-10% (Figure 1a), and CO 2 solubilities reported by Yang et al. [77] are larger than data measured by Servio and Englezos [79], Zhang et al. [83], Sun et al. [84], and Geng et al. [85] (Figure 1b). Data of Lee et al. [78] at 280 K are close to data from other studies, but data of Lee et al. [78] at 273 K and 276 K are much lower than data from other studies (Figure 1b). Calculations of this model for CO 2 solubility are consistent with measurements of Someya et al. [81] at 40 bar and 70 bar, but are lower than those measurements at 100 bar and 120 bar, since measurements of Someya et al. [81] show a positive pressure dependence of CO 2 solubility at HL W E.  Figure 1a shows that CO 2 solubility in water at HL W E increases with the increasing temperature whereas CO 2 solubility at VLE decreases with the increasing temperature at constant pressure. The prediction of this model for the temperature dependence of CO 2 solubility at HL W E is consistent with experimental data. As shown in Figures 2 and 3, the pressure dependence of CO 2 solubility at HL W E predicted by this model is negative, which is consistent with experimental data of Zhang [80], Zhang et al. [83], and Geng et al. [85]. It can also be seen from Figure 3 that the decreasing gradient of CO 2 solubility with pressure at constant temperature is small. For example, CO 2 solubility decreases by 10% when the pressure increases from 100 bar to 900 bar at 276.15 K according to the measurement of Geng et al. [85] and calculations of this model with Method 1. As mentioned before, Method 2 predicts larger pressure dependence of CO 2 solubility than Method 1. Calculations of Method 1 are more consistent with data of Geng et al. [85] than Method 2, whereas calculations of Method 2 are more consistent with data of Zhang [80] than Method 1. Experimental data of Yang et al. [77], Servio and Englezos [79], and Kim et al. [82] show very weak pressure dependence of CO 2 solubility at HL W E. Considering that the maximum pressure of these data is about 200 bar and the decrease of CO 2 solubility within 200 bar is close to the uncertainty of experimental data, it is indeed difficult to find the change tendency of CO 2 solubility with pressure from these data. On the other hand, experimental data measured by Zhang [80], Zhang et al. [83], and Geng et al. [85]), which cover a wider range of pressure than other experimental data, clearly show a negative pressure dependence of CO 2 solubility at HL W E. Anyway, more measurements on CO 2 solubility at HL W E are needed to verify this model.  Calculations of models of Hashemi et al. [86], Zhang et al. [83], and Yang et al. [77] were put into Figure 1, Figure 2, Figure 3, respectively. It can be seen that models of Zhang et al. [83] and Yang et al. [77] give a negative pressure dependence of CO 2 solubility at HL W E, which is consistent with our model. In contrast, the model of Hashemi et al. [86] gives a positive pressure dependence of CO 2 solubility at HL W E. The model of Yang et al. [77] gives a greater decrease gradient of CO 2 solubility with increasing pressure than our model using Method 2, while the decrease gradient of CO 2 solubility with increasing pressure predicted by the model of Zhang et al. [83] is comparable to that predicted by our model with Method 2. As shown in Figure 1, calculations of the model of Hashemi et al. [86] are lower than experimental data at 37 bar and 42 bar. Calculations of the model of Yang et al. [77] agree with experimental data reported in the same work, but are larger than those measured by other studies. The accuracy of the model of Zhang et al. [83] is comparable to our model since it adopted solubility models [101,113,114] which can accurately calculate CO 2 solubility at vapor-liquid equilibrium. However, the model of Zhang et al. [83] overestimates the occupancy fraction of CO 2 in small cage. Considering that composition of CO 2 hydrate phase at HL W VE calculated by Sun and Duan [18] is accurate and the variation of composition of CO 2 hydrate phase with pressure is small, we argue that calculations of this model for composition of CO 2 hydrate phase at HL W E are reliable. However, experimental data on composition of hydrate phase at HL W E are needed for verifying our model. In general, this model is generally more accurate than previous models in CO 2 solubility in aqueous solutions and in composition of hydrate phase. In addition, this model covers a larger range of pressure (up to 1000 bar) compared to previous models.  [77]. The dashed lines are CO 2 solubility at VL W E or liquid CO 2 -liquid water equilibrium (L C L W E). Symbols and represent experimental data of Geng et al. [85] and Yang et al. [77], respectively.
There are just a few experimental measurements of CO 2 solubility in aqueous electrolyte solutions at HL W E, including measurements of CO 2 solubility in aqueous NaCl solution made by Kim et al. [82] and Sun et al. [84], and measurement of CO 2 solubility in artificial seawater made by Zhang et al. [83]. Unfortunately, there are some discrepancies among these data. Data reported by Kim et al. [82] show a salting-in effect for CO 2 solubility at HL W E (in another word, CO 2 solubility in NaCl solution is greater than that in pure water at the same temperature and pressure). In contrast, data measured by Zhang et al. [83] and Sun et al. [84] show a salting-out effect. The absolute average deviation of this model from experimental data of Kim et al. [82] and Sun et al. [84] are given in Table 3. Calculations of this model for CO 2 solubility in aqueous NaCl solution at HL W E are lower than data of Kim et al. [82] because this model predicts a salting-out effect. Predictions of this model are consistent with experimental data of Sun et al. [84] in 1 wt% NaCl solution but are greater than their data in 3 wt% and 3.6 wt% NaCl solutions by 10-30%. Figure 4 shows that predictions of this model with Method 1 for CO 2 solubility in artificial seawater agree well with data of Zhang et al. [83] (the AAD is 1.6%). Considering that the effect of 35‰ seawater on CO 2 solubility is close to that of 3.6 wt% NaCl solution, the discrepancy between this model and data of Sun et al. [84] reflects the discrepancy between data of Zhang et al. [83] and Sun et al. [84].  For the effect of pressure on CO 2 solubility in aqueous electrolyte solutions at HL W E, experimental data measured by Zhang et al. [83] and Sun et al. [84] show a negative trends whereas data of Kim et al. [82] show very weak pressure dependence. As shown in Figure 4, this model predicts that CO 2 solubility in aqueous electrolyte solutions at HL W E decreases with the increasing pressure, similar to CO 2 solubility in pure water. More measurements on CO 2 solubility in aqueous NaCl solutions and in seawater at HL W E are needed to verify this model.

Conclusions
An accurate thermodynamic model was developed to predict CO 2 solubility in pure water and in seawater at hydrate-liquid water equilibrium. The van der Waals-Platteeuw basic hydrate model was used to describe the chemical potential of hydrate phase. Angledependent ab initio intermolecular potentials were used to calculate the Langmuir constants of CO 2 in hydrate cavities following our previous study. For the thermodynamic properties of the aqueous phase, two methods were adopted by this study. Method 1 is using the Pitzer model to calculate the activity of water and using the Poynting correction to calculate the fugacity of CO 2 dissolved in water, and Method 2 is using the SAFT-LJ EOS to calculate the activity of water and the fugacity of dissolved CO 2 . No parameters of this model were evaluated from experimental data of HL W E so that calculations of this model for HL W E are predictive.
The comparison of calculations of this model with experimental data shows that this model can predict CO 2 solubility in pure water and in seawater at HL W E with high accuracy. This model predicts that CO 2 solubility at HL W E increases with increasing temperature, which agrees with available experimental measurements. This model predicts that CO 2 solubility at HL W E decreases with increasing pressure and salinity. Most of the experimental data sets show negative pressure and salinity dependences of CO 2 solubility at HL W E except one experimental data set showing a positive pressure dependence of CO 2 solubility and one data set showing a positive salinity dependence of CO 2 solubility. The predictions of this model for pressure and salinity dependences of CO 2 solubility are consistent with most of experimental data. Two different methods to describe aqueous phase give similar change rate of CO 2 solubility with temperature and salinity whereas Method 2 gives larger decrease of CO 2 solubility with increasing pressure. This model covers a wider range of pressure (up to 1000 bar) than previous models and is generally more accurate than previous models in CO 2 solubility in aqueous solutions and in composition of hydrate phase. Anyway, more measurements on CO 2 solubility at HL W E are needed to verify this model. Extending our model to HL W E of CH 4 +CO 2 +N 2 hydrate will be done in future. A computer program for the calculation of CO 2 solubility in pure water and in seawater at hydrate-liquid water equilibrium can be obtained from the corresponding author via email.

Data Availability Statement:
There are no new data reported in this study and a computer program for the calculation of CO 2 solubility in pure water and in seawater at hydrate-liquid water equilibrium can be obtained from the corresponding author via email.

Acknowledgments:
The authors acknowledge three anonymous reviewers for their detailed, helpful, and pertinent comments which greatly improved the quality of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Appendix A. 1

. The Improved SAFT-LJ EOS
The SAFT-LJ EOS improved by Sun and Dubessy [55] is defined in terms of the dimensionless residual Helmholtz energy a res as follows: a res = a seg + a chain + a assoc + a polar + a ion , where a represents the dimensionless molar Helmholtz energy, the superscript seg, chain, assoc, polar, and ion represent the contributions from the repulsive and the attractive forces between Lennard-Jones segments, the chain formation, the intermolecular association, the multi-polar interaction, and ion-ion interaction, respectively. The following are formula for every Helmholtz energy term.

Appendix A.2. Lennard-Jones Segment Term
The equation for a seg is where m x is the segment number of mixtures, a seg 0 is the Helmholtz energy of one mole of non-associated spherical Lennard-Jones segments, and a HS is the molar Helmholtz energy of hard-sphere fluid. The Kolafa-Nezbeda equation [110] is used to calculate a seg 0 . All variables in Equation (A2) are calculated using the following equations: where the subscript x means the variable for mixture, the superscript * means the reduced variable, k B is the Boltzmann constant, T is temperature in K, V m is the molar volume in m 3 /mol, ρ is the density in mol/m 3 , η is the reduced packing density, m is the number of segments to form a non-spherical molecule, b is the molar volume of the segment, and ε is the Lennard-Jones energy parameter of the segment. ∆B 2,hBH is the residual second virial coefficient, d BH is the ratio of the hard-sphere diameter to the soft-sphere diameter of the segment. Coefficients C kl , γ, C k , D k , and D ln are component-independent, dimensionless, because they do not stand for any physical constants or variables. Their values were evaluated by Kolafa and Nezbeda [110] from computer simulation data for Lennard-Jones fluids. Parameters b, ε, and m are componentdependent. The value of the parameters b, ε, and m for pure CO 2 and pure H 2 O are evaluated from the saturated vapor pressure and liquid density data. They are given in Table A1. For mixtures, the calculation of parameters m x , b x , and ε x is based on van der Waals one-fluid mixing rule: where x i is the mole fraction of component i. Parameters for the i-j interactions, b ij and ε ij , are calculated by the following combining rules: The dimensionless parameters k 1,ij and k 2,ij are empirical corrections of the Lorentz-Berthelot rules and they are fitted to vapor-liquid equilibria data of the binary system. Appendix A. 3. Chain-Forming Term a chain takes into account the change of Helmholtz energy due to the presence of covalent chain-forming bonds among the segments belonging to one molecule. Including the chain-forming term is also the distinct character of SAFT-type EOSs from other EOSs based on perturbation theories. In SAFT, a chain is derived from Wertheim's theory in the limit of a covalent bonding between the segments. The resulted expression is as follows: where g LJ is the Lennard-Jones radial distribution function. ρ * i and T * i in Equation (A14) are defined as Appendix A. 4

. Association Term
SAFT EOS adopts the thermodynamic perturbation theory of Wertheim [46][47][48] to calculate the Helmholtz energy contribution due to the intermolecular association (hydrogenbonding force) as where the sum runs over the species i and all the associating sites of species i, and M i is the number of the associating sites of molecule i. X i A , the fraction of molecule i not bonded at site A, is obtained from the solution of the mass balance equation, where N Av is Avogrado constant, ∆ A i B j is the associating strength between site A of species i and site B of species j. The improved SAFT-LJ EOS only takes into account the selfassociation among water molecules and ignored the cross-associations between H 2 O and CO 2 , between H 2 O and ion, and between ion and ion (It was assumed that NaCl is fully dissociated in aqueous solution below 573 K) so that Equation (A18) can be solved analytically. Because H 2 O is modeled as a spherical molecule with four association sites, X A w can be calculated from the following equation: where the subscript w is the component H 2 O. The associating strength, ∆, depends on the radial distribution function and the interaction between two association sites. It is calculated using the following equation: where b w is the volume of one mole water molecules (not molar volume of fluid phase) in m 3 /mol, parameters ε assoc w and K assoc w are the association energy and association volume of water molecules, ρ * w and T * w are the reduced density and reduced temperature of water, respectively. The integral I is equal to 1.

. Multi-Polar Term
The improved SAFT-LJ EOS considered the dipole moment of H 2 O and the quadrupole moment of CO 2 but omit the quadrupole moment of H 2 O. The contribution of multi-polar interaction is calculated using the equation derived by Gubbins and Twu [111] based on the perturbation theory as follows: contains four two-body terms (a 3A1 , a 3A2 , a 3A3 , a 3A4 ) and four three-body terms (a 3B1 , a 3B2 , a 3B3 , a 3B4 ): a 3B = a 3B1 + 3a 3B2 + 3a 3B3 + a 3B4 .

. Charge-Charge Interaction Term
The restricted primitive model of mean spherical approximation [112] was adopted by the improved SAFT-LJ EOS to calculate a ion . In the restricted primitive model of mean spherical approximation (MSA), the solvent is regarded as a continuous medium with specified dielectric constant D and different ions are modeled as charged hard spheres of equal diameter. The following is the equation for a ion : where d is the diameter of hydrated ion, and χ is a dimensionless quantity defined by where κ is the Debye inverse screening length given by In Equation (A29), ε 0 is the vacuum permittivity, D w is the dielectric constant of water, z j is the valence of the charged ion j, e is the charge of an electron, and the summation is over all ions in the mixture.
d is the sole adjustable parameter of Equation (A27). The short-range interaction Lennard-Jones parameters for Na + and for Clare considered to be equal, i.e., b Na = b Cl , and ε Na = ε Cl . The parameters for the cross-interaction between ions and H 2 O are also assumed to be identical, i.e., b Na-H2O = b Cl-H2O , and ε Na-H2O = ε Cl-H2O . In the MSA theory, ions are assumed to be hard spheres and the dispersion force between ion and ion is not considered, thus ε Na is equal to 0. b Na , b Na-H2O , ε Na-H2O , and d were derived from the experimental data of the H 2 O-NaCl system and are listed in Table A2. Table A2. SAFT Similar to the cross-interaction parameters for ion-H 2 O, the values of the crossinteraction parameters for Na + -CO 2 are equals to those of Cl --CO 2 , i.e., b Na-CO 2 = b Cl-CO 2 , and ε Na-CO 2 = ε Cl-CO 2 . This study assumes that ε Na-CO 2 = 0. Equation (A30) gives the temperature-dependent binary interaction parameter k 2,ion-CO 2 for the calculation of b NaCl-CO 2 .