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

劉 鑫，陸現彩

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

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

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

1 模型和方法

1.1 假三元模型

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

x1+x2+x3= 1

(1)

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

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

PA=x3/2+x1

(2)

PB=x3/2+x2

(3)

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

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

根據亞規則固溶體的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能量參數，W可解釋為有序態的焓相對于兩個端員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)

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

Hmix=x1x2(x1W121+x2W122)+x1x3

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

Kx1x2x3W123

(17)

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

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

x2x3(x2W232+x3W233)+x3ΔHord

(18)

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

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

(19)

這樣，結合式(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)

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

1.2 Margules參數的計算

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

Hmix=x1x2(x1W21+x2W12)

(24)

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

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

在透輝石(diopside， CaMgSi2O6)-硬玉(jadeite， NaAlSi2O6)這種典型的耦合替代型的礦物固溶體中，存在兩種陽離子位置(圖2)，較小的M1位置可以被Mg2+和Al3+占據，較大的M2位置可以被Ca2+和Na+占據，具有C2/c空間群。透輝石的一對Ca2+和Mg2+可以被一對Na+和Al3+替代而成為硬玉。該體系中存在中間相礦物綠輝石(omphacite， Ca0.5Mg0.5Na0.5Al0.5Si2O6)，在低溫下具有典型的有序結構，也就是說 M1和M2可以被細分為M11、M12、M21和M22等4種亞位置，依次被Mg2+、Al3+、Ca2+和Na+完全占據，具有P2/n空間群?？紤]到發生了耦合替代，這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 透輝石-硬玉固溶體的晶胞結構(藍色為Si—O四面體，淺藍色為M1八面體位置，可以被Mg2+和Al3+占據，黃色為M2多面體位置，可以被Ca2+和Na+占據)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+)

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

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

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

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

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

2 結果和討論

2.1 混合熱力學性質

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

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

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

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

(25)

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

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

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

(26)

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

(27)

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

Svib=x1x2(x1S1+x2S2)

(28)

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

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

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

綠輝石的有序度隨溫度升高而發生顯著的變化(圖8)，有序-無序相變發生在1 148±25 K左右(虛線處)，左右兩側分別為P2/n相和C2/c相。結合相圖推斷此處相變應為一級相變。相變溫度與前人實驗結果1 138±10 K(Carpenter， 1982)、數學分析結果1 123 K (Greenetal.， 2007)和計算模擬結果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)兩種觀點。但是，由于本研究使用的模型相對簡單，僅用一個參數表征LRO和SRO。低溫下LRO和SRO共存，高溫下LRO消失，但是SRO還存在，這也是圖8中的有序度不變為零的原因。

圖 4 透輝石-硬玉固溶體的混合焓(a)， 構型熵(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 討論

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

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

3 結論

(1) 通過插入有序綠輝石作為虛擬第三端員，利用假三元模型揭示了透輝石-硬玉固溶體的熱力學性質，發現綠輝石的有序-無序相變為一級相變，相變溫度約為1 148± 25 K。研究獲得的有序度變化與前人的模擬和分析結果吻合較好。

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

致謝本文受南京大學優秀博士研究生創新能力提升計劃B資助。