幾個常用地球化學模擬軟件的比較

劉杰安，馮孝貴

清華大學 核能與新能源技術研究院，北京 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.