Modeling gas hydrate equilibria in electrolyte solutions

Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259 www.elsevier.com/locate/calphad

Modeling gas hydrate equilibria in electrolyte solutions Giles M. Marion a,∗ , David C. Catling b,c , Jeffrey S. Kargel d a Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512, USA b University of Washington, Seattle, Washington 98195, USA c University of Bristol, Bristol BS8 1RJ, United Kingdom d University of Arizona, Tucson, AZ 85721, USA

Received 31 January 2006; received in revised form 11 April 2006; accepted 22 April 2006 Available online 30 May 2006

Abstract Virtually all gas hydrate models use a statistical thermodynamic approach to describe solid phase gas hydrates. An alternative to this model is a “solid solution” model based on classical thermodynamics. Various models are used to describe associated gas phase and aqueous electrolyte properties. None of the electrolyte models, however, are state-of-the-art with respect to gas/electrolyte interactions. The objectives of this work were to (1) develop a classical thermodynamic approach for gas hydrate equilibria, (2) incorporate these hydrate chemistries into a state-of-the-art electrolyte model, and (3) validate the gas hydrate/electrolyte model. The model used the Pitzer approach to estimate as functions of temperature and pressure: (1) activity coefficients for cations, anions, and soluble gases, and (2) the activity of water. This model explicitly considers interactions between soluble gases and electrolytes. The Duan model was used to characterize the gas phase. Solubility products were calculated for CH4 ·6H2 O (T = 180–298 K) and CO2 ·6H2 O (T = 180–283 K). A “solid solution” model was used to define mixed CH4 –CO2 gas hydrates. Model fits to electrolyte/gas hydrate experimental data were an improvement over previously published models. For example, average errors for four salt–CO2 ·6H2 O systems were reduced from 7.2(±4.1)% to 3.0(±1.8)%. Similarly, average errors for 13 salt–CH4 ·6H2 O systems were reduced from 4.3(±1.1)% to 3.0(±0.6)%. The results of this work demonstrate the critical importance of the electrolyte component in gas hydrate models. c 2006 Elsevier Ltd. All rights reserved.

Keywords: Gas hydrate; Electrolyte; Pitzer approach; Modeling

1. Introduction Virtually all gas hydrate models use the statistical thermodynamic model of van der Waals and Platteeuw [1] as implemented by Parrish and Prausnitz [2] to describe solid phase gas hydrates [3–12]. An alternative to the statistical thermodynamic model for solid phase gas hydrates is a “solid solution” model based on classical thermodynamics [13]. Various models are used to describe associated gas phase and

∗ Corresponding address: Desert Research Institute, Division of Earth and Ecosystem Sciences, 2215 Raggio Parkway, Reno, NV 89512, USA. Tel.: +1 775 673 7349; fax: +1 775 673 7485. E-mail addresses: [email protected] (G.M. Marion), [email protected] (D.C. Catling), [email protected] (J.S. Kargel).

c 2006 Elsevier Ltd. All rights reserved. 0364-5916/$ - see front matter doi:10.1016/j.calphad.2006.04.002

aqueous phase properties. In a few gas hydrate models, the Pitzer approach is used to describe the activity of water in electrolyte solutions [5–8]. While the Pitzer approach is based solely on binary and ternary interaction parameters, this model extends easily to more complex solutions such as seawater, where there are seven independent components (Na+K+Mg+ Ca+Cl+SO4 +HCO3 +H2 O−1). However, none of the existing gas hydrate/electrolyte models are state-of-the-art with respect to gas/electrolyte interactions. The objectives of this work were to (1) develop a classical thermodynamic approach for gas hydrate equilibria, (2) incorporate these gas (carbon dioxide and methane) hydrate chemistries into a state-of-the-art electrolyte model parameterized for low temperatures and high pressures (the FREZCHEM model), and (3) validate the resulting gas hydrate/electrolyte model.

249

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

2. Theory 2.1. FREZCHEM model The FREZCHEM model is an equilibrium chemical thermodynamic model parameterized for concentrated electrolyte solutions using the Pitzer equations [14,15] for the temperature range from < −70 to 25 ◦ C and the pressure range from 1 to 1000 bar [16–22]. The model is currently parameterized for the Na–K–Mg–Ca–Fe–H–Cl–SO4 –NO3 –OH–HCO3 –CO3 –CO2 – O2 –H2 O system and includes 56 solid phases including ice, 11 chloride minerals, 14 sulfate minerals, 15 carbonate minerals, five solid-phase acids, three nitrate minerals, six acidsalts, and one iron oxide. An objective of this work was to develop a gas hydrate model based on classical chemical thermodynamic principles that can be incorporated seamlessly into the FREZCHEM model. A FORTRAN version of the resulting gas hydrate model is available from the senior author ([email protected]). 2.2. Pitzer approach In the Pitzer approach, the activity coefficients (γ ) as a function of temperature at 1.01 bar pressure for cations (M), anions (X), and neutral aqueous species (N), such as CO2 (aq) or CH4 (aq), are given by X 2 ln(γM ) = z M F+ m a (2BMa + Z CMa )   X X m a ΨMca + m c 2ΦMc + XX XX m c m a Cca + m a m a 0 ΨMaa 0 + z M X XX +2 m n λnM + m n m a ξnMa (1) X 2 m c (2BcX + Z CcX ) ln (γX ) = z X F +   X X m c ΨcXa + m a 2ΦXa + XX XX + m c m c0 Ψcc0 X + |z X | m c m a Cca X XX +2 m n λnX + m n m c ξncX (2) X X ln(γ N ) = m c (2λNc ) + m a (2λNa ) XX + m c m a ξNca (3) where B, C, Φ, Ψ , λ, and ξ are Pitzer-equation interaction parameters, m i is the molal concentration, and F and Z are equation functions. The activity of water (aw ) at 1.01 bar pressure is given by P   −φ m i (4) aw = exp 55.50844 where φ is the osmotic coefficient, which is given by: ( 3 XX
−Aφ I 2 2 φ + m c m a Bca + Z Cca (φ − 1) = P 1 m i 1 + bI 2   XX X φ + m c m c0 Φcc0 + m a Ψcc0 a   XX X φ + m a m a 0 Φaa 0 + m c Ψcaa 0

+

XX

+

XXX

m n m c λnc +

XX

m n m a λna )

m n m c m a ξn,c,a .

(5)

Note that Eqs. (1)–(3) and (5) include interaction terms for neutral species (n, e.g. CO2(aq) or CH4(aq) ) with cations (c) and anions (a), a factor missing from previous use of the Pitzer approach in gas hydrate models [5,7]. See Pitzer [14,15] or Marion and Farren [19] for a complete description of these equations that govern the temperature dependence of solution thermodynamic properties at 1.01 bar pressure. An important factor in gas hydrate models is that equilibria for all components and processes must be specified as functions of both temperature (Eqs. (1)–(5)) and pressure. In this work, we specified the pressure dependence of equilibrium constants (K ), activity coefficients (γ ), and the activity of water (aw ) as follows: The pressure dependence of equilibrium constants was estimated by

