1 引言

水不仅是许多地质流体和成矿流体中极为重要的组分，同时也是许多地外行星的重要成分。在地球和行星深部，水一般处于超临界状态。超临界水在密度、介电常数、粘度、扩散系数、电离平衡常数和溶解能力方面都与普通水相差悬殊(徐有生等，1995; 苏根利等，1998)。超临界水的热力学性质对于许多地质过程的研究都非常重要。例如，含水流体及水岩相互作用的热力学模拟可以用来解释有关矿物、岩石或矿床的成因所需的物理化学条件(Holland and Powell， 2011; Hu et al., 2000; 魏春景，2012; Zhang et al., 2013)。通过含水的流体包裹体热力学性质的测试、实验和计算研究可以获得相关地质过程的温度、压力和流体组成等物理化学条件或其变化范围(刘斌等，2000; 徐文刚等，2012; 席斌斌等，2010; Song et al., 2007; Mao et al., 2010，2013)。此外，超临界水对于工业过程的开发、环境工程及高温高压物理方面的研究也很重要(韩勇等，2010; Duan et al., 2008; Hu et al., 2007)，如以超临界水为介质的分解、氧化、催化、合成与转化过程(Akizuki et al., 2014; Guo et al., 2010; Marrone，2013)等。

鉴于超临界水的重要性，人们对它的压强(P)、摩尔体积(V)和温度(T)、相平衡等许多热力学性质进行了大量的实验测量(Mao et al., 2011; Wagner and Pruß，2002; Guo et al., 2015)，其中压强也常被简称为压力，而摩尔体积则常被简称为体积。20世纪90年代以来，国际上一些研究组采用新的高温高压实验技术(Brodholt and Wood， 1994; Frost and Wood， 1997b; Withers et al., 2000; Abramson and Brown， 2004)和原子-分子水平上的计算机模拟技术(Belonoshko and Saxena， 1991a; Brodholt and Wood， 1993; Wasserman et al., 1995; Destrigneville et al., 1996; Zhang and Duan， 2005)大大拓展了超临界水PVT数据的温压范围，其中部分数据列于表 1。

与此同时，人们也提出了许多热力学模型来计算超临界水的PVT性质(Gottschalk，1997; Bakker，2003)，其中一部分为Helmholtz自由能模型，而更多的则是状态方程。在众多的状态方程中，多数是半经验的维里型(virial-type)方程(Lee and Kesler， 1975; Benedict et al., 1940; Saxena and Fei， 1987; Brodholt and Wood， 1993; Pitzer and Sterner， 1994; Duan et al., 1992a，1996; Zhang and Duan， 2005; Duan and Zhang， 2006; Soave， 1995，1999)，或特定的统计力学理论的近似加上半经验处理的结果(Sun and Dubessy， 2012; Anderko and Pitzer， 1993; Duan et al., 1995; Churakov and Gottschalk， 2003; Ji et al., 2005)。这些方程一般都能用于高压，但它们的形式一般比较复杂，参数比较多，在给定温压下一般得不到体积的解析式，其中甚至有一部分方程无法通过积分得到Gibbs自由能(G)或Helmholtz自由能(A)、逸度系数()、剩余焓(H^R)、剩余熵(S^R)等重要热力学性质的解析式(Halbach and Chatterjee， 1982; Mäder and Berman， 1991; Grevel and Chatterjee， 1992)。这些方程在向混合物推广时一般需要很多混合物实验数据来拟合参数。相对而言，立方型方程形式简单，参数少，使用方便，外延性好，可以解析地求出体积根(具体方法见附录A)，而且在向混合物推广时所需的实验数据较少。因此，立方型方程在工程中的应用非常广泛(Valderrama，2003; Duan and Hu， 2004; Abdollahi-Demneh et al., 2010a，b)，但它们在多数情况下只用于中低压(一般在几十MPa以内)和中低温(多数在水的临界温度以下)(Duan and Hu， 2004; Hassanzadeh et al., 2008; Wong and Sandler， 1992; Springer et al., 2012; Wu et al., 2012)。

不过，特别值得一提的是Holloway(1977)所提出的Redlich-Kowng方程的一种简单修正形式：

