幾個常用地球化學(xué)模擬軟件的比較

劉杰安，馮孝貴

清華大學(xué) 核能與新能源技術(shù)研究院，北京 102201

放射性廢物尤其是高放廢物的安全處置問題已成為制約核能可持續(xù)發(fā)展的因素之一。由于高放廢物含有放射性強、發(fā)熱量大、毒性大、半衰期長的核素，因而對它的安全處置是一個世界性難題[1]。在眾多處置方案中，高放廢物地質(zhì)處置是開發(fā)時間最長、也是目前最有希望投入應(yīng)用的處置方案[2]。高放廢物地質(zhì)處置就是把高放廢物處置在深部的地質(zhì)處置庫中，使之與人類的生存環(huán)境永久隔離。它是一項以放射性核素的包容、阻滯為核心內(nèi)容，以多重屏障(地質(zhì)介質(zhì)屬于天然屏障，廢物體、包裝容器和緩沖回填材料等屬于工程屏障)為主要手段，以幾萬年以上公眾健康和環(huán)境保護(hù)為安全目標(biāo)的極其復(fù)雜的系統(tǒng)工程。其中，研究各關(guān)鍵核素在各種屏障材料中的吸附擴散規(guī)律(即核素遷移研究)，對于地質(zhì)處置庫的設(shè)計是必不可少的。由于地質(zhì)處置庫的設(shè)計必須考慮的時間尺度(10萬年至100萬年)遠(yuǎn)大于一般實驗研究所能達(dá)到的時間尺度，因此相關(guān)的研究不能僅僅依靠實驗，還必須借助于模型研究，即根據(jù)實驗測得的基本數(shù)據(jù)和基本的物理學(xué)和化學(xué)原理，建立描述核素擴散及水力學(xué)彌散的數(shù)學(xué)模型，用以預(yù)測一旦地質(zhì)處置庫受損后放射性核素的時間和空間分布。

上述的模型研究包括水力學(xué)和化學(xué)等多方面的內(nèi)容，現(xiàn)在已有相關(guān)的成熟軟件，其中關(guān)于化學(xué)方面的軟件就是所謂的地球化學(xué)模擬軟件。國際上被廣泛采用的地球化學(xué)模擬軟件有CHESS[3]、EQ3/6[4]、MINEQL+[5]、MINTEQA2[6]、PHREEQC[7]和WATEQ4F[8]等。國內(nèi)在這方面的研究工作起步比較晚，目前公開發(fā)表的僅有北京大學(xué)核環(huán)境化學(xué)組編寫的CHEMSPEC[9]。由于地球化學(xué)模擬軟件數(shù)量比較多，并且各軟件都有其自身的特點，因此，在應(yīng)用這些軟件之前，有必要了解這些軟件的概況及其主要的異同點。為此，本工作對其中常用的4個軟件CHEMSPEC(V2A)、EQ3/6(V8.0)、MINEQL+(V4.5)和VMINTEQ(V2.52，即MINTEQA2的可視版本)進(jìn)行比較研究，以期為相關(guān)軟件的選用提供參考。

1 各模擬軟件的原理與結(jié)構(gòu)比較

上述4個軟件CHEMSPEC、EQ3/6、MINEQL+和VMINTEQ均建立在熱力學(xué)平衡基礎(chǔ)上。各軟件將水溶液中的化學(xué)物質(zhì)區(qū)分為“組分(components)”和“物種(species)”兩大類，后者由前者通過化學(xué)反應(yīng)生成，后者為產(chǎn)物，前者為反應(yīng)物。同一元素有不同的化學(xué)形態(tài)時，參照地球化學(xué)的習(xí)慣選取其中的某一種作為組分(在選擇組分時不同軟件有時會存在差別)。組分和物種通過平衡常數(shù)K進(jìn)行關(guān)聯(lián)：對于組分，lgK=0；對于物種，取其生成反應(yīng)的平衡常數(shù)。各軟件中主要的控制方程有質(zhì)量(摩爾)平衡方程、電荷(電子)平衡方程、化學(xué)反應(yīng)平衡方程、溫度校正方程和活度系數(shù)校正方程等。盡管不同軟件采用的溫度校正方程和活度系數(shù)校正方程等各不相同，但它們的基本原理均相同。

