基于假三元模型的透輝石-硬玉固溶體熱力學(xué)性質(zhì)計算

劉 鑫，陸現(xiàn)彩

(內(nèi)生金屬礦床成礦機制研究國家重點實驗室，南京大學(xué) 地球科學(xué)與工程學(xué)院， 江蘇 南京 210023)

礦物的晶體化學(xué)特征反映了其形成的物理化學(xué)條件，礦物固溶體更是記錄了豐富的成因信息，綜合分析礦物的晶體化學(xué)變化和熱力學(xué)性質(zhì)，可重建成巖成礦過程的物理化學(xué)條件。在變質(zhì)巖石學(xué)研究中，目前有兩種主要方法估算變質(zhì)作用的溫度和壓力條件：傳統(tǒng)礦物溫壓計和視剖面圖溫壓計(Powell and Holland， 2008； 吳春明等， 2013； 吳佳林等， 2015)。與傳統(tǒng)礦物溫壓計僅僅利用實驗標定的礦物對或組合推測溫度壓力不同，視剖面圖方法以熱力學(xué)基本方程為基礎(chǔ)，利用礦物固溶體熱力學(xué)數(shù)據(jù)庫對特定全巖成分進行視剖面圖計算，可以得到礦物共生組合、礦物相豐度、成分以及溫度壓力條件等豐富的信息，是目前研究變質(zhì)作用物理化學(xué)條件和精細過程的最佳方法之一(Powelletal.， 1998； Powell and Holland， 2008； 魏春景等， 2009； 吳佳林等， 2015)。Powell等(1998)開發(fā)了用于計算視剖面圖的Thermocalc程序并建立了熱力學(xué)數(shù)據(jù)庫，前人先后提出了多種成分-活度模型，并用于變質(zhì)作用的精細研究。但是由于地質(zhì)過程的復(fù)雜性和固溶體礦物的多樣性，現(xiàn)有礦物熱力學(xué)數(shù)據(jù)庫遠遠不夠完善。研究礦物固溶體熱力學(xué)性質(zhì)的傳統(tǒng)方法主要依靠化學(xué)分析和量熱學(xué)，基于有限的不連續(xù)實驗數(shù)據(jù)擬合熱力學(xué)曲線，固有的實驗誤差和廣泛存在的動力學(xué)效應(yīng)會顯著降低熱力學(xué)參數(shù)的準確性。鑒于實驗方法的局限性，過去幾十年內(nèi)建立了多種礦物固溶體熱力學(xué)性質(zhì)研究的理論計算方法，然而這些方法建立的礦物固溶體模型很難同時考慮長程有序效應(yīng)(long-range order, LRO)和短程有序效應(yīng)(short-range order, SRO)。為了彌補這一不足，往往需要知道足夠多的原子間相互作用能，并且結(jié)合原子分布獲得固溶體的最低自由能，因而需要海量的構(gòu)型計算，而巨大的計算量和高昂的計算成本大大制約了這些理論方法的應(yīng)用和推廣。

本文介紹了一種全新的用于礦物固溶體熱力學(xué)性質(zhì)計算的假三元模型，該模型同時考慮了LRO和SRO效應(yīng)，僅僅通過少量構(gòu)型計算即可獲得精細的固溶體熱力學(xué)數(shù)據(jù)，能夠準確預(yù)測礦物相變邊界。本文將該方法用于研究榴輝巖中透輝石-硬玉固溶體系列的熱力學(xué)性質(zhì)。榴輝巖是一種標志性高級變質(zhì)巖，具有鮮明的構(gòu)造成因意義(陳意等， 2005； 魏春景等， 2009； 張麗娟等， 2018)，分析榴輝巖中的礦物固溶體體系的晶體化學(xué)特征，可利用視剖面圖方法重建榴輝巖形成的溫度壓力條件(Greenetal.， 2007； 魏春景等， 2009)，為地球動力學(xué)和構(gòu)造環(huán)境研究提供關(guān)鍵數(shù)據(jù)。

1 模型和方法

