什么控制了矿物中Ca同位素平衡分馏系数？

钙(Ca)

钙(Ca)，构成生物骨骼最重要的元素，虽然没能成为硅酸盐矿物的骨架，却也凭借着极高的丰度（地壳中丰度排第五），与活跃的地球化学性质，成为地球化学研究的关注对象。近二十年来，随着质谱技术的发展，我们可以精确测量Ca同位素的组成，并在火成岩中发现了近2‰ 的Ca同位素组成变化，这展现出Ca同位素在研究地球及类地行星岩浆演化过程的重要潜力。

前人提出了多种可能的机制来解释火成岩中巨大的钙同位素分馏，包括：(1)部分熔融与岩浆分异；(2)矿物间的同位素平衡分馏；(3)受外来物质的交代；(4)扩散或化学反应过程中的同位素动力学分馏。要在自然样品中区分这些过程，就必须获得矿物间的Ca同位素平衡分馏系数。但是，要如何获得这一参数呢？

平衡分馏系数的获得方法

一般来说，有三种方法来获得同位素平衡分馏系数，分别是：

1. 观测处于平衡状态下的自然样品。这一方法需要准确判断所测量的自然样品达到了平衡状态，而且很难获得同位素平衡分馏系数随温度、压力的变化关系。

2. 实验岩石学模拟。但这种方法主要受到实验条件的限制，一些高温、高压的条件在实验室中很难达到，且矿物间达到平衡往往需要花费较长的时间。

3. 理论计算。利用第一性原理计算可以精确地计算矿物的同位素平衡分馏系数。这种方法不引入经验参数，可以模拟较大范围内的温压条件，广泛应用于多种同位素体系的计算中。

那么，如何计算呢？

导致同位素平衡分馏的本质是同一位点上由不同同位素占据时能量存在细微差异。Urey公式给出了定量计算同位素平衡分馏系数的方法，即先计算“简约配分函数比”（RPFR），常被称为β值：

该公式针对晶体，即对任一矿物，只要获得了该矿物中同位素交换前后的振动频率v , 带入公式，即可获得该矿物的简约配分函数比，通常也将其作为矿物与自由原子之间的同位素平衡分馏系数。在获得简约配分函数比之后，两矿物之间的平衡分馏系数就是两者的β值之比，用对数形式表示即为：

 

Δ^44/40Ca_x-y (=δ^44/40Ca_x-δ^44/40Ca_y, x和y代表不同的矿物)就是处于平衡状态下矿物间的同位素组成差异。公式中所需要的振动频率可以通过第一性原理计算的方法计算得到。

“静态” 方法与“非静态”方法

获得频率之后还存在一个问题。对于任何一个给定的结构，都可以计算出相应的振动频率，并带入Urey公式得到其简约配分函数比。问题是，矿物的结构受到温度和压力的控制。一般而言，压力升高会使得矿物的体积压缩，而温度升高则会使得其体积增加。

文献中一般有两种方法给出矿物的简约配分函数比：

（1）第一种是计算矿物在0K和特定压力（一般是0GPa）下的结构，获得的频率直接带入Urey公式，得到简约配分函数比。这种方法忽略了温度对矿物结构的影响，我们称这种方法为“静态”方法，相应的压力和简约配分函数比为静态压力、静态简约配分函数比（static-RPFRs）。

（2）第二种方法则考虑温度对体积的影响。计算一系列体积下矿物的频率和自由能，从而得到矿物的状态方程P(V,T)，再据此将以体积和温度为自变量表示的简约配分函数比β(V,T)转化为压力和温度的函数β(P,T)。这种方法我们称为“非静态”方法。

显然，第二种方法更为严谨准确，但由于第二种方法需要更大的计算量，目前Ca同位素平衡分馏系数的计算多采用第一种方法。我们的工作也对比了这两种方法得出结果的差异，将在后面详细讨论。

矿物间的同位素平衡分馏系数

