几个常用地球化学模拟软件的比较

刘杰安，冯孝贵

清华大学 核能与新能源技术研究院，北京 102201

放射性废物尤其是高放废物的安全处置问题已成为制约核能可持续发展的因素之一。由于高放废物含有放射性强、发热量大、毒性大、半衰期长的核素，因而对它的安全处置是一个世界性难题[1]。在众多处置方案中，高放废物地质处置是开发时间最长、也是目前最有希望投入应用的处置方案[2]。高放废物地质处置就是把高放废物处置在深部的地质处置库中，使之与人类的生存环境永久隔离。它是一项以放射性核素的包容、阻滞为核心内容，以多重屏障(地质介质属于天然屏障，废物体、包装容器和缓冲回填材料等属于工程屏障)为主要手段，以几万年以上公众健康和环境保护为安全目标的极其复杂的系统工程。其中，研究各关键核素在各种屏障材料中的吸附扩散规律(即核素迁移研究)，对于地质处置库的设计是必不可少的。由于地质处置库的设计必须考虑的时间尺度(10万年至100万年)远大于一般实验研究所能达到的时间尺度，因此相关的研究不能仅仅依靠实验，还必须借助于模型研究，即根据实验测得的基本数据和基本的物理学和化学原理，建立描述核素扩散及水力学弥散的数学模型，用以预测一旦地质处置库受损后放射性核素的时间和空间分布。

上述的模型研究包括水力学和化学等多方面的内容，现在已有相关的成熟软件，其中关于化学方面的软件就是所谓的地球化学模拟软件。国际上被广泛采用的地球化学模拟软件有CHESS[3]、EQ3/6[4]、MINEQL+[5]、MINTEQA2[6]、PHREEQC[7]和WATEQ4F[8]等。国内在这方面的研究工作起步比较晚，目前公开发表的仅有北京大学核环境化学组编写的CHEMSPEC[9]。由于地球化学模拟软件数量比较多，并且各软件都有其自身的特点，因此，在应用这些软件之前，有必要了解这些软件的概况及其主要的异同点。为此，本工作对其中常用的4个软件CHEMSPEC(V2A)、EQ3/6(V8.0)、MINEQL+(V4.5)和VMINTEQ(V2.52，即MINTEQA2的可视版本)进行比较研究，以期为相关软件的选用提供参考。

1 各模拟软件的原理与结构比较

上述4个软件CHEMSPEC、EQ3/6、MINEQL+和VMINTEQ均建立在热力学平衡基础上。各软件将水溶液中的化学物质区分为“组分(components)”和“物种(species)”两大类，后者由前者通过化学反应生成，后者为产物，前者为反应物。同一元素有不同的化学形态时，参照地球化学的习惯选取其中的某一种作为组分(在选择组分时不同软件有时会存在差别)。组分和物种通过平衡常数K进行关联：对于组分，lgK=0；对于物种，取其生成反应的平衡常数。各软件中主要的控制方程有质量(摩尔)平衡方程、电荷(电子)平衡方程、化学反应平衡方程、温度校正方程和活度系数校正方程等。尽管不同软件采用的温度校正方程和活度系数校正方程等各不相同，但它们的基本原理均相同。

这些软件的结构都有一个共同的特点，即均由输入、输出、迭代算法及化学热力学数据库等几个模块组成。其中，输入模块一般有初始溶液的组分总浓度和参加反应的固相、气相等组成；输出模块则主要有各平衡组分/物种的浓度、溶液与固相和气相间的物质交换量等；迭代算法模块用于非线性方程组的求解，基本上都是采用改进的Newton-Raphson迭代法；化学热力学数据库模块包括标准状态(25 ℃，1 atm(1×105Pa))下的热力学平衡常数、平衡常数的温度系数或不同温度下的值、活度系数计算所需的参数、水溶液中各物种的化学计量因子等。

2 各模拟软件的功能比较