其中R是通用气体常数，b为协体积(covolume，14.6cm³·mol^-1)，a(T)为温度的多项式。该方程虽然非常简单，却能用于压力高达1GPa左右的超临界流体(H₂O和CO₂)。这表明立方型方程在高温、高压流体的热力学计算方面仍有很好的改进或发展潜力。此后陆续出现了许多关于超临界地质流体状态方程的改进或新方程开发方面的研究(Halbach and Chatterjee， 1982; Bowers and Helgeson， 1983; Holland and Powell， 1991; Kerrick and Jacobs， 1981; Mäder and Berman， 1991; Bottinga and Richet， 1981; Grevel and Chatterjee， 1992; Brodholt and Wood， 1993; Saxena and Fei， 1987; Belonoshko and Saxena， 1991a，b，1992; Pitzer and Sterner， 1994; Duan et al., 1992a，b，1996; Zhang and Duan， 2005; Duan and Zhang， 2006)，详见Gottschalk(1997)的综述。不过，除Bowers and Helgeson(1983)方程外，这些修正形式一般都不再是体积的三次方。

需要指出的是，许多PVT模型(状态方程或自由能模型)在参数拟合过程中使用了一些不准确或不可靠的数据(Bowers and Helgeson， 1983; Pitzer and Sterner， 1994; Brodholt and Wood， 1993; Frost and Wood， 1997a，b)，从而给模型的计算结果带来难以估计的偏差或不确定性。这主要是因为纯流体和混合物的PVT数据中都有相当一部分(其比例可能高达2/3左右)存在明显的系统误差或较大的不确定度(Hu et al., 2007; Wagner and Pruß，2002; Guo et al., 2015)。另外，和纯流体相比，流体混合物的PVT数据要少得多，其中高质量的数据更是缺乏，而且在温度-压力-组成(T-P-x)空间中的分布很不均匀，其中很多数据集中于一个有限的条件范围内。目前，像H₂O-NaCl和CO₂-H₂O这样拥有较多高温、高压PVT数据的地质流体混合物体系仍然非常少见(Hu et al., 2007; Mao and Duan， 2008)。即使是像CO₂-H₂O-NaCl非常重要的三元混合物体系，其高温高压PVT数据不仅数量少，误差通常也较大(Hu et al., 2007; Mao et al., 2010)，因而远不能满足建立一些具有宽广适用范围的可靠的热力学模型的需要。对于含有很多参数的状态方程，这种情况很可能会产生过度拟合(overfitting)现象，从而难以估计方程的计算或外延结果的可靠性。对于多组分体系的热力学模拟而言，发展立方型方程在当前的数据情况下是一种明智的选择。另外，近20年来，人们通过一些新的高温高压实验技术和计算机模拟技术所获得的超临界水的PVT数据覆盖了很宽的温压范围(表 1)。这些数据可以用来检验现有的状态方程，或者用来发展精度更高或适用范围更宽的状态方程。

考虑到立方型方程的优越性，本文作者在经过仔细评价的高精度PVT数据的基础上提出了一种新的多参数立方型状态方程。该方程可以在很宽的温压范围内精确地计算超临界水的PVT关系及相关的热力学性质。

2 状态方程 2.1 方程形式的选择

本文提出的多参数立方型状态方程采取如下形式：

其中R=8.31451J·K^-1·mol^-1，b为常量，a、c与d均为温度的函数，P、V和T的单位分别为MPa、cm³·mol^-1和K。由方程(2)可以得到第二维里系数B=b-a/RT。因此，a与温度之间的函数表达式的选择可参考第二维里系数的一些统计力学理论的推导结果。大量结果表明，许多真实流体的分子间势能与距离的关系非常接近于Lennard-Jones势能函数。由统计力学所导出的Lennard-Jones流体第二维里系数B与温度的多项式函数中温度的指数为1/4的整数倍。胡家文(2002)据此选择温度指数为1/4的整数倍或其相近值的简单多项式来拟合第二维里系数，结果较好。有鉴于此，a表达式中温度的指数也取1/4的整数倍。而c和d的表达式则在兼顾简洁性的同时根据拟合效果来确定。拟合结果表明，c和d随温度的变化比较剧烈。基于上述几点，经过大量的尝试，a、c和d最后采用了如下形式： 其中θ=T/T_C-1，其中T_C为水的临界温度(647.096K)。