利用第一性原理计算，已经陆续报道了多种矿物间的Ca同位素平衡分馏系数，如橄榄石、辉石、石榴石、长石与硫化钙等[2-8]，发现在大于1000K的高温下，矿物间也可以发生显著的Ca同位素平衡分馏。由于目前还没有关于同位素平衡分馏系数的实验报道，这些数据为解释岩浆与变质过程中的Ca同位素分馏行为提供了重要的指导。然而，目前还缺乏一些重要的含Ca矿物数据，同时影响矿物间Ca同位素平衡分馏系数的控制因素还不太清楚。

我们使用静态方法计算了一些矿物在静态零压下的简约配分函数比，包括角闪石(钾闪石 (KNaCaMg₅Si₈O₂₂F₂)和透闪石 

(Ca₂Mg₅[Si₄O₁₁]₂F₂))、黄长石 (镁黄长石 (Ca₂Mg[Si₂O₇])和铝黄长石 (Ca₂Al[AlSiO₇])), 含钠辉石 (硬玉 (Na, Ca, Mg, Al)Si₂O₆), 富钾碳酸盐 (butschliite (K₂Ca[CO₃]₂)) 和其他含钙矿物 (磷灰石 (Ca₅[PO₄]₃F), 硬石膏 (CaSO₄), 钙钛矿 (CaTiO₃) 及萤石 (CaF₂)). 此外，为探究压力对平衡分馏系数的控制，我们还利用非静态方法计算了透辉石和钙长石在0-10GPa范围内的结果。

我们将所有的静态计算结果（包括文献中的结果）汇总在了图1中，方便大家查询，图中使用的数据表格也作为附录材料上传到了论文中。在解读这些数据之前，我们必须提到一个第一性原理计算中的细节——“泛函”。

 

图1：矿物1000lnβ值汇总

不同泛函的计算结果

泛函，这里指的是第一性原理计算中使用的“交换关联泛函”，对于其具体的物理意义，大家可以自行查阅量子化学计算的书籍内容。在Ca同位素平衡分馏系数的计算中，经常使用两种泛函，分别是LDA和GGA。而GGA泛函也分为几种，其中PBE是最常用的。

根据第一性原理计算的一般经验，选取不同泛函的计算结果存在显著差异。尤其是对体积的预测，一般而言LDA预测的静态体积要比实验结果偏小，但是扣除热膨胀效应后这一偏差可以在一定程度上矫正。而GGA预测的静态体积要比实验结果偏大。这种差异也体现在简约配分函数比的计算上，LDA的计算结果显著不同于GGA的计算结果。因此，在计算平衡分馏系数时，不同方法得到的简约配分函数比不能直接进行比较。