2 !  P  0 0 −1V r P − P 0 1K r P − P 0 K = ln + (6) 0 RT 2RT KP where 0

X

0 V i + nV H2 O − VMX(cr)

0

X

0 K i + n K H2 O − K MX(cr)

1V r = 1K r =

0

0

0

0

(7) (8)

and K P is the equilibrium constant at pressure P (bar), K P0 is the equilibrium constant at standard pressure P 0 (1.01325 bar), R is the gas constant (83.1451 cm3 bar mol−1 deg−1 , 0 [15]), T is temperature (K), V i is the partial molar volume 0

at infinite dilution of the ith constituent, and K i is the molar compressibility at infinite dilution of the ith constituent [23]. 0 0 The values of V i and K i used in this paper are compiled in Marion et al. [20]. The pressure dependence of ion activity coefficients was estimated with the equation: 
!  0 V i − V i P − P0 γiP = . (9) ln 0 RT γiP In this equation, the partial molar volume for a cation (M) or anion (X) is given by: h X  i X 0 2 v v V M = V M + zM f + 2RT m a BMa + m c z c CMa h i XX
2 v 0 v m c m a zM Bca (10) + RT + z M Cca or h X  i X 0 2 v v V X = V X + zX f + 2RT m c BcX + m c z c CcX h i XX
2 v 0 v + RT m c m a zX Bca . (11) + |z X |Cca The equations for the volumetric Pitzer-equation parameters: f , 0 B v , B v , and C v are compiled in Marion et al. [20].

250

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

The pressure dependence of the activity of water is calculated by: ! ! 0 (V w − V w )(P − P 0 ) awP = ln . (12) 0 RT awP In this equation, the molar volume of water is calculated with P   Mw 1000 + m i Mi X Vw = − mi V i (13) 1000 ρ where Mw and Mi are the molecular masses of water and solution constituents, respectively. The solution density (ρ) in turn is calculated by P 1000 + m i Mi ρ= (14) P 0 1000 ex + m V + V i i ρ0 where m i is the molal concentration, Mi is the molecular mass, ρ 0 is the density of pure water at a given temperature and 0 pressure, V i is the partial molal volume at infinite dilution, and V ex is the excess volume of mixing given by XX  v  V ex v = fv +2 m c m a Bca + (Σ m c z c ) Cca . RT

(15) 0

Equations for the Pitzer-equation parameters ( f v , B v , C v ), V w , V i (Eqs. (12)–(15)), and ρ 0 are compiled in Marion et al. [20].

aqueous concentrations, in defining the solubility product is that is negates the necessity of having to extrapolate Henry’s law constants into the subzero temperature range to estimate CO2 (aq) (Eq. (16)). The compositions of even pure gas hydrates are variable. For the structure I of CO2 and CH4 gas hydrates, there are eight gas cavities (six large and two small) associated with 46 waters of hydration [7,8]. If all cavities are filled, then the formula would be M·5.75H2 O, where M is the gas molecule. If only the six large cavities are filled, then the formula would be M·7.67H2 O. Theoretical models, as a function of temperature and pressure, suggest CO2 hydration numbers between 5.9 and 6.3 (see Fig. 2 in Bakker et al. [7]). Because an objective of this paper was to develop a general-purpose model for gas hydrate equilibria, not a specific, compositionally dependent, model, we assigned the waters of hydration in our work a fixed value of 6.0 (Eqs. (17) and (18)). To estimate a gas hydrate solubility product requires knowing φg , Pg , and aw (Eq. (18)) as a function of temperature and pressure. The gas partial pressure, Pg , is experimentally measured. The activity of water, aw , is calculated by the FREZCHEM model (Eqs. (4), (5) and (12)), as is the gas fugacity coefficient (φg ) using a model developed by Duan et al. [24]. The equation used to calculate gas fugacity coefficients is given by: ln φ(T, P) = Z − 1 − ln Z +

2.3. Gas reactions + There are two types of reactions in this gas hydrate model that explicitly require estimation of the thermodynamic properties of gas components. The first type is the Henry’s law (K H ) reaction that controls equilibria between the gas and aqueous phases, which for a constituent “c” is given by: KH =

ac(aq) (γc ) (m c ) = f c(g) (φc ) (Pc )

(16)

where a is the activity in the aqueous phase, f is the fugacity in the gaseous phase, γ is the activity coefficient of the aqueous species, m is the molal concentration, φ is the fugacity coefficient, and P is the gas partial pressure. The second type of reaction is for gas hydrate equilibrium such as CO2 ·6H2 O, which is described by the reaction: CO2 ·6H2 O(cr) ↔ CO2 (g) + 6H2 O(l)

(17)

which is terms of the solubility product is given by:
f CO2 (g) (aw )6 K CO2 ·6H2 O = CO2 ·6H2 O(cr)

= φCO2 (g) PCO2 (g) (aw )6

(18)

where CO2 ·6H2 O in Eq. (18) is assumed to be a pure solid phase with an activity of unity. Gas concentrations were used to define gas hydrate equilibria because gas pressures are the experimental measurements in gas hydrate studies [8]. An additional benefit of using gas pressures directly, instead of

B C + Vr 2Vr2

D E + +G 4 4Vr 5Vr5

(19)

where Z =

D PV Pr Vr B C = =1+ + 2+ 4 RT Tr Vr Vr Vr     F −γ γ E . + 5 + 2 B + 2 exp Vr Vr Vr Vr2

(20)

Z is the compressibility factor that measures the departure of a real gas from an ideal gas (Z = 1). See Duan et al. [24] for the specification of the parameters for Eqs. (19) and (20). Given these parameters, it is possible to define all the terms in Eq. (20), except for Vr . In the FREZCHEM model, Brent’s method [25] is used to find the root of this equation (Vr ), which is then used to calculate Z (Eq. (20)), and finally φ (Eq. (19)), which is used in both Eqs. (16) and (18). Estimating solubility products for pure gas hydrates as outlined above is fairly straight-forward. Application of this model to real world situations is, however, generally more complicated because gases in natural environments are seldom pure gases. There are a number of models that can be used for gas mixtures (e.g., see [15,24,26]). In this work, we assumed that Amagat’s law of additive volumes holds for all pressures, which means that the ideal-solution law holds for the gaseous mixture, but not necessarily for the pure gases per se, and therefore the fugacity ( f ) is given by: f i = xi f i0 = xi φi0 Pi0

(21)

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