CHEMSPEC、EQ3/6、MINEQL+和VMINTEQ四个软件的功能列于表1。该4个软件主要的共同点是：(1) 4个软件都可以进行配合、氧化-还原、溶解-沉淀等3类基本化学反应的模拟计算；(2) 都可以采用Davies模型进行活度系数校正。主要的差别是：(1) MINEQL+和VMINTEQ具有处理吸附和表面配合这2类化学反应的能力，而CHEMSPEC和EQ3/6不能进行这2类化学反应的模拟计算；(2) CHEMSPEC和EQ3/6提供了多种活度系数校正方法，尤其是Pitzer模型使EQ3/6能够对高离子强度的体系进行模拟计算；(3) EQ3/6和MINEQL+可以处理离子交换反应，EQ3/6采用Vanselow和Gapon模型，而MINEQL+则定义了1个新的反应；(4) 4个软件中仅有EQ3/6可以进行反应动力学的模拟计算。

因为各软件不同的功能之间的差别是显而易见的，所以本工作的重点是对各软件的相同功能进行比较。由于各软件对于氧化-还原反应的处理方式都相同，即：在所研究的体系中，电子给体和受体间处于平衡状态，所有氧化-还原电对的电极电位都相等，由于氧化-还原电对的标准电极电位与氧化-还原反应的平衡常数可以通过Nernst方程关联在一起，因此各个软件在进行氧化-还原反应模拟计算中产生的差异，与它们在进行配合、溶解-沉淀反应模拟计算中产生的差异，两者在本质上并没有特别的不同。故下面的比较主要针对配合、溶解-沉淀反应进行。

3 4个软件的实例计算比较

一般来讲，各个软件计算结果的差别主要来源于其数据库的差别(各种校正被认为是数据库的一部分)。

表1 4个地球化学模拟软件的功能比较

注(Note)：△表示“有该项功能”(△ means “having the function”)

CHEMSPEC(V2A)采用瑞士 PSI的NAPSI_290502(260802)，共有66个组分和448个物种。

EQ3/6(V8.0)的数据库数量比较多，来源也有些不同。不同的数据库用3个字符的后缀进行区分，包括500、1kb、2kb、5kb、cmp、hmw、pit、shv、skb、sup、ymp和ypf等，其中ymp和ypf是Yucca Mountain Project专用数据库。各个数据库所包含的组分和物种数量有比较大的差异，cmp数据库所包含的组分和物种数量最多，如果没有特别说明，下面的计算均采用该数据库。能用Pitzer模型进行活度系数校正的数据库有3个：hmw、pit和ypf，其中hmw数据库只适用于25 ℃，pit数据库适用于0～100 ℃，ypf数据库适用于0～300 ℃。能进行压力校正的数据库有5个：500、1kb、2kb、5kb和shv。

MINEQL+(V4.5)采用MINTEQA2的数据库，也有原始的MINEQL的一些数据，其数据都是标准状态下无限稀释溶液的数据。

VMINTEQ(V2.52)是MINTEQA2的可视版，因此其数据库与MINTEQA2相同，也都是标准状态下无限稀释溶液的数据。

对于常压下的模拟计算而言，数据库的差别主要包括标准状态下平衡常数、温度校正方程以及活度系数校正方程等的差别。下面先采用单平衡点计算的方法，以比较这3种差别对计算结果的影响；然后再采用多平衡点计算的方法，以便于比较这4个软件对某一个实际体系进行模拟计算时，其结果随某一个参数在一定范围内变化的趋势。

3.1 单平衡点计算比较

1) 标准状态下平衡常数的差别

用4个软件进行模拟计算的结果列于表2。从表2可以看出：

表2 例1的计算结果

注(Note)： -表示“没有该项”(- means “none”)

表3 例1中的平衡常数

注(Note)： -表示“没有该项”(- means “none”)

(4) MINEQL+和VMINTEQ的计算结果的微小差别，一方面可能是2个软件中的次要物种有一点差别；另一方面也可能与计算结果的有效数字位数有关，MINEQL+中取4位有效数字，而VMINTEQ中取5位有效数字。

