海洋中水合物所处环境溶液(即海水)、当前普遍采用的水合物热力学抑制剂中的主要成分之一均为NaCl^[18]；并且对于NaCl这种溶质的相平衡研究也最为透彻。为使本研究更有针对性且便于对研究结果的验证，针对NaCl溶液环境进一步开展了相关的相平衡研究。

根据软件CSMHYD可以找到相平衡曲面上的点作为方程的拟合点。为使点的分布更均匀、适用性更强，应尽可能覆盖较大温度、压力和溶液盐浓度范围。目前水合物热力学抑制剂中NaCl的质量分数一般在20%以内^[19]，考虑到海水的稀释，工程中基本不会出现溶液环境浓度超出20%(w)的情况。同时其溶解度在25%左右，且受温度和压力影响较小，所以最终确定相平衡面中溶液NaCl质量分数X的范围是0%~25%。对于压力范围，目前已知的天然气水合物分布在水深500~4000 m，对应的压力范围是5~40 MPa^[19]。但高压范围内(80~200 MPa)的数据对于天然气水合物在油气安全开发和海水淡化中也有重要意义^[13]。因此，本研究相平衡面拟合中压力p范围为3~200 MPa。由此确定的温度T范围为257.8~316.8 K。据文献中拟合曲面的习惯，相平衡曲面(方程)的x、y、z轴分别为NaCl质量分数(%)、压力(MPa)、温度(K)。根据选定的范围确定拟合点的NaCl质量分数和温度，并由CSMHYD软件求出对应的压力p，共确定194个点作为拟合源点，见图 1。源点在拟合范围内均有分布，且在重要位置密度更大，以保证拟合结果在此区域的准确性。

图 1

图 1     拟合源点的分布 Figure 1     Distribution of fitting points

描述准确、合理且有物理含义的拟合函数形式对三维相平衡方程的拟合精度和适用性非常重要。由Clausius-Clapeyron方程可以看出，相平衡关系中临界温度T和压力p之间存在关系见式(1)。

式中：ΔH_d为天然气水合物分解焓；m为压缩因子；R为气体常数。

积分可以看出临界温度与压力的对数相关，见式(2)。

式中：a、b为常数。

考虑到离子浓度对压力和温度的影响，拟合函数中认为温度与压力的对数正相关，但具体形式与上式不同。对于NaCl质量分数X和温度T之间的关系，de Roo^[9]根据试验结果得到类似于温度和压力之间的关系见式(3)。

式中：c、d为常数；X′为摩尔分数，与质量分数正相关。

因此，拟合函数中认为温度与NaCl质量分数X的对数负相关。拟合源点对应的临界温度沿盐浓度方向单调递减，而自然对数函数是单调增函数，因此，拟合前将所有的点关于yOz平面对称，再向x轴正方向平移26个单位，得到一组新的拟合源点，从而方便使用自然对数函数进行拟合。拟合完毕后再进行逆向移动和对称，即可得到需要的显式函数表达式。根据上述分析，三维相平衡拟合方程形式采用分离变量形式(针对对称和平移后的拟合源点)，如式(4)所示。

式中: $\begin{array}{l} \mathit{A}(\mathit{X}) = {\mathit{a}_1} + {a_2}\ln X + {a_3}{[\ln\mathit{X}]^3}\\ {\kern 1pt} B(\mathit{X}) = {b_1} + {b_2}\ln X + {b_3}{[\ln\mathit{X}]^3}\\ C(\mathit{X}) = {c_1} + {c_2}\ln X + {c_3}{[\ln\mathit{X}]^3}\\ D(\mathit{X}) = {d_1} + {d_2}\ln X + {d_3}{[\ln\mathit{X}]^3} \end{array} $

a_i, b_i, c_i, d_i(i=1, 2, 3)为常数；T₀均为待拟合系数。

为提高拟合精度，假设函数中温度T与lnX和lnp的1~3次幂或者1~5次幂相关。为保证拟合函数在全域中的单调性得到满足，只保留了奇数次幂。试算发现，当幂次更高时，其对应的系数已经接近于零。拟合函数覆盖从零次到高次及其交叉项的所有可能项，可以使拟合方程对应的曲面与拟合源点更加接近。

1.2 拟合结果

使用最小化余差函数范式的方法确定拟合方程中的待定系数。首先给定系数初值，通过迭代方法找到函数某一组系数，该组系数使得函数在某一范围内整体的余差范式最小，最终得到了曲面方程的拟合结果见式(5)(已进行逆平移和对称)。