where xi is the mole fraction, f i0 is the fugacity of the pure gas, and φi0 is the fugacity coefficient at the temperature and total pressure (Pi0 ) of the gas phase [15]. See Fig. 14 in Duan et al. [24] for partial justification of this model for CH4 –CO2 gas mixtures. Another mixture problem is how to deal with mixed solidphase gas hydrates. Both methane and carbon dioxide form structure I gas hydrates. So CH4 –CO2 gas mixtures will ultimately lead to CH4 –CO2 gas hydrates. Equilibria for the simple gas hydrates can be represented by K CH4 ·6H2 O =

(CH4(g) )(H2 O)6 (CH4 ·6H2 O)

(22)

(CO2(g) )(H2 O)6 (CO2 ·6H2 O)

(23)

and K CO2 ·6H2 O =

251

3. Parameterization Gas fugacity coefficients (Eqs. (19) and (20)) were derived from the Duan et al. [24] model. The magnitudes of CH4 (g) and CO2 (g) fugacity coefficients as a function of pressure at 0.0 ◦ C are plotted in Fig. 1, which demonstrates that these gas fugacity coefficients are far removed from “ideal” values (φg = 1.0). There are an abundance of experimental gas partial pressures for gas hydrate equilibria across a broad range of temperatures and pressures (Fig. 2; [8]). The lower temperature limit in our model database for these systems is 180 K (Fig. 2) because this is the lower temperature limit of our model’s ability to estimate the activity of water (aw ) [16],which is needed to calculate the solubility product of gas hydrates (Eq. (18)). In our model, the upper temperature limit for methane hydrate is 298 K (25 ◦ C), which is the upper temperature limit for FREZCHEM; the upper temperature limit for carbon dioxide hydrate is 283 K

where parentheses represent activities (or fugacities). See Eq. (17) for the CO2 reaction. For pure gas hydrates, the terms in the denominators of Eqs. (22) and (23) are equal to unity. This is not the case, however, for mixed gas hydrates. In this work, we assumed that gas hydrate mixtures can be represented as solid solutions [13]. For the reaction: CO2 ·6H2 O + CH4(g) ⇔ CH4 ·6H2 O + CO2(g) ,

(24)

the equilibrium constant is given by (CH4 ·6H2 O)(CO2(g) ) K CO2 ·6H2 O = . (CO2 ·6H2 O)(CH4(g) ) K CH4 ·6H2 O

(25)

The activity of the solid phases can be represented as (X·6H2 O) = γX·6H2 O xX·6H2 O

(26)

where X is CH4 or CO2 , γ is the activity coefficient of the solid phase, and x is the mole fraction in the solid phase. Substituting Eq. (26) into Eq. (25) and rearranging terms leads to    xCH4·6H2 O γCH4 ·6H2 O K CO2 ·6H2 O (CH4(g) )  =   . (27) K CH4 ·6H2 O (CO2(g) ) xCO2 ·6H2 O γCO2 ·6H2 O

Fig. 1. The gas fugacity coefficients for methane and carbon dioxide at 0 ◦ C as a function of pressure using the Duan et al. [24] model.

The terms on the right-hand side of this equation are known (Fig. 3) or calculated by the model (gas fugacity = φg ·Pg ) (see previous discussions). If we assume an ideal solid solution (γCH4 ·6H2 O = γCO2 ·6H2 O ), then Eq. (27) simplifies to   xCH4 ·6H2 O K CO2 ·6H2 O (CH4(g) ) ≈  , (28) K CH4 ·6H2 O (CO2(g) ) xCO2 ·6H2 O which can be used directly to estimate the fraction of CH4 and CO2 in the gas hydrates since xCH4 ·6H2 O + xCO2 ·6H2 O = 1.0. Then these mole fractions can be substituted into the denominator of Eqs. (22) and (23) to convert these pure gas equations into mixed gas equations. Or alternatively, we could estimate γCH4 ·6H2 O and γCO2 ·6H2 O in Eq. (27) to get an improved model. Both of these approaches (Eqs. (27) and (28)) will be presented in the following text.

Fig. 2. Experimental data for methane and carbon dioxide hydrate equilibria. Data are from Sloan [8].

252

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

with Eq. E5.2.1 in Sloan [8]. The resulting methane and carbon dioxide hydrate densities for M·6H2 O were 0.9144 and 1.1204 g/cm3 , respectively. The molar volumes at infinite 0 dilution (V i ), needed in Eqs. (6) and (7), for these gas hydrates were calculated with the equation: 0

Vi =

Fig. 3. The solubility products of CH4 ·6H2 O and CO2 ·6H2 O at a hypothetical 1 atm total pressure. Symbols are data; lines are model fits. Two lines for Ln (KCH4 ·6H2 O ) are to improve the fits at T > 285 K and ≤285 K.

(10 ◦ C), which is the temperature above which liquid CO2 (l) becomes the thermodynamically stable phase. Given gas partial pressures (Fig. 2), a model to calculate φ (Eqs. (19) and (20)), and a model to calculate aw (Eqs. (4), (5) and (12)), then the solubility product (Eq. (18)) can be calculated for gas hydrates (Fig. 3). The actual solubility product calculations are made at the experimental gas partial pressures, which vary widely (Fig. 2). Eq. (6) was used to adjust all these pressure-dependent estimates (K P ) to a hypothetical 1.0 atm total pressure (K P0 ), which is what is presented in Fig. 3. Eq. (6) requires a specification of molar volumes and compressibilities for all constituents of chemical reactions. The latter are summarized in Marion et al. [20], except for solidphase gas hydrates. Gas hydrate densities (ρ) were estimated

M Wi ρi

(29)

where M Wi is the molecular mass for M·H2 O. The calculated molar volumes for methane and carbon dioxide hydrates were both 135.75 cm3 mol−1 , because these two gas hydrates are both structure I, and the differences in densities are exactly matched by differences in molecular mass (Eq. (29)). Molar compressibilities for solid phases in Eqs. (6) and (8) are ignored in the FREZCHEM model, except for ice; see Marion et al. [20] for the consequences of this simplification, which are minor. In the model validation and the application of this model, it is necessary to specify the activity coefficients of soluble gases (Eq. (3)). Table 1 summarizes the binary and ternary soluble methane Pitzer-equation parameters used in the FREZCHEM model; all of these parameters were taken from Duan et al. [26]. Only NaCl and CaCl2 salts are parameterized for temperature and pressure (Table 1); these salts cover brine compositions of 0–6 m, temperatures of 0–250 ◦ C, and pressures of 0–1600 bar (Duan et al. [26]). Fortunately, variations in λ and ξ with temperature and pressure are small (see Fig. 1a, b in Duan et al. [26]). Furthermore, given that most gas hydrates form in seafloor sediments dominated by NaCl, then this model should be accurate for seafloor sediment gas hydrates. For goodnessof-fit of these methane-salt parameterizations, see Figs. 2–4, 6 and 7 and Tables 6, 7 in Duan et al. [26]. An earlier paper [17] summarized the corresponding parameters for soluble carbon dioxide; all of those parameters were taken from He and Morse [27]. The latter database for