幸运的是，不同泛函产生的差偏是系统性的。如图2所示，两种不同方法得到的1000lnβ值之间的偏差接近于常数，即GGA=LDA-0.20 (1000K)。因此，在计算平衡分馏系数时，这两种方法带来的偏差可以被经验性扣除。例如，1000K下，两种方法得到的1000lnα_{apatite-fluorite}(-0.13 ‰ vs -0.12‰）和1000lnα_{anorthite-grossular}(-0.68 ‰ vs -0.71‰)是一致的。

 

图2 两种泛函计算结果对比

（1000lnβ at 1000K）

为了便于对比，我们将GGA方法得到的结果进行矫正后汇入了图1中，矫正公式为GGA = LDA - (0.20 ± 0.05)×(10⁶/T²), 这一公式是拟合现有的少数几个数据得到的，如果后续有新的研究报道GGA 结果，建议重新拟合后再使用。矫正之后，不同方法得到的结果便可以很方便地进行对比。我们发现，不同类型矿物的简约配分函数比具有很大的差异，1000K下，1000lnβ最小为0.96 (oldhamite)，最大可达2.36 (forsterite)。此外，我们还可以探讨控制Ca同位素平衡分馏系数的因素。

力常数——最直接的控制因素

 

在高温下(v in cm⁻¹ < 1.39 T in K)，Urey公式可以简化为高温近似形式。这一高温近似式告诉我们，简约配分函数比取决于两同位素的相对质量差异，与温度成平方反比，且受原子平均力常数控制。从图3也可以看出Ca的简约配分函数比与原子力常数呈很好的线性关系。所以对于给定的同位素体系，力常数是简约配分函数比最直接的控制因素。

 

这一高温近似式告诉我们，简约配分函数比取决于两同位素的相对质量差异，与温度成平方反比，且受原子平均力常数控制。从图3也可以看出Ca的简约配分函数比与原子力常数呈很好的线性关系。所以对于给定的同位素体系，力常数是简约配分函数比最直接的控制因素。

 

图3. 1000lnβ 随力常数的变化

然而，力常数这一物理量无法通过简单的测量和计算获取。在实际的地质示踪过程中，如果没有平衡分馏系数时，我们往往需要容易获取的晶体化学特征来粗略判断平衡分馏的尺度与方向。因此我们需要了解晶体化学是如何通过影响力常数来控制简约配分函数比的。

键长与配位数

键强是力常数的直接体现，受键长和配位数的影响。更短的键长、更高的价态以及更少的配位数，意味着更强的键强，更易富集重的同位素。自然界中，Ca只有一个化合价，因此，在没有平衡分馏系数报道情况下，常通过键长和配位数来粗略判断同位素的相对富集顺序。我们发现了RPFR对键长具有很好的负相关关系，但是其对配位数（CN）的相关性却很差（图4 a,b）例如，镁黄长石和钙铝榴石都是8配位的，单从配位数判断两者应该具有相似的钙同位素组成。然而，根据我们的计算结果，Δ^44/40Ca_{grossular-akermanite}在1000K时可以达到0.9‰。此外，根据依配位数分组的1000lnβ-键长图中（图4 c），在特定的键长下，具有较高配位数的矿物反而具有较高的RPFR值（除硬石膏和萤石外，可能与矿物中Ca还与F配位有关）。

 

图4 键长和配位数对简约配分函数比的控制

当然，简约配分函数比与配位数的相关性较差，并不意味着配位数对同位素平衡分馏没有影响。这是因为矿物中钙周围的配位环境多样，受阴离子类型和压力等多种因素的影响。如果这些因素不相等，则很难单独通过配位数(一个离散值的整数，特别是大多数矿物的配位数集中在6和8上，图4b)来评估矿物中的配位环境。因此，如果不考虑键长、阴离子类型等因素，仅根据配位数，无法确定不同矿物间重同位素的富集顺序。

总之，相比于配位数，以键长作为判断平衡分馏方向的半经验参数更为可靠。如果两种矿物硅酸盐具有相似的Ca-O键长，那么具有较高配位数的矿物易富集更重的Ca同位素。然而，当Ca与不同的阴离子结合时，这一规则将不再适用（图4 c）。例如，萤石的配位数为8，橄榄石的配位数为6，萤石与Ca/(Mg+Ca)_mol=1/4的橄榄石具有相似的键长(2.342 vs. 2.333 Å)，但该橄榄石(RPFR=1.78 at 1000K)要比萤石(RPFR=1.54 at 1000K)更容易富集重的Ca同位素。这体现了阴离子对Ca同位素平衡分馏的控制。

阴离子的控制

与Ca结合的阴离子可以显著影响简约配分函数比。S与O属于同一主族，陨硫钙石(CaS)与石灰石(CaO)具有相似的晶体结构，但由于S2-的离子半径更大，因此Ca-S键长(2.826 Å) 要比Ca-O键长(2.416 Å)更长，因此，石灰石更易富集重的Ca同位素。

此外，尽管大多数Ca的有氧酸盐中Ca直接与O结合，但有氧酸根的中心原子也会影响Ca同位素平衡分馏。例如方解石 （CaCO₃）具有比硬石膏 （CaSO₄）更高的1000lnβ 值。这可能反映了阴离子结构和中心原子电负性的差异。一方面C的电负性比S略微偏小（2.55 vs 2.58），另一方面在CO₃^2-中，C原子与O原子成sp²杂化的σ键和离域大π键，而S原子与O原子成sp²杂化的σ键。因此，当与O原子成键时，C原子要比O原子丢失更多的电子云，使得CO₃^2-中的O原子含有更多的负电荷，使得碳酸盐中的Ca-O键比硫酸盐中的Ca-O键结合更紧密。

固溶体的组成

理论计算的对象多为成分简单的晶体，但自然界中的含钙矿物大多以复杂的固溶体形式存在。因此，理解固溶体成分对矿物简约配分函数比的影响对联系理论计算结果与实际自然样品的观测结果至关重要。前人已在多种矿物中发现了成分效应：

1. 在斜方辉石中^{[2, 3]}，随着Ca/(Mg+Ca)的比例减小，RPFR显著增加，但当Ca/(Mg+Ca)小于1/32或大于2/16时，成分对RPFR的影响消失。而Fe替换Ca时，含Fe 斜方辉石比相同成分的无铁斜方辉石略微偏轻 (0.07‰ at 1400K)；

2. 橄榄石中^[5-7]，Ca/(Mg+Ca)比例1/8与1/16之间均发现了显著的成分效应，但比例在1/16至1/64之间的成分效应仍存在争议。Huang^[5]和Antonelli^[6]的结果中未发现该区间中的成分效应，但Song^[7]的结果中发现RPFR在此成分区间内仍随Ca/(Mg+Ca) 比例变化。这可能是不同研究中构建的超胞结构的差异，Song 测试了所有可能的超胞结构，并选取了能量最低（最能代表自然矿物）的构型进行计算，而其他研究中选用的超胞可能具有更高的能量。因此Song 的结果应更为可靠;

3. 在石榴石中^[8]，1000lnβ 随着Ca/(Ca+Fe+Mg) 从1/24 增至 12/12 时降低了0.41;

4. 在长石中^[6]，拉长石 (An50) 的RPFR 比钙长石更低;

5. 在方解石或白云石中^[3,4]，当晶格中的Ca/(Mg+Ca)比例小于1/2时，RPFR随着Ca含量的增加而减小，而此规律在Ca/(Mg+Ca)大于1/2时反转。这可能是因为当Ca/(Mg+Ca)>1/2时晶格中主要由Ca控制，而Ca/(Mg+Ca)<1/2时转为Mg;

在此基础上，我们进一步测试了透辉石中Na-Al替换Ca-Mg对Ca同位素分馏的影响。我们构建了4种可能的替换结构：1) Na → Mg; Al → Ca; 2) Na, Al → Ca; 3) Na → Ca; Al → Mg; 4) Na, Al → Mg. 结果显示构型1和2（Al → Ca）具有比构型3和4 (Al → Mg)更高的RPFR。然而，Na⁺的离子半径与Ca²⁺相似，而Al³⁺的离子半径与Mg²⁺相似，因此构型3（RPFR=1.50）更能代表实际硬玉（含Na透辉石）的构型。而1000K下无Na透辉石的RPFR在1.45-1.56范围内，因此并未发现前人提出的辉石的“硬玉效应”。

