幔源岩浆在地壳中分层侵位的控制因素：二维热-力学模拟.pdf

书书书 （）：岩石学报 ：崔晓娜，陈林 幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟 岩石学报，（）：，：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟崔晓娜，陈林 ，中国科学院地质与地球物理研究所，岩石圈演化国家重点实验室，北京 中国科学院大学地球与行星科学学院，北京 ，收稿，改回 ：，（）：，：，（，）：（）（），（），（），（）（），；摘要幔源岩浆在地壳内的上升和聚集样式不仅依赖于岩浆自身的性质，还取决于围岩的强度和热状态。已有数值和物理模型大多关注岩浆自身物性对其上升过程的影响，而对围岩流变强度或热状态如何影响岩浆的上升和聚集过程的研究相对薄弱，尤其是周期性的幔源岩浆在壳内分层侵位的受控因素仍然不清楚。本文利用二维热 力学数值模拟方法，通过发展本文受国家自然科学基金项目（）和国家重点研发计划项目（）联合资助第一作者简介：崔晓娜，女，年生，博士生，地球动力学专业，：通讯作者：陈林，男，年生，研究员，主要从事地球动力学数值模拟研究，：多期岩浆脉冲和岩墙生成算法，研究了岩浆从深部地幔上升至地壳内部侵位的动力学过程，系统测试了地壳（围岩）地温梯度（上、下地壳的地温梯度分别以 和 表示）和地壳强度对岩浆上升过程和聚集样式的影响。模拟结果表明：（）地壳地温梯度对岩浆的侵位深度有重要影响，岩浆侵入冷地壳（），岩浆主体在岩石圈深度聚集，地表的相对高差小于 ；岩浆侵入温地壳（）在下地壳底部聚集形成岩浆房，上升至上、下地壳界面，岩浆房上方的地表地形呈现中心拗陷两翼隆起的形态，地表最大高程可达 ；岩浆侵入热地壳（）仅在下地壳底部聚集形成岩浆房，地形演化特征与温地壳背景的情况类似，但最大高程小于 。（）在同等地壳地温梯度条件下，上、下地壳的相对强度决定了岩浆的聚集样式：下地壳的强度越弱，岩浆更易在下地壳聚集形成岩浆房；上地壳的强度越弱，岩浆更易在地表喷发。进一步分析表明：岩浆的聚集样式受地壳地温梯度与地壳流变分层性的共同控制，地壳越热且流变强度分层性越显著则越有利于岩浆在地壳中的多层级侵位；每一期岩浆脉冲的供给均会导致岩浆房内部的超压值骤增。我们的模拟结果对理解火山喷发前壳内岩浆的赋存状态及岩浆活动区的地形演化具有启示意义。关键词岩浆；地壳；地温梯度；粘度；数值模拟中图法分类号 上地幔部分熔融形成的低粘度、低密度的基性岩浆在地壳中聚集和侵位的过程通常会伴随岩浆房、岩墙或岩席的形成（，；，）。岩浆的侵位和喷发往往会导致地震、海啸等自然灾害的发生（，；，），也会形成铜矿、金矿等金属矿床（，；杨林等，）。因此，理解幔源岩浆从深部上升至地表的动力学过程具有重要的科学意义。高分辨率地球物理成像结果揭示，我国东北长白山火山下方的地壳内部存在多层级的岩浆聚集区（，；，）。岩石地球化学数据表明，地表出露的玄武岩呈现不同期次的岩浆补给（郭文峰等，；，）、岩浆房多次聚集的特征（，；，），且在地表喷发的火山数量远小于地下赋存岩浆的体积（，）。由此可见，幔源岩浆在地壳中的侵位具有多期次岩浆供给和分层聚集的特点。然而，幔源岩浆侵位过程中产生的分层聚集样式受何种因素控制仍不清楚。针对幔源岩浆的上涌过程，前人已开展了大量数值和物理模拟研究。早期数值模拟研究主要关注岩浆侵位对围岩带来的热影响（，；，），或是探讨单期幔源岩浆在均质地壳的聚集过程（，；，）。近期的数值模拟研究主要关注岩浆的传输通道（，；，）。围岩的流变强度及其分层性决定了岩浆的传输方式，包括熔融底辟、减压通道或张性裂隙（，），而岩墙和底辟可以共同促进岩浆向上传播（，）。伴随岩浆侵位，侵入体周围发育的裂隙由初期呈径向展布的剪裂隙转换为向上扩展的主干张性断裂（卢靖雯等，），同时岩浆侵位会对围岩施加附加变形，导致近场与远场的小尺度变形（常成和罗纲，）。（）利用三维热 力学模拟方法探究地壳流变和温度结构对基性岩侵位形态的影响，最近 （）考虑了岩浆熔融抽取过程中的岩浆成分演化，模拟壳内多期次玄武岩席聚集形成岩浆房的机制。物理模拟研究早期主要探究岩墙向岩席的转换机制，通过将流体注入固体明胶的模拟实验，验证了当固体明胶存在刚度对比层、层间弱接触面或外部应力时，可以导致垂向传播的裂隙运动轨迹转为水平向（，；，）。最近，（）利用物理模拟方法探究了岩浆沿岩席长距离水平传播的形成机制，发现岩浆的流变学性质具有重要影响。然而，物理模拟受相似比的限制，难以直接运用于理解实际的地质过程（，）。综上，前人的相关研究要么没有考虑岩浆的多期补给，要么没有考虑岩墙作为不同层级岩浆房之间的连接通道，尤其是针对多期幔源岩浆在地壳中分层侵位过程及其控制因素的研究仍然很薄弱。在本研究中，我们聚焦幔源岩浆从深部到地壳的聚集运移过程，通过在二维热 力学数值模拟程序中引入多期岩浆脉冲和岩墙生成算法，模拟幔源岩浆从深部地幔上升到地壳的动力学过程，系统测试地壳地温梯度和流变强度的影响，探究幔源岩浆在地壳内部分层侵位的控制因素。我们将首先介绍数值模拟方法与模型设计；然后描述地壳地温梯度和地壳强度对壳内岩浆赋存状态的影响；最后，讨论幔源岩浆在地壳内聚集样式的控制因素。数值模拟方法我们使用开源的二维热 力学 程序（，）模拟幔源岩浆多期次侵位的动力学过程。本程序在求解动量守恒方程、质量守恒方程与能量守恒方程时，联合交错网格有限差分方法和粒子标记（，）技术。在节点上解方程，在 上传递物性参数，可以有效避免利用有限差分方法求解连续性方程时可能引起的数值发散问题。同时，该方法考虑了岩石的粘 弹 塑性流变本构关系和岩石的部分熔融，可以在时间与空间尺度上精细模拟岩浆侵位的动力学过程。控制方程（）质量守恒 连续性方程 岩石学报 ，（）：（）式中，和 分别代表速度矢量的水平和垂直分量。（）动量守恒 斯托克斯方程 （）（，）（）式中，和 分别代表水平和垂直坐标；是密度，为压力 、温度 、岩石成分 和部分熔融百分数 的函数；是重力加速度；、分别为偏应力张量，其本构关系为：（）()（）式中，、为应变率张量；代表有效粘度。（）能量守恒 热传导方程 （，）（，）（）（）（）式中，是等压热容；代表温度对时间的物质导数；和分别代表水平和垂直方向的热流；是热导率，与压力 、温度 和岩石成分 有关（，）；是热膨胀系数；代表放射性生热率，依赖于岩石成分；代表剪切生热，表示形变过程中由机械能转化的热能；代表绝热生热，表示随着压力变化而产生的热量；代表潜热，表示岩石发生相变过程中释放或消耗的热。岩石流变本程序考虑了岩石的粘 弹 塑性（）流变，体偏应变率 包括三个分量：（）（）（）（）（）（），（）式中，为剪切模量；为偏应力的目标共旋转时间导数；为岩石的屈服应力；为第二偏应力不变量 ()；为塑性势；为塑性因子。以上方程中的偏应力分量 是基于粘 弹 塑性本构关系以显式的一阶时间有限差分形式表示为粘 弹性应力的目标时间导数：（），（）式中，为弹性时间步长；是上一个时间步长由 技术求解并经过平流和旋转后的偏应力张量（，；，）；为粘 弹因子；为有效的粘塑粘度。()（）其中，式中，是气体常数；、和 是岩石物理实验约束的流变参数，分别代表物质常数、活化能、活化体积和应力指数；为 粘 度 因 子，用 于 线 性 改 变 物 质 的 有 效 粘 度（，）。由于实验室条件下物质成分的变化，所测得的流变参数也具有很大的不确定性。为了将这些结果应用到实际地质过程，我们使用粘度因子对物质在实验室条件下测得的有效粘度进行线性调节（刘丹红和陈林，）。在这里，我们通过改变粘度因子调控地壳的强度。熔融岩石（，为岩石的部分熔融百分数）的粘度为 。一方面，我们假定岩浆上升是一个十分快速的地质过程，在该过程中没有明显的热损失，因此粘度基本没有变化；另一方面，将部分熔融岩石的粘度取为常数，可以简化一些复杂因素的影响，便于我们重点考察围岩的粘 弹 塑性流变对幔源岩浆上升过程带来的影响。此外，在模拟过程中，我们设定粘度的上限为 ，下限为 ，采用的流变参数见表 。物质的塑性变形遵从 屈服准则（，）：（）（）其中 是岩石的内聚力；是有效内摩擦角。崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟表 数值模型中用到的主要物理参数 物质流变）（）（）（）（）（）（）（）（）（）（）（）空气 上地壳湿石英 （）下地壳斜长石 （）岩石圈干橄榄石 （）软流圈干橄榄石 （）熔体 （）基性岩浆 （）参考文献），注：其他通用物质参数：热膨胀系数 ，压缩系数 ；）（）（），式中 代表岩石在 和 条件下的密度；）（）；（）；（）；）当 ，（）（）；当 ，；当 ，（）（）；当 ，；当 ，；当 ，（）；当 ，（）；当 ，；）（）；（）；（）；（）；（）；（）；（）；（）；（）；（）部分熔融根据岩石物理实验数据的约束，程序考虑了岩石的部分熔融行为（，；，）。假设部分熔融百分数 与温度为简单的线性关系：，（）式中，为岩石的固相线，为岩石的液相线，各种岩石的固相线和液相线见表 。部分熔融岩石的有效密度 取决于部分熔融百分数：（）（）式中，和 分别代表岩石的固相密度和熔融密度，它们都是关于温度和压力的函数，计算表达式见表 。其中，我们对幔源岩浆密度的选取参照了 （）的结果，并且与该作者汇编的结果进行了对比，在 的压力范围内，我们的结果中岩浆的密度与前人的实验结果和理论计算结果均一致（电子版附图 ），说明岩浆密度的选取具有合理性。部分熔融对有效热容量（）和热扩散系数（）的影响由下式给出（，）：()()()()（）多期岩浆脉冲的实现火山喷发需要体积足够大、低粘度、低密度、部分熔融且持续供应的岩浆热源，以提供浮力、维持高温（，）。我们通过预设岩石圈之下的岩浆房来模拟幔源岩浆上升的动力学过程，在已有的二维粘 弹 塑性数值模拟程序（，）的基础上加入多期岩浆脉冲算法，在岩石圈以下多期次、定量的稳定提供岩浆热源，尽量保证岩浆通道不变形。算法实现细节为：（）我们添加新一期岩浆脉冲的判断条件为：深度 以下，部分岩浆通道和岩浆源区中熔体的总含量小于 ；（）当熔体含量满足上述条件时，我们用性质相同但标号不同的 补充新一期次的幔源岩浆；（）在模拟过程中，我们总共提供六期岩浆脉冲，为岩浆的持续上涌提供动力。岩墙的生成算法岩墙作为拉伸作用形成的薄弱带，它的形成与围岩的应力场有关（，；，）。在本文中，拉伸应力标记为正、挤压应力标记为负。岩墙作为岩浆房之间的连接通道（，；，），岩浆沿岩墙传播（，；，），其宽度随深度变化（，）。目前对于岩墙的热 力学模拟主要探究岩墙对岩浆上升过程的影响，例如，（）认为花岗岩岩浆房会形成辐射状岩墙，为了探究岩浆上升过程中岩墙和底辟的相互作用，因此设定临界的应力和第二应力不变量作为垂直状岩墙形成的条件。（）通过岩石热 力学模拟方法，考虑岩墙 岩石学报 ，（）：图 岩墙形成的示意图 作为玄武岩熔融提取的通道。我们通过探究幔源岩浆侵位对地壳应力状态带来的影响，认为垂直状岩墙是岩浆上升的通道，辐射状岩墙是岩浆房体积扩张时围岩的同期破裂，表征岩浆房的发育。因此，本文中考虑岩墙是在岩浆侵位对上覆围岩造成拉伸环境中形成的（，），作为岩浆向浅部运移的通道。我们在程序中加入岩墙生成算法，岩墙的形成和扩展由岩浆房超压值和地壳应力共同控制。岩浆房超压值为岩浆房内部的总压强与静岩压力的差值（岩浆房内部的超压值相差无几），岩浆房内部的超压值会随着每一期次的岩浆补给而增加、随着岩浆房的体积增大而减小，超压值演化表征了岩浆房的发育过程；当岩浆底侵时，岩浆聚集使得岩浆房之上产生应力集中的区域，即岩浆房的发育对上覆的上地壳与下地壳形成了水平的拉伸应力，这种拉伸应力即为地壳应力，地壳的应力场决定了岩墙垂直向发育。以下地壳岩墙的形成过程为例，具体算法如下：（）首先，由于多期次的岩浆脉冲供应，壳内岩浆房发育至一定体积，将岩浆房内部超压值小于 作为岩浆房发育结束的时刻与判断是否形成岩墙的初始时刻。（）随后，当岩浆房之上的地壳应力大于 时，发育长度为岩浆房顶部至上、下地壳界面、宽度为 的岩墙（见图）。经过多次测试，宽度的岩墙是保证岩浆上升的最窄宽度，过窄的岩墙会导致通道快速冷却、无法输运岩浆，因此我们在模拟中设定岩墙的最小宽度为 。模型设计为了模拟幔源岩浆上升的动力学过程，我们设计了如图所示的初始模型。模型在水平方向上的宽度为 ，垂向深度为 ，包含 上地壳、下地壳和 岩石圈地幔，代表了一个典型的大陆岩石圈，具体物性参数见表 。我们初始在岩石圈底部预设一个半径为 、温度为 的圆形幔源岩浆区域，在顶部设置了一个宽度为 的岩浆通道，连接到地壳底部。整个模型空间被 个节点组成的规则网格离散化，空间分辨率为 ，总共 万个粒子（）随机分布在网格空间中，粒子分辨率高达 。模型的所有力学边界均为自由滑边界。在模型的上部，我们设计一个厚 的伪空气层，其密度和粘度分别为 、，用于模拟岩浆上升过程中产生的地表起伏（，；，）。模型的热边界条件：顶边界固定为 ，侧边界为绝热边界，底边界为 固 定 为 （，；，），整个空气层温度固定为 ，岩石圈底边界以下为绝热地温梯度 。为了区分地壳的不同热状态，我们采用分段函数赋予了地壳两种初始温度结构：第一种温度结构初始固定莫霍面的温度，设定地表和莫霍底界面的初始温度分别为 和 ，通过改变上、下地壳界面的初始温度（、）来调节上、下地壳的地温梯度（、），让温度随深度线性增加；第二种温度结构初始固定岩石圈底部的温度，设定地表和岩石圈底界面的初始温度分别为 和 ，通过改变莫霍面的初始温度（、）来调节地壳和岩石圈的地温梯度，同样让温度随深度线性增加。模拟结果为了探究来自地幔深部的多期岩浆在地壳中分层聚集的动力学过程，我们基于多期岩浆脉冲和岩墙算法的 程序，开展了 组数值模拟实验（表 ）。在模拟中，系统测试了地壳地温梯度和地壳强度对岩浆上涌过程的影响。在这里，我们选择上地壳的粘度因子 、下地壳的粘度因子 、的模型作为参考模型（即表 中的 ）。参考模型图 为参考模型（，）的演化结果。来自岩石圈底部的幔源岩浆穿过岩石圈到达莫霍面，初始阶段（；图 ），岩浆在莫霍面之上聚集形成椭圆状岩浆房，伴随同时的垂向增厚和横向扩展，粘度结果显示岩浆房周围出现局部的低粘度区域（图 ），表征岩浆房的扩张趋势。随后，由于岩浆侵位对上覆地壳造成的拉伸力，下地壳形成岩浆通道使得岩浆由此上升（图 ），水平拉伸应力瞬时消散（图 ）。时，岩浆上涌至上、下地壳界面，多期次岩浆脉冲带来的高超压导致岩浆通道不断变宽，在岩浆房内部，二期岩浆及多期岩浆以韧性底辟形式对流，促进多期岩浆混合及岩浆崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟图 模型设计（）为模型全空间；（）为（）中黑色方框的放大区域岩性色标见标注，不同的紫色表示六期基性岩浆，不同的绿色表示六期熔体，蓝色为冷却后的基性岩浆，红色为冷却后的熔体 （）；（），表 数值模拟实验 模型 熔体含量（）（）（）（）相关图件 图 、附图 附图 附图 附图 图 、附图 附图 附图 附图 附图 附图 附图 岩石学报 ，（）：续表 模型 熔体含量（）（）（）（）相关图件 附图 附图 图 、附图 附图 图 、附图 附图 附图 附图 附图 图 图 图 房上升（图 ），粘度结果显示岩浆房之上局部的低粘度区域（图 ），表征岩浆的上升趋势。最终 时，上、下地壳界面之上岩浆聚集形成扁平状岩浆房，但岩浆房顶部的应力累积不足（图 ），达到稳态至形成双层聚集的岩浆房，此时下地壳岩浆房体积扩张，并且内部卷入极少量的地壳和岩墙等物质（图 ），上地壳与下地壳减薄程度不同的原因在于两者对于岩浆侵位带来的空间问题响应不同。在 的时间尺度内，岩浆房内部保持活跃，进行持续的岩浆对流。图 中 为参考模型的水平偏应力随时间的变化，下地壳与上地壳岩浆侵位均对岩浆房之上的局部区域形成水平拉伸应力。上地壳累积的拉伸应力在脆韧性转换深度处最大，向地表逐渐降低，向深部为韧性区应力累积很弱；下地壳累积的拉伸应力在上、下地壳界面处最大，向深部减小，同样韧性区应力累积很弱，并且下地壳岩浆侵位导致上、下地壳界面隆起，可以被上地壳韧性区容纳（图 ）。同时，我们选取参考模型（）中三个固定点的水平偏应力结果来描述岩浆上涌对下地壳、上地壳、地表带来的影响，在图 中用红色的圆标记。下地壳岩浆房之上红色固定点拉伸力的初始波峰表示下地壳岩浆聚集的初始阶段，第二个 小波峰表示 添加的第二期岩浆脉冲，时拉伸应力由 骤降表征岩浆房破裂，表征该深度已经转变为岩墙或岩浆成分（见图 ）。与之相似，上地壳脆韧性转换深度处，应力的第一个波峰表示下地壳岩浆聚集过程中上地壳的同期应力累积，第二个 小波峰表示添加的第二期岩浆脉冲，第三个 的波峰表示上地壳岩浆的聚集过程，随后应力逐渐消散。岩浆初始聚集时地表应力的累积滞后，随后在岩浆开始沿岩墙向上传播至源区六期岩浆脉冲结束的过程中，浅表一直处于拉伸环境，当岩浆脉冲停止后，应力出现负值表示地壳浅部转为挤压状态。我们选择下地壳岩浆房底部与上地壳岩浆房底部的两个固定点从超压演化分析岩浆房的发育过程，在图 中用白色的圆标记。下地壳岩浆房内的超压因一期岩浆初始聚集而骤增（图 ），后因岩浆房的不断发育，超压值随着岩浆房体积增大而持续减小，后五期的岩浆脉冲均可引起超压值骤增，使岩浆房维持高超压状态而体积不断扩张，时超压值骤降表征岩浆房的破裂（见图 ），此时为下地壳岩浆房发育结束的时刻与形成下地壳岩墙的时刻，最终下地壳岩浆房内部的超压值趋于稳定，当岩浆脉冲停止后超压负值表征岩浆上升停滞。上地壳岩浆房内部的超压初始值对应下地壳岩浆房破裂的时刻（图 ），超压值由负转正表征上、下地壳界面深度的初始岩浆聚集，后四期岩浆脉冲引起的超压最大值小于 ，当岩浆脉冲停止后超压负值亦表征岩浆上升停滞。地壳地温梯度的影响在参考模型的基础上，我们测试了地壳地温梯度的影响。如图 ，设定莫霍面的初始温度分别为 、（、；见表 ），根据单一变量测试原则，其余参数设置与参考模型相同。首先，我们对比 时冷地壳、温地壳与热地壳背景下岩浆的赋存状态。图 中，冷地壳的地温梯度为 ，低的地温梯度使得地壳的粘度过大（），即使岩浆高超压也无法克服地壳的强流变实现岩浆聚集，此时岩浆主要在岩石圈底部聚集。图 中，温地壳的地温梯度为 ，在地壳具有较高的地温梯度时，岩浆在莫霍面之上聚集形成椭圆状岩浆房，并且岩浆可以上升至上、下地壳界面。图 中，热地壳的地温梯度为 ，高的地温梯度使得地壳的粘度降低，此时岩浆仅在莫霍面之上聚集形成席状岩浆房，产生以上现象是基于岩崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟图 参考模型（）（）表示的是成分场随时间的演化（；见表 ），图中岩石色标同图 ，标记的白色圆为两个超压记录点（）（）表示的是粘度场随时间的演化（）（）为水平偏应力随时间的变化，标记的红色圆为三个应力记录点 图（）为相应三个记录点的水平偏应力随时间的演化 （）（）（；），（）（）（）（）（），（）（）（），（）浆的正浮力与相应地壳地温梯度下的地壳强度共同作用的结果。我们进一步分析了岩浆侵位对地形带来的影响，图 分别为冷地壳、温地壳、热地壳背景下地形随时间的变化，温地壳与热地壳的地形演化呈现相同的趋势：侵位均导致地表抬升，初始中心隆起后转变为中心凹陷两侧隆起形态（图 ，），而冷地壳时，地表的地形无明显起伏（图 ）。速度骤增均对应多期岩浆脉冲的岩浆供应时刻（图 ）。图 为 时的地形垂直剖面，地表最大隆起值发生在温地壳下，地表最高高程约为 ，热地壳次之，地表最高高程约为 ，冷地壳时，地表的相对高差小于 。综上，地壳地温梯度对岩浆的侵位深度具有重大影响，岩石学报 ，（）：图 岩浆房内超压值随时间的演化（）为下地壳岩浆房的超压值随时间的演化；（）为上地壳岩浆房的超压值随时间的演化 选取的超压点位置见图 （）；（）图 地壳地温梯度对岩浆侵位深度的影响（）（）分别对应冷地壳、温地壳、热地壳（、；见表 ）背景下的模拟结果（）（）为相应的粘度结果（）为初始的地温梯度 岩石色标同图 ，粘度色标同图 （）（）、（，；），（）（）（），温地壳有利于岩浆上升到更浅的深度，地壳的地温梯度过高或过低都不利于岩浆的上升。并且，在岩浆多期次侵位的过程中，地形多期次抬升，壳内岩浆房的高度对地表的隆起高度贡献更大。地壳强度的影响为了探究不同的地壳强度对岩浆聚集的影响，我们引入了粘度因子方法调节地壳强度，在参考模型的基础上，系统测试了 、，、时的地壳强度对岩浆上升过程的影响（；见表 ）。在这里，我们展示 、三个端元模型的演化结果，其他模拟结果见电子版附图 。图 为 、（）的模拟结果，在地壳强流变的情况下，岩浆主体在岩石圈底部聚集，仅有小体积的岩浆在下地壳底部侵位。当固定上地壳强度、仅变弱下地壳强度时（，图 ），岩浆在莫霍面之上形成了扁平状岩浆房，外部为老期次岩浆内部为新期次岩浆，上、下地壳界面发崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟图 地壳地温梯度对地形的影响（）（）分别为冷地壳、温地壳、热地壳背景下地形随时间的演化（）实线与虚线分别为左翼隆起的最高高程与中心拗陷处的地表高程随时间的变化（）为 时地形的垂直剖面（）为地表的抬升速率 颜色标记与图 相同 （）（），（），（）（）生宽缓的隆升。当固定下地壳强度、仅变弱上地壳强度时（，图 ），岩浆倾向于向上传播，破裂前的下地壳发育球状岩浆房，岩浆由岩墙穿过下地壳后（见图 ），高超压驱动下的岩浆在上地壳自发上升至地表，上、下地壳界面处的韧性区仅有少量岩浆聚集，地表喷发主要为新期次的岩浆，包裹极少量的下地壳岩墙、上地壳岩石等物质。综上，我们基于不同的上、下地壳强度背景下模拟岩浆上涌的过程，表明上、下地壳的相对强度对岩浆的聚集样式具有重大影响，地壳强度对岩浆房的横向宽度影响更大，地壳越弱，岩浆房的宽度越宽。岩石学报 ，（）：图 地壳强度对岩浆上升的影响（）（）为 、（），、（），、（）三个模型的模拟结果（）（）为相应的粘度结果 岩石色标同图 ，粘度色标同图 （）（），（）；，（）；，（）（）（），讨论 岩浆聚集样式的控制因素根据模型测试，我们发现有两种因素会影响壳内岩浆的聚集样式。首先是地壳的地温梯度，地壳越冷，岩浆无法在壳内聚集，地壳越热，岩浆越易在莫霍面之上聚集；其次是上、下地壳的相对强度，岩浆主体在流变较弱的深度聚集。在模拟过程中，我们设定来自地幔深部的岩浆热源提供六期岩浆脉冲，添加到岩石圈内部的岩浆体积是相同的，排除了由于岩浆体积不同带来的影响，因此岩浆的聚集样式主要取决于地壳的强度和热状态。对于岩浆聚集样式的控制因素，（）通过模拟单期幔源岩浆在均质下地壳聚集的过程，强调围岩内聚力和内摩擦系数在控制侵入体的高度和形态上起着重要作用，韧性的下地壳有利于岩浆的横向扩张，塑性下地壳使得岩浆沿断层上涌至地表。（）通过模拟镁铁质岩浆注入预设长英质岩浆房的过程，发现地壳的塑性强度决定了岩浆上升的高度，岩浆侵位低强度的地壳仅会停留在预设的下地壳岩浆房，岩浆侵位中等强度的地壳可以沿断层上升至上地壳，岩浆侵位高强度的地壳仅会上升至上、下地壳界面。（）利用三维热 力学模拟的方法探究地壳温度和流变对岩浆侵位形态的影响，发现地壳的温度越高，岩浆更易在莫霍面聚集形成岩浆房，地壳的塑性强度越小，岩浆侵位的深度越浅。但他们的结果均没有形成分层聚集的岩浆房。我们的模拟结果与前人的研究结果基本一致，但我们通过引入粘度因子来调节地壳的粘性流变性质，发现地壳粘性流变分层性在控制岩浆多层级的聚集过程中具有重要影响。以参考模型的地温梯度为基础，我们系统地测试了不同上地壳、下地壳强度时岩浆的上升过程，并将模拟结果总结在图 中。结果显示，当上地壳与下地壳均具有较强的流变时（、），幔源岩浆主体在岩石圈底部聚集；当下地壳的强度较弱（、）时，在莫霍界面处可以形成下地壳岩浆房；进一步，当上地壳的强度较弱时（、），在上、下地壳界面和莫霍面可以形成双层聚集的上地壳岩浆房和下地壳岩浆房。岩浆在岩石圈、莫霍面、上、下地壳界面等不同深度聚集与上、下地壳的强度呈现了一定的趋势，强度越弱岩浆房越扁平，上地壳岩浆房更难形成，并且上地壳岩浆房的厚度总是小于下地壳岩浆房的厚度。对长白山火山活动的启示长白山火山历史上曾发生多次剧烈喷发（，崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟图 上、下地壳强度对岩浆聚集样式的影响在这些模式中，上地壳与下地壳的粘度因子是不同的，观测到岩浆聚集的三种模式：仅在岩石圈聚集、仅在下地壳聚集与双层聚集 ：，），目前该火山仍存在着较高的潜在喷发危险（刘嘉麒等，；，）。前人运用了不同的地震学探测方法探究了长白山火山区下方的地壳结构，结果显示，该火山下方的地壳内存在多层聚集的岩浆房，在下地壳表现为球状岩浆房、在上地壳表现为扁平状岩浆房（张先康等，；，；，；陈棋福等，；，）。大地电磁结果显示该地区火山下的地壳存在显著的电导率不均一的特征，对应为岩石成分或熔流体含量的不均一（汤吉等，；仇根根等，；，），预示着壳内存在多层级的流变性质分界面。来自地幔深部的低波速、高温物质为长白山地区的岩浆活动持续提供热源（，；，；，），维持了长白山火山现今较热的地温状态（平均热流值大于 ，）。我们的模拟结果表明，当地壳围岩流变分层性越显著、地温梯度越高时，更容易形成多层级聚集的岩浆房，长白山壳内结构特征与附图 所示的结果具有很好的一致性。因此，我们推测长白山地区壳内岩浆的聚集样式主要受围岩地温梯度和地壳内强烈的流变分层性的影响。幔源岩浆侵位过程中地壳热状态对岩浆活动区地形演化的影响在这里，我们通过设定不同的背景地壳地温梯度模拟了幔源岩浆的上升过程，进一步分析地表的地形演化结果发现：在不同的背景地壳地温梯度下，岩浆活动区的地形演化具有明显的分类特征，冷地壳背景下岩浆底侵对壳内及地表基本没有影响，而温地壳和热地壳背景下，由于壳内岩浆多期次聚集形成岩浆房，岩浆房之上的地表地形出现了多期次的抬升。这种地壳热状态与地表高程之间的对应关系与许多地质实例相吻合。例如，我国大同火山的背景热流值约为 （，），属于中低热流区，对应为我们的模拟中的冷地壳背景，该地区的地形相对高差较小（电子版附图 ），与我们针对冷地壳的地形模拟结果具有高度的一致性。长白山火山区的背景热流值约为 （，；附图 ），对应为我们的模拟中的温地壳背景，相对于周围的地形，该火山区地形的相对高差约 （附图 ）；黄石超级火山是美国西部的一个地幔柱成因火山（，），该地区的平均热流值高于 （，），为典型的热地壳，黄石超级火山区作为岩浆底侵区域（，），相对于斯内克河平原，相对高差约 （，附图 ）。这两个火山区均由于岩浆底侵导致了地表抬升，与我们针对温地壳和热地壳的地形模拟结果高度吻合。因此，在幔源岩浆侵位的过程中，地壳的热状态与地表的地形演化之间具有极大的相关关系，结合我们对不同地壳地温梯度背景下岩浆侵位过程中岩浆活动区地形的演化结果，可以得出：岩浆活动区的地形演化能够反映岩浆侵位时地壳的热状态，地壳越冷，地形的相对变化越小，地壳越热，壳内聚集的岩浆房会在地表形成相应的地形抬升，抬升幅度与岩浆房的高度正相关。模型局限性我们采用高分辨率的二维热 力学模型模拟了幔源岩浆上升的精细动力学过程。但是，我们的模型也存在一些局限性，首先是岩浆的粘度，由于岩浆熔流体含量的复杂性，我们考虑部分熔融岩石的流变时对粘度做了一定的简化，这与实际地质过程中熔融岩石的粘度范围存在差别（，），熔融岩石的粘度与温度和压力有关，模拟过程中综合考虑以上因素会对模拟结果带来影响，但是不会影响我们对壳内岩浆的聚集样式的探究。其次是岩墙的宽度，实际地质过程中岩墙的宽度为厘米或者米级尺度（，），受计算条件的限制，我们模拟的岩墙宽度超过了实际的分辨率。虽然我们的模型不能解释自然界所有的复杂性，但是我们可以在极短时间内模拟更为接近实际的过程，并给定了岩浆上升的一般性模型。未来我们将添加两相流的模拟（，；，），考虑挥发分出溶带来的影响（，）。结论我们通过二维热 力学方法模拟了幔源岩浆上涌的动力学过程，探究了地壳地温梯度和地壳强度对岩浆上升及其聚集过程的影响，主要获得以下结论：（）地壳地温梯度对岩浆的侵位深度有重要影响。岩浆侵位冷地壳，岩浆主体在岩石圈深度聚集；岩浆侵位温地壳 岩石学报 ，（）：形成下地壳岩浆房，并且岩浆可以侵位到上、下地壳界面；岩浆侵位热地壳仅形成下地壳席状岩浆房。同时，伴随岩浆脉冲的多期次侵位，地表呈现多期次快速抬升的特征。在温地壳与热地壳条件下，岩浆房上方的地表地形呈现中心凹陷两侧隆起的形态；在冷地壳条件下，地表基本无起伏。（）地壳流变性质的分层性决定了岩浆侵位的样式。同一地壳地温梯度下，强的上地壳与强的下地壳（、）使得岩浆主体在岩石圈聚集，无法上升至地壳；当上地壳的强度强、下地壳的强度弱时（、），岩浆在下地壳聚集形成岩浆房；当上地壳的强度弱、下地壳的强度强时（、），岩浆倾向于在地表喷发；在较弱的上地壳与较弱的下地壳的情况下（、），岩浆在壳内更容易形成双层聚集的岩浆房。此外，我们还发现在壳内岩浆多层聚集过程中，周期性岩浆脉冲的补给会导致岩浆房内部超压值呈现瞬时骤增随后下降的特征。（）根据我们的分析，长白山火山壳内岩浆的分层聚集样式主要受控于地壳的热状态和流变性质分层性；火山区地形的相对变化能够反映岩浆侵位时地壳的热状态，地壳越冷，地形的相对变化越小，地壳越热，壳内聚集的岩浆房会在地表形成显著的地形抬升。我们的模拟结果对理解岩浆在壳内的聚集样式和岩浆活动区的地形演化具有启示意义。致谢 本 研 究 使 用 的 源 程 序 来 自 ，所有的图片均使用 绘制（：；，）。感谢三位审稿专家和编辑提出的建设性意见，这些意见极大地帮助我们提高了文章的质量。感谢北京大学的李扬教授、中国科学院地质与地球物理研究所的刘博达副研究员对本文前期的研究工作提供的建议，感谢课题组颜智勇、唐嘉萱、谢仁先、司润港、向宵对本文的帮助。感谢国家超级计算天津中心提供“天河一号”计算平台，感谢北京超级云计算中心（）提供的高性能计算资源。，（）：，：，（），：，：，：，（）：，（）：，（）：，（）：，：，（）：，（）：，（）：，（）：（），（），（）：，：，（）：（）：，：，（）：，（）：，（）：，（）：，（）：，（）：，（）：，（）：，（）：，崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟 ，（）：，（）：，（）：，（）：，（）：（）：，（）：，（）：，（）：，（）：，：，：（）：，：，（）：，（）：，：，：，（）：，：，（）：，（），（）：（），（）：，（）：（），：，（）：（），（）：，：，（）：，（）：，（）：，（）：（），（）：，：，（）：（）：，：，（）：，：，（）：，：，：，（）：，（）：，：，（）：，：，（）：，（）：，岩石学报 ，（）：，（）：（），（）：，（）：，：，（）：，（）：，：，：，：，：，（）：（），（）：，（）：，：，（）：（），：，（）：附中文参考文献常成，罗纲 斑岩矿床侵入体顶部破裂系统形成的力学机制：多场耦合数值模拟的启示地球物理学报，（）：陈棋福，艾印双，陈?长白山火山区深部结构探测的研究进展与展望 中国科学（地球科学），（）：卢靖雯，王勤，刘春 花岗质岩浆侵位对围岩裂隙发育和热结构影响的数值模拟 地质学报，（）：郭文峰，刘嘉麒，郭正府 长白山上新世以来玄武岩成分演变规律及其成因 岩石学报，（）：刘丹红，陈林 边界驱动对流与克拉通岩石圈减薄：二维热力学模拟 中国科学（地球科学），（）：刘嘉麒，陈双双，郭文峰，孙春青，张茂亮，郭正府 长白山火山研究进展 矿物岩石地球化学通报，（）：仇根根，裴发根，方慧，杜炳锐，张小博，张鹏辉，袁永真，何梅兴，白大为 长白山天池火山岩浆系统分析地球物理学报，（）：汤吉，邓前辉，赵国泽，李文军，宣飞，晋光文，白登海，詹艳，梁竞阁，蒲兴华，王继军，李国深，洪飞，马明志，陈风学 长白山天池火山区电性结构和岩浆系统 地震地质，（）：杨林，王庆飞，赵世宇，李华健，赵鹤森，董超一，刘学飞，邓军 造山型金矿构造控矿作用 岩石学报，（）：张先康，张成科，赵金仁，杨卓欣，李松林，张建狮，刘宝峰，成双喜，孙国伟，潘素珍 长白山天池火山区岩浆系统深部结构的深地震测深研究 地震学报，（）：崔晓娜等：幔源岩浆在地壳中分层侵位的控制因素：二维热 力学模拟