Table 1 Binary and ternary soluble methane Pitzer-equation parameters used in the FREZCHEM model taken or derived from Duan et al. [26] [Numbers are in computer scientific notation where e±x x stands for 10±x x ] Pitzer parameter

Equation parametera a1

λCH4 ,Na λCH4 ,K λCH4 ,Mg λCH4 ,Ca c λCH4 ,Fe d λCH4 ,Cl λCH4 ,SO4 λCH4 ,HCO3 λCH4 ,CO3 ξCH4 ,Na,Cl ξCH4 ,K,Cl ξCH4 ,Mg,Cl ξCH4 ,Ca,Cl ξCH4 ,Fe,Cl d

9.92230792e−2 0.13909 0.24678 −5.64278808e0 0.24678 0.00 0.03041 0.00669 0.16596 −0.00624 −0.00382 −0.01323 −0.00468 −0.01323

b

a b c d

a2

a3

a4

a5

2.57906811e–5

8.51392725e–3

Parameters are for the equation: Pitzer Parameter = a1 + a2 T + a3 T 2 + a4 T 3 + a5 /T , where T is temperature (K). This equation also contains the terms: +1.83451402e−2 ∗ P/T–8.07196716e−6 ∗ P2 /T, where P is pressure (bar) and T is temperature (K). This equation also contains the term: +5.27816886e−5 ∗ P, where P is pressure (bar). Assumed the same as the magnesium analogues.

1.00057752e3

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

253

carbon dioxide parameters is much broader in scope than the corresponding database for soluble methane parameters (Table 1). For example, Pitzer ion-interaction parameters (λ and ξ ) for CO2 (aq) with H+ , Na+ , K+ , Ca2+ , Mg2+ , Cl− , − SO2− 4 , and HSO4 are defined for the temperature range from 0 to 90 ◦ C and the concentration range from 0.01 to 6 m (He and Morse [27]). Because these CO2 parameterizations lack a pressure component, Eq. (9) was used in this paper to adjust these activity coefficients for pressure. For goodness-of-fit of these carbon dioxide-salt parameterizations, see Figs. 1 and 3 in He and Morse [27]. It is necessary to specify the Henry’s law constants (Eq. (16)) for both methane and carbon dioxide, which are given by: ln(K H,CH4 ) = −4.30210345e1 + 6.83277221e–2 · T 5.68718730e3 + − 3.56636281e–5 · T 2 T 5.79133791e1 − 6.11616662e–3 · P + (680 − T ) + 7.85528103e–4 · P · ln(T ) P + 9.42540759e–2 · T 1.92132040e–2 · P − (680 − T ) 9.17186899e–6 · P 2 + T

Fig. 4. A comparison of two models for estimating the total gas pressure of mixed CH4 –CO2 gas hydrates. Table 2 The equations relating the activity coefficient of CH4 ·6H2 O to the mole fraction of CH4 in the gas phase and temperature (K) (see Fig. 5).

(30)

[26], and ln(K H,CO2 ) = 2.495691e2 + 4.570806e–2 · T 1.593281e4 − 4.045154 · ln(T ) − T 1.541270e6 + T2

(31)

[17,28]. Eq. (30) is written as a function of both pressure (bar) and temperature (K), while Eq. (31) is written only as a function of temperature (pressure = 1.01 bar). Eq. (6) was used to adjust Eq. (31) for variable pressures. Earlier, we presented two models for calculating the stability of mixed methane/carbon dioxide gas hydrates (Eqs. (27) and (28)). We tested these concepts with experimental data from Adisasmito et al. [29], Dholabhai and Bishnoi [30], and Dholabhai et al. [31] for mixed CH4 –CO2 gas hydrates. Eq. (28) gives estimates of mole fractions that are used with Eq. (22) or (23) to estimate equilibrium gas pressures for mixed CH4 –CO2 solutions assuming an ideal solid solution (γCH4 ·6H2 O = γCO2 ·6H2 O ). This model overestimates the experimental measurements by 1% to 43% (Fig. 4). On the other hand, Eq. (27) in Fig. 4 is a fit of our model to the experimental data. In this case, we set γCO2 ·6H2 O = 1.00 and adjusted γCH4 ·6H2 O to optimize the fit to the experimental data. In Fig. 5 are estimates of γCH4 ·6H2 O as a function of xCH4 (g) and temperature. Within each of the eight datasets, the xCH4 (g) values are similar, and the numerical designation in the legend of Fig. 5 represents the average value. For the top four datasets [high xCH4 (g) ], it was possible to derive independent equations

XCH4 (g)

Equation

1.000 0.914 0.870 0.829 0.770 ≤ 0.595

γCH4 ·6H2 O γCH4 ·6H2 O γCH4 ·6H2 O γCH4 ·6H2 O γCH4 ·6H2 O γCH4 ·6H2 O

r2 = 1.000 = 2.07740 − 3.977471 × 10−3 T = 1.70079 − 2.722018 × 10−3 T = 2.04624 − 4.044674 × 10−3 T = 2.83506 − 6.937219 × 10−3 T = 4.22891 − 1.2173844 × 10−2 T

0.920 0.634 0.869 0.868 0.688

(Fig. 5, Table 2); for the bottom four datasets [low xCH4 (g) ] the data showed such scatter (Fig. 5) that we fit a single equation to these four datasets (Table 2). Note that as xCH4 (g) approaches 1.00 [pure CH4 (g) system], γCH4 ·6H2 O approaches 1.00 as it should theoretically for a pure solid phase. To estimate γCH4 ·6H2 O at an intermediate xCH4 (g) from the equations in Table 2, one can generate estimates above and below the desired xCH4 (g) and linearly interpolate between them. Applying these equations to estimate γCH4 ·6H2 O (with γCO2 ·6H2 O = 1.0) and using Eq. (27) led to a great improvement in the model estimate of total gas pressures for mixed CH4 –CO2 hydrates (Fig. 4); only 5 of 48 comparisons are off by more then 3%. Using Eq. (28) with γCH4 ·6H2 O = γCO2 ·6H2 O led to errors as large as 43% (Fig. 4). The system of equations in Table 2 can be simplified because both the intercepts and slopes of these equations are highly correlated to xCH4 (g) . The linear equation γCH4 ·6H2 O = a + bT with a = 1.00 + 7.998(1.00 − xCH4 (g) ) and b −0.03021(1.00 − xCH4 (g) ) can be written as

(32) =

γCH4 ·6H2 O = 1.00 + (7.998 − 0.03021T )(1.00 − xCH4 (g) ).(33) We can test these two alternative models (Table 2 equations or Eq. (33)) for estimating γCH4 ·6H2 O by rewriting Eq. (27) as

254

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