2.2 方程参数的拟合

为了保证方程的可靠性，新方程参数拟合所用的PVT数据均来自International Association for the Properties of Water and Steam(IAPWS)推荐的高精度模型IAPWS-95(273.16~1275K，0~1GPa)(Wagner and Pruß，2002)和IAPWS-IF97方程的5区(1073.15~2273.15K，P≤10MPa)(Wagner et al., 2000)。IAPWS-95是在1990年以前的PVT、热容等多种热力学数据基础上建立的，其密度偏差如下：在P≤10MPa、T≤423K范围内液态密度的不确定度为0.0001%~0.001%；气相密度的不确定度为0.05%或更小。密度的不确定度随温度或压力的升高而增大，但一般不会超过0.1%(除了极端的条件)。临界区压力的不确定度为0.1%。IAPWS-95方程适用于273.16~1275K、0~1000MPa和0~1.332g/cm³，其精度和可靠性得到了国际同行的认可，并被广泛应用。不过，上述模型是为水而精心设计的，不容易向混合物推广。它的主要目的是为科学研究和一般性的应用提供尽可能精确的数据。根据上述数据，利用最小二乘法可以得到方程中的各个参数(表 2)。

本文作者发现，在临界点附近的超临界区，PVT性质的拟合十分困难，这是因为流体在临界点附近具有奇异性。在平均场理论框架内所建立的解析形式的状态方程在原则上无法正确地描述这种奇异性。虽然一些非经典(因而也是非解析的)的理论模型或方法可以准确地描述这种奇异性，如重整化群和跨接函数法，但它们通常都极为复杂，计算量也非常大(Kiselev and Ely， 2004，2007; Kiselev et al., 2006)。这一缺陷大大限 制了它们的实际应用。目前有不少半经验方法可以近似地表达流体的临界奇异性(Thiéry and Dubessy， 1996; Kedge and Trebble， 2002，2004; Chen and Mi， 2001; Mi et al., 2003)，但相应的方程一般都是针对工业流体(主要是纯流体，混合物比较少见)，因而一般只能用于比较有限的温压范围。也有一些为特定流体精心设计的状态方程能够精确地适用于临界区，但它们一般都比较复杂，或者含有很多非线性参数，难以向混合物推广，如Wagner and Pruß(2002)、Span and Wagner(1996)、Kim(2007)和Pitzer and Sterner(1994)方程。考虑上述情况，新方程的参数优化未考虑临界点附近的区域。经过反复尝试，新方程参数优化的温度下限取为723.15K。

2.3 新方程的检验以及与其它方程的比较

为了检验新方程的精度和外延性，本文不仅采用了高精度模型IAPWS-95和IAPWS-IF97的PVT计算值，同时还采用了一些高温高压PVT数据。这些数据包括两类：一是通过人工合成包裹体技术所得的PVT数据，如Brodholt and Wood(1994)、Frost and Wood(1997b)和Withers et al.(2000)；二是通过分子动力学(MD)和Monte Carlo(MC)模拟所得的TIP4P和SPCE模型水的PVT数据，如Brodholt and Wood(1993)(TIP4P，MD)、Wasserman et al.(1995)(SPCE，MD)、Destrigneville et al.(1996)(TIP4P，MC)和Zhang and Duan(2005)(SPCE，MD)。这些数据覆盖了很宽的温压范围(表 1)。

Guo et al.(2015)的评价结果表明，Zhang and Duan(2005)的数据不仅温压范围很宽，精度也较高，与IAPWS-95模型的结果及其外延值吻合得非常好。Wasserman et al.(1995)的高压数据(1GPa以上)与Zhang and Duan(2005)的结果相当吻合。另外两组高压模拟数据(1GPa以上)以及三组实验数据的不确定度较小，但均有一定的系统偏差。这些系统偏差很有规律，因而易于校正。本文利用IAPWS-95模型或其适当外延结果对几组数据的系统偏差进行了修正。