$ \begin{align} T(X,p)=247.268+3.7167\text{ln}(\text{26-}{X})\text{+} \\ 0.0793175{{[\text{ln}(\text{26-}{X})]}^{3}}+ \\ \text{ln}p\left\{ \text{9}.\text{57275-0}.\text{547938ln}(\text{26-}{X})\text{+} \right. \\ \left. 0.104963[\text{ln}{{(\text{26-}{X})}^{\text{3}}}] \right\}+ \\ {{[\text{ln}p]}^{3}}\left\{ -0.118983+ \right. \\ 0.0123896\text{ln}(\text{26-}{X}\text{)-} \\ \left. 0.00260223{{[\text{ln}(\text{26-}{X})]}^{3}} \right\}+ \\ {{[\text{ln}p]}^{5}}\left\{ 0.00388941- \right. \\ 0.0000814427\text{ln}(\text{26-}{X})\text{+} \\ \left. 0.0000307206{{[\text{ln}(\text{26-}{X})]}^{3}} \right\} \\ \\ \end{align} $

(5)

拟合优度已非常接近1。平移的单位可以由26改为其他正数，该数字决定了方程适用的NaCl的质量分数范围(本文是0%~26%，但在临近26%时与数据源点的误差较大，因此，建议在质量分数为25%以下时使用)。

图 2(a)给出由拟合方程确定的曲面及其拟合源点，图 2(b)为3%、9%、17%、25% 4种质量分数下的截线与拟合源点的对比。由图 2可以看出，曲面/线与拟合点在整个拟合范围差别都很小。拟合方程反映了工程应用范围内天然气水合物临界相平衡处温度、压力、化学浓度的基本规律：当固定化学浓度时，临界温度随压力的增加按照自然对数的幂函数形式增加；当固定压力时，临界温度随质量分数的增加按照自然对数的幂函数形式减小。

图 2

图 2     拟合结果和拟合源点的对比 Figure 2     Comparison between fitting resuits and fitting points

2 相平衡方程的验证和比较

将本研究提出的相平衡方程与多篇文献中覆盖面较广的二维和三维试验数据进行对比，从而验证方程在描述实际水合物分解中的适用性。

2.1 与固定盐浓度下试验数据及拟合曲线的比较

Dickens和Quinby-Hunt^[10]及Liu^[20]通过试验测定了纯水和海水环境中天然气水合物的相平衡点，并给出了海水环境下的拟合曲线，见式(6)(单位已换算成与本研究一致，下同)。

$ T(\mathit{p}) = \frac{1}{{3.79 \times {{10}^{ - 3}} - 2.83 \times {{10}^{ - 4}} \times \log p}} $

(6)

此后，Brown等^[21]、Peltzer和Brewer^[22]根据Dickens海水中相平衡试验数据分别得到了另外的经验公式。

图 3给出本研究结果与纯水(X=0)和海水环境(X=3.35%)试验数据以及上述3个拟合方程的对比。可见，本研究方程能够很好地满足试验结果，与Dickens等^[10]及其他文献^[21-22]的拟合曲线也十分接近。但除纯水和海水的化学环境，本研究结果还可适用于0%~25%的任何盐浓度，适用范围更广。

图 3

图 3     方程对应T-p曲线与文献中试验点和曲线的对比 Figure 3     Comparison between T-p curves by phase equilibrium equation in this paper and test data/curves in literature

2.2 与三维相平衡试验点的比较

de Roo等^[9]在温度为261.85~285.98 K、最大压力为11.0 MPa的条件下，通过试验研究了NaCl不饱和溶液中，温度和压力对甲烷水合物形成的影响。Maekawa等^[11]得到了在18 MPa以下压力范围内水合物相平衡的试验数据。这两组中低压力下实验得出的结果和本文的相平衡曲面之间的对比如图 4所示。由图 4可以看出，本文曲面和实验点在中低压条件时高度吻合。

图 4

图 4     相平衡曲面与deRoo等和Maekawa等试验点的对比 Figure 4     Comparison between phase equilibrium surface in this paper and de Roo's and Maekawa's test points

Jager和Sloan^[12]测量了在纯水和4种不同浓度NaCl溶液中，压力不超过70 MPa的情况下天然气水合物分解的相平衡数据。Hu等^[13]的研究揭示了在高盐度和高压的极端条件下，天然气水合物的不寻常的形成过程，压力提高至约200 MPa。将这两组高压和高盐度条件下测得的数据与本文曲面对比，如图 5所示，在压力较高的区域(大于50 MPa)，预测值开始出现了稍大的偏差，温度预测值均略高于试验点(15%＜X＜20%，50 MPa＜p＜150 MPa时的温度误差在3 K以内)。这是CSMHYD软件基于的理论不完善导致的^[8]。由于该软件基于理想固溶假设，因此只能用于描述此假设下的水合物相平衡，而在高压时不再严格满足这一条件。另外，Hu等^[13]认为在过饱和盐溶液中是液体-水合物-固体NaCl-气体水合物四相平衡，形成和盐的析出是相互竞争的关系，在高盐度时也会造成与实验曲面的误差。