2) 温度校正方程的差别

CHEMSPEC、MINEQL+和VMINTEQ等3个软件的数据库中只有25 ℃下的数据，对于其它温度下的模拟计算，需要进行温度校正。

CHEMSPEC对平衡常数的温度校正采用如下公式(A、B、C、D、E为经验常数)：

lgKT=A+BT+C/T+DlgT+E/T2

(1)

若数据库中没有提供相应的温度校正系数，而提供了焓变值，则用Van’t Hoff公式校正。

MINEQL+采用Van’t Hoff公式对反应平衡常数进行校正。

VMINTEQ一般也采用Van’t Hoff公式对反应平衡常数进行温度校正。另外，对其中25种组分，还可以采用以下经验公式进行校正(A、B、C、D、E、F、G为经验常数)：

lgKT=A+BT+C/T+DlgT+

ET2+F/T2+GT1/2

(2)

此经验公式比Van’t Hoff公式准确，所以对这25种组分，VMINTEQ比MINEQL+的温度修正更准确。

由于焓变是温度的函数，所以以上3个软件的温度应用范围比较窄，均为0～100 ℃。

EQ3/6数据库中收集了除25 ℃以外的数据，所以其温度适用范围更广，但各个数据库对应的温度不同。常用的cmp数据库适用于0～300 ℃。

例2：计算不同温度(25～90 ℃)时Na2HPO4-KH2PO4标准缓冲溶液的pH值。

图1 不同软件计算结果相对误差随温度的变化

以文献[10-11]中的标准pH值为基准，对4种软件计算出的pH值进行比较，其相对误差示于图1。从图1可以看出，EQ3/6计算出的pH值的相对误差最小，其绝对值始终小于1%。这是由于其数据库中含有0～300 ℃时的经验常数，所以其温度校正更准确。CHEMSPEC和VMINTEQ的校正结果基本一致，在温度较低时误差较小，随着温度的增加，误差也越来越大，如在90 ℃时，误差达到了3.5%。上述3个软件在25 ℃时的误差均比较小，但MINEQL+在25 ℃时的误差却比较大，这是由于计算时没有进行活度系数校正引起的(本工作曾尝试进行活度系数校正，但结果没有收敛)。

3) 活度系数校正方程的差别

对于水溶液中的带电组分，CHEMSPEC采用Davies公式、扩展的Debye-Hückel公式和B-dot公式计算离子活度系数(根据溶液离子强度I和数据库是否有某离子的Debye-Hückel参数a和b进行选择)。EQ3/6数据库较多，各数据库对应的活度系数计算公式有所区别，包括Davies公式、B-dot公式和Pitzer公式。MINEQL+采用Davies公式计算活度系数。VMINTEQ采用Davies公式和修正的Debye-Hückel公式计算活度系数。Davies公式、扩展的Debye-Hückel公式和B-dot公式都只适用于离子强度I小于0.5 mol/L的情况。当离子强度大于0.5 mol/L时，只能用Pitzer公式。

对于水溶液中的电中性组分，除MINEQL+没有进行活度系数校正以外，其余3个软件都有相应的校正公式。

CHEMSPEC采用B-dot公式计算活度系数；MINEQL+采用Davies公式计算活度系数；VMINTEQ主要采用扩展的Debye-Hückel公式计算活度系数，但CaCl+和MgCl+采用Davies公式；EQ3/6采用3种方法计算活度系数，分别为Davies公式(cmp数据库)、B-dot公式(cmp数据库)和Pitzer公式(hmw数据库)。计算结果列于表4和表5。从表4和表5可以看出：