另外，Brodholt and Wood(1994)、Frost and Wood(1997b)和Withers et al.(2000)在相同条件合成的同一批包裹体样品的摩尔体积或密度会有显著的差异。Abramson and Brown(2004)发现，如果同一批包裹体样品的平均摩尔体积用最小摩尔体积(对应于最大密度)来代替，则结果与状态方程计算值的吻合情况将会有显著的改善。他们认为，实验中的包裹体或多或少都有泄漏。因此，最小摩尔体积应该代表没有泄漏或泄漏最少时的结果。Guo et al.(2015)的计算表明，同一批样品摩尔体积的差异明显超出了样品合成条件波动的影响范围。也就是说，包裹体的泄漏很可能是产生这种差异的主要因素。这与Abramson and Brown(2004)的结论的评价结果是一致的。另外，Guo et al.(2015)还发现，Withers et al.(2000)用来计算包裹体密度的基本方程有系统偏差，因为确定该方程参数所用的PVT数据有明显的系统偏差。本文利用IAPWS-95模型(或其适当外延结果)以及建立在Zhang and Duan(2005)分子模拟结果基础上的状态方程Duan and Zhang(2006)对几组数据的系统偏差进行了校正。这些修正的数据连同Wasserman et al.(1995)的高压数据(1GPa以上)和Zhang and Duan(2005)的结果一起用于检验新方程。

为了更全面地了解新方程，有必要将新方程与一些现有的地质流体状态方程在相同条件下进行系统的比较。考虑到地质流体状态方程为数较多，本文只能选择一些比较知名且有代表性的状态方程与新方程进行比较。这些方程包括：

(1) Bowers and Helgeson(1983)(BH83)方程：

它是Holloway(1977)(H77)的一个变种，其中a(T) 被改成了温度的指数函数。

(2) Kerrick and Jacobs(1981)(KJ81)方程：

其中y=b/4V。

(3) Pitzer and Sterner(1994)(PS94)方程：

其中ρ=1/V，每个c_i表达式中的非零参数实际上只有1~4个。方程共有27个参数。

(4) Duan et al.(1996)(DMW96)方程(一种基于对应态原理的状态方程)：

其中P_m=3.0626σ³P/ε，T_m=154T/ε，σ和ε为Lennard-Jones流体分子间势能参数，m是指对比量的参考流体为methane。

上述方程原作者所声明的适用范围见表 3。不过，此处应该说明的是，其中多数方程真正的适用范围与表 3中的情况有显著的差别。例如，DMW96方程可以适用的最低温度并不是1.3T_C(841.2K)(段振豪等，2000)，而是约900K(1.39T_C)。H77方程可以适用的最高压力实际上约为1GPa，最低压力实际上也不止500bar，因为在此压力附近方程的最大体积偏差达10%以上。BH83方程可以适用的最高压力约为1GPa左右，远低于原作者所声称的4GPa；最低压力约为1500bar，而不是500bar。由于上述两个方程不适于低压，导致它们对逸度系数的计算值不准确。由BH83方程计算出的临界温度约为702K，比实验值(647.096K)高出约55K。这样，在大约702K以下的超临界流体区(单相)，BH83方程却给出一个两相区(汽-液平衡区)，从而导致状态方程在上述区域及其附近几十度的范围内出现严重的偏差(图 1a)。KJ81方程和许多其它的地质流体状态方程在临界温度的实验值和计算值上下几十度的范围内也都有类似的情况。由KJ81方程计算出的临界温度约为686K，比实验值高出约39K，由此导致方程在590~740K之间的体积偏差均较大(图 1b)。尤其是在647.096K附近，体积偏差非常大，其中最大偏差接近70%。因此，KJ81方程实际的适用范围要明显小于表 3中的结果。

本文对新方程和上述几个方程做了多种情况下的检验与比较。首先是在新方程的适用范围内用IAPWS-95和IAPWS-IF97来检验新方程和表 3中的方程。新方程对IAPWS-95(723.15~1273.15K，0~1GPa)模型的数据的再现结果见图 2和图 3，其中平均和最大偏差分别为0.263%和3.13%。这个结果明显优于H77、BH83、KJ81和DMW96方程，而与PS94方程的结果很接近。从图 2中可以看出，PS94方程的整体精度(平均偏差)不如新方程，只是在723K附近略优于新方程；而且，PS94方程有明显的系统偏差。