1.1 假三元模型

將一種固溶體礦物的分子式寫為MR，其晶格中有數(shù)量相同的兩種陽離子位點(或者人為地分為兩種)，記為α和β位點，都可被M1和M2兩種陽離子占據(jù)。假設(shè)一個固溶體具有這種結(jié)構(gòu)，并且包含A和B兩種陽離子，則兩個端員礦物為AR和BR，其中A和B都可以占據(jù) α和β兩種位點。可在二元固溶體端員(AαAβ)R2和(BαBβ )R2之間考慮一個(AαBβ)R2端員作為中間有序態(tài)(圖1)，這個相可以是自然存在的礦物，也可以僅僅是一種有序態(tài)。假如3個端員組分的濃度分別標記為x1、x2和x3，則存在：

x1+x2+x3= 1

(1)

圖 1 二元固溶體的假三元模型Fig. 1 Fictive ternary model for binary solid solution

由于僅存在2個獨立端員，引入的(AαBβ)R2端員是為了在規(guī)則固溶體模型中引入LRO參數(shù)Q，該組分也可以使用兩端員摩爾分數(shù)PA和PB和分子式APABPBR 表示。那么，(x1,x2,x3)和(PA,PB)的關(guān)系如下：

PA=x3/2+x1

(2)

PB=x3/2+x2

(3)

如果PAα和PAβ分別是兩種陽離子位點A原子出現(xiàn)的概率(B原子同理)，應(yīng)同時滿足：

PAα=x1+x3

(4)

PAβ=x1

(5)

PBα=1-PAα=x2

(6)

PBβ=1-PAβ=x2+x3

(7)

Q=(PAα-PAβ)/(PAα+PAβ) =x3/2PA(PB≥0.5)

(8)

Q=(PBβ-PBα)/(PBα+PBβ)=x3/2PB(PB≤0.5)

(9)

有序度Q的變化范圍是[0, 1]。Q=1意味著某個組分固溶體的最有序狀態(tài)，其組分組成只能沿著圖1所示三角形的2個腰邊[分別為(AαAβ)R2-(AαBβ)R2連接線和(BαBβ )R2-( AαBβ )R2連接線]變化，A和B原子在α和β位置的占據(jù)上具有明顯的選擇性，總是盡可能先占據(jù)滿其中一種位置。Q=0意味著PAα=PAβ和PBβ=PBα，也就是說，A和B原子在α和β位置完全無序分布，組分沿著圖1三角形的底邊[(AαAβ )R2-(BαBβ )R2兩個端員的連接線]變化。x=0.5和Q=1則對應(yīng)完全有序中間態(tài)。

根據(jù)亞規(guī)則固溶體的Bragg-Williams(BW)模型(Bragg and Williams， 1934， 1935)，固溶體的混合焓可以表示如下：

Hmix=PAαPBα(PAαWAα+PBαWBα)+PAβPBβ(PAβWAβ+
PBβWBβ)+PAαPBβWαβ+PAβPBαWβα

(10)

這里的W都是Margules能量參數(shù)，W可解釋為有序態(tài)的焓相對于兩個端員AR和BR的機械混合焓的增量，即Wαβ=ΔHord。將(4)～(7)式代入(10)式后化簡可得：

Hmix=x1x2[x1(WAα+WAβ+Wαβ+Wβα)+x2(WBα+
WBα+Wαβ+Wβα)]+x1x3(x1WAβ+x3WBβ)+x2x3
(x2WBα+x3WAα)+x3Wαβ+x1x2x3(2WAα+2WBβ+
Wαβ+Wβα)

(11)

同時，定義以下等式：

WAα+WAβ+Wαβ+Wβα=W121

(12)

WBα+WBβ+Wαβ+Wβα=W121

(13)

WAB=W131，WBβ=W133，Wβα=W232，WAα=W233

(14)

2WAα+2WBβ+Wαβ+Wβα=W123

(15)

這樣，式(11)可以改寫為：