這些軟件的結(jié)構(gòu)都有一個共同的特點，即均由輸入、輸出、迭代算法及化學(xué)熱力學(xué)數(shù)據(jù)庫等幾個模塊組成。其中，輸入模塊一般有初始溶液的組分總濃度和參加反應(yīng)的固相、氣相等組成；輸出模塊則主要有各平衡組分/物種的濃度、溶液與固相和氣相間的物質(zhì)交換量等；迭代算法模塊用于非線性方程組的求解，基本上都是采用改進(jìn)的Newton-Raphson迭代法；化學(xué)熱力學(xué)數(shù)據(jù)庫模塊包括標(biāo)準(zhǔn)狀態(tài)(25 ℃，1 atm(1×105Pa))下的熱力學(xué)平衡常數(shù)、平衡常數(shù)的溫度系數(shù)或不同溫度下的值、活度系數(shù)計算所需的參數(shù)、水溶液中各物種的化學(xué)計量因子等。

2 各模擬軟件的功能比較

CHEMSPEC、EQ3/6、MINEQL+和VMINTEQ四個軟件的功能列于表1。該4個軟件主要的共同點是：(1) 4個軟件都可以進(jìn)行配合、氧化-還原、溶解-沉淀等3類基本化學(xué)反應(yīng)的模擬計算；(2) 都可以采用Davies模型進(jìn)行活度系數(shù)校正。主要的差別是：(1) MINEQL+和VMINTEQ具有處理吸附和表面配合這2類化學(xué)反應(yīng)的能力，而CHEMSPEC和EQ3/6不能進(jìn)行這2類化學(xué)反應(yīng)的模擬計算；(2) CHEMSPEC和EQ3/6提供了多種活度系數(shù)校正方法，尤其是Pitzer模型使EQ3/6能夠?qū)Ω唠x子強度的體系進(jìn)行模擬計算；(3) EQ3/6和MINEQL+可以處理離子交換反應(yīng)，EQ3/6采用Vanselow和Gapon模型，而MINEQL+則定義了1個新的反應(yīng)；(4) 4個軟件中僅有EQ3/6可以進(jìn)行反應(yīng)動力學(xué)的模擬計算。

因為各軟件不同的功能之間的差別是顯而易見的，所以本工作的重點是對各軟件的相同功能進(jìn)行比較。由于各軟件對于氧化-還原反應(yīng)的處理方式都相同，即：在所研究的體系中，電子給體和受體間處于平衡狀態(tài)，所有氧化-還原電對的電極電位都相等，由于氧化-還原電對的標(biāo)準(zhǔn)電極電位與氧化-還原反應(yīng)的平衡常數(shù)可以通過Nernst方程關(guān)聯(lián)在一起，因此各個軟件在進(jìn)行氧化-還原反應(yīng)模擬計算中產(chǎn)生的差異，與它們在進(jìn)行配合、溶解-沉淀反應(yīng)模擬計算中產(chǎn)生的差異，兩者在本質(zhì)上并沒有特別的不同。故下面的比較主要針對配合、溶解-沉淀反應(yīng)進(jìn)行。

3 4個軟件的實例計算比較

一般來講，各個軟件計算結(jié)果的差別主要來源于其數(shù)據(jù)庫的差別(各種校正被認(rèn)為是數(shù)據(jù)庫的一部分)。

表1 4個地球化學(xué)模擬軟件的功能比較

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

CHEMSPEC(V2A)采用瑞士 PSI的NAPSI_290502(260802)，共有66個組分和448個物種。

EQ3/6(V8.0)的數(shù)據(jù)庫數(shù)量比較多，來源也有些不同。不同的數(shù)據(jù)庫用3個字符的后綴進(jìn)行區(qū)分，包括500、1kb、2kb、5kb、cmp、hmw、pit、shv、skb、sup、ymp和ypf等，其中ymp和ypf是Yucca Mountain Project專用數(shù)據(jù)庫。各個數(shù)據(jù)庫所包含的組分和物種數(shù)量有比較大的差異，cmp數(shù)據(jù)庫所包含的組分和物種數(shù)量最多，如果沒有特別說明，下面的計算均采用該數(shù)據(jù)庫。能用Pitzer模型進(jìn)行活度系數(shù)校正的數(shù)據(jù)庫有3個：hmw、pit和ypf，其中hmw數(shù)據(jù)庫只適用于25 ℃，pit數(shù)據(jù)庫適用于0～100 ℃，ypf數(shù)據(jù)庫適用于0～300 ℃。能進(jìn)行壓力校正的數(shù)據(jù)庫有5個：500、1kb、2kb、5kb和shv。

MINEQL+(V4.5)采用MINTEQA2的數(shù)據(jù)庫，也有原始的MINEQL的一些數(shù)據(jù)，其數(shù)據(jù)都是標(biāo)準(zhǔn)狀態(tài)下無限稀釋溶液的數(shù)據(jù)。