(1) 对于H+，其活度相同(pH=5)，但各软件计算出的浓度各不相同；EQ3/6采用相同的数据库(cmp数据库)和不同的活度系数计算方法(Davies公式和B-dot公式)，其结果有较大差异(1.391 5×10-5和1.281 0×10-5)；采用Davies公式的结果更接近于MINEQL+的结果(1.391 5×10-5和1.410×10-5)，采用B-dot公式的结果更接近于CHEMSPEC的结果(1.281 0×10-5和1.241 1×10-5)；这说明在离子强度较大时，计算活度系数的不同方法之间的差别与稀溶液相比是显著增加的；

(2) 物种CaCl+和MgCl+在有的数据库有，而在另一些数据库中却没有；有这两个物种时，对应的计算结果里游离的Mg2+和Ca2+浓度都较小；

表4 例3的计算结果(CHEMSPEC,MINEQL+和VMINTEQ)

注(Note)： -表示“没有该项”(- means “none”)

表5 例3的计算结果(EQ3/6)

注(Note)： -表示“没有该项”(- means “none”)

(3) 与例1相比，MINEQL+和VMINTEQ两个软件计算的结果之间的差距明显增加，这是因为：一方面，例1的体系是稀溶液，而本例的离子强度较大(各软件计算的平均值约为0.25 mol/L)，这种情况下计算活度系数的不同方法之间的差别也较大；另一方面，MINEQL+和VMINTEQ这两个软件的数据库中的物种有差别，在VMINTEQ中有CaCl+和MgCl+这两个物种，而在MINEQL+中没有。

3.2 多平衡点计算比较

上面的单平衡点计算不能反映计算结果随某一个参数变化的趋势，为了从变化趋势上对4个软件进行比较，因此进行了下面的计算。

例4：计算不同pH值(pH=5～9)时总浓度为10-6mol/L的Am(Ⅲ)在北山51号井地下水(地下水组成见文献[9])中的化学形态分布。

首先需要说明两点：(1) MINEQL+的数据库中没有Am(Ⅲ)的数据，有关Am(Ⅲ)的数据均采用CHEMSPEC中的数据。(2) 在MINEQL+的计算过程中，pH=8.5和pH=9两个点的计算没有进行活度系数的校正(本工作曾尝试进行活度系数校正，但结果没有收敛)。

图2—图5列出了Am(Ⅲ)的几种主要的化学形态随体系pH值的变化趋势，从中可以看出：

(3) 在图2、图3和图5中，当pH>6时，各物种的浓度随着pH值的增加而逐渐减少，这主要是由于生成了AmCO3OH(cr)沉淀的缘故；另外，当pH值达到9.5(对CHEMSPEC，pH=9)时，还会产生Am(OH)3(cr)沉淀；

图2 Am3+计算浓度的比较

图计算浓度的比较

图计算浓度的比较

图计算浓度的比较

3.3 计算中的其它差别

1) 数据库中数据的增补

4个软件中，MINEQL+增补数据最方便，可以直接在软件用户界面中增加组分和物种，也可以直接对物种的平衡常数进行修改。VMINTEQ可以直接在软件用户界面上增加物种，并对物种的平衡常数进行修改；但对于基本组分的增加，需要修改基本组分的数据文件comp.mdb。CHEMSPEC和EQ3/6增补数据都比较麻烦。CHEMSPEC需要对组分文件components.dtb、物种文件species.dtb和化学计量矩阵文件stoichiometry.dtb进行修改，并需要使用相关的数据库生成和格式转换程序。而EQ3/6为了增补组分和物种，需要修改数据文件data0，然后通过运行EQPT程序产生相应的data1文件；但如果只需要临时改变某些物种的平衡常数，也可以直接在编辑输入文件时进行修改。

2) 收敛性

迭代计算的收敛性与初值的选取和迭代算法有关。通过具体实例的计算发现(收敛判据均采用各软件的默认值)，MINEQL+在计算过程中最容易发生不收敛的情况(参见例2和例4)。VMINTEQ中允许产生的沉淀不能超过20种；然而在例4中，若将所有饱和度大于0的物种都设为沉淀的话(没有超过20)，计算就不能收敛，这就需要将部分最不可能生成沉淀的物种从沉淀列表中剔除掉。在本工作计算例子中，EQ3/6和CHEMSPEC没有出现不收敛的情况。