与IAPWS-IF97(1073.15~2273.15K，0~10MPa)相比，新方程的平均和最大偏差分别为0.0235%和0.211%(图 2和图 3)。在同样条件下，PS94与DMW96方程的平均和最大偏差都与新方程非常接近；而H77、BH83和KJ81方程的平均和最大偏差分别达到了新方程的2.2~6.6和2.6~4.8倍。

本文在下面三种条件下对新方程和表 3中的方程做了更全面的检验与比较：

(ⅰ)上述方程所声明的适用范围与新方程的适用范围的重叠区，结果如表 4和图 4所示。可以看出，在此条件下，绝大多数方程均明显不如新方程，其中只有PS94方程与新方程的精度最为接近。与IAPWS-95模型的结果相比，新方程在723K附近不如PS94方程，但在总体上(平均偏差)优于PS94方程；与高温、高压下的实验或分子模拟结果相比，新方程的最大体积偏差(2.25%)明显优于PS94方程的结果(2.88%)，但其平均体积偏差(1.05%)略逊于PS94方程(0.945%)。

(ⅱ)上述重叠区以外、1.4GPa以下的高温区，结果如表 5和图 5a所示。显然，在这种情况下，新方程的最大偏差和平均偏差均显著地优于其它各个方程。

(ⅲ)上述重叠区以外、1.4~2.0GPa的区域，结果如表 6和图 5b所示。在这种情况下，新方程显著地优于H77、BH83、KJ81和PS94方程。只有DMW96方程的精度与新方程最为接近，其最大偏差略优于新方程，但其总体精度(平均偏差)则略逊于新方程。

新方程与几组高温、高压PVT数据的直接比较见图 6。很明显，新方程的计算值与实验及分子模拟结果在很宽的温压范围内都吻合得很好。

3 由状态方程计算一些重要的热力学性质

除密度或体积外，在涉及流体的地球物理或地球化学模拟中，人们经常需要计算流体的膨胀系数α、压缩系数β、逸度系数、剩余焓H^R与剩余熵S^R(即同温同压下实际流体与理想气体的焓与熵之差)、Gibbs或Helmholtz自由能(G或A)、热力学能(U)等许多热力学性质。这些性质对于人们了解地球或行星内部过程中力和能量的转化或传递非常重要。例如，许多地质过程常常伴随显著的热效应，而热物理性质的计算需要用到H^R和S^R。流体的膨胀系数和压缩系数则与流体的许多物理性质(如声速和粘度)密切相关。逸度系数是实际流体的热力学性质对理想气体偏离程度的一种量度。它可以用来计算流体的摩尔Gibbs自由能G(即化学势)：

其中f=P，f为逸度，p⁰为标准压力，G⁰(T)为标准态下的摩尔Gibbs自由能，其值可以由标准态下的焓、熵和热容计算而得到。(11)式表明，逸度系数也是剩余Gibbs自由能G^R的一种量度，即G^R=RTln。从Gibbs自由能出发，经过适当的微积分可得到许多重要的热力学性质。许多涉及流体的地质过程(如水-岩相互作用)的热力学模拟或计算都离不开Gibbs自由能。对上述性质的准确计算对于人们正确理解、分析或模拟有关的地质过程非常重要。考虑到Gibbs或Helmholtz自由能等热力学性质很容易由前几种性质的简单运算而得到，此处只需给出前几种性质的算式。

3.1 表达式的推导

(1)膨胀系数和压缩系数

由于方程(2)不是V的显函数，(12)式中的求导难以直接实现。为此可利用 将(12)式转变化为 其中

(2)逸度系数、剩余焓和剩余熵

这些性质需要通过状态方程的积分才能得到。新方程第二项中体积的二次式(V²+cV+d)能否被分解成V的两个一次因式完全取决于c²-4d的值。如令Δ=c²-4d，则由于新方程所适用的温压范围很宽，Δ的值可能会出现正、负和零三种情况。相应地，a/(V²+cV+d)的积分式也随会Δ而变化，从而导致逸度系数、剩余焓、剩余熵的表达式也随会Δ而变化。