Hmix=x1x2(x1W121+x2W122)+x1x3(x1W131+

x3W133)+x2x3(x2W232+x3W233)+x3ΔHord+

x1x2x3W123

(16)

這就是三元體系的亞規(guī)則固溶體模型混合焓表達式，適用于熱力學(xué)性質(zhì)不對稱的體系。研究表明(Liuetal., 2019)，該表達式會顯著高估混合焓和有序-無序相變邊界，可對最后一項引入校正因子K解決：

Hmix=x1x2(x1W121+x2W122)+x1x3

(x1W131+x3W133)+x2x3(x2W232+x3W233)+x3ΔHord+

Kx1x2x3W123

(17)

值得說明的是，K因子的引入，并不會影響任意二元體系(x1、x2和x3其中一個為0，對應(yīng)于圖1三角形的3條邊)的焓效應(yīng)，這也是引入該因子的前提。當K從1變化到0時，基于該展開式的計算結(jié)果與實驗結(jié)果吻合度越來越高，根本原因在于BW模型沒有考慮短程有序效應(yīng)，因子K的引入和最后一項的消失，相當于對SRO效應(yīng)校正，混合焓展開式為：

Hmix=x1x2(x1W121+x2W122)+x1x3(x1W131+x3W133)

x2x3(x2W232+x3W233)+x3ΔHord

(18)

值得注意的是，原本設(shè)定的LRO參數(shù)Q在經(jīng)過校正之后攜帶了SRO效應(yīng)信息，因此本方法的Q綜合反映了LRO和SRO效應(yīng)，很難將兩者完全區(qū)分開。最后，BW模型下的每摩爾可交換原子的構(gòu)型熵應(yīng)該采取以下形式：

Sconf=-R(PAαlnPAα+PBαlnPBα)/2-
R(PAβlnPAβ+PBβlnPBβ)/2

(19)

這樣，結(jié)合式(18)和(19)，可以寫出自由能的表達式：

Gmix(Q,x,T) =Hmix-TSconf

(20)

式(20)表明，當溫度和組分固定時，一定存在有序度Qeq，使得自由能最小。這個最小自由能即為固溶體的混合自由能。因此用Qeq表示的平衡自由能、熵和焓可用以下方程表示：

Gmix(x,T)=Geq(Qeq,x,T) =min[G(Q,x,T)]

(21)

Smix(x,T)=Seq(x,T)

(22)

Hmix(x,T)=Gmix(x,T)+TSmix(x,T)

(23)

至此，礦物固溶體的所有混合熱力學(xué)性質(zhì)均可通過組分、溫度、有序度Q和Margules參數(shù)計算，其中Margules參數(shù)能夠使用經(jīng)驗勢或第一性原理方法計算。

1.2 Margules參數(shù)的計算

在計算中需要6個Margules參數(shù)：W121、W122、W131、W133、W232和W233，可以通過單缺陷構(gòu)型方法(single-defect method, SDM)和準隨機構(gòu)型計算(quasi-random Structure, QRS)兩種策略(Sluiter and Kawazoe， 2002； Vinogradetal.， 2006； Vinogradetal.， 2013)獲得。所謂的隨機結(jié)構(gòu)就是離子自由排布的構(gòu)型，這對于非常大的晶胞很容易實現(xiàn)。然而，在超胞大小受限的計算中，獲得完全隨機的構(gòu)型異常困難，只能用盡可能隨機的構(gòu)型代替，也就是準隨機構(gòu)型，這些準隨機構(gòu)型中可能存在短程有序現(xiàn)象。針對特定組分，計算若干種準隨機構(gòu)型的能量，求其平均后即可用亞規(guī)則固溶體模型(式24)的焓擬合Margules參數(shù)。

Hmix=x1x2(x1W21+x2W12)

(24)

此外，還需要知道ΔHord，這僅需要計算中間組分的有序結(jié)構(gòu)相對于端員礦物的機械混合焓的增量就可以獲得。

1.3 透輝石-硬玉固溶體系列的計算