Fig. 5. The relationship between the activity coefficient of CH4 ·6H2 O and the mole fraction of CH4 (g) and temperature. Symbols are data; lines are model fits (see Table 2 for equations).

xCH4 ·6H2 O K CO2 ·6H2 O (CH4 (g))(γCO2 ·6H2 O ) = = K D. xCO2 ·6H2 O K CH4 ·6H2O (CO2 (g))(γCH4 ·6H2 O )

(34)

The gas hydrate equilibrium constants (K ) are known (Fig. 3), the gas fugacities are calculated with Eq. (21), γCO2 ·6H2 O = 1.0 (assumed), and γCH4 ·6H2 O is calculated with the equations in Table 2 (or Eq. (33)). Given that xCH4 ·6H2 O + xCO2 ·6H2 O = 1.0, then the mole fractions of gas hydrates can be calculated by xCO2 ·6H2 O = 1.0/(1.0 + K D )

(35)

and xCH4 ·6H2 O = K D /(1.0 + K D ).

(36)

For the case where xCH4 (g) = 0.770. (xCO2 (g) = 0.230), T = 280 K, P = 41.1 bar (equilibrium total gas pressure for a mixed CO2 –CH4 gas hydrate), the calculated xCH4 ·6H2 O = 0.663 for both Table 2 equations and Eq. (33), and xCO2 ·6H2 O = 0.337. For another case where xCH4 (g) = 0.265 (xCO2 (g) = 0.735), T = 280 K, P = 31.3 bar (equilibrium total gas pressure for a mixed CO2 –CH4 gas hydrate), the calculated xCH4 ·6H2 O = 0.179 (Table 2 equations) and 0.180 (Eq. (33)) and the corresponding xCO2 ·6H2 O = 0.821 and 0.820. Both Table 2 equations and the simpler Eq. (33) give, within roundoff error, identical results. In both of the above hypothetical cases, xCH4 ·6H2 O /xCH4 (g) <1.0 and xCO2 ·6H2 O /xCO2 (g) > 1.0, which means that CO2 is preferentially incorporated into mixed gas hydrates. This is probably largely due to the greater solubility of CO2 (aq) compared to CH4 (aq) (cf., Eqs. (30) and (31)). Note that both the equations in Table 2 and Eq. (33) will produce γCH4 ·6H2 O values >1 at subzero temperatures (extrapolations of equations in Fig. 5). As equilibrium gas pressure drops (and temperature drops), the system approaches ideal behaviour (see Eq. (28) in Fig. 4). In programming the FREZCHEM model for these relationships, we set an upper bound of γCH4 ·6H2 O = 1.0 at subzero temperatures. Also all xCH4 (g) values ≤ 0.595 are assumed to fit a single equation (Table 2, Fig. 5) or are assigned a value of 0.595 in Eq. (33). In this model, replacing the unit activities of pure simple hydrates in Eqs. (22) and (23) with the activities of the mixed

hydrates based on either Eq. (27) or (28) produces identical equilibrium gas pressures with either Eq. (22) or (23). 4. Validation The fit of the model to the experimental data (Figs. 3–5) is a test of the internal consistency of the experimental data and the derived model. While these fits are encouraging, they are not validation. To validate the model requires a comparison of the model output to independent data. Of special importance in this paper was to demonstrate that the FREZCHEM electrolyte model can accurately simulate the role of electrolytes in gas hydrate equilibria. We will examine several cases dealing primarily with the role of NaCl in gas hydrate equilibria. In Fig. 6, we compare the % difference between model and experimental measurements of methane gas concentrations for NaCl systems in equilibrium with CH4 ·6H2 O. What is most striking about these plots is the large variability in % difference

Fig. 6. The percent difference between model calculated and experimental methane gas concentrations for NaCl systems in equilibrium with CH4 ·6H2 O, where % Diff. = [(model − expt.)/expt.] × 100.

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

Fig. 7. The percent difference between model calculated and experimental carbon dioxide gas concentrations for NaCl systems in equilibrium with CO2 ·6H2 O, where % Diff. = [(model − expt.)/expt.] × 100.

both within and between datasets. Approximately 34% of the model estimates are greater than 10% in “apparent” error. Whether or not this error in model calculations is real or not will be discussed below. Most of the comparisons to the Kobayashi et al. [32] dataset are positive (i.e., our model overestimates the experimental measurements); on the other hand, most of the comparisons to the de Roo et al. [33] dataset are negative. The best fit of our model is to the Dholabhai et al. [34] dataset where % differences are 2%–3% (Fig. 6); but this dataset only contains comparisons at a single low concentration (3.0 wt.%, 0.529 m). Much of this variability, both within and between datasets, is clearly attributable to variation in measured equilibrium methane gas concentrations. For example, two measurements in the Kobayshi dataset at 273.3 K and NaCl = 20 wt.% produced equilibrium methane concentrations of 61.8 and 71.9 bar; our model predicted 74.1 bar at this temperature and NaCl concentration, which leads to “apparent” errors of 19.9% and 3.1%, respectively. In another comparison in the

255

Kobayashi dataset at two similar temperatures (276.3 and 276.4 K) and 20 wt.% NaCl, the measured equilibrium methane gas concentrations were 136.6 and 108.2 bar, respectively. Our model predicts values of 115.7 and 117.0 bar at these temperatures and NaCl concentrations, which leads to “apparent” errors of −15.3% and 8.1%, respectively. It is clear from these two cases and examples to be discussed below for CO2 ·6H2 O that is is difficult to accurately establish the equilibrium gas concentrations in these experiments. In Fig. 7, we compare the % difference between model and experimental estimates of carbon dioxide gas concentrations for NaCl solutions in equilibrium with CO2 ·6H2 O. Only 9% of the carbon dioxide comparisons are off by >10% (Fig. 7), which is an improvement over the 34% of Fig. 6 for methane. Again, a large part of this variability is clearly attributable to variation in measured gas concentrations. For example in the Dholabhai et al. [37] dataset at 278 K and NaCl = 5 wt.%, measurements in two cases gave carbon dioxide gas concentrations of 30.0 and 37.7 bar; our model predicts 28.0 bar at this temperature and NaCl concentration, which leads to “apparent” errors of −6.6% and −25.5%. The best fit of our model is to the Englezos and Hall [6] dataset where the % differences, with one exception, are within ±5% (Fig. 7 [35,36]). In Fig. 8, we compare equilibrium carbon dioxide gas partial pressures versus temperature for four datasets from Englezos and Hall [6] that include pure water, 10.0 wt.% NaCl (1.901 m), 15.2 wt.% NaCl (3.067 m), and 10.57 wt.% CaCl2 (1.065 m). The average absolute value of the error between our model and experimental measurements for these four datasets are 1.4%, 2.4%, 5.6%, and 2.7%, respectively [average overall error = 3.0(±1.8)%]. These errors in model fits to this specific experimental dataset are an improvement over the model errors of Englezos and Hall [6] of 2.4%, 6.2%, 12.2%, and 8.0%, respectively [average overall error = 7.2(±4.1)%], and the model errors of Bakker et al. [7] of 1.2%, 4.8%, 9.2%, and 5.7%, respectively [average overall error = 5.2(±3.3)%]. In general, our model smoothes the experimental bumps, with the possible exception of one datapoint at 15.2 wt.% NaCl and