考虑到方程(2)是显压型状态方程，即压力是其它变量(体积和温度)的显函数，逸度系数的推导可以使用分部积分法，即将积分自变量P转换成V(Hu et al.，2008)：

其中Z=PV/RT，V^id为同温、同压下理想气体的摩尔体积。

根据焓和剩余焓的定义以及理想气体状态方程(PV=RT)可知H^R=U^R+PV^R=U^R+PV-RT。据此，剩余焓的推导可以使用下式：

当Δ>0时 当Δ=0时， 当Δ＜0时， 当Δ→0时(即u→∞时)，上式中的反正切函数值趋于π/2。本文作者发现，此时π的有效数字的多少对剩余焓计算值的准确性影响很大。因此，π的取值必须有足够多的有效数字，比如π=3.1415926535897932384626。

在ln和H^R均已知的情况下，可由G^R=H^R-TS^R得到剩余熵，其中G^R=RTln。

3.2 表达式的检验

文献中一般很难见到高温、高压流体的逸度系数、剩余焓或剩余熵数据。为了方便，本文用IAPWS-95模型计算出这些性质，然后用它们来检验新方程。考虑到极低压力下的真实气体可视为理想气体，根据理想气体的H或S与压力的关系不难由IAPWS-95模型算出真实流体的H^R或S^R，进而算出G^R。然后由G^R=RTln可得ln。

另外，文献中也很难见到膨胀系数α和压缩系数β的实验数据。新方程对α和β计算值的准确性可用等压热容(C_P)和等容热容(C_V)之差来间接地检验。这是因为理论上已证明

其中C_P-C_V可用IAPWS-95模型计算出来。

新方程与IAPWS-95模型对ln、C_P-C_V、H^R和S^R计算值的比较见图 7和图 8。由图 7可以看出，方程(2)对ln的计算值与IAPWS-95模型的结果吻合得相当好。虽然在不同的温度范围内，由于Δ的符号不同而导致ln的表达式有所不同，但在同一压力下由方程(2)算出的ln-T曲线是非常光滑的，如图 7b所示(相应的Δ符号有变化)。由图 8a可以看出，方程(2)对C_P-C_V的计算值与IAPWS-95模型的结果也吻合得非常好。其中C_P-C_V与压力的关系曲线上尖峰的的位置和峰高都预测得非常好。此外，方程(2)对H^R和S^R的计算结果也与IAPWS-95模型的结果吻合得较好，如图 8b和8c所示。

4 结论

有鉴于高质量的混合物PVT数据的数量及其在温度-压力-组成空间中的分布情况远不如纯流体数据，本文利用2273.15K和1GPa以下水的高精度PVT数据建立了超临界水的一种新的立方型状态方程。新方程可以精确地再现723.15~2273.15K、0~1.4GPa范围内的PVT数据，同时也能较好地适用于2273.15~4273.15K、1.4~2.0GPa。这说明新方程具有良好的外延性。新方程不仅显著地优于同类的立方型方程，而且也优于许多常用的多参数高次维里型方程。除了PVT性质以外，通过对新方程的适当微积分运算可以得到膨胀系数、压缩系数、逸度系数、剩余焓、剩余熵等许多其它重要的热力学性质的解析表达式。这些表达式的计算结果与高精度模型IAPWS-95的结果均很吻合。

致谢 衷心地感谢两位匿名评审专家对本文的评论及所提出的建设性修改意见。

附录A 解析法求三次方程的根

任何一种立方型状态方程均可转化成如下形式的三次标准方程：

其中X可以是体积V或密度(ρ=1/V)。令 则有 其中 令 代入(A3)式有 从(A4)及(A5)式可知，p³与q³为下列二次方程的解： 即 令

(1)若Δ>0，则p、q各有一个实根，即p₁和q₁：

p、q分别还有两个共轭复根p₂和p₃、q₂和q₃： 方程(A3)的三个根为

(2)若Δ=0，则从(A10)式可得

方程式(A3)的三个实根为 其中两个为重根。

(3)若Δ＜0，则

其中 令 则θ=θ₀、θ₀+2π、θ₀+4π三个值，即方程(A3)有三个不同的实根：