在透輝石(diopside， CaMgSi2O6)-硬玉(jadeite， NaAlSi2O6)這種典型的耦合替代型的礦物固溶體中，存在兩種陽離子位置(圖2)，較小的M1位置可以被Mg2+和Al3+占據(jù)，較大的M2位置可以被Ca2+和Na+占據(jù)，具有C2/c空間群。透輝石的一對Ca2+和Mg2+可以被一對Na+和Al3+替代而成為硬玉。該體系中存在中間相礦物綠輝石(omphacite， Ca0.5Mg0.5Na0.5Al0.5Si2O6)，在低溫下具有典型的有序結(jié)構(gòu)，也就是說 M1和M2可以被細分為M11、M12、M21和M22等4種亞位置，依次被Mg2+、Al3+、Ca2+和Na+完全占據(jù)，具有P2/n空間群。考慮到發(fā)生了耦合替代，這4個亞位點可以被重新分為兩類，一類包括M11和M21，另一類包括M12和M22，前者可以組合標記為α位置，后者可以組合標記為β位置。這樣，綠輝石的分子式可寫作[Ca0.5Mg0.5]α[Na0.5Al0.5]βSi2O6，作為透輝石(端員1)和硬玉(端員2)的中間相(端員3)，從而建立了該體系的假三元模型。在高溫下，綠輝石晶胞的M11、M12、M21和M22等4種亞位置兩兩等價，可以看做只有M1和M2兩種陽離子位置，具有C2/c空間群。

圖 2 透輝石-硬玉固溶體的晶胞結(jié)構(gòu)(藍色為Si—O四面體，淺藍色為M1八面體位置，可以被Mg2+和Al3+占據(jù)，黃色為M2多面體位置，可以被Ca2+和Na+占據(jù))Fig. 2 Unit cell of diopside-jadeite solid solution (Si—O tetrahedrons are blue, M1 octahedrons are light blue which can be occupied by Mg2+ and Al3+; M2 polyhedrons are yellow, which can be occupied by Ca2+ and Na+)

由于α和β位置為組合所得，只能利用準隨機結(jié)構(gòu)來獲得Margules參數(shù)。因此，除純的硬玉、透輝石和綠輝石晶胞外，還需建立一系列的準隨機構(gòu)型。在硬玉-透輝石、硬玉-綠輝石和綠輝石-透輝石3個二元固溶體路徑上分別選擇7、3 和3 個組分點進行計算，計算使用6 × 3 × 3硬玉結(jié)構(gòu)超胞。計算使用GULP 軟件包(Gale， 1997, 2005； Gale and Rohl， 2003)，使用的經(jīng)驗勢如表1所示(Vinogradetal.， 2007)。為了降低誤差，每個組分均選擇10個準隨機構(gòu)型進行計算，所有構(gòu)型能量減去機械混合能后求平均值，再通過擬合亞規(guī)則固溶體模型的焓獲得Margules參數(shù)(表2)，Margules參數(shù)擬合效果如圖3。

表1 透輝石-硬玉的經(jīng)驗勢模型Table 1 Force field parameters for diopside-jadeite system

注：*所有Buckingham勢的截斷半徑都是12 ?。

表 2 透輝石-硬玉固溶體的Margules參數(shù)Table 2 Margules parameters for diopside-jadeite solid solution

圖 3 透輝石-硬玉體系的Margules參數(shù)圖解Fig. 3 Margules parameters for diopside-jadeite solid solution

2 結(jié)果和討論

2.1 混合熱力學(xué)性質(zhì)

透輝石-硬玉體系混合焓等溫線(圖4a)具有微弱的不對稱性，中間組分的巨大低谷是由于結(jié)構(gòu)中部分有序造成的，而且主要體現(xiàn)為LRO。該體系的構(gòu)型熵等溫線(圖4b)顯示，273 K下的綠輝石的熵為0，隨著溫度的升高，熵逐漸增加直到完全無序態(tài)的最大值，中間組分的構(gòu)型熵低谷在1 148 K左右基本消失，對應(yīng)混合焓的迅速增大，指示了有序-無序相變的發(fā)生。在較低溫度下，混合自由能等溫線(圖4c)被有序綠輝石組分分隔成兩部分，此時綠輝石的自由能為負值；隨溫度升高，等溫線逐漸變平滑，且所有組分的自由能均具負值?；旌献杂赡艿淖兓欢ǔ潭壬戏从沉似浞€(wěn)定性，但對穩(wěn)定性的準確判斷還需公切線分析的相圖來說明。