除Na⁺外，我们还测试了K⁺替换Ca对碳酸盐和角闪石中Ca同位素分馏的影响。我们发现钾闪石（richterite）比透闪石（tremolite）的Ca-O键长更长（2.507 Å vs 2/474 Å），因而钾闪石具有比透闪石更易富集轻的Ca同位素。反映了当一个K⁺和一个Na⁺替换一个Ca²⁺时会增加未替换Ca的Ca-O键长。与之相反，当碳酸盐中K⁺替换Ca²⁺时（butschliite vs calcite），富K碳酸盐比方解石更易富集重的Ca同位素。这是由于方解石中的阳离子被CO₃^2-层分隔，当一层Ca²⁺被一层K⁺替换后，由于K⁺比Ca²⁺体积更大，这一替换将挤压Ca²⁺层，从而使得Ca-O键长降低。

压力的影响

 

压力对RPFR的影响是晶体晶格随压力增加而压缩的表现。我们的结果表明，无论钙长石或透辉石，其RPFR均与压力呈正相关关系（图5），这是因为随着压力的增加，Ca-O键长会相应变短。 

 

图5. 简约配分函数随压力的变化

值得注意的是，我们发现随着压力的增加，透辉石与钙长石间的Ca同位素分馏会发生倒转。尽管钙长石仅在低压下稳定，这种倒转难以发生在现代地壳中，但倒转压力随着温度升高而降低，因此有可能发生在具有较高热梯度的环境中。