3) 运行时间

对于简单的体系，MINEQL+和VMINTEQ的运行时间在1 s以上，EQ3/6的运行时间为0.1～1 s，而CHEMSPEC的运行时间小于0.1 s(电脑配置：Intel(R) Core(TM)2 Duo CPU T5450 @ 1.66 GHz；1 GB内存)。CHEMSPEC的运行时间明显少于其它的几个软件，这可能与其在算法方面的改进有关[9]。

4) 单位换算

EQ3/6的输入与输出文件中以mol/g或mg/kg溶液为单位，而其它3个软件都以mol/L或mg/L为单位，这就需要进行换算。如果已知溶液的密度，则可以很容易地进行换算；如果溶液的密度是未知数，则可以假定其密度为1 kg/L。对于稀溶液来说这种假定误差较小；但当溶液不是足够稀时，会产生一定的误差。

4 小 结

各个地球化学模拟软件具有相同的基本功能，均可进行配合、溶解-沉淀、氧化-还原等反应的模拟计算。各软件有各自的特点，CHEMSPEC在算法上有所改进，更有利于计算的收敛；EQ3/6有众多的数据库，能适用于高温、高压及高离子强度下的计算；MINEQL+很容易对数据库中的数据进行增补；VMINTEQ具有强大的处理吸附和表面配合反应的能力。

计算结果的差别，主要来源于数据库中数据的不同。对于常压下的模拟计算而言，数据库中标准状态下平衡常数、温度校正方程以及活度系数校正方程等的差别，都会引起计算结果的偏差。对标准状态下的稀溶液，各个软件的计算结果基本一致。

[1] 王 驹,徐国庆,郑华铃，等.中国高放废物地质处置研究进展:1985-2004[J].世界核地质科学,2005,22(1):5-16.

[2] Laurence S C.Site Selection and Characterization Processes for Deep Geological Disposal of High Level Nuclear Waste[R].Albuquerque: Sandia National Laboratories, 1997．

[3] van der Lee J.Thermodynamic and Mathematical Concepts of CHESS, Technical Report Nr.LHM/RD/98/39[R].Fontainebleau, France: École des Mines de Paris, 1998.

[4] Wolery T J.EQ3/6: Software for Geochemical Modeling, Version 8.0, UCRL-CODE-2003-009[R].Livermore, California: Lawrence Livermore National Laboratory, 2002.

[5] Schecher W D, McAoy D C．MINTEQL+: A Chemical Equilibrium Modeling System, Version 4.5 for Windows, User’s Manual[R].Hallowell, Environmental Research Software, 2003.

[6] Allison J D, Brown D S, Novo-Gradac K J.MINTEQA2/PRODEFA2, a Geochemical Assessment Model for Environmental Systems, Version 3.0 User’s Manual, EPA/600/3-91/021[R].Athens, Georgia: Environmental Research Laboratory, 1991.

[7] Parkhurst D L, Appelo C A J.User’s Guide to PHREEQC (Version 2): A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations, Water-Resources Investigations Report 99-4259[R].Colorado, Denver: 1999.

[8] Ball J W, Nordstrom D K.User’s Manual for WATEQ4F, With Revised Thermodynamic Data Base and Test Cases for Calculating Speciation of Major, Trace, and Redox Elements in Natural Waters, U.S.Geological Survey Open-File Report 91-183[R].Menlo, 1991.

[9] 王祥云,陈 涛,刘春立.化学形态分析软件CHEMSPEC及其应用[J].中国科学B辑-化学,2009,39(11):1-12.

[10] 顾庆超主编.新编化学用表[M].南京:江苏教育出版社,1998:903.

[11] 常文保,李克安.简明分析化学手册[M].北京:北京大学出版社,1981:259.