VMINTEQ(V2.52)是MINTEQA2的可視版，因此其數(shù)據(jù)庫與MINTEQA2相同，也都是標(biāo)準(zhǔn)狀態(tài)下無限稀釋溶液的數(shù)據(jù)。

對于常壓下的模擬計算而言，數(shù)據(jù)庫的差別主要包括標(biāo)準(zhǔn)狀態(tài)下平衡常數(shù)、溫度校正方程以及活度系數(shù)校正方程等的差別。下面先采用單平衡點計算的方法，以比較這3種差別對計算結(jié)果的影響；然后再采用多平衡點計算的方法，以便于比較這4個軟件對某一個實際體系進(jìn)行模擬計算時，其結(jié)果隨某一個參數(shù)在一定范圍內(nèi)變化的趨勢。

3.1 單平衡點計算比較

1) 標(biāo)準(zhǔn)狀態(tài)下平衡常數(shù)的差別

用4個軟件進(jìn)行模擬計算的結(jié)果列于表2。從表2可以看出：

表2 例1的計算結(jié)果

注(Note)： -表示“沒有該項”(- means “none”)

表3 例1中的平衡常數(shù)

注(Note)： -表示“沒有該項”(- means “none”)

(4) MINEQL+和VMINTEQ的計算結(jié)果的微小差別，一方面可能是2個軟件中的次要物種有一點差別；另一方面也可能與計算結(jié)果的有效數(shù)字位數(shù)有關(guān)，MINEQL+中取4位有效數(shù)字，而VMINTEQ中取5位有效數(shù)字。

2) 溫度校正方程的差別

CHEMSPEC、MINEQL+和VMINTEQ等3個軟件的數(shù)據(jù)庫中只有25 ℃下的數(shù)據(jù)，對于其它溫度下的模擬計算，需要進(jìn)行溫度校正。

CHEMSPEC對平衡常數(shù)的溫度校正采用如下公式(A、B、C、D、E為經(jīng)驗常數(shù))：

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

(1)

若數(shù)據(jù)庫中沒有提供相應(yīng)的溫度校正系數(shù)，而提供了焓變值，則用Van’t Hoff公式校正。

MINEQL+采用Van’t Hoff公式對反應(yīng)平衡常數(shù)進(jìn)行校正。

VMINTEQ一般也采用Van’t Hoff公式對反應(yīng)平衡常數(shù)進(jìn)行溫度校正。另外，對其中25種組分，還可以采用以下經(jīng)驗公式進(jìn)行校正(A、B、C、D、E、F、G為經(jīng)驗常數(shù))：

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

ET2+F/T2+GT1/2

(2)

此經(jīng)驗公式比Van’t Hoff公式準(zhǔn)確，所以對這25種組分，VMINTEQ比MINEQL+的溫度修正更準(zhǔn)確。

由于焓變是溫度的函數(shù)，所以以上3個軟件的溫度應(yīng)用范圍比較窄，均為0～100 ℃。

EQ3/6數(shù)據(jù)庫中收集了除25 ℃以外的數(shù)據(jù)，所以其溫度適用范圍更廣，但各個數(shù)據(jù)庫對應(yīng)的溫度不同。常用的cmp數(shù)據(jù)庫適用于0～300 ℃。

例2：計算不同溫度(25～90 ℃)時Na2HPO4-KH2PO4標(biāo)準(zhǔn)緩沖溶液的pH值。

圖1 不同軟件計算結(jié)果相對誤差隨溫度的變化

以文獻(xiàn)[10-11]中的標(biāo)準(zhǔn)pH值為基準(zhǔn)，對4種軟件計算出的pH值進(jìn)行比較，其相對誤差示于圖1。從圖1可以看出，EQ3/6計算出的pH值的相對誤差最小，其絕對值始終小于1%。這是由于其數(shù)據(jù)庫中含有0～300 ℃時的經(jīng)驗常數(shù)，所以其溫度校正更準(zhǔn)確。CHEMSPEC和VMINTEQ的校正結(jié)果基本一致，在溫度較低時誤差較小，隨著溫度的增加，誤差也越來越大，如在90 ℃時，誤差達(dá)到了3.5%。上述3個軟件在25 ℃時的誤差均比較小，但MINEQL+在25 ℃時的誤差卻比較大，這是由于計算時沒有進(jìn)行活度系數(shù)校正引起的(本工作曾嘗試進(jìn)行活度系數(shù)校正，但結(jié)果沒有收斂)。

3) 活度系數(shù)校正方程的差別