Fig. 8. Model and experimental gas partial pressures of CO2 in equilibrium with CO2 ·6H2 O as a function of salt concentrations and temperature. Filled symbols and solid lines are data from Englezos and Hall [6]; open symbols and dashed lines are model fits.

256

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

274.45 K (Fig. 8), which is also an outlier in Fig. 7. Both the Englezos and Hall [6] and the Bakker et al. [7] models use a simpler version of the Pitzer equation for defining the activity of water in electrolyte solutions than our model that incorporates interactions between neutral species (dissolved gases) and salts (Eqs. (4) and (5)) and corrects for pressure (Eq. (12)). We also fit our model to experimental data for CO2 ·6H2 O in the presence of KCl and CaCl2 salt concentrations [8]. But these salt patterns were similar to those of NaCl (Figs. 6–8); so, with one exception (CaCl2 in Fig. 8), these KCl and CaCl2 data were not presented. We made similar comparisons to Dholabhai et al. [34] datasets for CH4 ·6H2 O in the presence of NaCl (one dataset), mixtures of NaCl and KCl (six datasets), and mixtures of NaCl and CaCl2 (six datasets). Using a model developed by Englezos and Bishnoi [5], Dholabhai et al. [34] found average errors for NaCl, mixtures of NaCl + KCl, and mixtures of NaCl and CaCl2 , of 3.74%, 5.54%, and 3.50%, respectively. Applying our model to these 13 datasets led to reduced errors of 2.60%, 3.68%, and 2.75%, respectively. As was the case for our earlier comparisons to salt–CO2 ·6H2 O systems, the Englezos and Bishnoi [5] model for salt–CH4 ·6H2 O systems uses a simpler version of the Pitzer model for calculating the activity of water than our model (Eqs. (4) and (5)). As was also the case for the previously discussed CO2 ·6H2 O systems, the CH4 ·6H2 O systems showed an increasing error with increasing salt concentration. For example, the average error for the 3% NaCl–3% KCl dataset was 0.93%; on the other hand, the average error for the 10% NaCl–12% KCl dataset was 5.28%. Increasing errors with increasing salinities (see patterns in the discussion above) have been observed in earlier work (e.g., [6,7,34]). Eq. (18) is the model used in this work to estimate equilibrium gas pressures. Since the errors in estimating gas partial pressures for pure water systems are minor (1.4% in Fig. 8), we can assume that the equilibrium constants are well-defined in Eq. (18). The one property in Eq. (18) that is most directly related to salinity is the activity of water (aw ), which in our model is estimated from the osmotic coefficient (φ) (Eqs. (4) and (5)). From the latter equations, as solution concentrations (m i ) increase, the potential for error in aw increases. However, a direct test of the FREZCHEM model to accurately simulate aw between 0.1 and 6.0 m NaCl at 298 K compared to an experimental database from Robinson and Stokes [38] led to a maximum error of 0.040% at 1.01 bar pressure. The calculated aw at 6.0 NaCl and 100 bar of pressure is only 0.13% smaller than at 1.01 bar. An error in aw of 1% raised to the 6th power (Eq. (18)) would cause an error of 6% in the calculated equilibrium gas concentration. It is possible that errors in calculated aw may contribute to the salinity effect on errors. However, the problem is more likely due to difficulties in establishing gas equilibria at high salt concentrations. There could be interactions between gas hydrates and salts that are not directly captured by Eqs. (16) and (18). For example, this could be due to “salting out” of dissolved gases from solution (removal) as salinity increases. This phenomenon would affect the connection between soluble and atmospheric gas concentrations (Eqs. (16) and (18)), which could affect both

φ and Pg in Eq. (18). The fact that experimental measurements appear to get bumpier with increasing salinity (Fig. 8) is also an argument supporting the difficulty of establishing gas equilibria with increasing salinity. Earlier in Figs. 6 and 7, we pointed out the high variability in equilibrium gas concentrations. Fig. 8 documents the role of salinity as an inhibitor of CO2 ·6H2 O formation. At 275 K, the equilibrium carbon dioxide gas pressures for pure water, 10.0 wt. % NaCl (1.901 m), and 15.2 wt.% NaCl (3.067 m) are approximately 15, 25, and 39 bar. As salinity increases, aw decreases, and Pg must increase to satisfy the solubility product (Eq. (18)). The gas hydrate model presented in this work is based on an existing electrolyte model for calculating aw in Eq. (18) (FREZCHEM model). The gas fugacity coefficient (φ) of Eq. (18) is based on the published Duan et al. [24] model. And finally, the equilibrium constants for CH4 ·6H2 O and CO2 ·6H2 O are based on model fits to pure gas/water data (Figs. 2 and 3). None of the salt data displayed in Figs. 6–8 or discussed in the text were used to parameterize our gas hydrate model. The fact that our model can accurately reproduce these independent experimental data is validation for the approach used in developing this gas hydrate model. 5. Limitations We pointed our earlier that the Pitzer-equation interaction parameters needed to calculate the activity coefficient of CH4 (aq) (Eq. (3)) have only a limited temperature range, except for NaCl and CaCl2 (Table 1). Fortunately, our solubility product equations for gas hydrates (Eqs. (18), (22) and (23)) are written in terms of gas pressures and not gas solution concentrations. So this temperature limitation for estimating activity coefficients of CH4 (aq) is not a major limitation of our model. The comparable databases for CO2 (aq) are much broader and even less of a limitation (see Marion [17]). At a temperature of 0 ◦ C, the fugacity coefficients of CH4 (g) and CO2 (g) vary significantly from ideality (φ = 1.0) as a function of pressure (Fig. 1). One potential limiting factor in using the Duan et al. [24] model to define fugacity coefficients is that the lower temperature limit for this model is 0 ◦ C (273 K). In Fig. 9A, we depict how the fugacity coefficient for methane changes according to the Duan model between 273 and 1473 K at 1 and 20 bar of pressure. Also included in Fig. 9A are FREZCHEM model estimates of the methane gas fugacity coefficient at subzero temperatures for systems in equilibrium with CH4 ·6H2 O at 1, 5, 10, 15, 20, and 26.45 bar of gas pressure; these calculations are an extrapolation of the Duan model into the subzero temperature range. Fig. 9B depicts similar calculations for CO2 fugacity coefficients in equilibrium with CO2 ·6H2 O. The equilibrium temperatures of methane and carbon dioxide hydrate at 1 bar pressure in these two cases are 193 and 218 K, respectively. The point of these diagrams is that they demonstrate that these subzero extrapolations are consistent with the trends at higher temperatures. Furthermore, the magnitude of the fugacity coefficients at subzero temperatures range from 0.916 to 0.994 for systems in equilibrium with gas hydrates (Fig. 9). Compared

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

