This paper presents a radically new approach to cubic equations of state (EOS) in which the Gibbs-Helmholtz 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 Gibbs-Helmholtz 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 (CO2) emissions from the consumption of carbon-based 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 CO2 emissions from 1000 coal-fired power plants over 100 years would cause only a 1-millimeter 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 CO2 decreases more rapidly than pore seawater with depth due to heating. Thus there is a neutral buoyancy zone (NBZ). Additionally, pumping liquid CO2 into the sediments at temperatures near 275.15 K and pressures around 30 MPa produces locally supersaturated conditions of CO2 in seawater, creating an environment that favors the formation and sustainability of a liquid carbon dioxide and CO2 hydrates thermodynamically. Thus various combinations of a liquid CO2 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 long-term impact of storing carbon dioxide in deep ocean sediments, accurate quantitative descriptions of the growth and possible dissolution of liquid CO2 and CO2 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 CO2 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 Gibbs-Helmholtz equation as a constraint;(2)both closed-form 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 Gibbs-Helmholtz 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 Gibbs-Helmholtz 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., charge-charge, charge-dipole, dipole-dipole, charge-quadrupole, 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 CO2 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 CO2 and seawater are equal). Note that this occurs at 2600 m where and  MPa and the mass density of seawater is 1.045 g/cm3. For comparison, mass densities of liquid CO2 at experimental conditions were calculated in two different ways for a temperature of  : () using the classical Soave-Redlich-Kwong (SRK) equation and () by the novel multiscale Gibbs-Helmholtz 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 CO2 predicted by SRK is 0.9627 g/cm3 (just outside the figure). On the other hand, the novel Gibbs-Helmholtz constrained equation for pure components that will be described in this paper predicts liquid CO2 densities (shown as the red curve) much closer to extrapolated experimental values (shown by - - -). More specifically, at and 26.5 MPa, the proposed Gibbs-Helmholtz constrained equation of state predicts a liquid CO2 density of 1.0378 g/cm3, which compares favorably with the extrapolated experimental value of 1.045 g/cm3. 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 Gibbs-Helmholtz 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 CO2 and seawater.

1.2. Current Needs

To be successful in modeling deep ocean sedimentary storage of CO2, any equation of state must:(1)accurately model multiple liquid phases involving both nonelectrolyte (liquid CO2) and strong electrolyte (seawater with dissolved CO2) phase behavior,(2)accurately predict the mass density inversions that occur between seawater and liquid CO2 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 liquid-liquid-hydrate (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 segment-based NRTL activity coefficient (SAC) methods developed by Chen and Song [5, 6], the predictive Soave-Redlich-Kwong (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. EOS-Based 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 (charge-dipole), and long (charge-charge) 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 Huron-Vidal type or -based mixing rule (Huron and Vidal [12]), where activity coefficients are first predicted by LIFAC and then modified by iso-activity constraints to obtain salt-free 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 pressure-composition 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 liquid-liquid equilibrium at conditions relevant to CO2 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 Waals-Platteeuw theory to predict various equilibria between vapor, liquid, ice, and hydrates. However, while Yoon and Yamamoto report liquid-liquid-hydrate 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 long-range Coulombic interactions and calculate solvent-ion interactions using a potential functions (e.g., square-well, Yukawa potentials). See , for example, Galindo et al. [8]. and Bezhadi et al. [9]. who use the SAFT-Variable Range (SAFT-VR) equation to study vapor-liquid 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 high-order transcendental density functionality makes it difficult to guarantee that all density roots can be calculated reliably.

2.3. Activity Coefficient Models

The recently developed segment-based 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 low-pressure applications. Thus, in our opinion, the SAC approach is not applicable to deep ocean sedimentary CO2 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 high-pressure environment of carbon storage and the lack of model parameters for light gases eliminate the use of the segment-based 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 high-pressure liquid-liquid 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 Gibbs-Helmholtz 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 Gibbs-Helmholtz equation, we have that Therefore . Using the Soave form of the Redlich-Kwong (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 first-order, 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 one-sided 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

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 Gibbs-Helmholtz 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 isothermal-isobaric (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 high-pressure 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 closed-form 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 Gibbs-Helmholtz 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 right-hand 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

The proposed mixing rule for is thermodynamically rigorous since it has been derived using the Gibbs-Helmholtz 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 CO2-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 Lahey-Fijitsu (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  cm3/mol and , where is the Soave alpha function. The numerical results for the Gibbs-Helmholtz constrained cubic EOS shown in Figure 1 were determined using  cm3/mol and (13) with  cm3bar/mol, which was determined by performing a set of 25 NTP Monte Carlo simulations at low pressure (). The potential energy model for CO2 consisted of Lennard-Jones forces and electrostatic forces. Table 1 shows the calculated densities for both SRK and the G-H constrained (GHC) EOS.

The numerical values shown in Table 1 are those that have been plotted in Figure 1.

4.2. Other Validated Liquid CO2 Density Predictions

Table 2 shows the liquid CO2 experimental data of Magee and Ely [20].

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 32-term extended Benedict-Rubin-Webb 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 Gibbs-Helmholtz 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 CO2 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 cm3/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,  cm3/mol, gives a much better estimate of the excluded volume. The calculated values of ranged from to  cm3bar/mol. Comparisons of specific volumes at 273.15 K are shown in Table 3 and Figure 4 in cm3/g, which is the form in which Kell and Whalley [21] reported their experimental results.

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 cm3/g or 6.374 cm3/mol, overcompensates for polarity, and produces specific volumes that are roughly 1 cm3/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 ( cm3/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 CO2-H2O Mixtures

To validate the mixing rule given by (19) we compared numerical density calculations for CO2-H2O mixtures with the high-pressure 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:  cm3bar/mol,  cm3bar/mol, and  cm3bar/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.

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 Gibbs-Helmholtz 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 Gibbs-Helmholtz constrained (GHC) equation and is truly predictive.

The GHC equation was compared to the classical SRK equation and validated using experimental high-pressure liquid density data for CO2, water, and CO2-H2O 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.