2011, Vol. 54
2. 北京大学地球与空间科学学院地球物理学系,北京 100871
2. School of Earth and Space Sciences, Peking University, Beijing 100871, China
作为上地幔的主要矿物，橄榄石及其高压相物质的性质在很大程度上决定了上地幔的物质结构、物质性质及动力学行为.对橄榄石相变动力学的深入研究，是研究俯冲带与周围地幔的相互作用[1~4]以及深源地震形成机制[5]等重大科学问题的关键.橄榄石及其相似矿物的实验研究表明，地幔中具有橄榄石结构的矿物到其高压相相变的主要机制是以不连贯的晶面成核和晶界长大为主的成核长大型相变[6, 7].
成核长大型相变的相变动力学理论计算包括三个方面:成核理论、长大理论和相变体积分数计算[8].成核率是单位时间、单位体积内新相晶核形成的数目；长大率是新相晶核尺度的增加速度；相变体积分数计算是考虑成核率、长大率随时间变化的情况下对新相物质所占总体积的比率随时间变化的计算.
Cahn[9]严格地推导了成核长大型相变的相变体积分数计算公式.应用Cahn[9]的理论，Rubie和 Ross[10]基于已有橄榄石相似矿物的相变动力学实验数据，计算了成核率和长大率对相变动力学计算结果的影响.结果表明在冷的俯冲环境中，成核过程对橄榄石的相变动力学过程影响很小.此后，长大动力学成了橄榄石相变动力学的研究重心[11~14].王曙光等[15]考虑到晶核形成时表面能对其形状的影响，对Cahn[9]的相变体积分数计算公式进行了修正.结果表明，在应用修正的相变体积分数计算公式对实验数据进行拟合并得到成核率、长大率数据时，所得的长大率相同，但成核率存在很大差异.具体影响取决于生成相颗粒与反应相颗粒之间以及生成相颗粒之间的表面能.结果同时表明，晶核的形状和成核率虽然基本上不影响对橄榄石的相变动力学过程的估计，但是会影响对相变后矿物颗粒粒度大小的估计，进而影响对俯冲带流变结构的估计.因此，对于俯冲带应力场及流变结构的研究需要首先基于实验数据对橄榄石的晶核形状及成核率进行准确估计.
目前橄榄石的成核率数据总体来说是缺乏的，仅有的成核率数据质量也需要进一步的检验.Rubie等[16]的实验给出了地幔橄榄石相似矿物Ni2SiO4的成核率及成核率参数.最近，Kubo等[13]的实验给出了镁橄榄石Mg2SiO4 的成核率数据.王曙光等[17]基于目前有限的地幔橄榄石(Mg, Fe)2SiO4 的相变动力学数据，对其成核率范围进行了约束.
本文从成核理论出发，考虑新相晶核表面能对其形状的影响，修正了晶面成核情形的成核动力学方程.基于新的成核率方程和王曙光等[15]修正的相变体积分数计算公式提出了利用相变动力学实验数据对成核率参数进行初步约束及判断成核率实验数据质量的方法.利用目前已有的成核率资料对成核率参数进行了初步确定，并对这些成核率资料的质量进行了判定.
2 非均匀成核的成核动力学方程 2.1 均匀成核情形下的Gibbs自由能变化橄榄石相变的非均匀成核理论是以经典的均匀成核理论为基础的.均匀成核理论是假定在反应相矿物颗粒内的任何位置都能成核(没有优势成核位置)，当新相晶粒的Gibbs自由能超过成核势垒时，新相的晶粒成为可以进一步长大成新相物质的晶核.下面我们以Christian[8]的工作为例，简要介绍传统的成核理论.
假设生成相物质以球的形状集合在一起的时候，会产生如下的Gibbs自由能变化:
|
(1) |
其中，r是新相集合体的半径，vβ 是新相的分子体积，gβ，gα 分别是生成相和反应相分子的Gibbs自由能，σαβ 是新相集合体和反应相物质之间单位面积的表面能.
方程(1)等号右边第一项代表了生成相集合体与体积相关的Gibbs自由能变化，第二项代表了与表面相关的Gibbs自由能变化.当Gibbs自由能的变化为极大值时，与Gibbs自由能变化的极大值相对应的临界半径为
|
(2) |
Gibbs自由能变化的极大值(临界Gibbs自由能)为
|
(4) |
临界晶核表面所包含的分子数为
|
(5) |
其中，sβ 是生成相分子的表面积.
2.2 传统的成核方程均匀成核的成核动力学方程为
|
(6) |
其中，T是温度，ΔGa 是摩尔成核活化能，R是气体常数，k是玻耳兹曼(Boltzmann)常数，Δgc 是临界成核功(临界Gibbs自由能).I0 是指前系数，
|
(7) |
其中，Nv 是单位体积内反应相物质的分子数目，h是亥姆霍兹(Helmholtz)常数.在成核动力学计算中一般把I0 取为待定常数，其数值由实验数据来约束.
新相晶核形状为透镜体.晶界成核的临界 Gibbs自由能Δg cB =ΦΔgc, 其中，Φ 是形状因子，表征了透镜体与同半径的球体的体积比，其表达为
|
(8) |
其中，θ 是反应相晶粒与生成相晶粒之间的接触角度.形状因子的取值范围是(0，1).
在均匀成核方程(6)基础上，传统的非均匀成核方程表示为
|
(9) |
I0B为待定常数.方程等号右面第二个指数项可以表示为
这种直接对均匀成核率方程(6)中临界成核功项做修正而得到非均匀成核率方程(9)，加上将系数I0 和I0B 取为待定常数的传统成核理论[8]方法虽然简化了方程，但忽略了均匀成核和非均匀成核情况下I0 和I0B 这两个系数的变化，无法根据实验数据同时拟合耦合在一起的形状因子和表面能:Φ(σαβ)3.下面，我们将在传统的成核理论[8]基础上，从新相晶核的形状出发，参考均匀成核率方程(6)的推导过程，推导出非均匀成核率方程，并进一步用以分离形状因子和表面能.
2.3 修正的非均匀成核的成核动力学方程颗粒边界成核的成核率方程可以表达为和均匀成核情形的成核率方程完全类似的形式:
|
(10) |
其中，
|
(11) |
NvB 是颗粒边界成核情形单位体积内反应相物质的分子数.对于橄榄石的晶体结构，由于其单位体积内的晶面的表面积为3.35/d，d为反应相物质的切角八面体晶粒中相互平行的两个正方形表面之间的距离，这时，
|
(12) |
其中，sα 为反应相物质的分子表面积.
下面，我们给出方程(10)中与新相物质晶核形状有关的参数:临界半径rcB、临界Gibbs 自由能(Gibbs自由能变化的极大值)Δg cB、晶核体内所包含的临界分子数ncB 以及表面所包含的临界分子数ocB.
对于颗粒边界成核情形，生成相晶粒的Gibbs自由能变化表示为
|
(13) |
其中，Φ 为形状因子，rB 为与球冠密合的球的半径，θ 为新相晶粒与反应相晶粒的接触角度，σαα 为反应相晶粒之间单位面积的表面能.由于反应相晶粒与生成相晶粒在接触角上的平衡，所以
|
(14) |
方程(13)可进一步化简为
|
(15) |
晶面成核情形下的临界半径:
|
(16) |
晶面成核情形下的临界Gibbs自由能变化
|
(17) |
晶面成核情形下的晶核体积内的所包含的临界分子数
|
(18) |
晶面成核情形下的晶核表面上的所包含的临界分子数:
|
(19) |
将以上参量的表达式代入方程(11)，得
|
(20) |
进一步地，假设r0 为生成相物质的分子半径，忽略反应相物质和生成相物质分子半径的微小差别时，我们可以近似地把方程(20)表示为
|
(21) |
王曙光等[15]考虑到新相晶核表面能对形状的影响，假设新相晶核形状为透镜体，对边界成核晶界长大的橄榄石相变体积分数计算公式进行了修正.应用相变体积分数计算公式，拟合高温高压条件下 X 光原位观测的相变动力学数据，可以得到对应的温度压力条件下的成核率和长大率.应用Cahn[9]的传统相变体积分数计算公式和新的计算公式，对相同的实验数据进行拟合时，得到的成核率之间存在如下的关系[15]:
|
(22) |
其中，IcB 和ImB 分别为应用传统相变体积分数计算公式和修正的相变体积分数计算公式拟合得到的成核率.
将方程(22)，(3)和(17)代入到方程(10)，我们可以得到用Cahn[9]的相变体积分数计算公式时的成核率与物质的物理参数间的关系:
|
(23) |
利用方程(21)，方程(23)可整理为如下形式:
|
(24) |
应用方程(24)拟合成核率数据，可以确定新相晶粒与反应相晶粒的接触角度θ，以及新相物质与反应相物质之间的表面能σαβ.
定义
|
(25) |
|
(26) |
并定义如下变量:
|
(27) |
|
(28) |
可以发现方程(24)表达的是lnIcB′ 和Δg cB′ 之间的线性关系:
|
(29) |
通过拟合方程(29)得到系数x1，x2 之后，可以进一步通过方程(25)和(26)计算接触角度θ 以及反应相物质与生成相物质之间的表面能σαβ.
由方程(25)，(26)，可以得到接触角度θ 与拟合参数x1，x2 的关系:
|
(30) |
其中，形状因子Φ 也是接触角度θ 的函数(见方程(8)).
θ 的取值范围为(0，π/2)，方程(30)等式左侧为θ 的单调递增函数，取值范围是(0，1)，拟合系数x1，x2 满足关系:
|
(31) |
下面我们将尝试利用现有的成核率数据确定新相晶粒与反应相晶粒的接触角度θ 以及新相物质与反应相物质之间的表面能σαβ.
3.2 对橄榄石及其相似矿物成核实验数据的拟合及成核率参数的确定目前对橄榄石及其相似矿物相变过程的实验研究主要分为两类，一类是利用X 光谱直接对相变体积分数随时间的变化进行原位测量，并通过对其拟合同时确定成核率和长大率；另一类是对经历不同反应时间的样品淬火卸载后进行切片观测，得到不同反应时间后新相的厚度，从而确定其长大率.目前仅有的利用橄榄石相似矿物的X 光原位观测相变动力学数据得到的成核率计算结果来自于Rubie等[16]和Kubo等[13].另外，王曙光等[17]通过对挤碰物理图像的分析，对边界成核情形的相变动力学理论中表征挤碰程度的无量纲参数进行了修正.在此基础上，王曙光等[17]利用第二类观测资料(经过淬火切片观测得到长大率)对橄榄石的成核率进行了估算(见表 1).
|
表 1 非均匀成核率方程中热参数拟合所使用的晶界成核率数据 Table 1 Data in parameter refinement of heterogeneous nucleation rate equation |
本文尝试性地根据表 1中的成核率数据对方程(29)中的成核率参数进行拟合.拟合过程中所使用的其他参数以及拟合结果也列于表 1中.同时，我们在表 1中给出了根据方程(30)得出的关于成核率参数的理论约束.可以看出，除了Rubie等[16]的成核率数据能够满足理论约束外，其他两组数据满足不了理论约束.
硅酸盐矿物的不连贯晶面的表面能在0.5~1J·m-2(Tomoaki Kubo, 个人交流).Cooper 和 Kohlstedt[18]在室温条件下观测的橄榄石颗粒之间的表面能为0.9J·m-2，目前对于橄榄石和尖晶石颗粒之间的表面能只有Rubie等[16]做了估计，他们应用传统成核率方程给出的σαβΦ1/3 为5×10-2J·m-2.当设表面能为0.6J· m-2 时，得到形状因子为6×10-4.
我们将根据实验结果得到的关于Ni2SiO4 和 Mg2SiO4 的点(lnIcB′ ，Δg cB′)以及按照约束条件计算的参考拟合线表示在图 1 中.只有斜率比参考线更小(绝对值更大)的直线才是符合条件的拟合线.Ni2SiO4 的拟合直线(图 1 中实线)符合这样的约束，但是Mg2SiO4 的拟合直线(图 1中的虚线)不符合，表现为温度高(对应横坐标值更大)的数据点，成核率偏大(对应纵坐标更大)，而且温度越高，成核率的偏离越大.
|
图 1
Ni2SiO4 和Mg2SiO4 成核率参数函数关系图
图中横坐标为 ，其中，ΔgcB是边界成核情形的临界Gibbs自由能，k是玻耳兹曼常数，T是温度，Φ 是形状因子，表征同样底面积的球冠和半球的体积比，σαβ 是反应相颗粒和生成相颗粒接触的表面能，Vβ 是生成相物质的摩尔体积，ΔGd 是摩尔驱动能.纵坐标是lnIcB′ =lnIcB - ，其中IcB 是按照传统的Cahn的相变体积分数对实验数据拟合得到的成核率[9], ΔGa 是摩尔活化能，d是颗粒大小，h是普朗克常数，r0 是橄榄石分子大小.图中方块和菱形分别是Ni2SiO4 和Mg2SiO4 的成核率所对应的数据点.实线和虚线分别是Ni2SiO4 和Mg2SiO4 的成核率数据点所对应的拟合直线.阴影线是符合接触角度在(0，π/2)之间的参考线，只有比相同截距的参考线斜率更小的拟合直线才是符合条件的.
Fig. 1
Nucleation parameter relation in Ni2SiO4 and Mg2SiO4
The abscissa is ，whereΔgcB is the critical Gibbs free energy of grain boundary nucleation.k is Boltzmann constant, Tis temperature, Φ is shape factor which represents the volume ratio of lens to the sphere with the same radius, a is interfacial energybetween the product nuclei and react nuclei, V is molar volume of the product material, AGd is driving energy.The ordinate is lnIcB′ =lnIcB -, where IcB is the nucleation rate from fitting the traditional Cahn's formula to the experimental data[9]，△Ga ismolar activation energy, d is grain size, h is Planck constant, r0 is molecule size of olivine.Square and diamond are the nucleation rate data points of Ni2SiO4 and Mg2SiO4 respectively.Solid line and dash line are the nucleation rate fitting lines of Ni2SiO4 and Mg2SiO4 respectively.Shadowed lines are the reference lines according to the constraint of contact angle within (0，π/2). 0nly refined lines with the same intercept but smaller slope satisfy the constraint of contact angle.
|
对于(Mg, Fe)2SiO4 的实验数据，由于Mosenfelder等[11]的长大率拟合直线与对应的观测数据在900℃时存在很大差别，而且目前无从判断在900℃下，应用拟合直线计算的长大率与观测的数据哪个更准确，我们在此分别讨论这两种情况.
图 2 为(Mg, Fe)2SiO4 成核率参数函数关系图.应用Mosenfelder等[11](见表 1)的数据(在图 2中表示为上三角、方块和圆圈)计算成核率.由这些实验点得到的拟合直线不能够拟合出合理的接触角度.这与图 1中Mg2SiO4 的拟合结果类似:拟合直线的斜率比参考线的大(绝对值小)，且高温条件下的成核率偏大.应用Mosenfelder等[11]的长大率拟合直线计算18GPa, 900 ℃条件下的长大率所对应的成核率(表 1中标下划线的数据)时，对应的拟合直线比参考线的斜率更小(绝对值大)，可以拟合出合理的接触角度.
|
图 2 (Mg, Fe)2SiO4 成核率参数函数关系图 图中三角、方块、圆圈分别对应于18GPa, 900 ℃、1000 ℃、1100 ℃下的成核率数据，空心和实心分别对应于相应温压条件下成核率数据的最小值和最大值.倒三角为应用Mosenfelder等[11]的拟合直线计算的长大率数值得到的成核率参数.正三角为应用 Mosenfelder等[11]的表 1中直接读取的长大率数据得到的成核率参数.横纵坐标以及阴影线的表达和图 1中相同. Fig. 2 Nucleation parameter relation in (Mg?Fe)2SO4 Triangle, square, circle represent nucleation rate under 18 GPa and 900 °C，1000 °C，1100 °C respectively.0pen and solid symbols respond to the mimmum and maximum value of the nucleation rate under the certain temperature and pressure condition.Reversed triangle scalculated by the fitting line of Mosenfelder et al.[11], triangles are directly read from Table 1 in Mosenfelder et al.[11] .The abscissa and ordinate and the shadow grey line are same as that in Fig.1. |
本文考虑到表面能对晶核形状的影响，修正了非均匀成核方程.新的成核率方程还原了传统成核率方程中简化为常数并由实验数据确定的指前系数.由此可以用来确定出成核率方程中新相晶粒与反应相晶粒的接触角度及新相物质与反应相物质之间的表面能.而传统的非均匀成核率方程只能确定指前系数及耦合在一起的接触角度和表面能.
结合本文修正的非均匀成核率方程，以及王曙光等[15]修正的相变体积分数计算公式，本文提出了一种新的方法，可以基于以往的实验数据确定成核率方程中的物理量:接触角度和表面能.由于接触角度的取值范围是(0，π/2)，这种方法可以用来确定实验数据的质量.
针对目前仅有的成核率实验数据进行的成核率方程中物理量的确定表明，Ni2SiO4 的实验数据可以给出接触角度为0.57°，表面能为32.33J·m-2.而Mg2SiO4 的实验数据不能拟合出范围在(0，π/2)的接触角度，可能的原因是高温条件下的成核率数据严重偏离真值[24].(Mg, Fe)2SiO4 的实验数据是否能拟合出合理的接触角度与900℃下的实验观测有直接关系.应用实验观测的900 ℃的长大率数据所计算的成核率，不能够拟合出合理的接触角度，当应用拟合直线计算出的长大率数据得到的成核率替换该观测值时，可以拟合出合理的接触角度.这可能和实验加热过程中相变的影响有关.
应用Ni2SiO4 的成核率数据，以及(Mg, Fe)2SiO4经过拟合曲线计算的成核率数据，得到的接触角度几乎相同，都是0.57°.这种一致性是偶然还是真实的物理情形，需要进一步研究.
由于成核率对于俯冲带的颗粒粒度结构以及流变结构的决定性影响，对于成核率的研究是俯冲带动力学研究的关键之一.进一步的研究需要补充更多的低温或者淬火后先加压再加温的动力学实验，以精确地确定橄榄石成核的动力学特征.
致谢本文的研究过程中和Kubo Tomoaki教授及景志成博士进行了有益的讨论，匿名审稿人提出的中肯的修改意见亦对本文有所帮助.在此一并表示感谢.
| [1] | Bina C R, Stein S, Marton F C, et al. Implications of slab mineralogy for subduction dynamics. Phys. Earth Planet. Inter. , 2001, 127(1-4): 51-66. DOI:10.1016/S0031-9201(01)00221-7 |
| [2] | Rubie D C. The olivine(spinel transformation and the rheology of subducting lithosphere. Nature , 1984, 308(5959): 505-508. DOI:10.1038/308505a0 |
| [3] | Karato S I, Riedel M R, Yuen D A. Rheological structure and deformation of subducted slabs in the mantle transition zone: implications for mantle circulation and deep earthquakes. Phys. Earth Planet. Inter. , 2001, 127(1-4): 83-108. DOI:10.1016/S0031-9201(01)00223-0 |
| [4] | Schelming H, Monz R, Rubie D C. The influence of olivine metastability on the dynamics of subduction. Earth Planet. Sci. Lett. , 1999, 165(1): 55-66. DOI:10.1016/S0012-821X(98)00249-0 |
| [5] | Sung C M, Burns R G. Kinetics of high-pressure phase transformations: implications to the evolution of the olivine(spinel transition in the downgoing lithosphere and its consequences on the dynamics of the mantle. Tectonophysics , 1976, 31(1-2): 1-32. DOI:10.1016/0040-1951(76)90165-7 |
| [6] | Burnley P C. The effect of nonhydrostatic stress on the olivine-spinel transformation in magnesium germinate. Davis: University of California, 1990 : 187 . |
| [7] | Brearley A J, Rubie D C, Ito E. Mechanisms of the transformations between the alpha, beta and gamma polymorphs of Mg2SiO4 at 15 GPa. Phys. Chem. Miner. , 1992, 8: 343-358. |
| [8] | Christian J W. The Theory of Transformations in Metals and Alloys Part I. Equilibrium and General Kinetic Theory. Oxford: Pergamon, 1975 : 449 . |
| [9] | Cahn J W. The kinetics of grain boundary nucleated reactions. Acta Metall. , 1956, 4(5): 449-459. DOI:10.1016/0001-6160(56)90041-4 |
| [10] | Rubie D C, Ross C R II. Kinetics of the olivine-spinel transformation in subducting lithosphere: experimental constraints and implications for deep slab processes. Phys. Earth Planet. Inter. , 1994, 86(1-3): 223-241. DOI:10.1016/0031-9201(94)05070-8 |
| [11] | Mosenfelder J L, Marton F C, Ross C R, et al. Experimental constraints on the depth of olivine metastability in subducting lithosphere. Phys. Earth Planet. Inter. , 2001, 127(1-4): 165-180. DOI:10.1016/S0031-9201(01)00226-6 |
| [12] | Kubo T, Ohtani E, Kato T, et al. Effects of water on the α-β transformation kinetics in San Carlos olivine. Science , 1998, 281(5373): 85-87. DOI:10.1126/science.281.5373.85 |
| [13] | Kubo T, Ohtani E, Funakoshi K. Nucleation and growth kinetics of the α-β transformation in Mg2SiO4 determined by in situ synchrotron powder X-ray diffraction. Am. Miner. , 2004, 89: 285-293. DOI:10.2138/am-2004-2-305 |
| [14] | Hosoya T, Kubo T, Ohtani E, et al. Water controls the fields of metastable olivine in cold subducting slabs. Geophys. Res. Lett. , 2005, 32: L17305. DOI:10.1029/2005GL023398 |
| [15] | 王曙光, 叶国扬, 宁杰远, 等. 对一类相变体积分数计算公式的修正. 地球物理学报 , 2007, 50(6): 1794–1801. Wang S G, Ye G Y, Ning J Y, et al. Modification on the volume fraction formulae of a sort of phase transformations. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1794-1801. |
| [16] | Rubie D C, Tsuchida Y, Yagi T, et al. An in situ X ray diffraction study of the kinetics of the Ni2SiO4 olivine-spinel transformation. J. Geophys. Res. , 1990, 95(B10): 15829-15844. DOI:10.1029/JB095iB10p15829 |
| [17] | 王曙光, 宁杰远, 叶国扬, 等. 对橄榄石相变形核率的估算. 地球物理学报 , 2007, 50(1): 122–125. Wang S G, Ning J Y, Ye G Y, et al. Evaluation of the nucleation rate of olivine. Chinese J. Geophys. (in Chinese) , 2007, 50(1): 122-125. |
| [18] | Cooper R F, Kohlstedt D L. Interfacial energies in the olivine-basalt system. In Akimoto S, Manghnami M H, eds. High Pressure Research in Geophysics. Adv. Earth Planet. Sci., 1982, 12: 217~228 |
| [19] | Akimoto S I, Fujisawa H, Kalsura T. The olivine-spinel transition in Fe2SiO4 and Ni2SiO4. J. Geophys. Res. , 1965, 70(8): 1969-1977. DOI:10.1029/JZ070i008p01969 |
| [20] | Morishima H, Kato T, Suto M, et al. The phase boundary between α- and β- Mg2SiO4 determined by in situ X-ray observation. Science , 1994, 265(5176): 1202-1203. DOI:10.1126/science.265.5176.1202 |
| [21] | Akaogi M, Ito E, Navrotsky A. Olivine-modified spinel-spinel transitions in the system Mg2SiO4-Fe2SiO4: calorimetric measurements, thermochemical calculation, and geophysical application. J. Geophys. Res. , 1989, 94(B11): 15671-15685. DOI:10.1029/JB094iB11p15671 |
| [22] | Ringwoodite A E. Composition and Petrology of the Earth's Mantle. New York: McGraw-Hill, 1975 . |
| [23] | 叶国扬. 关于俯冲带内部410-km间断面起伏的地震学探测. 北京: 北京大学地球与空间科学学院地球物理系, 2007 . Ye G Y. Seismological study on the undulation of the 410-km discontinuity in subduction zones (in Chinese). Beijing: Department of Geophysics, School of Earth and Space Sciences, Peking University, 2007 . |
| [24] | 王曙光, 宁杰远. 加载、加热路径上橄榄石的相变及其对相变动力学参数计算的影响. 地球物理学报 , 2010, 53(1): 111–118. Wang S G, Ning J Y. Nucleation on heating/loading path and its effects on the evaluation of nucleation rate. Chinese J. Geophys. (in Chinese) , 2010, 53(1): 111-118. DOI:10.3969/j.issn.0001-5733.2010.01.012 |