257

Fig. 9. The Duan et al. [24] model estimates of gas fugacity coefficients at 1 and 20 bar between 273 and 1473 K compared to extrapolations below 273 K generated by the FREZCHEM model for equilibria of (A) methane hydrate and (B) carbon dioxide hydrate.

to the broader range of fugacity coefficients (Fig. 1), these subzero systems are close to ideal (φ = 1.0). For these reasons, we feel that the subzero temperature extrapolations of the Duan model are unlikely to cause significant error in equilibrium gas hydrate calculations. Carbon dioxide is important in governing bicarbonate/carbonate solution and mineral equilibria as well as controlling CO2 ·6H2 O formation. Many of the equations used for bicarbonate/carbonate equilibria are based on Plummer and Busenberg [28]; see Marion [17] for a complete description of bicarbonate/carbonate equilibria and model parameterizations. The maximum CO2 gas pressure for these bicarbonate/carbonate relationships is 1.01 bar (1 atm). Introducing CO2 gas hydrates into the FREZCHEM model necessitates increasing CO2 gas pressures up to 45 bar at 283 K (Fig. 2). For a hypothetical 1 × 10−3 m NaHCO3 solution, the FREZCHEM model predicts pH values of 8.16, 4.72, and 3.17 at 3.6 × 10−4 , 1.01, and 45 bar of CO2 pressure, respectively. While the first two points are within the range of our model, 45 bar of CO2 pressure is well beyond the range of the bicarbonate/carbonate model. So caution is necessary in interpreting mixed bicarbonate/carbonate and gas hydrate chemistries at high CO2 gas pressures. It may even be necessary, at times, when working with alkaline solutions such as seawater to remove alkalinity from the model when simulating high CO2 (g) pressures for gas hydrates. Finally, the model outlined in this paper is limited to methane hydrate, carbon dioxide hydrate, and their mixtures. There are innumerable other gases that can form gas hydrates, such as ethane, propane, nitrogen, and sulfur dioxide [8]. Incorporation of other gases into this model would require similar parameterizations to those of this study. 6. Discussion The main purposes of this study were to develop a gas hydrate model that could be seamlessly incorporated into a state-ot-the-art electrolyte model, and then validate the model.

In this section, we will first demonstrate a simple application that has relevance for all the Earth’s oceans. Fig. 10 depicts the stability of CH4 ·6H2 O (x = 1.0) and water–ice in seawater as a function of temperature and pressure (depth). This figure is based on “standard” seawater [39], which has the composition: Na = 0.48610 m, K = 0.01058 m, Mg = 0.05475 m, Ca = 0.01065 m, Cl = 0.56664 m, SO4 = 0.02927, and alkalinity = 0.00230 m. The depth (Y -axis) is calculated with P (bar) = ρ (g/cm3 )·g (cm/s2 )·D (km)/10

(37)

where ρ is density and g is gravitational acceleration (980.6 cm/s2 ). Densities of seawater were estimated using the FREZCHEM model. The equilibrium temperature of water–ice stability decreases almost linearly with increasing depth (pressure) (Fig. 10), which is in contrast to earlier work [40] that suggested that the temperature of water–ice stability actually increased over part of the pressure range. In our model, the equilibrium temperature of water–ice is −1.92 ◦ C at the surface (1.01 bar) and approximately −6 ◦ C at 5.0 km, while the earlier work [40] suggested a temperature of 0 ◦ C at at the surface (1.01 bar) and approximately −2 ◦ C at 5.0 km. The equilibrium temperature of water–ice in seawater is −1.92 ◦ C at 1.01 bar [41]. Experimental evidence clearly shows that the stability of water–ice decreases in temperature with increasing pressure, at least up to 2000 bar (≈20 km) [42,43]. Fig. 10 also includes the stability line for equilibrium between CO2 ·6H2 O and CO2 (g) (x = 1.0). Because CO2 gas pressures at equilibrium with CO2 ·6H2 O are less than for analogous CH4 systems (Fig. 2), CO2 ·6H2 O can occur closer to the surface (Fig. 10). Formation or dissolution of gas hydrates in natural settings requires an accurate gas hydrate/electrolyte model. Gas hydrates typically form in semi-closed environments such as in oceanic sediments and permafrost [44]. As gas hydrates form, water is removed from solution (Eq. (17)), which leads to concentration of salts in the residual water. Similarly, as gas hydrates dissolve, released water dilutes

258

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259

Fig. 10. The stability of water–ice and CH4 ·6H2 O in “standard” seawater (solid lines). The dashed line separates the analogous equilibrium between CO2 (g) and CO2 ·6H2 O. In both cases, the assumption was made that the mole fraction in the gaseous phase (x) is 1.0 [i.e., pure CH4 (g) or CO2 (g)].

salt concentrations. Changing salt concentrations changes the chemical environment and properties such as ionic strength that affect activity coefficients and the activity of water (Eqs. (1)–(5)). Accurately modeling gas hydrate chemistry in natural settings requires a state-of-the-art electrolyte model in order to optimize the predictive power of the model. Gas hydrates are hypothesized to have played a role in the geochemical evolution of Earth, Mars, and Europa [44– 46]. Application of a gas hydrate model to other-worldly environments where the chemistries may be quite different from seawater on Earth is another strong argument for a state-of-theart electrolyte model as a component of a gas hydrate model. 7. Conclusions We demonstrated that a classical thermodynamic approach to modeling gas hydrate chemistry is a viable alternative to the “standard” statistical thermodynamic model. Integrating these gas hydrate chemistries into a state-of-the-art electrolyte model generally led to an improvement in model accuracy for salt systems. In a comparison of four salt–CO2 ·6H2 O systems, the average error decreased from 7.2(±4.1%) to 3.0(±1.8)%. In a comparison of 13 salt–CH4 ·6H2 O systems, the average error decreased from 4.3(±1.1)% to 3.0(±0.6)%. Integrating gas hydrate chemistries into a state-of-the-art electrolyte model is the antithesis of the usual approach to developing gas hydrate models. The latter models are usually based on state-of-theart gas hydrate models to which salt effects are added [8]. The results of this work clearly demonstrate the critical importance of the electrolyte component in gas hydrate models. Understanding gas hydrates as a potential energy source, global climate driver, greenhouse gas repository, factor in oceanic slope stabilities, or factor in planetary geochemistry requires state-of-the-art knowledge and the ability to model controlling processes accurately. Integrating gas hydrate