對于水溶液中的帶電組分，CHEMSPEC采用Davies公式、擴展的Debye-Hückel公式和B-dot公式計算離子活度系數(shù)(根據(jù)溶液離子強度I和數(shù)據(jù)庫是否有某離子的Debye-Hückel參數(shù)a和b進(jìn)行選擇)。EQ3/6數(shù)據(jù)庫較多，各數(shù)據(jù)庫對應(yīng)的活度系數(shù)計算公式有所區(qū)別，包括Davies公式、B-dot公式和Pitzer公式。MINEQL+采用Davies公式計算活度系數(shù)。VMINTEQ采用Davies公式和修正的Debye-Hückel公式計算活度系數(shù)。Davies公式、擴展的Debye-Hückel公式和B-dot公式都只適用于離子強度I小于0.5 mol/L的情況。當(dāng)離子強度大于0.5 mol/L時，只能用Pitzer公式。

對于水溶液中的電中性組分，除MINEQL+沒有進(jìn)行活度系數(shù)校正以外，其余3個軟件都有相應(yīng)的校正公式。

CHEMSPEC采用B-dot公式計算活度系數(shù)；MINEQL+采用Davies公式計算活度系數(shù)；VMINTEQ主要采用擴展的Debye-Hückel公式計算活度系數(shù)，但CaCl+和MgCl+采用Davies公式；EQ3/6采用3種方法計算活度系數(shù)，分別為Davies公式(cmp數(shù)據(jù)庫)、B-dot公式(cmp數(shù)據(jù)庫)和Pitzer公式(hmw數(shù)據(jù)庫)。計算結(jié)果列于表4和表5。從表4和表5可以看出：

(1) 對于H+，其活度相同(pH=5)，但各軟件計算出的濃度各不相同；EQ3/6采用相同的數(shù)據(jù)庫(cmp數(shù)據(jù)庫)和不同的活度系數(shù)計算方法(Davies公式和B-dot公式)，其結(jié)果有較大差異(1.391 5×10-5和1.281 0×10-5)；采用Davies公式的結(jié)果更接近于MINEQL+的結(jié)果(1.391 5×10-5和1.410×10-5)，采用B-dot公式的結(jié)果更接近于CHEMSPEC的結(jié)果(1.281 0×10-5和1.241 1×10-5)；這說明在離子強度較大時，計算活度系數(shù)的不同方法之間的差別與稀溶液相比是顯著增加的；

(2) 物種CaCl+和MgCl+在有的數(shù)據(jù)庫有，而在另一些數(shù)據(jù)庫中卻沒有；有這兩個物種時，對應(yīng)的計算結(jié)果里游離的Mg2+和Ca2+濃度都較小；

表4 例3的計算結(jié)果(CHEMSPEC,MINEQL+和VMINTEQ)

注(Note)： -表示“沒有該項”(- means “none”)

表5 例3的計算結(jié)果(EQ3/6)

注(Note)： -表示“沒有該項”(- means “none”)

(3) 與例1相比，MINEQL+和VMINTEQ兩個軟件計算的結(jié)果之間的差距明顯增加，這是因為：一方面，例1的體系是稀溶液，而本例的離子強度較大(各軟件計算的平均值約為0.25 mol/L)，這種情況下計算活度系數(shù)的不同方法之間的差別也較大；另一方面，MINEQL+和VMINTEQ這兩個軟件的數(shù)據(jù)庫中的物種有差別，在VMINTEQ中有CaCl+和MgCl+這兩個物種，而在MINEQL+中沒有。

3.2 多平衡點計算比較

上面的單平衡點計算不能反映計算結(jié)果隨某一個參數(shù)變化的趨勢，為了從變化趨勢上對4個軟件進(jìn)行比較，因此進(jìn)行了下面的計算。

例4：計算不同pH值(pH=5～9)時總濃度為10-6mol/L的Am(Ⅲ)在北山51號井地下水(地下水組成見文獻(xiàn)[9])中的化學(xué)形態(tài)分布。

首先需要說明兩點：(1) MINEQL+的數(shù)據(jù)庫中沒有Am(Ⅲ)的數(shù)據(jù)，有關(guān)Am(Ⅲ)的數(shù)據(jù)均采用CHEMSPEC中的數(shù)據(jù)。(2) 在MINEQL+的計算過程中，pH=8.5和pH=9兩個點的計算沒有進(jìn)行活度系數(shù)的校正(本工作曾嘗試進(jìn)行活度系數(shù)校正，但結(jié)果沒有收斂)。

圖2—圖5列出了Am(Ⅲ)的幾種主要的化學(xué)形態(tài)隨體系pH值的變化趨勢，從中可以看出：

(3) 在圖2、圖3和圖5中，當(dāng)pH>6時，各物種的濃度隨著pH值的增加而逐漸減少，這主要是由于生成了AmCO3OH(cr)沉淀的緣故；另外，當(dāng)pH值達(dá)到9.5(對CHEMSPEC，pH=9)時，還會產(chǎn)生Am(OH)3(cr)沉淀；