图 5

图 5     相平衡曲面与Jager4及H u等对比试验点的对比 Figure 5     Comparison between phase equilibrium surface in this paper and Jager's and Hu's test points

总体来说，本研究相平衡曲面方程的拟合函数选取更为合理，高盐度、高压区域的预测值比现有曲面精确^{[9, 11]}。

3 三维相平衡曲面部分特点的验证

3.1 化学因素对T-p曲线的影响

基于相平衡方程，图 6给出不同盐浓度下的T-p曲线。由图 6可见，压力较大时，各曲线之间基本满足平移关系；但在压力较小时很难直观判断。因此，下面对相平衡方程两边对X求偏导，给出温度沿盐浓度方向的变化梯度。若该偏导数在压力p方向变化很小，则意味着盐浓度改变时，其温度差值在不同压力下基本相同，表示平移关系近似成立。求偏导后的结果见式(7)。

图 6

图 6     不同盐质量分数下的盐溶液T-p曲线 Figure 6     T-p curves for various mass fractions of salt solutions

$ \begin{array}{l} T{(\mathit{X}, \mathit{p})_x} = \left\{ {3.7167 - 0.547938\ln p + } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.0123896{[\ln\mathit{p}]^3} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.0000814427{[\ln\mathit{p}]^5} + \left\{ {0.237952 + } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0.314889\ln p - 0.00780668{[\ln\mathit{p}]^3} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {0.0000921617{{[\ln\mathit{p}]}^5}} \right\}\left. {\ln {{(26 - \mathit{X})}^2}} \right\}/\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (26 + \mathit{X}) \end{array} $

(7)

图 7为式(7)对应曲面图形。可见当盐浓度较低、压力高于10 MPa时，其偏导数沿压力方向几乎不发生变化，可以认为T-p曲线满足平移关系；但是在压力较低时，其偏导沿压力方向变化较为剧烈，T-p曲线在这一区域并不能视为平移。

图 7

图 7     相平衡方程对盐浓度求偏导后对应的曲面 Figure 7     Surface corresponding to the partial derivative of the phase equilibrium equation with respect to the salt concentration

3.2 Clausius-Clapeyron线性行为的验证

在不考虑盐浓度影响时，天然气水合物相平衡关系满足Clausius-Clapeyron方程。Dickens等^[10]、孙始财等^[23]以及其他研究人员的试验都证明了这一点。为验证本文方程是否具有这一特点，绘制研究得到的1/T-X-lnp关系的曲面见图 8，曲面整体接近一个平面。在不同的盐浓度下，1/T和lnp之间基本上保持线性关系，局部线性程度很好，整体有一定非线性，与前人的研究结果一致。

图 8

图 8     Clausius-Clapeyron线性行为的验证 Figure 8     Verification for Clausius-Clapeyron linear behavior

4 结论与展望

(1) 基于天然气水合物相平衡的统计热力学理论和CSMHYD软件，拟合得到温度-压力-NaCl盐浓度三维相平衡曲面方程。方程形式简洁、物理关系明确，适用的压力和NaCl溶液浓度区间几乎涵盖了目前工程中可能遇到的所有情况，且在描述已有的试验数据时具有很高的吻合度，特别在高压、高浓度区域比现有的拟合曲线或曲面更优，适用范围更广。

(2) 基于方程揭示了在压力高于10 MPa时，化学因素将近似导致T-p曲线平移，但压力较小时非平移所能描述；Clausius-Clapeyron线性行为在整个曲面上基本成立。

(3) 后续研究将基于更加准确的理论模型和软件(如CSMGem)以类似的方法拟合出更精确的曲面，并进一步考虑对水合物赋存于土体中时，更复杂的孔隙环境下相平衡曲面方程的修正。

参考文献