chemistries into a state-of-the-art electrolyte model as was done in this work is a step in the right direction. Acknowledgements Funding was provided by a NASA Planetary Geology and Geophysics Project, “An Aqueous Geochemical Model for Cold Planets”, and a NASA EPSCoR Project, “Building Expertise and Collaborative Infrastructure for Successful Astrobiology Research, Technology, and Education in Nevada”. We thank Annette Risley and Lisa Wable for assistance in preparing the manuscript. The FORTRAN code for this gas hydrate model is available from GMM. References [1] J.H. van der Waals, J.C. Platteeuw, Adv. Chem. Phys. 2 (1959) 1–57. [2] W.R. Parrish, J.M. Prausnitz, Ind. Eng. Chem. Process Des. Dev. 11 (1972) 26–34. [3] G.D. Holder, G. Corbin, K.D. Papadopoulos, Ind. Eng. Chem. Fundam. 19 (1980) 282–286. [4] P. Englezos, Ind. Eng. Chem. Res. 31 (1992) 2232–2237. [5] P. Englezos, P.R. Bishnoi, AIChE J. 34 (1988) 1718–1721. [6] P. Englezos, S. Hall, Can. J. Chem. Eng. 72 (1994) 887–893. [7] R.J. Bakker, J. Dubessy, M. Cathelineau, Geochim. Cosmochim. Acta 60 (1996) 1657–1681. [8] E.D. Sloan Jr., Clathrate Hydrates of Natural Gases, second edn, Marcel Dekker, New York, 1998. [9] O.Ye. Zatsepina, B.A. Buffett, J. Geophys. Res. 103 (1998) 24,127–24,139. [10] K. Nasrifar, M. Moshfeghian, J. Chem. Thermodyn. 33 (2001) 999–1014. [11] M.A. Clarke, P.R. Bishnoi, Fluid Phase Equilib. 220 (2004) 21–35. [12] R. Masoudi, B. Tohidi, A. Danesh, A.C. Todd, Fluid Phase Equilib. 215 (2004) 163–174. [13] W. Stumm, J.J. Morgan, Aquatic Chemistry: An Introduction Emphasizing Chemical Equilibria in Natural Waters, Wiley-Interscience, New York, 1970. [14] K.S. Pitzer, in: K.S. Pitzer (Ed.), Activity Coefficients in Electrolyte Solutions, second edn, CRC Press, Boca Raton, 1991, pp. 75–153. [15] K.S. Pitzer, Thermodynamics, third edn, McGraw-Hill, New York, 1995.

G.M. Marion et al. / Computer Coupling of Phase Diagrams and Thermochemistry 30 (2006) 248–259 [16] G.M. Marion, Geochim. Cosmochim. Acta 66 (2002) 2499–2516. [17] G.M. Marion, Geochim. Cosmochim. Acta 65 (2001) 1883–1896. [18] G.M. Marion, D.C. Catling, J.S. Kargel, Geochim. Cosmochim. Acta 67 (2003) 4251–4266. [19] G.M. Marion, R.E. Farren, Geochim. Cosmochim. Acta 63 (1999) 1305–1318. [20] G.M. Marion, J.S. Kargel, D.C. Catling, S.D. Jakubowski, Geochim. Cosmochim. Acta 69 (2005) 259–274. [21] G.M. Marion, S.A. Grant, CRREL Spec. Rept. 94-18. USACRREL, Hanover, New Hampshire, 1994. [22] M.V. Mironenko, S.A. Grant, G.M. Marion, R.E. Farren, CRREL Spec. Rept. 97-5. USACRREL, Hanover, New Hampshire, 1997. [23] B.S. Krumgalz, A. Starinsky, K.S. Pitzer, J. Soln. Chem. 28 (1999) 667–692. [24] Z. Duan, N. Møller, J.H. Weare, Geochim. Cosmochim. Acta 56 (1992) 2605–2617. [25] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing, second edn, Cambridge University Press, Cambridge, 1992. [26] Z. Duan, N. Møller, J. Greenberg, J.H. Weare, Geochim. Cosmochim. Acta 56 (1992) 1451–1460. [27] S. He, J.W. Morse, Geochim. Cosmochim. Acta 57 (1993) 3533–3554. [28] L.N. Plummer, E. Busenberg, Geochim. Cosmochim. Acta 46 (1982) 1011–1040. [29] S. Adisamito, R.J. Frank, E.D. Sloan, J. Chem. Eng. Data 36 (1991) 68–71. [30] P.D. Dholabhai, P.R. Bishnoi, J. Chem. Eng. Data 39 (1994) 191–194. [31] P.D. Dholabhai, J.S. Parent, P.R. Bishnoi, Fluid Phase Equilib. 141 (1997) 235–246.

259

[32] R. Kobayashi, H.J. Withrow, G.B. Williams, D.L. Katz, Proc. of the 30th Annual Convention, Natural Gasoline Association of America, 1951, pp. 27–31. [33] J.L. de Roo, C.J. Peters, R.N. Lichtenthaler, G.A.M. Diepen, AIChE J. 29 (1983) 651–657. [34] P.D. Dholabhai, P. Englezos, N. Kalogerakis, P.R. Bishnoi, Can. J. Chem. Eng. 69 (1991) 800–805. [35] S.D. Larson, Phase studies of the two-component carbon dioxide-water system involving the carbon dioxide hydrate. Ph.D. Dissertation, Univ. Illinois, 1955. [36] J.G. Vlahakis, H.-S. Chen, M.S. Suwandi, A.J. Barduhn, Syracuse University Research and Development Report 830, 1972. [37] P.D. Dholabhai, N. Kalogerakis, P.R. Bishnoi, J. Chem. Eng. Data 38 (1993) 650–654. [38] R.A. Robinson, R.H. Stokes, Electrolyte Solutions, second edn (revised), Butterworths, London, 1970. [39] F.J. Millero, M.L. Sohn, Chemical Oceanography, CRC Press, Boca Raton, Florida, 1992. [40] K.A. Kvenvolden, Rev. Geophys. 31 (1993) 173–187. [41] F.J. Millero, Physical Chemisty of Natural Waters, Wiley-Interscience, New York, 2001. [42] G.M. Marion, S.D. Jakubowski, Cold Reg. Sci. Tech. 38 (2004) 211–218. [43] W. Wagner, A. Saul, A. PruB, J. Phys. Chem. Ref. Data 23 (1994) 515–525. [44] M.D. Max (Ed.), Natural Gas Hydrate in Oceanic and Permafrost Environments, Kluwer Academic Publ., Dordrecht, 2000. [45] M.D. Max, S.M. Clifford, Geophys. Res. Lett. 28 (2001) 1787–1790. [46] J.S. Kargel, J. Kaye, J.W. Head, III, G.M. Marion, R. Sassen, R. Crowley, O. Prieto, S.A. Grant, D. Hogenboom, Icarus 148 (2000) 226–265.