平衡分馏系数受压力影响是因为两种矿物的1000lnβ随压力变化的斜率不同。在Mg同位素中也发现了类似的现象^[12]，其中辉石和橄榄石的Mg同位素1000lnβ的斜率比镁铝榴石和镁铁榴石更陡。这种现象可能反映了与配位数有关：钙在钙长石中是6配位和6.75配位(取决于统计方式)，而在透辉石中是8配位的；Mg在辉石和橄榄石中均为6配位，而在镁铝榴石和镁铁榴石中均为8配位。这可能是因为具有更低配位数的Ca或Mg与周围离子的堆积更加紧密，因此随着压力增加，键长变短的速度更快。相反，如果两矿物的配位数相同，那么1000lnβ随压力变化的斜率相似，计算1000lnα时可以抵消，因而压力对其矿物间同位素平衡分馏的影响可以忽略。例如石榴石和单斜辉石的Ca原子配位数都是8，因而在0GPa和10GPa两种压力下的Ca同位素平衡分馏系数几乎没有变化。

静态结果与非静态结果的对比

无论是静态RPFR还是非静态RPFR中都发现了Ca同位素平衡分馏系数在5GPa以上发生倒转的情况（图6），但在同样的温度和压力下，两种RPFR值具有显著差异。这一差异来源于所计算的结构是否考虑了热膨胀的影响。因为在高温下的热膨胀效应更为显著，所以结构差异更为明显。然而，随着温度升高，RPFR的绝对值变小，因此两种方法的差异随着温度升高也缓慢变小。而两种方法计算得到的0GPa下1000lnα_{diopside-anorthite}在1000K下比较相似（0.08 vs 0.15）。因为非静态方法较静态方法需要更大的计算量（5-10倍），所以大多数矿物仅计算静态结果足以。

 

图6. 0GPa和5GPa下diopside和anorthite的静态RPFR与非静态RPFR

 小结

我们计算了在0GPa下角闪石(钾闪石 (KNaCaMg₅Si₈O₂₂F₂)和透闪石 (Ca₂Mg₅[Si₄O₁₁]₂F₂)) 黄长石(镁黄长石 (Ca₂Mg[Si₂O₇])和铝黄长石 (Ca₂Al[AlSiO₇])), 含钠辉石 (i.e., 硬玉 (Na, Ca, Mg, Al)Si₂O₆), 富钾碳酸盐 (K₂Ca[CO₃]₂) 和其他含钙矿物 (磷灰石 (Ca₅[PO₄]₃F), 硬石膏 (CaSO₄), 钙钛矿 (CaTiO₃)和萤石 (CaF₂))和0-5GPa下透辉石和钙长石的Ca同位素简约配分函数比。通过汇总目前所有的简约配分函数比数据，我们发现RPFR直接受力常数控制。但RPFR与配位数关系较弱，因此在定性估计同位素分馏方向时，用键长估计要比用配位数更为可靠。为探究成分效应，我们发现Na掺杂对辉石中的Ca同位素分馏影响较弱，而富K角闪石（richterite）比无K角闪石（tremolite）更易富集轻的Ca同位素。anorthite 和 diopside 的 RPFR 与压力呈正相关，且理论上两者的同位素平衡分馏在可以在高压下反转。

文章信息

本工作以“The factors controlling equilibrium inter-mineral Ca isotope fractionation: Insights from first-principles calculations”为题，发表在于地球科学领域国际著名学术期刊Geochimica et Cosmochimica Acta上（DOI: 10.1016/j.gca.2022.07.021）, 通讯作者为中国科学技术大学地球和空间科学学院特任副研究员康晋霆博士，第一作者为博士研究生肖子聪。共同作者包括周陈硕士，吴忠庆教授，黄方教授。所有计算在中国科学技术大学超级计算中心完成。本工作受到科技部重点研发计划(2018YFA0702600)和国家自然科学基金(Grant No. 41903002)资助。