[1]	魏厚振, 颜荣涛, 韦昌富, 等. 含天然气水合物沉积物相平衡问题研究综述[J]. 岩土力学, 2011, 32(8): 2287-2294. DOI:10.3969/j.issn.1000-7598.2011.08.007
[2]	孙雷, 李士伦, 刘健仪, 等. 凝析气藏开采过程中天然气水合物生成条件预测方法研究[J]. 天然气工业, 2002, 22(1): 61-66. DOI:10.3321/j.issn:1000-0976.2002.01.018
[3]	BARRER R M, STUART W I. Non-stoicheiometric clathrate compounds of water[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1957, 243(1233): 172-189. DOI:10.1098/rspa.1957.0213
[4]	VAN DER WAALS J H, PLATTEEUW J C. Clathrate solutions[J]. Advances in Chemical Physics, 1959, 2: 1-57.
[5]	PARRISH W R, PRAUSNITZ J M. Dissociation pressures of gas hydrates formed by gas mixtures[J]. Industrial & Engineering Chemistry Process Design and Development, 1972, 11(1): 26-35.
[6]	PITZER K S, MAYORGA G. Thermodynamics of electrolytes. II. Activity and osmotic coefficients for strong electrolytes with one or both ions univalent[J]. The Journal of Physical Chemistry, 1973, 77(19): 2300-2308. DOI:10.1021/j100638a009
[7]	HOLDER G D, HAND J H. Multiple-phase equilibria in hydrates from methane, ethane, propane and water mixtures[J]. AIChE Journal, 1982, 28(3): 440-447. DOI:10.1002/aic.690280312
[8]	BALLARD A L, SLOAN J R E D. The next generation of hydrate prediction I. Hydrate standard states and incorporation of spectroscopy[J]. Fluid Phase Equilibria, 2002, 194-197: 371-383. DOI:10.1016/S0378-3812(01)00697-5
[9]	DE ROO J L, PETERS C, LICHTENTHALER R N, et al. Occurrence of methane hydrate in saturated and unsaturated solutions of sodium chloride and water in dependence of temperature and pressure[J]. AIChE Journal, 1983, 29(4): 651-657. DOI:10.1002/aic.690290420
[10]	DICKENS G R, QUINBY-HUNT M S. Methane hydrate stability in seawater[J]. Geophysical Research Letter, 1994, 21(19): 2115-2118. DOI:10.1029/94GL01858
[11]	MAEKAWA T, ITOH S, SAKATA S, et al. Pressure and temperature conditions for methane hydrate dissociation in sodium chloride solutions[J]. Geochemical Journal, 1995, 29(5): 325-329.
[12]	JAGER M D, SLOAN E D. The effect of pressure on methane hydration in pure water and sodium chloride solutions[J]. Fluid Phase Equilibria, 2001, 185(1/2): 89-99.
[13]	HU Y, MAKOGON T Y, KARANJKAR P, et al. Gas hydrates phase equilibria and formation from high concentration NaCl brines up to 200 MPa[J]. Journal of Chemical & Engineering data, 2017, 62(6): 1910-1918.
[14]	SLOAN E D, KOH C A. Clathrate hydrates of natural gases[M]. 3rd ed. NewYork: CRC Press, 2008.
[15]	SLOAN E D, KOH C A. CSMHYD. [CP/DK]. 2019. http://hydrates.mines.edu/CHR/Software.html#opennewwindow.
[16]	CLENNELL M B, HOVLAND M, BOOTH J S, et al. Formation of natural gas hydrates in marine sediments: 1. Conceptual model of gas hydrate growth conditioned by host sediment properties[J]. Journal of Geophysical Research, 1999, 104(B10): 22985-23003. DOI:10.1029/1999JB900175
[17]	HENRY P, THOMAS M, CLENNELL M B. Formation of natural gas hydrates in marine sediments: 2. Thermodynamic calculations of stability conditions in porous sediments[J]. Journal of Geophysical Research, 1999, 104(B10): 23005-23022. DOI:10.1029/1999JB900167
[18]	张洋, 李广雪, 刘芳. 天然气水合物开采技术现状[J]. 海洋地质前沿, 2016, 32(4): 63-68.
[19]	盛杰.深水钻井液中天然气水合物抑制技术研究[D].青岛: 中国石油大学(华东), 2009.
[20]	LIU C L, YE Y C, SUN S C, et al. Experimental studies on the p-T stability conditions and influencing factors of gas hydrate in different systems[J]. Science China Earth Sciences, 2013, 56(4): 594-600. DOI:10.1007/s11430-012-4564-3
[21]	BROWN K M, BANGS N L, FROELICH P N, et al. The nature, distribution, and origin of gas hydrate in the Chile Triple Junction region[J]. Earth and Planetary Science Letters, 1996, 139(3/4): 471-483.
[22]	PELTZER E T, BREWER P G. Practical physical chemistry and empirical predictions of methane hydrate stability[M]//MAX M D. Natural Gas Hydrate: In Oceanic and Permafrost Environments. Dordrecht: Springer, 2000: 17-28.
[23]	孙始财, 刘昌岭, 业渝光, 等. 氯盐溶液中甲烷水合物高压分解条件及影响因素[J]. 物理化学学报, 2011, 27(12): 2773-2778. DOI:10.3866/PKU.WHXB20112773

X