圖2 Am3+計算濃度的比較

圖計算濃度的比較

圖計算濃度的比較

圖計算濃度的比較

3.3 計算中的其它差別

1) 數(shù)據(jù)庫中數(shù)據(jù)的增補

4個軟件中，MINEQL+增補數(shù)據(jù)最方便，可以直接在軟件用戶界面中增加組分和物種，也可以直接對物種的平衡常數(shù)進(jìn)行修改。VMINTEQ可以直接在軟件用戶界面上增加物種，并對物種的平衡常數(shù)進(jìn)行修改；但對于基本組分的增加，需要修改基本組分的數(shù)據(jù)文件comp.mdb。CHEMSPEC和EQ3/6增補數(shù)據(jù)都比較麻煩。CHEMSPEC需要對組分文件components.dtb、物種文件species.dtb和化學(xué)計量矩陣文件stoichiometry.dtb進(jìn)行修改，并需要使用相關(guān)的數(shù)據(jù)庫生成和格式轉(zhuǎn)換程序。而EQ3/6為了增補組分和物種，需要修改數(shù)據(jù)文件data0，然后通過運行EQPT程序產(chǎn)生相應(yīng)的data1文件；但如果只需要臨時改變某些物種的平衡常數(shù)，也可以直接在編輯輸入文件時進(jìn)行修改。

2) 收斂性

迭代計算的收斂性與初值的選取和迭代算法有關(guān)。通過具體實例的計算發(fā)現(xiàn)(收斂判據(jù)均采用各軟件的默認(rèn)值)，MINEQL+在計算過程中最容易發(fā)生不收斂的情況(參見例2和例4)。VMINTEQ中允許產(chǎn)生的沉淀不能超過20種；然而在例4中，若將所有飽和度大于0的物種都設(shè)為沉淀的話(沒有超過20)，計算就不能收斂，這就需要將部分最不可能生成沉淀的物種從沉淀列表中剔除掉。在本工作計算例子中，EQ3/6和CHEMSPEC沒有出現(xiàn)不收斂的情況。

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內(nèi)存)。CHEMSPEC的運行時間明顯少于其它的幾個軟件，這可能與其在算法方面的改進(jìn)有關(guān)[9]。

4) 單位換算

EQ3/6的輸入與輸出文件中以mol/g或mg/kg溶液為單位，而其它3個軟件都以mol/L或mg/L為單位，這就需要進(jìn)行換算。如果已知溶液的密度，則可以很容易地進(jìn)行換算；如果溶液的密度是未知數(shù)，則可以假定其密度為1 kg/L。對于稀溶液來說這種假定誤差較??；但當(dāng)溶液不是足夠稀時，會產(chǎn)生一定的誤差。

4 小 結(jié)

各個地球化學(xué)模擬軟件具有相同的基本功能，均可進(jìn)行配合、溶解-沉淀、氧化-還原等反應(yīng)的模擬計算。各軟件有各自的特點，CHEMSPEC在算法上有所改進(jìn)，更有利于計算的收斂；EQ3/6有眾多的數(shù)據(jù)庫，能適用于高溫、高壓及高離子強度下的計算；MINEQL+很容易對數(shù)據(jù)庫中的數(shù)據(jù)進(jìn)行增補；VMINTEQ具有強大的處理吸附和表面配合反應(yīng)的能力。

計算結(jié)果的差別，主要來源于數(shù)據(jù)庫中數(shù)據(jù)的不同。對于常壓下的模擬計算而言，數(shù)據(jù)庫中標(biāo)準(zhǔn)狀態(tài)下平衡常數(shù)、溫度校正方程以及活度系數(shù)校正方程等的差別，都會引起計算結(jié)果的偏差。對標(biāo)準(zhǔn)狀態(tài)下的稀溶液，各個軟件的計算結(jié)果基本一致。

[1] 王 駒,徐國慶,鄭華鈴，等.中國高放廢物地質(zhì)處置研究進(jìn)展:1985-2004[J].世界核地質(zhì)科學(xué),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] 王祥云,陳 濤,劉春立.化學(xué)形態(tài)分析軟件CHEMSPEC及其應(yīng)用[J].中國科學(xué)B輯-化學(xué),2009,39(11):1-12.

[10] 顧慶超主編.新編化學(xué)用表[M].南京:江蘇教育出版社,1998:903.

[11] 常文保,李克安.簡明分析化學(xué)手冊[M].北京:北京大學(xué)出版社,1981:259.