2018, Vol. 61
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
Wilson(1963)根据夏威夷岛链的年龄递增现象，提出“热点构造”的概念后，Morgan(1971)在此基础上首次系统性的提出了“地幔柱”的概念，他认为地球内部存在起源于地幔的缓慢上升的细长柱状热物质流，这是地幔柱形成的初始过程中所需要的地幔柱尾部通道的雏形.经过几十年的研究，国内外学者对地幔柱的认识逐步深入，“地幔柱”理论也逐渐成为一种新的大地构造学说.地幔柱是地球内部热对流的一种方式，起源于670 km的热/化学边界层(Christensen and Harder, 1991; Burov and Poliakov, 2001; d′Acremont et al., 2003)或者2900 km的核幔边界(Farnetani and Richards, 1994; Lavelle, 1997; Kellogg and King, 1997).温度较高的地幔热物质由软流圈或下地幔一路上升到地幔的顶部，并与岩石圈相互作用.在此过程中，地幔柱的头部和岩石圈底部会发生部分熔融，形成玄武岩岩浆，同时可能于短时间内(少于一百万年)沿地壳薄弱带或断裂带大量喷出地表.喷发地点若位于洋壳处，则形成海底高原，产生热点；如果岩浆于大陆地壳内部喷发，大量的溢流玄武岩覆盖于地表之上，冷却后形成大火成岩省(LIPs).目前一般认为，地幔柱是热点轨迹或者由大量溢流玄武岩组成的大火成岩省最有可能的产生机制(Morgan, 1981；Richards et al., 1989；Campbell and Griffiths, 1990).
目前针对地幔柱的研究范围涉及地幔柱的形成条件、形成机制以及退化消亡等方面，主要集中于通过地幔柱与岩石圈的相互作用分析证明地幔柱的存在(李建康和王登红，2005).地质学、地球物理和地球化学(卢记仁，1996；张本仁，2001；徐义刚，2002)等常规方法着眼于地幔柱与岩石圈相互作用后在地表产生的静态信息，反演推测出其在相互作用过程中的动态信息.而实验方法(Campbell and Griffiths, 1990)和数值模拟方法(d′Acremont et al., 2003; Burov and Guillou-Frottier, 2005)则是通过模型设置，正演直接得到地幔柱与岩石圈相互作用过程中的动力学信息，同时对反演所得模型进行检验.数值模拟方法相比较实验方法而言，在调节各种物理化学参数方面更为方便，极大地缩短了模拟周期，同时可直接得到温度场、压力场、速度矢量场等更多信息.因此，数值模拟方法对地幔柱与岩石圈相互作用的研究起着非常重要的作用.
Campbell和Griffiths(1990)通过实验室模拟的方法，针对不同密度、不同黏度的物质的对比，模拟了地幔柱的形态、生成和演化过程，给出了与热点轨迹和大火成岩省的形成相关的第一个地幔柱物理模型，形象地模拟出地幔柱的结构：巨大的球状顶冠和狭窄的尾部结构.计算机技术的发展使数值模拟逐渐成为研究地幔柱的重要手段之一，国内外学者利用各种地质模型和数值方法对地幔柱动力学进行了大量的研究.Olson(1990)通过简单的数值模拟指出地幔柱在上升过程中，其头部会热熔蚀岩石圈的底部；Christensen和Harder(1991)在不区分不同圈层物质组分的前提下，提出了三维的地幔对流模型；Farnetani和Richards(1994)将地幔柱设置为与温度弱相关的均一黏度流体，利用数值模拟得出海底高原和大陆溢流玄武岩与地幔柱活动可能存在的关系，并建立了地球圈层中2600 km以上的物质黏度随深度的变化关系；Manglik和Christensen(1997)将地幔熔融后的残余物质所产生的浮力作为熔融量的影响因素之一，解释了在地幔柱的上升过程中岩石圈中岩石熔融速率的变化规律.d′Acremont等(2003)和Burov等(2007)讨论了源于上地幔的地幔柱与岩石圈相互作用的过程.前者以加勒比大火成岩省为例，模拟得出源于400 km的地幔柱与大洋岩石圈相互作用进而形成海洋大火成岩省的动力学过程；后者将初始地幔柱深度设置为650 km，利用地幔柱与岩石圈相互作用后模拟得到的温度场、黏度场和不同地表地形特征来讨论地幔柱的存在性.国内学者(李建康和王登红，2004；Leng和Gurnis，2012)主要模拟源于核幔边界处的地幔柱演化过程.前者模拟得到在地幔柱动力演化过程中的温度场时空变化曲线；后者研究了2870 km厚度的地幔中，地幔柱的形状随黏度变化的关系.蒙伟娟等(2005)以峨眉山大火成岩省为例，模拟了源于670 km的地幔柱与岩石圈相互作用的过程中，模型区域内的温度、黏度和地表地形的变化过程；并通过比较不同黏度模型的模拟结果，得出选择非牛顿流体的必要性.
由上述可以发现，对地幔柱的模拟方法的重心逐渐从实验室模拟向计算机数值模拟转移，模型的设置由均一介质变为各圈层介质物性相异，针对地幔柱的流变性设置也由牛顿流体转变为非牛顿流体.目前针对地幔柱与岩石圈相互作用中的模拟，主要研究的是温度、黏度和地表地形的变化.但对于大火成岩省的形成机制而言，地幔柱与岩石圈相互作用过程中熔融程度的变化是研究大火成岩省溢流玄武岩的起源及演化过程的关键.
本文利用开源有限元程序Ellipsis3D，集中模拟地幔柱与岩石圈的相互作用过程中，不同圈层物质在不同的作用阶段产生的熔融程度的变化.在非牛顿流体近似的基础上，采用与温度、压力、应变率强相关的黏度模型，对地幔柱在上地幔中上升以及与岩石圈接触后发生相互作用的过程中，不同圈层物质的熔融程度变化进行二维的数值模拟.通过数值模拟计算，得出不同初始温度条件和不同的固液相线的设置对物质初始熔融深度及熔融程度的影响.比较对于不同初始温度的地幔柱在相同相线设置下的熔融程度的差异，得出对于不同初始温度的地幔柱，采用不同固液相线模型进行计算的合理性.分析不同圈层物质在动态演化过程中的相图、熔融程度图，讨论不同圈层物质的熔融机制.最后，展示了熔融过程的数组模拟在探讨峨眉山大火成岩省的形成机制上的应用.
1 物理模型本文涉及的数值模拟是利用Ellipsis3D程序(O′Neill et al., 2006)对计算所涉及的，由质量守恒方程、动量守恒方程和能量守恒方程组成的偏微分方程组进行求解，其中流变关系采用非牛顿流体近似.Ellipsis3D是在Ellipsis(Moresi et al., 2003)的基础上进行修正，并实现了三维直角坐标的计算后，得到的适合模拟地球介质大应变的有限元程序.与普通的有限元程序不同，Ellipsis3D采用Lagrange移动积分点和欧拉网格进行混合编码，其主要特点是采用欧拉网格，利用粒子追踪岩石的流变历史，并将粒子作为有限单元的积分点，处理伴有应变弱化或损耗以及复杂的发展性边界条件的黏弹性-塑性流变问题.与处理地幔对流问题的有限元程序相比，Ellipsis3D更适合模拟岩石圈的动力学行为.
Ellipsis3D中的关键算法主要有Uzawa算法、完全多重网格方法(Full Multigrid Method)及Lagrange积分点有限元法.联立质量守恒方程和Navier-Stokes矩阵方程，利用Uzawa算法求解瞬时速度场，避免了由于计算区域节点较多导致的无法直接对矩阵求逆的问题；完全多重网格法弥补了利用共轭梯度算法求解高自由度的线性方程组时，迭代收敛速度较低的问题，加快了迭代收敛速度；Lagrange积分点有限元方法的中心思想是在保持欧拉网格通用性和鲁棒性的同时，通过粒子追踪物质的变形过程和流变历史，并将相互独立的网格和粒子动态耦合.该方法虽然计算效率比高斯积分点的有限元方法低，但是在处理粒子与网络耦合问题上占有绝对优势.
1.1 控制方程
假设流体不可压缩，根据Boussinesq近似，描述质量守恒的连续性方程
|
(1) |
在Ellipsis3D中，应力张量σ分成两部分，压力P和偏应力τ：
|
(2) |
其中δij为克罗内克函数；从而得到忽略惯性项后的描述动量守恒的运动方程：
|
(3) |
其中f是体力；描述能量守恒的能量方程为
|
(4) |
其中，CP是等压热容，
|
(5) |
根据Ranalli(1995)高温高压岩石物性实验，地壳和地幔岩石有如下的流变定律关系式：
|
(6) |
其中
|
(7) |
从上式可以看出，岩石的流变性质与应变率，活化能，温度强相关.为了便于处理，数值模拟中一般用有效黏滞系数ηeff来表征岩石的流变性质.根据有效黏滞性系数的定义：
|
(8) |
将(7)式代入(8)式中可得岩石的有效黏滞性系数可表示为
|
(9) |
除了上述流变定律之外，岩石在受力条件下的屈服机制及强度弱化机制也是岩石的流变性质非常重要的一个方面.数值模拟中一般采用如下屈服机制(O′Neill et al., 2006)：
|
(10) |
其中，屈服应力遵循Mohr-Coulomb准则：
|
(11) |
其中，C0为内聚力，μ为摩擦系数.
1.3 熔融损耗关系地球深部高温的地幔岩石在地幔中快速上升，到达某一深度时，上升的地幔岩石的温度与在该深度对应的压力下，物质的固相线温度相等，地幔热物质开始熔融.如图 1所示.
|
图 1 熔融原理图 Fig. 1 Melting principle |
假定地幔残余物分批次进行熔融析出和熔融损耗，在Ellipsis3D中对熔融程度进行了参数化.首先定义一个超固相线温度Tss：
|
(12) |
其中，T是追踪粒子的温度，T(P)s和T(P)l是在压力P下橄榄岩的固相线温度(Takahashi and Kushiro, 1983)和液相线温度(McKenzie and Bickle, 1988)，如图 2所示.McKenzie和Bickle(1988)证明了熔融比例F和超固相线温度Tss间的关系：
|
图 2 橄榄岩的相线图 虚线表示线性表达(Takahashi and Kushiro, 1983), 实线表示多项式表达(Gasparik, 1990; Ohtani, 1986)；红色表示液相线，黑色表示固相线；蓝色和品红色分别代表熔融程度为30%和60%. Fig. 2 Phase diagram for peridotites The thick lines are the solidus (black) and liquidus (red) for the linear (dashed) and polynomial (solid) parameterization. The thinner lines correspond with phase equilibrium lines for which a degree of depletion is reached of 30% (blue) and 60% (magenta). |
|
(13) |
熔融产物导致了在能量方程中需要增加一个潜能项，能量方程变为
|
(14) |
其中，ΔS是熔融过程中的熵变，F是熔融比例.
1.4 模型地幔柱和大陆岩石圈的相互作用发生在上地幔，所以地幔柱的起源不需要太深.因此，本文参考Burov和Guillou-Frottier(2005)，模型的垂向尺度设置限制在670 km；本文中模型采用笛卡儿坐标系，计算区域为1340 km×670 km，初始圆形地幔柱设置在区域底部，低黏度、高温度(+300K)、半径为70 km；研究区域网格划分为160×80.初步设置莫霍面深度为45 km(黄汲清等，1977)，其中上地壳深度20 km, 代表性岩石为花岗岩；下地壳厚度25 km, 代表性岩石为镁铁质麻粒岩；上地幔代表性岩石为橄榄岩.详细参数见表 1和表 2.
|
表 1 数值模拟过程中共同的符号标记和物性参数(Turcotte和Schubert，2002) Table 1 The common notations and parameters in the numerical modeling |
|
|
表 2 数值模拟过程中不同圈层的物性参数 Table 2 The physical parameters of different layers in the numerical modeling |
本文采用与d′Acremont等(2003)相似的温度条件.假定在岩石圈内部，热传导为热量传递的主要方式，在上地幔则近似绝热温度梯度来确定初始温度场.岩石圈底部温度为1603 K，区域底部670 km处温度为1873 K.另外，边界条件设置如下：左右两侧为固定边界，水平速度vxx=0, 垂直方向为滑动边界；区域底部为岩石静压力，上表面为自由表面.初始模型图及边界条件剖面如图 3所示.
|
图 3 实验设置和边界条件 Fig. 3 Setup and boundary conditions |
图 4描述的是在初始温度为2173 K，线性固、液相线的条件下，模型区域物质熔融程度变化的模拟结果；由下至上(a—d)描述的是地幔柱在上地幔中上升及其与岩石圈相互作用时的熔融程度的变化过程.该过程大致可分为三个阶段，分别是：地幔柱的初融阶段、地幔柱与岩石圈的纵向作用阶段、地幔柱的横向展平阶段.
|
图 4 不同阶段下地幔柱与岩石圈相互作用图 第一列：熔融比例图；第二列：相图-青色代表岩石圈, 灰绿色代表上地幔, 红色代表地幔柱. Fig. 4 The interaction of plume and lithosphere in different stages The first column: depletion fields; the second column: phase fields-cyan represents lithosphere, grayish-green represents upper mantle, red represents plume. |
第一阶段：如图 4a-4b，初始位置位于转换带的高温低黏的地幔柱异常体，由于与周围物质存在着温度和密度等方面的差异而产生扰动，由于重力不稳定产生向上的速度，进而在地幔中开始快速上升.由于地幔柱在上地幔中上升速度较快，其与周围的地幔岩石热交换量较小，因此它的温度基本保持恒定.但在此过程中，地幔柱所受压力迅速减小，变化较大，其对应的固相线温度随之快速减小.到达某一深度，在该深度的压力下的固相线温度恰好与地幔柱固有温度相等，即地幔柱的温度线与地幔柱内物质的固相线相交，地幔柱开始熔融.
第二阶段：如图 4b-4c，地幔柱在上地幔中继续上升，到达岩石圈底部，并开始不断纵向侵蚀岩石圈.当地幔柱到达岩石圈底部时，地幔柱的柱头部分已经完成了一定程度的熔融，而柱尾则大部分在上升过程中，由于与周围相对低温的物质进行热交换而未能达到其固相线温度，仅靠近柱头中心的极少部分柱尾发生少量熔融.地幔柱与岩石圈纵向作用的过程中，压力改变较小，温度梯度较大，地幔柱和与其接触的岩石圈地幔发生热交换，地幔柱周围的岩石圈地幔温度升高，超过其在该深度下的固相线温度，开始熔融.此时地幔柱本身仍以减压熔融为主，但对于岩石圈地幔，升温熔融占其熔融的主导地位.
第三阶段：如图 4c-4d，地幔柱以横向展平为主，不断水平向侵蚀岩石圈底部.在地幔柱的展平过程中，其中心温度变化较小.随着展平范围扩大，与高温地幔柱接触的相对低温的岩石圈地幔范围增大，地幔柱边缘物质与周围冷物质发生热交换，大量的岩石圈地幔温度达到其固相线温度，开始大规模熔融.熔融地幔物质随熔融地幔柱物质继续向两侧水平运动，直至熔融物质在岩石圈范围内冷却.
2.2 熔融程度的影响因素 2.2.1 不同固、液相线的选取对熔融模拟结果的影响为了更加清楚的对比不同固、液相线的选取对熔融模拟结果的影响，地幔柱的初始温度设置为高于模型区域底部边界300 K，即2173 K.通过数值模拟计算，我们得到了在相同的初始温度条件，不同的固、液相线模型下地幔柱与岩石圈相互作用过程中的熔融程度的变化，如图 5所示.图 5中，第1、2列表示线性固、液相线的条件下，地幔柱与岩石圈相互作用过程中的熔融程度和相图的模拟结果；第3、4列表示多项式固、液相线的条件下，地幔柱与岩石圈相互作用过程中的熔融程度和相图的模拟结果.从下至上分别表示(a)地幔柱的初融位置，(b)地幔柱到达岩石圈底部的熔融程度，(c)地幔柱与岩石圈纵向作用过程中的熔融程度和(d)地幔柱在岩石圈地幔中横向展平过程中的熔融程度.
|
图 5 初始温度为2173K时，不同相线模型下地幔柱与岩石圈相互作用图 第1、2列：利用线性固相线模拟; 第3、4列：利用多项式固相线模拟; 第1、3列：不同相线模型下的熔融比例图；第2、4列：不同相线模型下的相图-青色代表岩石圈, 灰绿色代表上地幔, 红色代表地幔柱.由下至上(a—d)：地幔柱与岩石圈相互作用过程随时间变化图. Fig. 5 The depletion fields and phase fields for both linear and polynomial parameters of the phase-diagram when the initial temperature of plume is 2173K The first and second column: linear solidus and liquidus; the third and fourth column: polynomial solidus and liquidus. The first and third column: depletion fields; the second and four column: phase fields-cyan represents lithosphere, grayish-green represents upper mantle, red represents plume.Bottom to up (a—d): the interaction variation of plume and lithosphere over time. |
首先是地幔柱在上地幔中上升过程中的初始熔融深度(图 5a).我们可以看出，在地幔柱上升的0.2 Ma内，线性模型模拟的地幔柱在深度约为209 km处首次出现了小程度的熔融；多项式模型模拟的地幔柱则在上升约0.14 Ma，位于302 km左右开始熔融，与线性模型的模拟结果相差近100 km.造成这样的差距的原因是地幔柱在上地幔中快速上升时，温度变化小，但其所受压力迅速减小，变化较大，因此地幔柱在上地幔上升过程中发生的熔融应以减压熔融为主，熔融最早出现在地幔柱柱头中心.多项式固、液相线的模型可以很好的刻画在地幔深部，由于压力的变化造成固相线温度的变化；而线性模型则适用于在浅部温度变化较大的物质的熔融程度的模拟，但物质在深部由于快速减压而产生的熔融无法得到很好的反映.地幔柱继续上升至岩石圈底部(图 5a-5b)，此时地幔柱内部最大熔融程度发生在深度为190 km处的地幔柱中心，利用线性模型所模拟的地幔柱最大熔融程度约为11%；而利用多项式模型(模型B)所模拟的地幔柱其自身大部分已经具有相当程度的熔融，熔融程度由与上地幔物质接触的边缘向地幔柱中心逐渐增大，最大熔融程度达到68%.地幔柱纵向侵蚀岩石圈底部发生在0.2~0.26 Ma(线性模型)以及0.19~0.28 Ma(多项式模型)之间.地幔柱和与其接触的岩石圈地幔进行热交换，部分岩石圈地幔温度达到其固相线温度，发生熔融.线性模型模拟的地幔柱与岩石圈由于发生纵向作用使岩石圈地幔产生熔融最大深度为58 km，对于多项式模型模拟结果而言，最大深度约为50 km；线性模型模拟结果显示地幔柱最大熔融程度为44%，而多项式模型的模拟结果则高达89%，如图 5c所示.对于横向展平熔融结果而言，如图 5d，线性模型的模拟结果熔融深度范围为59~126 km，多项式模型模拟的深度范围为50~160 km.同时，在地幔柱的横向运动过程中，地幔柱和岩石圈地幔整体的平均熔融程度减小，即地幔柱在横向展平作用过程中，与周围相对冷的岩石圈地幔接触面积逐渐增大，导致其自身热量耗损较快.因此，地幔柱的横向运动过程是以地幔柱自身物质熔融程度减小，其周围岩石圈地幔物质熔融程度增加为结果，在扩张过程中，地幔柱内部逐渐开始冷却.
2.2.2 不同初始温度的地幔柱的选取对熔融模拟结果的影响将地幔柱初始温度设置调整为1903 K，利用线性模型和多项式模型模拟的熔融结果如图 6所示.其中第1、2列表示线性模型的模拟结果对比，第3、4列表示多项式模型模拟结果对比.第1、3列的模拟结果中，地幔柱初始温度为2173 K，第2、4列地幔柱初始温度为1903 K.
|
图 6 不同初始温度的地幔柱与岩石圈相互作用的熔融程度 第1列：地幔柱初始温度2173K，固相线为线性表达；第2列：地幔柱初始温度1903 K，固相线为线性表达；第3列：地幔柱初始温度2173K，固相线为多项式表达；第4列：地幔柱初始温度1903 K，固相线为多项式表达. Fig. 6 The depletion fields of variable initial temperatures for both linear and polynomial parameters of the phase-diagram The first and the third column: the initial temperature of plume is 2173 K; The second and the fourth column: the initial temperature of plume is 1903 K; The first and second column: linear solidus and liquidus; The third and fourth column: polynomial solidus and liquidus. |
首先是线性模型下不同地幔柱初始温度的熔融程度对比，如图 6中第1列和第2列所示.在第2列，即地幔柱初始温度为1903 K，线性固、液相线的模拟结果中，地幔柱的初始熔融深度约159 km，纵向作用深度为67 km，最大熔融程度变为32%.相比较第1列，即初始温度为2173 K中所得结果——初融深度209 km，纵向作用深度58 km-地幔柱的初始熔融深度变浅，地幔柱与岩石圈相互作用产生熔融的纵向作用深度变深，也就是说地幔柱与岩石圈纵向作用的范围减小.地幔柱内部最大熔融比例从44%变为32%.物质熔融程度减小，则地幔柱与岩石圈相互作用的熔融混合物熔融程度，即平均熔融比例减小.
其次是多项式模型下不同地幔柱初始温度的熔融程度对比，如图 6中第3列和第4列所示.在第四列，即地幔柱初始温度为1903 K的多项式固、液相线的模拟结果中，地幔柱的初始熔融深度为176 km，纵向作用深度为75 km，最大熔融程度约为39%.相比较初始温度为2173 K的模拟结果——初融深度302 km，纵向作用深度50 km，最大熔融程度89%-存在大幅度的变化.由此可以看出，对于地幔柱小范围的温度变化，多项式模型模拟的熔融结果响应更加明显.而当地幔柱到达岩石圈底部，与岩石圈地幔发生相互作用时，由于地幔柱自身温度较低，其与周围的岩石圈地幔的温度差减小，通过热传递达到热平衡时的温度相对较低，岩石圈地幔的熔融程度较低，岩石圈熔融物质与地幔柱熔融物质混合后的平均熔融量因此减小.同时，从地幔柱发生初始熔融时间及地幔柱到达岩石圈并与其发生相互作用的时间可以得出，相对高温的地幔柱在上地幔中的上升速度比低温的地幔柱上升速度快，地幔柱与岩石圈相互作用的整个过程的时间较短.
3 对峨眉山大火成岩省的应用地幔柱研究是判断大火成岩省的重要手段之一.大火成岩省(Large igneous province, LIPs)的主要岩体一般为镁铁质火山岩及伴生侵入岩，其在空间分布特征上呈现大尺度连续(Coffin and Eldholm, 1994).根据其组成和结构的差异，一般包括大陆溢流玄武岩(CFB)及其伴生侵入岩、大洋盆地溢流玄武岩、海底高原、海岭、海山群、被动火山边缘和火山链等.大火成岩省作为地球上目前所知的最大的火山作用，其形成过程记录了地球在某一特定历史时期，物质和能量由地球内部向地表迁移，在很短时间内完成巨大的岩浆喷出量；其成因明显不同于洋底扩张过程，这个过程要求在地幔深部存在相应巨大的热异常.传统的板块构造理论主要涉及的是板块间的横向运动的动力学过程，其理论核心不适用于大火成岩省的形成机制，无法对产生大规模的镁铁质火成岩成因和动力学机制给出很好的解释.地幔柱理论的应运而生为其提供了可行的动力学机制，地幔柱的研究逐渐成为研究和探索大火成岩省相关的一系列地质科学关键问题的重要手段之一.地幔柱学说的出现弥补了板块构造理论强调水平运动而忽略垂向运动的缺点，它与板块构造理论共同构成了地球构造理论中两个不可或缺的相互补充完善的内容.
峨眉山大火成岩省(ELIP)主要分布于我国的扬子克拉通西苑(云南、四川和贵州三省境内)，是存在于我国境内的唯一得到广泛认可的大火成岩省.峨眉山大火成岩省以溢流的晚二叠纪玄武岩为岩性主体，其在空间分布上南北纵贯1000 km，东西横跨超900 km(宋谢炎等，2002)，出露面积约为2.5×105~3.0×105 km2，体积为0.3×106~0.6×106 km3(Thompson et al., 2001; Wignall，2011)，何斌等(2003)考察了峨嵋山玄武岩下覆茅口灰岩，根据其剥蚀差异及类圆形的等厚线分布，将峨眉山大火成岩省自西向东分为深度剥蚀带(内带)、部分剥蚀带(中带)、短暂沉积间断带或古风化壳(外带)以及连续沉积带等.内带是直径达400 km的圆形区域，中带为内带边缘以外的宽度300 km的弧形环带.李宏博(2010)对峨眉山大火成岩省进行区域调研和野外探测，在区调资料的基础上同时结合野外观察数据，推测其中心收敛于永仁一带.峨眉山玄武岩具有与多数大陆大火成岩省相似的结构分布，并未形成由于柱尾的大规模熔融而在地表形成线状火山链结构，而是地幔热物质在上升过程中大量熔融，一部分沿着岩石圈软弱带或破裂带上升至地表，对原始岩体进行覆盖，形成大规模的上覆溢流玄武岩，这可能是地幔柱的柱头熔融为熔体主导的结果.这与本文模拟结果符合，地幔柱上升过程过程中，柱头熔融程度较大，柱尾则随着地幔柱上升的动力学过程未形成大规模熔融.另一部分地幔柱柱头处的熔体则通过底部侵蚀的方式对下地壳底部进行侵蚀作用，在其冷却后成为地壳的一部分.
徐义刚等(2002)等根据分异熔融反演模型，采用线性的橄榄岩的固、液相线设置.地幔柱从地球深部开始近似绝热上升，从地球化学的角度初步估算出形成峨眉山玄武岩所需的地幔温度大于1550 ℃(1823 K)，根据其物质成分判断其熔融始于约140 km处，终于约60 km处，最大平均熔融程度为16.5%左右.本文中以上下地壳分别为花岗岩和镁铁质麻粒岩，地幔为橄榄岩的模型对峨眉山地幔柱的熔融情况进行模拟.模型设置地幔柱初始温度1903 K时，线性固相线模拟所得结果初融深度159 km，纵向作用深度67 km，局部最大熔融程度32%，这与地球化学反演结果并不完全一致，其原因是地球化学反演中地幔柱初始温度设置较低，且采用线性固液相线进行反演导致其结果具有一定的局限性.与此同时，本文中熔融情况的数值模拟仅展现的是不含断裂带情况下，地幔柱与岩石圈相互作用过程中的熔融程度的变化.而在实际情况中，由于构造复杂性导致的物质熔融和喷发后在地表产生的物质的组成等方面，都将与本文的不含断裂无喷发的模拟结果存在着一定程度的差异.对含断裂带模型设置的熔融情况的讨论，有待进一步的研究.但总体而言，对壳幔作用中涉及的各圈层物质的熔融过程进行数值模拟，能够帮助我们更好地确定，或至少在一定程度上限定不同大火成岩省的原始岩浆的起源位置的深度范围及其所可能存在的物质成分.
4 结论数值模拟能够动态的反映出在不同的地球动力学过程中，不同圈层的物质熔融程度的动态变化规律.通过对不同圈层物质熔融程度的计算以及对熔融相关的模拟结果的分析和综合讨论，可以得出以下几点：
(1) 地幔柱与岩石圈相互作用过程中的熔融情况可以分为三个阶段：地幔柱的初融阶段，地幔柱自身熔融占主导，减压熔融为主因；地幔柱纵向侵蚀岩石圈阶段，岩石圈地幔发生了由升温主导的部分熔融，地幔柱由于压力的持续减小继续熔融，物质整体的熔融程度增加；地幔柱的横向展平阶段，岩石圈地幔随着地幔柱的扩展其熔融范围增加，以升温熔融为主，地幔柱自身熔融程度减小.
(2) 线性固、液相线温度跨度较大，但覆盖的压力范围较小，因此它能够较好的模拟浅部物质的熔融情况，而对于较深的地幔物质的熔融，线性模型存在一定的局限；多项式固、液相线弥补了线性模型对深部模拟的缺陷，温压关系的表达较为全面，对于地幔柱物质在深部上升过程中所产生的熔融能够很好的展现.
(3) 对于地幔柱物质温度的小程度变化，利用多项式模型模拟的熔融程度结果的响应比线性模型的模拟结果更加明显.利用多项式模型模拟的熔融结果中，低温地幔柱的初始熔融深度、纵向作用深度以及最大熔融程度相比高温地幔柱均改变了50%左右；线性模型中改变量在20%~30%之间.高温地幔柱的模拟结果中，初始熔融位置较深，纵向作用深度较浅，发生在地幔柱中心的熔融较为彻底，最大熔融程度较大.
(4) 模拟得到峨眉山大火成岩省的初始熔融深度在159~176 km之间，纵向作用深度范围在70 km左右，最大熔融程度最高可达39%.但此地区区域构造复杂，深大断裂带的分布导致本文中的结果还存在一定的局限.对于通过在岩石圈内预设断裂带或软弱带的方式，模拟在地幔柱与含断裂的岩石圈相互作用以及沿破裂带上升喷出地表的动力学过程中，不同物质的熔融情况有待进一步研究.
Burov E, Poliakov A.
2001. Erosion and rheology controls on synrift and postrift evolution:verifying old and new ideas using a fully coupled numerical mode. Journal of Geophysical Research, 106(B8): 16461-16481.
DOI:10.1029/2001JB000433 |
|
Burov E, Guillou-Frottier L.
2005. The plume head-continental lithosphere interaction using a tectonically realistic formulation for the lithosphere. Geophysical Journal International, 161(2): 469-490.
DOI:10.1111/gji.2005.161.issue-2 |
|
Burov E, Guillou-Frottier L.
2007. Plume head-lithosphere interactions near intra-continental plate boundaries. Tectonophysics, 434(1-4): 15-38.
DOI:10.1016/j.tecto.2007.01.002 |
|
Campbell I H, Griffiths R W.
1990. Implications of mantle plume structure for the evolution of flood basalts. Earth and Planetary Science Letters, 99(1-2): 79-93.
DOI:10.1016/0012-821X(90)90072-6 |
|
Coffin M F, Eldholm O.
1994. Large igneous provinces:Crustal structure, dimensions, and external consequences. Reviews of Geophysics, 32(1): 1-36.
DOI:10.1029/93RG02508 |
|
Christensen U, Harder H.
1991. 3D convection with variable viscosity. Geophysical Journal International, 104(1): 213-226.
DOI:10.1111/gji.1991.104.issue-1 |
|
d'Acremont E, Leroy S, Burov E B.
2003. Numerical modelling of a mantle plume:The plume head-lithosphere interaction in the formation of an oceanic large igneous province. Earth and Planetary Science Letters, 206(3-4): 379-396.
DOI:10.1016/S0012-821X(02)01058-0 |
|
Farnetani C G, Richards M A.
1994. Numerical investigations of the mantle plume initiation model for flood basalt events. Journal of Geophysical Research:Solid Earth, 99(B7): 13813-13833.
DOI:10.1029/94JB00649 |
|
Gasparik T.
1990. Phase relations in the transition zone. Journal of Geophysical Research, 95(B10): 15751-15769.
DOI:10.1029/JB095iB10p15751 |
|
Griffiths R W, Campbell I H.
1990. Stirring and structure in mantle starting plumes. Earth and Planetary Science Letters, 99(1-2): 66-78.
DOI:10.1016/0012-821X(90)90071-5 |
|
He B, Xu Y G, Xiao L, et al.
2003. Generation and spatial distribution of the Emeishan large igneous province:New evidence from stratigraphic records. Acta Geologica Sinica (in Chinese), 77(2): 194-202.
|
|
Huang J Q, Ren J S, Jiang C F.
1977. An outline of the Tectonic chacarteristics of China. Acta Geologica Sinica (in Chinese), 51(2): 117-135.
|
|
Kellogg L H, King S D.
1997. The effect of temperature dependent viscosity on the structure of new plumes in the mantle:Results of a finite element model in a spherical, axisymmetric shell. Earth and Planetary Science Letters, 148(1): 13-26.
|
|
Kirby S H.
1983. Rheology of the lithosphere. Reviews of Geophysics, 21(6): 1458-1487.
DOI:10.1029/RG021i006p01458 |
|
Lavelle J W.
1997. Buoyancy-driven plumes in rotating, stratified cross flows:Plume dependence on rotation, turbulent mixing, and cross-flow strength. Journal of Geophysical Research:Oceans, 102(C2): 3405-3420.
DOI:10.1029/96JC03601 |
|
Leng W, Gurnis M.
2012. Shape of thermal plumes in a compressible mantle with depth-dependent viscosity. Geophysical Research Letters, 39(5): L05310.
|
|
Li H B, Zhang Z C, Lu L S, et al.
2010. Geometry of the mafic dykes warms in Emeishan large igneous province:Implications for mantle plume. Acta Petrologica Sinica (in Chinese), 26(10): 3143-3152.
|
|
Li J K, Wang D H.
2004. A preliminary research on the numerical simulation of the emeiMantle plume. Acta Geoscientica Sinica (in Chinese), 25(5): 509-514.
|
|
Li J K, Wang D H.
2005. Advances in the numerical simulation of the mantle plume. Geological Science and Technology Information (in Chinese), 24(4): 13-20.
|
|
Lu J R.
1996. Dynamical characteristics of the Emei mantle plume. Acta Geoscientica Sinica (in Chinese), 17(4): 424-438.
|
|
Manglik A, Christensen U R.
1997. Effect of mantle depletion buoyancy on plume flow and melting beneath a stationary plate. Journal of Geophysical Research:Solid Earth, 102(B3): 5019-5028.
DOI:10.1029/96JB03623 |
|
McKenzie D, Bickle M J.
1988. The volume and composition of melt generatted by extension of the lithosphere. Journal of Petrology, 29(3): 625-679.
DOI:10.1093/petrology/29.3.625 |
|
Meng W J, Chen Z A, Bai W M.
2015. Numerical simulation on process of the plume-lithosphere interaction. Chinese Journal of Geophysics (in Chinese), 58(2): 495-503.
|
|
Moresi L, Dufour F, Mühlhaus H B, et al.
2003. A lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. Journal of Computational Physics, 184(2): 476-497.
DOI:10.1016/S0021-9991(02)00031-1 |
|
Morgan W J.
1971. Convection plumes in the lower mantle. Nature, 230(5288): 42-43.
DOI:10.1038/230042a0 |
|
Morgan W J.
1981. Hotspot tracks and the opening of the Atlantic and Indian Oceans. The Sea, ideas and observations on progress in the study of the seas, 7: 443-487.
|
|
O'Neill C, Moresi L, Müller D, et al.
2006. Ellipsis 3D:A particle-in-cell finite-element hybrid code for modelling mantle convection and lithospheric deformation. Computers & Geosciences, 32(10): 1769-1779.
|
|
Ohtani E, Kato T, Sawamoto H.
1986. Melting of a model chondritic mantle to 20 GPa. Nature, 322: 352-353.
DOI:10.1038/322352a0 |
|
Ranalli G. 1995. Rheology of the Earth, 2nd edition, Chapman and Hall, London, p. 413.
|
|
Richards M A, Duncan R A, Courtillot V E.
1989. Flood basalts and hot-spot tracks:Plume heads and tails. Science, 246(4926): 103-107.
DOI:10.1126/science.246.4926.103 |
|
Song X Y, Hou Z Q, Wang Y L.
2002. The mantle plume features of Emeishan basalts. Journal of Mineralogy and Petrology (in Chinese), 22(4): 27-32.
|
|
Sun Y J, Dong S W, Fan T Y.
2013. 3D rheological structure of the continental lithosphere beneath China and adjacent regions. Chinese Journal of Geophysics, 56(9): 2936-2946.
|
|
Takahashi E, Kushiro I.
1983. Melting of a dry peridotite at high pressures and basalt magma genesis. American Mineralogist, 68(9-10): 859-879.
|
|
Thompson G M, Ali J R, Song X Y, et al.
2001. Emeishan basalts, SW China:Reappraisal of the formation's type area stratigraphy and a discussion of its significance as a large igneous province. Journal of the Geological Society, 158(4): 593-599.
DOI:10.1144/jgs.158.4.593 |
|
Turcotte D L, Schubert G. 2002.
Geodynamics, applications of Continuum Physics to Geological Problems, 2nd edition, Cambridge University Press, Cambridge. (2nd edition). Cambridge: Cambridge University Press.
|
|
Weertman J.
1970. The creep strength of the Earth's mantle. Reviews of Geophysics, 8(1): 145-168.
DOI:10.1029/RG008i001p00145 |
|
Wignall P B.
2001. Large igneous provinces and mass extinctions. Earth-Science Reviews, 53(1-2): 1-33.
DOI:10.1016/S0012-8252(00)00037-4 |
|
Wilson J T.
1963. A possible origin of the Hawaiian Islands. Canadian Journal of Physics, 41(6): 863-870.
DOI:10.1139/p63-094 |
|
Xu Y G.
2002. Mantle plumes, large igneous provinces and their geologic consequences. Earth Science Frontiers (in Chinese), 9(4): 341-353.
|
|
Zhang B R.
2001. Magmatic activities from plume-source in the Qinling orogenic belt and its dynamic significance. Earth Science Frontiers (in Chinese), 8(3): 57-66.
|
|
何斌, 徐义刚, 肖龙, 等.
2003. 峨眉山大火成岩省的形成机制及空间分布:来自沉积地层学的新证据. 地质学报, 77(2): 194–202.
|
|
黄汲清, 任纪舜, 姜春发, 等.
1977. 中国大地构造基本轮廓. 地质学报, 51(2): 117–135.
|
|
李宏博, 张招崇, 吕林素.
2010. 峨眉山大火成岩省基性墙群几何学研究及对地幔柱中心的指示意义. 岩石学报, 26(10): 3143–3152.
|
|
李建康, 王登红.
2004. 峨眉地幔柱动力学数值模拟的初步研究. 地球学报, 25(5): 509–514.
|
|
李建康, 王登红.
2005. 地幔柱数值模拟研究进展. 地质科技情报, 24(4): 13–20.
|
|
卢记仁.
1996. 峨眉地幔柱的动力学特征. 地球学报, 17(4): 424–438.
|
|
蒙伟娟, 陈祖安, 白武明.
2015. 地幔柱与岩石圈相互作用过程的数值模拟. 地球物理学报, 58(2): 495–503.
DOI:10.6038/cjg20150213 |
|
宋谢炎, 侯增谦, 汪云亮, 等.
2002. 峨眉山玄武岩的地幔热柱成因. 矿物岩石, 22(4): 27–32.
|
|
孙玉军, 董树文, 范桃园.
2013. 中国大陆及邻区岩石圈三维流变结构. 地球物理学报, 56(9): 2936–2946.
DOI:10.6038/cjg20130908 |
|
徐义刚.
2002. 地幔柱构造、大火成岩省及其地质效应. 地学前缘, 9(4): 341–353.
|
|
张本仁.
2001. 秦岭地幔柱源岩浆活动及其动力学意义. 地学前缘, 8(3): 57–66.
|
|