Research Article  Open Access
A MultiScale GibbsHelmholtz Constrained Cubic Equation of State
Abstract
This paper presents a radically new approach to cubic equations of state (EOS) in which the GibbsHelmholtz equation is used to constrain the attraction or energy parameter, a. The resulting expressions for for pure components and for mixtures contain internal energy departure functions and completely avoid the need to use empirical expressions like the Soave alpha function. Our approach also provides a novel and thermodynamically rigorous mixing rule for . When the internal energy departure function is computed using Monte Carlo or molecular dynamics simulations as a function of current bulk phase conditions, the resulting EOS is a multiscale equation of state. The proposed new GibbsHelmholtz constrained (GHC) cubic equation of state is used to predict liquid densities at high pressure and validated using experimental data from literature. Numerical results clearly show that the GHC EOS provides fast and accurate computation of liquid densities at high pressure, which are needed in the determination of gas hydrate equilibria.
1. Introduction and Motivation
Carbon storage in deep ocean sediments is one of the many technologies being explored to mitigate carbon dioxide (CO_{2}) emissions from the consumption of carbonbased fuels (e.g., oil, coal, and even biofuels) and the associated greenhouse gas (GHG) effects that have impacted global climate change. Carbon storage has been identified by the National Academy of Engineering as one of fourteen grand challenge problems of the 21st century [1]. The feasibility of permanent deep ocean storage rests on two factors—neutral buoyancy (Brewer et al. [2], Dornan et al. [3]) and hydrate formation—and recent impact estimates indicate that the associated displacement of seawater in porous sediments by CO_{2} emissions from 1000 coalfired power plants over 100 years would cause only a 1millimeter rise in sea levels (Combs [4]).
Pumping carbon dioxide in deep ocean sediments is challenging. Ocean depth gives rise to a point of neutral buoyancy and geothermal heating in the sediments produces a second, deeper point of neutral buoyancy because the density of CO_{2} decreases more rapidly than pore seawater with depth due to heating. Thus there is a neutral buoyancy zone (NBZ). Additionally, pumping liquid CO_{2} into the sediments at temperatures near 275.15 K and pressures around 30 MPa produces locally supersaturated conditions of CO_{2} in seawater, creating an environment that favors the formation and sustainability of a liquid carbon dioxide and CO_{2} hydrates thermodynamically. Thus various combinations of a liquid CO_{2} phase, liquid seawater and one or more hydrate phases will coexist. There are also chemical reactions in seawater that produce carbonate and bicarbonate ions as well as the presence of sodium (Na) and chloride (Cl) ions in seawater; thus strong electrolyte solution behavior must be considered. To assess and predict the short and longterm impact of storing carbon dioxide in deep ocean sediments, accurate quantitative descriptions of the growth and possible dissolution of liquid CO_{2} and CO_{2} hydrate reservoirs as well as the interactions between these reservoirs and the ocean over multiple time and length scales are needed. This, in turn, requires a concerted effort of computer modeling, experimental work, and observational data to build understanding and to validate computer models.
The focus of this paper is on one aspect of these complex simulations—the rapid and accurate computation of liquid properties using equations of state (EOS) for modeling the growth and dissolution of CO_{2} liquid and gas hydrate reservoirs. This class of reservoir simulations can require hundreds of thousands of property, phase split, and phase equilibrium computations, and it has long been recognized that cubic equations of state generally do a poor job predicting accurate liquid densities at high pressure. Moreover, more complicated equations of state like the statistical associating fluid theory (SAFT) equation and its many variants require nested iteration for hydrogen bonding compounds and thus will not provide rapid enough computations for reservoir simulations.
The novel contributions of this paper are the following:(1)the development of a new thermodynamically rigorous framework for pure component parameters in cubic equations of state (EOS) based on using the GibbsHelmholtz equation as a constraint;(2)both closedform and integral multiscale expressions for the energy parameters in cubic EOS that directly incorporate molecular level information (specifically internal energy departure functions) obtained from Monte Carlo or molecular dynamics simulations;(3)the development of a new mixing rule based on the GibbsHelmholtz equation for the prediction of phase properties and phase equilibria for mixtures;(4)numerical testing and validation with experimental data from the open literature.
The significance of the modeling effort contained in this article is the development and validation of a radically new, fast, reliable, and truly predictive cubic EOS framework using a novel multiscale modeling approach that makes combined use of molecular simulations and classical thermodynamics. Specifically, the GibbsHelmholtz equation is used to constrain the energy parameter, a(T, ) in the van der Waals family of cubic EOS through the internal energy departure function. The proposed framework is very general, uses Monte Carlo simulations to evaluate the internal energy departure function, and can readily accommodate the usual molecular interactions for nonelectrolyte species (e.g., van der Waals forces) as well as electrostatic interactions (i.e., chargecharge, chargedipole, dipoledipole, chargequadrupole, etc.). Electrostatic interactions for weak and strong electrolyte systems (or other interactions) can be directly taken into account by using an appropriate potential energy function at the molecular level. Moreover, the proposed modeling framework, when coupled with excess Gibbs free energy () models for solid and hydrate phases, has the potential to rapidly, reliably, and accurately determine all types of phase equilibrium involving vapor (), liquid (), solid ice (), and gas hydrate () phases and be readily incorporated into a variety of (reservoir) simulation programs.
1.1. Preview
Figure 1 is reproduction of Figure 4 in Brewer et al. [2] who report experimental mass densities for liquid CO_{2} at seven distinct pressures and two temperatures (273.15 K and 283.15 K). The o’s and ’s in Figure 1 represent experimental data points while the curves shown by and    in Figure 1 are extrapolated from least squares fits of the experimental data. Also shown in Figure 1 is an estimate of the upper point of neutral buoyancy (i.e., where the mass density of CO_{2} and seawater are equal). Note that this occurs at 2600 m where and MPa and the mass density of seawater is 1.045 g/cm^{3}. For comparison, mass densities of liquid CO_{2} at experimental conditions were calculated in two different ways for a temperature of : () using the classical SoaveRedlichKwong (SRK) equation and () by the novel multiscale GibbsHelmholtz constrained (GHC) equation of state for pure components proposed in this work. The classical SRK equation (shown in green) does extremely poorly, predicting densities that are way too low (e.g., at a depth of 2600 m, the density of liquid CO_{2} predicted by SRK is 0.9627 g/cm^{3} (just outside the figure). On the other hand, the novel GibbsHelmholtz constrained equation for pure components that will be described in this paper predicts liquid CO_{2} densities (shown as the red curve) much closer to extrapolated experimental values (shown by   ). More specifically, at and 26.5 MPa, the proposed GibbsHelmholtz constrained equation of state predicts a liquid CO_{2} density of 1.0378 g/cm^{3}, which compares favorably with the extrapolated experimental value of 1.045 g/cm^{3}. Not surprisingly, the proposed approach also predicts an upper point of neutral buoyancy at approximately 2845 m (29.0 MPa), which is in good agreement with the experimental value of 2600 m reported by Brewer et al. [2]. Details of these computations are given in Section 4.
The key conclusions that can be drawn from the numerical comparisons shown in Figure 1 are the following.(1)The classical SRK equation does a very poor job of matching liquid carbon dioxide densities at high pressure.(2)The novel multiscale GibbsHelmholtz constrained (GHC) cubic equation of state for pure components proposed in this paper predicts a good approximate upper point of neutral buoyancy between pure liquid CO_{2} and seawater.
1.2. Current Needs
To be successful in modeling deep ocean sedimentary storage of CO_{2}, any equation of state must:(1)accurately model multiple liquid phases involving both nonelectrolyte (liquid CO_{2}) and strong electrolyte (seawater with dissolved CO_{2}) phase behavior,(2)accurately predict the mass density inversions that occur between seawater and liquid CO_{2} in the neutral buoyancy zone (i.e., both the upper and lower points of neutral buoyancy),(3)be computationally tractable and provide a high level of computational reliability since phase equilibrium computations need to be performed over many finite elements of any reservoir simulation.(4)be able to be coupled to models for ice and/or gas hydrate phases in order to successfully model liquidliquidhydrate (LLH) or other equilibrium.
2. Literature Survey
A brief survey of literature relevant to phase equilibrium in electrolyte and gas hydrate systems is presented.
2.1. Existing Models
Current models for multiphase equilibrium capable of modeling strong electrolyte systems include the segmentbased NRTL activity coefficient (SAC) methods developed by Chen and Song [5, 6], the predictive SoaveRedlichKwong (PSRK) equation of Kiepe et al. [7], the electrolyte SAFT equation (Galindo et al. [8], Bezhadi et al. [9]), and a variety of very specific models (e.g., Chen and Millero [10]). Each of these approaches has advantages and disadvantages.
2.2. EOSBased Models
The PSRK equation of Kiepe et al. [7] is a cubic equation of state that uses based mixing rule for the attraction parameter in the SRK equation coupled to the Liquid Functional Activity Coefficient (LIFAC) group contribution method of Yan et al. [11]. LIFAC accounts for the influence of short (van der Waals), middle (chargedipole), and long (chargecharge) effects of strong electrolytes on the activity coefficients of the nonelectrolyte species. That is, phase behavior for nonelectrolyte components is directly modeled using the SRK equation with a HuronVidal type or based mixing rule (Huron and Vidal [12]), where activity coefficients are first predicted by LIFAC and then modified by isoactivity constraints to obtain saltfree estimates of activity coefficients for use in mixing rules. Also because of the presence of ions, a large number of parameters are needed to capture the interactions between gases and ionic groups in the PSRK electrolyte approach. Many results for gas solubility, Henry’s law constants, and pressurecomposition diagrams are presented and compared to experimental data, and while the work of Kiepe et al. [7]. does address strong electrolyte systems, no studies of liquidliquid equilibrium at conditions relevant to CO_{2} storage are reported. Additionally, recent work by Yoon et al. [13] shows that the PSRK model (not the PSRK electrolyte model) can be used in conjunction with van der WaalsPlatteeuw theory to predict various equilibria between vapor, liquid, ice, and hydrates. However, while Yoon and Yamamoto report liquidliquidhydrate equilibrium for cyclopropane hydrates, they only consider nonelectrolyte systems with water.
The Statistical Associating Fluid Theory (SAFT) equation has also been extended to electrolyte systems. All approaches thus far use a mean spherical approximation (MSA) for longrange Coulombic interactions and calculate solvention interactions using a potential functions (e.g., squarewell, Yukawa potentials). See , for example, Galindo et al. [8]. and Bezhadi et al. [9]. who use the SAFTVariable Range (SAFTVR) equation to study vaporliquid equilibrium in strong electrolyte mixtures. While SAFT has proven to be a very popular approach, the iterative computation of densities has embedded in it another level of iteration for finding mole fractions of unbonded sites in mixtures that exhibit hydrogen bonding (i.e., strong electrolytes). This coupled with the fact that SAFT tends to yield highorder transcendental density functionality makes it difficult to guarantee that all density roots can be calculated reliably.
2.3. Activity Coefficient Models
The recently developed segmentbased NRTL activity coefficient (SAC) models of Chen and Song [5, 6] have been tested on a wide array of organic electrolyte mixtures common to the pharmaceutical industry and very good results have been obtained. However, no hydrate forming gases are included in the SAC solute data, and because this approach is an activity coefficient approach, it is more suitable for lowpressure applications. Thus, in our opinion, the SAC approach is not applicable to deep ocean sedimentary CO_{2} storage.
2.4. Perspective
The complex functionality and nested iteration structure for hydrogen bonding systems makes SAFT unsuitable for use in reservoir simulation where fast and reliable phase equilibrium computations are required. On the other hand, the highpressure environment of carbon storage and the lack of model parameters for light gases eliminate the use of the segmentbased activity coefficient (SAC) approach. This leaves cubic equations of state from the van der Waals family, which we believe represent a reasonable starting point for the development of an EOS for use in reservoir simulation. Cubic EOS are a reasonable approach because it is straightforward to guarantee that all density roots at specified conditions of temperature, pressure, and composition can be found using a reliable equation solver coupled with polynomial deflation. However, while this makes them computationally tractable and provides a high level of computational reliability, the main difficulty in using cubic EOS has always been their inability to accurately predict liquid density and their limited development for strong electrolyte solutions. Inaccurate density calculations have a domino effect since liquid densities are used to compute fugacity coefficients, which in turn are used to compute other thermodynamic properties . Additionally, the PSRK approach requires a very large number of parameters and has not been tested on highpressure liquidliquid equilibrium. Therefore one of the main focuses of the computer modeling aspects of this paper is to improve the ability of cubic EOS to predict accurate liquid densities for nonelectrolyte and electrolyte solutions. In this paper, we focus specifically on nonelectrolyte systems but do consider systems with electrostatic effects.
3. A New Multiscale Modeling Approach to PVT Behavior of Fluids
To improve the ability of cubic EOS to predict accurate liquid densities for nonelectrolyte, we have developed a new multiscale methodology for cubic equations of state by using the GibbsHelmholtz equations to constrain the attraction parameter, . The resulting multiscale modeling framework(1)is thermodynamically rigorous,(2)easily allows molecular level information to be included,(3)has been developed for pure components and mixtures,
3.1. Pure Components
Any thermodynamic departure function is defined as the difference between a real property and an ideal gas property at the same temperature and pressure. In the material that follows, we denote departure functions using the superscript D. From the differential form of the Gibbs free energy and the definition of fugacity, , it follows that and that where and denote the Gibbs free energy of departure and ideal gas Gibbs free energy respectively, evaluated at the same temperature and pressure , is a fugacity coefficient, and is the gas constant. From the GibbsHelmholtz equation, we have that Therefore . Using the Soave form of the RedlichKwong (RK) equation for a pure component given by and setting , where is the molar volume of the solid phase, we have that At high pressure, is a reasonable assumption. To see this consider any pair of isotherms in the compressed liquid region of a generic diagram like the one shown in Figure 2. Note that the isotherms are packed closely together and therefore is small. Thus we let which gives . Subsequent differentiation of with respect to temperature at constant pressure yields Use of the expressions: () , () , and () in the expression for and some algebraic rearrangement gives Equation (5) is a firstorder, inhomogeneous differential equation that defines at fixed pressure and is easily solved. The general solution to the homogeneous form of (5) is given by . The particular solution to (5) can be calculated using variation of parameters or simple onesided Green’s functions. To do this, we assume that there is a solution of the form , where is a weighting parameter that must be determined. Straightforward application of variation of parameters applied to (5) gives the condition , which yields Thus . Differentiation of and substitution of and easily show that this form of the particular solution satisfies (5). Substituting the expression for from (5) gives Therefore the complete solution to (5) is where the constant can be evaluated by applying the boundary condition at . That is, The foregoing analysis yields the following thermodynamically constrained integral expression for
Observations
Approximating , the excluded volume, by , where is the solid molar volume is reasonable since the solid phase represents a “close packed” phase. The standard expression often overestimates the excluded volume. Moreover, solid density data, , for many pure components are readily available.
Equation (10) involves . The effect of pressure on is incorporated implicitly through the evaluation of the internal energy departure function and thus (10) is really an expression for .
Because was developed using the GibbsHelmholtz equation as a constraint, it is thermodynamically rigorous. We do not require the acentric factor or empirical expressions for such as the alpha function given by Soave [14]
Other similar empirical expressions (e.g., Mathias [15], Mathias and Copeman [16], Melham [17], Twu et al. [18]) are not required either. Also, the boundary condition used to evaluate the constant in (8) is the critical point and often represents an experimentally validated condition for most pure compounds.
can be readily determined using isothermalisobaric (NTP) Monte Carlo (MC) or molecular dynamics (MD) simulations. While this requires an empirical potential energy model for the component(s) of interest, it has the flexibility to readily accommodate a variety of interactions that take place between particles (van der Waals forces, Coulombic forces, etc.). Moreover, any temperature and pressure dependence of can be automatically included in the molecular simulations. We use NTP Monte Carlo simulations to measure .
Note that the expression for involves V and therefore (10) is coupled to the EOS. To decouple the equation of state from (10) we invoke the highpressure limit that . Under this assumption, , in which case (10) becomes
The assumption is only approximate at finite pressure but, as we show, only introduces a small error in calculated liquid density.
If one further assumes that temperature functionality of is a weak enough to not affect the integrals in (12), then the following closedform solution for results
Note that the functionality of involves a linear term in temperature and a term. Also note that (13) still allows to be determined from molecular simulations and thus still includes implicit pressure effects on . Equation (13) represents a radically new expression for that relates energies measured at the molecular scale, , to the energy parameter in cubic equations of state. There is nothing even close to (13) in the open literature!.
3.2. Mixtures
Modeling and computations involving mixtures are always handled using mixing rules. Our analysis for mixtures begins with the fundamental thermodynamic relationship: where is the Gibbs free energy departure function of the mixture, ϕ_{M} is the fugacity coefficient for the mixture, denotes the mole fraction of the th component in the mixture, and is the fugacity coefficient for pure component . Differentiation and application of the GibbsHelmholtz equation for mixtures give where is the enthalpy departure function of the mixture.
To develop a closed form mixing rule for we do the following.(1)Differentiate with respect to at constant .(2)Let , where each pure component excluded volume , where is the density of pure solid component .(3)Differentiate with respect to at constant and .(4)Use the approximations and , which imply and , respectively.
These assumptions lead to the differential equation where is the internal energy departure function of the mixture, and where the functions and their first derivatives are known and given by (13) applied to pure component . Moreover, assumptions () and (5) show that , and since must be computed to calculate , given and , the entire right hand side of (16) is well defined.
Equation (16) can be solved in exactly the same way that (5) was solved. This yields the following solution or mixing rule for Note that there are three terms on the righthand side of (17). The last two terms in the mixing rule are quadratic in composition (since they involve the product of ), the last term explicitly involves the pressure, and the second term involves the quantity , which accounts for temperature, pressure and composition dependence in an implicit manner. Note also that there is no need to use empirical relationships such as the geometric combining rule to define the coefficients of the cross composition terms (i.e., the terms). Rather all terms including the cross composition terms in the proposed mixing rule are defined in terms of physically meaning quantities (i.e., ’s, ’s, ’s, etc.). Moreover, the cross composition terms are not assumed to be symmetric.
All that remains is to establish a boundary condition for to evaluate the constant in (17). To do this we use Kay’s rules [19] to calculate mixture critical properties and and assume that for any mixture of composition . This gives and the following expression for
Observations
The proposed mixing rule for is thermodynamically rigorous since it has been derived using the GibbsHelmholtz equation as a constraint for both the mixture and for each pure component.
has the same temperature dependence for mixtures as for pure components. It involves a linear term in as well as dependence.
There is explicit pressure dependence in the mixing rule through the term .
This new mixing rule contains a number of terms that are quadratic in composition as well as internal energy departure functions for the mixture and pure components. The departure functions account for temperature, pressure, and composition dependence of in an implicit manner.
The proposed mixing rule is a multiscale mixing rule since it can easily incorporate and from molecular simulations.
4. Numerical Results and Comparisons of Modeling and Experiments
In this section, numerical results for liquid densities using the new proposed framework are presented and compared with experimental data. Specifically, liquid densities for pure carbon dioxide, pure water, and CO_{2}water mixtures are calculated using (13) and (19) and compared to high pressure experimental density data. All bulk phase property computations and molecular simulations were performed in double precision arithmetic on a Dell High Precision 670 workstation using the LaheyFijitsu (LF95) compiler.
4.1. Numerical Results in Figure 1
The results for the classical SRK equation shown in Figure 1 were computed using , bar, and , which yield cm^{3}/mol and , where is the Soave alpha function. The numerical results for the GibbsHelmholtz constrained cubic EOS shown in Figure 1 were determined using cm^{3}/mol and (13) with cm^{3}bar/mol, which was determined by performing a set of 25 NTP Monte Carlo simulations at low pressure (). The potential energy model for CO_{2} consisted of LennardJones forces and electrostatic forces. Table 1 shows the calculated densities for both SRK and the GH constrained (GHC) EOS.

The numerical values shown in Table 1 are those that have been plotted in Figure 1.
4.2. Other Validated Liquid CO_{2} Density Predictions
Table 2 shows the liquid CO_{2} experimental data of Magee and Ely [20].
 
*Reproduced from Magee and Ely [20, runs 900 through 1100 in Table V, page 1178–1179]. 
It is important to note that runs 900 through 1100 are the only runs that measure sample mass using a gravimetric method for fixed experimental apparatus volume; all other runs determine molar density using a 32term extended BenedictRubinWebb equation. See Magee and Ely [20, page 1167].
Figure 3 shows a comparison of the liquid molar densities calculated using the SRK equation and the proposed multiscale GibbsHelmholtz constrained equation of state with the experimental data of Magee and Ely [20]. Note that the proposed multiscale method using the expression for given by (13) does an exceptionally good job of matching experimental liquid CO_{2} density data over a wide range of temperatures and pressures while the Soave expression of consistently underestimates the liquid molar density.
4.3. Validated Numerical Results for Densities of Compressed Water
We have calculated liquid densities for water at pressures between 5.30 and 1026.39 bars at various temperatures using the SRK and the proposed GHC equation and compared those numerical results to the experimental data of Kell and Whalley [21]. Since water is a polar compound, we have included the Peneloux [22] volume translation given by the expression in the computations using the SRK equation, where is the volume predicted by the traditional SRK equation and is the SRK predicted volume modified by the Peneloux volume translation. The critical properties for water used in these calculations were , , , and . Note that the value of calculated from critical properties is 21.082 cm^{3}/mol and represents a very poor estimate of the excluded volume.
For the GHC equation, the solid volume, , of very high density amorphous ice was used to estimate the excluded volume. This value, cm^{3}/mol, gives a much better estimate of the excluded volume. The calculated values of ranged from to cm^{3}bar/mol. Comparisons of specific volumes at 273.15 K are shown in Table 3 and Figure 4 in cm^{3}/g, which is the form in which Kell and Whalley [21] reported their experimental results.
 
*Experimental data reproduced from Kell and Whalley [21, Table $4$, page 586]. ${V}^{\prime \text{SRK}}$ includes Peneloux volume translation $({V}^{\prime \text{SRK}}={V}^{\text{SRK}}0.40768((R{T}_{c})/{p}_{c})[0.29441({p}_{c}{V}_{c})/R{T}_{c}])$. 
From Table 3, it is easily seen that the SRK equation with the Peneloux [22] volume translation does poorly at matching the experimental data. However, this is not surprising since the Peneloux volume translation term, , is an empirical correction. For water, the volume translation has a value of 0.354111 cm^{3}/g or 6.374 cm^{3}/mol, overcompensates for polarity, and produces specific volumes that are roughly 1 cm^{3}/g smaller than they should be. On the other hand, without the Peneloux volume translation the specific volumes calculated by the SRK equation are extremely poor ( cm^{3}/g). In contrast, the GHC equation does an excellent job of matching experimental specific volumes. Moreover, we believe that the excellent match provided by the GHC equation is due to the fact that it makes combined use of molecular scale information and a rigorous classical relationship to determine .
4.4. Validated Density Calculations for CO_{2}H_{2}O Mixtures
To validate the mixing rule given by (19) we compared numerical density calculations for CO_{2}H_{2}O mixtures with the highpressure experimental data of Teng et al. [23] at 278 K. For the SRK calculations we included the Peneloux volume translation for water. For the GHC equation, we used one fixed set of values of pure component internal energy departure functions: cm^{3}bar/mol, cm^{3}bar/mol, and cm^{3}bar/mol in order to illustrate that repeated and costly molecular simulations are not required for good phase property predictions. The comparison is shown in Table 4.
 
*
Experimental data reproduced from Teng et al. [23, Table $2$, page 1306].${\rho}^{\prime \text{SRK}}$ includes the Peneloux correction for water. 
Again, the GHC equation clearly outperforms the SRK equation with the Peneloux translation.
4.5. Reliability and Computational Speed
At the equation of state level, the liquid density calculations presented in Section 4 are 100 reliable when the global terrain methodology of Lucia and Feng [24] and polynomial deflation are used. Moreover, phase density computations using the terrain approach for either pure components or mixtures require roughly 0.01 s to find all three roots. At the molecular scale, NTP Monte Carlo calculations require substantial computational resources. If MC simulations are to be performed each time the temperature, pressure, and/or composition change, then significant computational costs may result. On the other hand, if these molecular scale calculations are performed at nominal conditions of temperature, pressure and composition and interpolated or extrapolated, then considerable reductions in computational overhead will result with little error at the bulk phase length scale.
5. Conclusions
A radically new cubic equation of state approach was presented in which the excluded volume parameter, , was approximated using solid molar volumes and new expressions for the energy parameter, , were derived by using the GibbsHelmholtz equation to constrain the value of . The resulting expressions for and a necessarily (1)are thermodynamically rigorous,(2)avoid the need for empirical correlations like the alpha function of Soave,(3)have temperature functionality with a linear term, and a term. (4)Involve internal energy departure functions, and thus make the new approach a multiscale cubic equation of state approach.
Moreover, the mixing rule or expression for has a pressure explicit term clearly showing that the energy parameter is pressure dependent, albeit weak. The resulting new cubic equation approach is called the GibbsHelmholtz constrained (GHC) equation and is truly predictive.
The GHC equation was compared to the classical SRK equation and validated using experimental highpressure liquid density data for CO_{2}, water, and CO_{2}H_{2}O mixtures. Numerical testing clearly shows that the GHC equation clearly uniformly outperforms the SRK equation and matches experimental liquid density data exceptionally well.
Future work with the GHC equation will involve further numerical testing and experimental validation of liquid densities for a wide variety of mixtures as well as predictions of phase equilibrium for systems involving multiple liquid phases and gas hydrates.
References
 http://www.engineeringchallenges.org/.
 P. G. Brewer, G. Friederich, E. T. Peltzer, and F. M. Orr Jr., “Direct experiments on the ocean disposal of fossil fuel ${\text{CO}}_{2}$,” Science, vol. 284, no. 5416, pp. 943–945, 1999. View at: Publisher Site  Google Scholar
 P. Dornan, S. Alavi, and T. K. Woo, “Free energies of carbon dioxide sequestration and methane recovery in clathrate hydrates,” Journal of Chemical Physics, vol. 127, no. 12, Article ID 124510, 2007. View at: Publisher Site  Google Scholar
 A. Combs, “An ocean trap for carbon dioxide,” MIT Technology Review, May 2009. View at: Google Scholar
 C.C. Chen and Y. Song, “Solubility modeling with a nonrandom two liquid segment activity coefficient model,” Industrial & Engineering Chemistry Research, vol. 43, pp. 8354–8362, 2004. View at: Google Scholar
 C.C. Chen and Y. Song, “Extension of nonrandom twoliquid segment activity coefficient model for electrolytes,” Industrial & Engineering Chemistry Research, vol. 44, no. 23, pp. 8909–8921, 2005. View at: Publisher Site  Google Scholar
 J. Kiepe, S. Horstmann, K. Fischer, and J. Gmehling, “Application of the PSRK model for systems containing strong electrolytes,” Industrial & Engineering Chemistry Research, vol. 43, no. 20, pp. 6607–6615, 2004. View at: Google Scholar
 A. Galindo, A. GilVillegas, P. J. Whitehead, G. Jackson, and A. N. Burgess, “SAFTVRE: phase behavior of electrolyte solutions with the statistical associating fluid theory for potentials of variable range,” Journal of Physical Chemistry B, vol. 103, no. 46, pp. 10272–10281, 1999. View at: Google Scholar
 B. Behzadi, B. H. Patel, A. Galindo, and C. Ghotbi, “Modeling electrolyte solutions with the SAFTVR equation using Yukawa potentials and the meanspherical approximation,” Fluid Phase Equilibria, vol. 236, no. 12, pp. 241–255, 2005. View at: Publisher Site  Google Scholar
 C.T. Chen and F. J. Millero, “Precise equation of state for seawater in oceanic ranges of salinity, temperature and pressure,” Deep Sea Research, vol. 24, pp. 365–369, 1977. View at: Google Scholar
 W. Yan, M. Topphoff, C. Rose, and J. Gmehling, “Prediction of vaporliquid equilibria in mixedsolvent electrolyte systems using the group contribution concept,” Fluid Phase Equilibria, vol. 162, no. 12, pp. 97–113, 1999. View at: Publisher Site  Google Scholar
 M.J. Huron and J. Vidal, “New mixing rules in simple equations of state for representing vapourliquid equilibria of strongly nonideal mixtures,” Fluid Phase Equilibria, vol. 3, no. 4, pp. 255–271, 1979. View at: Google Scholar
 J.H. Yoon, Y. Yamamoto, T. Komai, and T. Kawamura, “PSRK method for gas hydrate equilibria: I. simple and mixed hydrates,” AIChE Journal, vol. 50, no. 1, pp. 203–214, 2004. View at: Publisher Site  Google Scholar
 G. Soave, “Equilibrium constants from a modified RedlichKwong equation of state,” Chemical Engineering Science, vol. 27, no. 6, pp. 1197–1203, 1972. View at: Google Scholar
 P. M. Mathias, “A versatile phase equilibrium equation of state,” Industrial and Engineering Chemistry Process Design and Development, vol. 22, no. 3, pp. 385–391, 1983. View at: Google Scholar
 P. M. Mathias and T. W. Copeman, “Extension of the PengRobinson equation of state for polar fluids and fluid mixtures,” Fluid Phase Equilibria, vol. 13, p. 91, 1983. View at: Google Scholar
 G. A. Melhem, “A modified PengRobinson equation of state,” Fluid Phase Equilibria, vol. 47, p. 189, 1989. View at: Google Scholar
 C. H. Twu, D. Bluck, J. R. Cunningham, and J. E. Coon, “A cubic equation of state with a new alpha function and a new mixing rule,” Fluid Phase Equilibria, vol. 69, pp. 33–50, 1991. View at: Google Scholar
 W. B. Kay, “Density of hydrocarbon gases and vapors at high pressure and temperature,” Industrial & Engineering Chemistry Research, vol. 28, p. 1014, 1936. View at: Google Scholar
 J. W. Magee and J. F. Ely, “Specific heats (Cv) of saturated and compressed liquid and vapor carbon dioxide,” International Journal of Thermophysics, vol. 7, no. 6, pp. 1163–1182, 1986. View at: Publisher Site  Google Scholar
 G. S. Kell and E. Whalley, “The PVT properties of water: I. Liquid water in the temperature range 0 to 150 degrees C and at pressures up to 1 kb,” Philosophical Transactions of the Royal Society A, vol. 258, no. 1094, pp. 565–614, 1965. View at: Google Scholar
 A. Peneloux, E. Rauzy, and R. Freze, “A consistent correction for RedlichKwongSoave volumes,” Fluid Phase Equilibria, vol. 8, no. 1, pp. 7–23, 1982. View at: Google Scholar
 H. Teng, A. Yamasaki, M.K. Chun, and H. Lee, “Solubility of liquid ${\text{CO}}_{2}$ in water at temperatures from 278 K to 293 K and pressures from 6.44 MPa to 29.49 MPa and densities of the corresponding aqueous solutions,” Journal of Chemical Thermodynamics, vol. 29, no. 11, pp. 1301–1310, 1997. View at: Google Scholar
 A. Lucia and Y. Feng, “Global terrain methods,” Computers & Chemical Engineering, vol. 26, no. 45, pp. 529–546, 2002. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2010 Angelo Lucia. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.