2.2 超額自由能和溫度組分相圖

在物理化學(xué)研究中，通常用超額熱力學(xué)性質(zhì)表示實際固溶體和理想固溶體的熱力學(xué)性質(zhì)的差值，超額自由能(GE)是應(yīng)用最廣的一個熱力學(xué)參數(shù)，其表達式為：

+(1-x) ln(1-x)]

(25)

透輝石-硬玉固溶體的超額自由能(圖5)顯示：在1 123 K以下，在綠輝石相的組分附近存在超額自由能的極小值，但是當溫度高于1 123 K時，極小值迅速消失。

除超額自由能外，礦物固溶體(A,B)R某端員組分的固相活度系數(shù)γ是另外一種常用的熱力學(xué)參數(shù)，與固溶體超額自由能存在如下關(guān)系(Redlich and Kister， 1948； Plummer and Busenberg， 1987)：

lnγA=(GE+xB?GE/?xA)/RT

(26)

lnγB=(GE+xA?GE/?xB)/RT

(27)

盡管有研究認為，晶格內(nèi)的原子振動對自由能的貢獻可以忽略(Bennyetal.， 2009； Ruiz-Hernandezetal.， 2010； Wangetal.， 2011)，但亦有相反的觀點(Vinogradetal.， 2007； Liuetal.， 2016)。本研究采用了Vinograd 等(2007)對振動效應(yīng)的校正公式：

Svib=x1x2(x1S1+x2S2)

(28)

其中，參數(shù)S1= 2.5 J/mol·K，S2= 4.5 J/mol·K。經(jīng)過修正后的該體系的固相活度見圖6。

基于自由能等溫線，通過公切線分析可得透輝石-硬玉固溶體的溫度-組分相圖(圖7)。對自由能考慮振動效應(yīng)修正之后(虛線)，相圖的溫度邊界有所降低，透輝石一側(cè)降低100 K左右，而硬玉一側(cè)降低200～300 K。同時可見兩個較大的固溶隙，由中間相綠輝石分隔開。富透輝石一側(cè)的固溶隙要比富硬玉一側(cè)略大，相邊界溫度(～1 623 K)也更高(～1 473 K)，這一認識與前人研究結(jié)果一致(Greenetal.， 2007； Vinogradetal.， 2007)。

2.3 綠輝石的有序-無序相變

綠輝石的有序度隨溫度升高而發(fā)生顯著的變化(圖8)，有序-無序相變發(fā)生在1 148±25 K左右(虛線處)，左右兩側(cè)分別為P2/n相和C2/c相。結(jié)合相圖推斷此處相變應(yīng)為一級相變。相變溫度與前人實驗結(jié)果1 138±10 K(Carpenter， 1982)、數(shù)學(xué)分析結(jié)果1 123 K (Greenetal.， 2007)和計算模擬結(jié)果1 150±20 K(Vinogradetal.， 2007) 非常相似。需要說明的是，很多研究分別討論M1和M2兩個位置的有序度QM1和QM2， 分別提出QM1=QM2(Burton and Davidson，1988； Vinograd， 2002a， 2002b)和QM1=2QM2(Rossietal.， 1983； Vinogradetal.， 2007)兩種觀點。但是，由于本研究使用的模型相對簡單，僅用一個參數(shù)表征LRO和SRO。低溫下LRO和SRO共存，高溫下LRO消失，但是SRO還存在，這也是圖8中的有序度不變?yōu)榱愕脑颉?/p>

圖 4 透輝石-硬玉固溶體的混合焓(a)， 構(gòu)型熵(b)和混合自由能(c)等溫線(溫度范圍為273 K到2 073 K，間隔100 K)Fig. 4 Isotherms of enthalpy (a), entropy (b) and free energy of mixing (c) for diopside-jadeite solid solution (The temperature is from 273 K to 2 073 K with a interval of 100 K)

圖 5 透輝石-硬玉體系的超額自由能Fig. 5 The excess free energy for diopside-jadeite solid solution

圖 6 透輝石-硬玉體系的固相活度(實線和虛線分別表示硬玉和透輝石的活度，間隔100 K)Fig. 6 Activity plots for diopside-jadeite solid solution (The solid lines and dash lines show the activity of jadeite and diopside respectively with an interval of 100 K)

圖 7 透輝石-硬玉固溶體的溫度-組分相圖(實線和虛線分別表示考慮振動貢獻修正前后的相邊界)Fig. 7 Temperature-composition diagram for diopside-jadeite solid solution(The solid lines and dash lines show the bounda-ries without and with vibrational correction respectively)

圖 8 綠輝石的有序度隨溫度的變化Fig. 8 Illustration of LRO parameter varying with temperature

2.4 討論

本研究獲得的透輝石-硬玉固溶體的熱力學(xué)性質(zhì)數(shù)據(jù)可以用于相關(guān)變質(zhì)作用的熱力學(xué)研究，尤其適合MORB成分的榴輝巖相變質(zhì)作用研究。事實上，研究榴輝巖相變質(zhì)作用需要多種礦物固溶體的熱力學(xué)數(shù)據(jù)，前人已經(jīng)通過實驗和計算的方法得到了石榴子石(Weietal.， 2004)、綠泥石(Hollandetal.， 1998)、單斜輝石(Greenetal.， 2007)、角閃石(Dieneretal.， 2007)、斜長石(Holland and Powell， 1996)等礦物固溶體的熱力學(xué)參數(shù)，但準確性、可靠性和系統(tǒng)性存在顯著差異，本文得到的透輝石-硬玉體系的熱力學(xué)數(shù)據(jù)較為系統(tǒng)，并且充分考慮了兩種有序度，端員組分的參數(shù)與前人的數(shù)據(jù)高度一致，表明可有效補充當前只存在端員組分熱力學(xué)數(shù)據(jù)的不足，適用于成分更為復(fù)雜含輝石-硬玉體系變質(zhì)巖的成因研究。

同時，本研究還驗證了假三元模型的有效性，該方法具有計算規(guī)模小、理論模型簡單的特點，可適用于其他相似替代型固溶體系列礦物的研究。這一方法的應(yīng)用，可以解決當前熱力學(xué)數(shù)據(jù)不夠系統(tǒng)的現(xiàn)實問題。由于可以避免動力學(xué)因素的干擾和實驗條件和技術(shù)的限制，這一方法至少可以部分替代實驗方法的運用，作為實驗數(shù)據(jù)的有效補充，為巖石成因和礦物成因研究提供數(shù)據(jù)支撐。

3 結(jié)論

(1) 通過插入有序綠輝石作為虛擬第三端員，利用假三元模型揭示了透輝石-硬玉固溶體的熱力學(xué)性質(zhì)，發(fā)現(xiàn)綠輝石的有序-無序相變?yōu)橐患壪嘧?，相變溫度約為1 148± 25 K。研究獲得的有序度變化與前人的模擬和分析結(jié)果吻合較好。

(2) 基于計算模擬結(jié)果和實驗結(jié)果的對比，證實了假三元模型用于二元固溶體熱力學(xué)性質(zhì)研究的有效性，并從一定程度上解決固溶體短程有序效應(yīng)無法衡量的問題。利用準隨機結(jié)構(gòu)計算Margules參數(shù)和假三元模型可以快速計算二元固溶體和耦合置換固溶體的熱力學(xué)性質(zhì)。

致謝本文受南京大學(xué)優(yōu)秀博士研究生創(chuàng)新能力提升計劃B資助。