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

崔晓娜 陈林

上地幔部分熔融形成的低粘度、低密度的基性岩浆在地壳中聚集和侵位的过程通常会伴随岩浆房、岩墙或岩席的形成(Blacketal., 2021; Caricchietal., 2021)。岩浆的侵位和喷发往往会导致地震、海啸等自然灾害的发生(Carvajaletal., 2022; Journeauetal., 2022),也会形成铜矿、金矿等金属矿床(Parketal., 2021; 杨林等, 2023)。因此,理解幔源岩浆从深部上升至地表的动力学过程具有重要的科学意义。高分辨率地球物理成像结果揭示,我国东北长白山火山下方的地壳内部存在多层级的岩浆聚集区(Zhuetal., 2019; Yangetal., 2021)。岩石地球化学数据表明,地表出露的玄武岩呈现不同期次的岩浆补给(郭文峰等, 2014; Liuetal., 2021)、岩浆房多次聚集的特征(Leeetal., 2021; Yietal., 2021),且在地表喷发的火山数量远小于地下赋存岩浆的体积(Crisp, 1984)。由此可见,幔源岩浆在地壳中的侵位具有多期次岩浆供给和分层聚集的特点。然而,幔源岩浆侵位过程中产生的分层聚集样式受何种因素控制仍不清楚。

针对幔源岩浆的上涌过程,前人已开展了大量数值和物理模拟研究。早期数值模拟研究主要关注岩浆侵位对围岩带来的热影响(Annenetal., 2006; Petford and Gallagher, 2001),或是探讨单期幔源岩浆在均质地壳的聚集过程(Gerya and Burg, 2007; Schubertetal., 2013)。近期的数值模拟研究主要关注岩浆的传输通道(Gudmundsson, 2006; Rummeletal., 2020)。围岩的流变强度及其分层性决定了岩浆的传输方式,包括熔融底辟、减压通道或张性裂隙(Kelleretal., 2013),而岩墙和底辟可以共同促进岩浆向上传播(Caoetal., 2016)。伴随岩浆侵位,侵入体周围发育的裂隙由初期呈径向展布的剪裂隙转换为向上扩展的主干张性断裂(卢靖雯等, 2022),同时岩浆侵位会对围岩施加附加变形,导致近场与远场的小尺度变形(常成和罗纲, 2022)。Gorczyk and Vogt (2018)利用三维热-力学模拟方法探究地壳流变和温度结构对基性岩侵位形态的影响,最近Rummeletal.(2020)考虑了岩浆熔融抽取过程中的岩浆成分演化,模拟壳内多期次玄武岩席聚集形成岩浆房的机制。物理模拟研究早期主要探究岩墙向岩席的转换机制,通过将流体注入固体明胶的模拟实验,验证了当固体明胶存在刚度对比层、层间弱接触面或外部应力时,可以导致垂向传播的裂隙运动轨迹转为水平向(Kavanaghetal., 2006,2015; Menandetal., 2010)。最近,Williamsetal.(2022)利用物理模拟方法探究了岩浆沿岩席长距离水平传播的形成机制,发现岩浆的流变学性质具有重要影响。然而,物理模拟受相似比的限制,难以直接运用于理解实际的地质过程(Burovetal., 2014)。综上,前人的相关研究要么没有考虑岩浆的多期补给,要么没有考虑岩墙作为不同层级岩浆房之间的连接通道,尤其是针对多期幔源岩浆在地壳中分层侵位过程及其控制因素的研究仍然很薄弱。

在本研究中,我们聚焦幔源岩浆从深部到地壳的聚集运移过程,通过在二维热-力学数值模拟程序中引入多期岩浆脉冲和岩墙生成算法,模拟幔源岩浆从深部地幔上升到地壳的动力学过程,系统测试地壳地温梯度和流变强度的影响,探究幔源岩浆在地壳内部分层侵位的控制因素。我们将首先介绍数值模拟方法与模型设计;然后描述地壳地温梯度和地壳强度对壳内岩浆赋存状态的影响;最后,讨论幔源岩浆在地壳内聚集样式的控制因素。

1 数值模拟方法

我们使用开源的二维热-力学MATLAB程序(Gerya, 2019)模拟幔源岩浆多期次侵位的动力学过程。本程序在求解动量守恒方程、质量守恒方程与能量守恒方程时,联合交错网格有限差分方法和粒子标记(marker-in-cell, MIC)技术。在节点上解方程,在marker上传递物性参数,可以有效避免利用有限差分方法求解连续性方程时可能引起的数值发散问题。同时,该方法考虑了岩石的粘-弹-塑性流变本构关系和岩石的部分熔融,可以在时间与空间尺度上精细模拟岩浆侵位的动力学过程。

1.1 控制方程

(i)质量守恒-连续性方程

(1)

式中,Vx和Vy分别代表速度矢量的水平和垂直分量。

(ii)动量守恒-斯托克斯方程

(2)

(3)

(4)

(5)

(iii)能量守恒-热传导方程

(6)

1.2 岩石流变

(7)

(8)

(9)

其中,式中,R是气体常数;AD、Ea、Va和n是岩石物理实验约束的流变参数,分别代表物质常数、活化能、活化体积和应力指数;f为粘度因子,用于线性改变物质的有效粘度(Beaumontetal., 2006)。由于实验室条件下物质成分的变化,所测得的流变参数也具有很大的不确定性。为了将这些结果应用到实际地质过程,我们使用粘度因子对物质在实验室条件下测得的有效粘度进行线性调节(刘丹红和陈林,2020)。在这里,我们通过改变粘度因子调控地壳的强度。

熔融岩石(5%

表1 数值模型中用到的主要物理参数

物质的塑性变形遵从Drucker-Prager屈服准则(Ranalli,1995):

σyield=C0+Psin(φ)

(10)

其中C0是岩石的内聚力;φ是有效内摩擦角。

1.3 部分熔融

根据岩石物理实验数据的约束,程序考虑了岩石的部分熔融行为(Gerya and Yuen, 2003a; Burg and Gerya, 2005)。假设部分熔融百分数M与温度为简单的线性关系:

(11)

式中,Tsolidus为岩石的固相线,Tliquidus为岩石的液相线,各种岩石的固相线和液相线见表1。

部分熔融岩石的有效密度ρeff取决于部分熔融百分数M:

ρeff=ρsolid-M(ρsolid-ρmolten)

(12)

式中,ρsolid和ρmolten分别代表岩石的固相密度和熔融密度,它们都是关于温度和压力的函数,计算表达式见表1。其中,我们对幔源岩浆密度的选取参照了Sanloup (2016)的结果,并且与该作者汇编的结果进行了对比,在0～5GPa的压力范围内,我们的结果中岩浆的密度与前人的实验结果和理论计算结果均一致(电子版附图1),说明岩浆密度的选取具有合理性。

图1 岩墙形成的示意图

部分熔融对有效热容量(Cpe)和热扩散系数(αe)的影响由下式给出(Burg and Gerya, 2005):

(13)

1.4 多期岩浆脉冲的实现

火山喷发需要体积足够大、低粘度、低密度、部分熔融且持续供应的岩浆热源,以提供浮力、维持高温(Caricchietal., 2021)。我们通过预设岩石圈之下的岩浆房来模拟幔源岩浆上升的动力学过程,在已有的二维粘-弹-塑性数值模拟程序(Gerya and Burg, 2007)的基础上加入多期岩浆脉冲算法,在岩石圈以下多期次、定量的稳定提供岩浆热源,尽量保证岩浆通道不变形。算法实现细节为:(1)我们添加新一期岩浆脉冲的判断条件为:深度120km以下,部分岩浆通道和岩浆源区中熔体的总含量小于5%;(2)当熔体含量满足上述条件时,我们用性质相同但标号不同的marker补充新一期次的幔源岩浆;(3)在模拟过程中,我们总共提供六期岩浆脉冲,为岩浆的持续上涌提供动力。

1.5 岩墙的生成算法

岩墙作为拉伸作用形成的薄弱带,它的形成与围岩的应力场有关(Gudmundsson, 2006; Bertelsenetal., 2018)。在本文中,拉伸应力标记为正、挤压应力标记为负。岩墙作为岩浆房之间的连接通道(Cox, 1980; Reuberetal., 2018),岩浆沿岩墙传播(Becerriletal., 2013; Scuderoetal., 2019),其宽度随深度变化(Geshietal., 2010, 2020)。目前对于岩墙的热-力学模拟主要探究岩墙对岩浆上升过程的影响,例如,Caoetal.(2016)认为花岗岩岩浆房会形成辐射状岩墙,为了探究岩浆上升过程中岩墙和底辟的相互作用,因此设定临界的应力和第二应力不变量作为垂直状岩墙形成的条件。Rummeletal.(2020)通过岩石热-力学模拟方法,考虑岩墙作为玄武岩熔融提取的通道。我们通过探究幔源岩浆侵位对地壳应力状态带来的影响,认为垂直状岩墙是岩浆上升的通道,辐射状岩墙是岩浆房体积扩张时围岩的同期破裂,表征岩浆房的发育。因此,本文中考虑岩墙是在岩浆侵位对上覆围岩造成拉伸环境中形成的(Gudmundsson, 1988),作为岩浆向浅部运移的通道。

我们在程序中加入岩墙生成算法,岩墙的形成和扩展由岩浆房超压值和地壳应力共同控制。岩浆房超压值为岩浆房内部的总压强与静岩压力的差值(岩浆房内部的超压值相差无几),岩浆房内部的超压值会随着每一期次的岩浆补给而增加、随着岩浆房的体积增大而减小,超压值演化表征了岩浆房的发育过程;当岩浆底侵时,岩浆聚集使得岩浆房之上产生应力集中的区域,即岩浆房的发育对上覆的上地壳与下地壳形成了水平的拉伸应力,这种拉伸应力即为地壳应力,地壳的应力场决定了岩墙垂直向发育。以下地壳岩墙的形成过程为例,具体算法如下:(1)首先,由于多期次的岩浆脉冲供应,壳内岩浆房发育至一定体积,将岩浆房内部超压值小于120MPa作为岩浆房发育结束的时刻与判断是否形成岩墙的初始时刻。(2)随后,当岩浆房之上的地壳应力大于30MPa时,发育长度为岩浆房顶部至上、下地壳界面、宽度为2km的岩墙(见图1)。经过多次测试,2km宽度的岩墙是保证岩浆上升的最窄宽度,过窄的岩墙会导致通道快速冷却、无法输运岩浆,因此我们在模拟中设定岩墙的最小宽度为2km。

2 模型设计

为了模拟幔源岩浆上升的动力学过程,我们设计了如图2所示的初始模型。模型在水平方向上的宽度为160km,垂向深度为180km,包含20km上地壳、20km下地壳和90km岩石圈地幔,代表了一个典型的大陆岩石圈,具体物性参数见表1。我们初始在岩石圈底部预设一个半径为10km、温度为1700K的圆形幔源岩浆区域,在顶部设置了一个宽度为3km的岩浆通道,连接到地壳底部。整个模型空间被321×361个节点组成的规则网格离散化,空间分辨率为500m,总共120万个粒子(markers)随机分布在网格空间中,粒子分辨率高达150m。

图2 模型设计

模型的所有力学边界均为自由滑边界。在模型的上部,我们设计一个厚20km的伪空气层,其密度和粘度分别为1kg/m3、1013Pa·s,用于模拟岩浆上升过程中产生的地表起伏(Gerya and Burg, 2007; Schmelingetal., 2008)。模型的热边界条件:顶边界固定为273K,侧边界为绝热边界,底边界为固定为1600K (Burg and Gerya, 2005; Huetal., 2018),整个空气层温度固定为273K,岩石圈底边界以下为绝热地温梯度0.5K/km。为了区分地壳的不同热状态,我们采用分段函数赋予了地壳两种初始温度结构:第一种温度结构初始固定莫霍面的温度,设定地表和莫霍底界面的初始温度分别为273K和TMoho=900K,通过改变上、下地壳界面的初始温度(TUC=586K、686K、786K)来调节上、下地壳的地温梯度(GUC、GLC),让温度随深度线性增加;第二种温度结构初始固定岩石圈底部的温度,设定地表和岩石圈底界面的初始温度分别为273K和1573K,通过改变莫霍面的初始温度(TMoho=673K、773K、873K、973K、1073K)来调节地壳和岩石圈的地温梯度,同样让温度随深度线性增加。

3 模拟结果

为了探究来自地幔深部的多期岩浆在地壳中分层聚集的动力学过程,我们基于多期岩浆脉冲和岩墙算法的MATLAB程序,开展了27组数值模拟实验(表2)。在模拟中,系统测试了地壳地温梯度和地壳强度对岩浆上涌过程的影响。在这里,我们选择上地壳的粘度因子fUC=0.05、下地壳的粘度因子fLC=0.2、TUC=686K、TMoho=900K的模型作为参考模型(即表2中的model15)。

表2 数值模拟实验

3.1 参考模型

图3为参考模型(model15,fUC=0.05,fLC=0.2,GUC=20K/km,GLC=10K/km)的演化结果。来自岩石圈底部的幔源岩浆穿过岩石圈到达莫霍面,初始阶段(4.8kyr;图3ai),岩浆在莫霍面之上聚集形成椭圆状岩浆房,伴随同时的垂向增厚和横向扩展,粘度结果显示岩浆房周围出现局部的低粘度区域(图3aii),表征岩浆房的扩张趋势。随后,由于岩浆侵位对上覆地壳造成的拉伸力,下地壳形成岩浆通道使得岩浆由此上升(图3bi),水平拉伸应力瞬时消散(图3biii)。8.6kyr时,岩浆上涌至上、下地壳界面,多期次岩浆脉冲带来的高超压导致岩浆通道不断变宽,在岩浆房内部,二期岩浆及多期岩浆以韧性底辟形式对流, 促进多期岩浆混合及岩浆房上升(图3ci),粘度结果显示岩浆房之上局部的低粘度区域(图3cii),表征岩浆的上升趋势。最终13.1kyr时,上、下地壳界面之上岩浆聚集形成扁平状岩浆房,但岩浆房顶部的应力累积不足(图3ciii),达到稳态至形成双层聚集的岩浆房,此时下地壳岩浆房体积扩张,并且内部卷入极少量的地壳和岩墙等物质(图3di),上地壳与下地壳减薄程度不同的原因在于两者对于岩浆侵位带来的空间问题响应不同。在13.1kyr的时间尺度内,岩浆房内部保持活跃,进行持续的岩浆对流。

图3 参考模型

图3中aiii～diii为参考模型的水平偏应力随时间的变化,下地壳与上地壳岩浆侵位均对岩浆房之上的局部区域形成水平拉伸应力。上地壳累积的拉伸应力在脆韧性转换深度处最大,向地表逐渐降低,向深部为韧性区应力累积很弱;下地壳累积的拉伸应力在上、下地壳界面处最大,向深部减小,同样韧性区应力累积很弱,并且下地壳岩浆侵位导致上、下地壳界面隆起,可以被上地壳韧性区容纳(图3bii)。同时,我们选取参考模型(model15)中三个固定点的水平偏应力结果来描述岩浆上涌对下地壳、上地壳、地表带来的影响,在图3aiii中用红色的圆标记。下地壳岩浆房之上红色固定点拉伸力的初始波峰表示下地壳岩浆聚集的初始阶段,第二个40MPa小波峰表示4.5kyr添加的第二期岩浆脉冲,4.8kyr时拉伸应力由38MPa骤降表征岩浆房破裂,0MPa表征该深度已经转变为岩墙或岩浆成分(见图1)。与之相似,上地壳脆韧性转换深度处,应力的第一个波峰表示下地壳岩浆聚集过程中上地壳的同期应力累积,第二个34MPa小波峰表示添加的第二期岩浆脉冲,第三个41MPa的波峰表示上地壳岩浆的聚集过程,随后应力逐渐消散。岩浆初始聚集时地表应力的累积滞后,随后在岩浆开始沿岩墙向上传播至源区六期岩浆脉冲结束的过程中,浅表一直处于拉伸环境,当岩浆脉冲停止后,应力出现负值表示地壳浅部转为挤压状态。

我们选择下地壳岩浆房底部与上地壳岩浆房底部的两个固定点从超压演化分析岩浆房的发育过程,在图3di中用白色的圆标记。下地壳岩浆房内的超压因一期岩浆初始聚集而骤增(图4a),后因岩浆房的不断发育,超压值随着岩浆房体积增大而持续减小,后五期的岩浆脉冲均可引起超压值骤增,使岩浆房维持高超压状态而体积不断扩张,4.8kyr时超压值骤降表征岩浆房的破裂(见图1),此时为下地壳岩浆房发育结束的时刻与形成下地壳岩墙的时刻,最终下地壳岩浆房内部的超压值趋于稳定,当岩浆脉冲停止后超压负值表征岩浆上升停滞。上地壳岩浆房内部的超压初始值对应下地壳岩浆房破裂的时刻(图4b),超压值由负转正表征上、下地壳界面深度的初始岩浆聚集,后四期岩浆脉冲引起的超压最大值小于30MPa,当岩浆脉冲停止后超压负值亦表征岩浆上升停滞。

图4 岩浆房内超压值随时间的演化

3.2 地壳地温梯度的影响

在参考模型的基础上,我们测试了地壳地温梯度的影响。如图5d,设定莫霍面的初始温度分别为TMoho=783K、873K、973K (model23、model24、model25;见表2),根据单一变量测试原则,其余参数设置与参考模型相同。首先,我们对比13.1kyr时冷地壳、温地壳与热地壳背景下岩浆的赋存状态。图5ai中,冷地壳的地温梯度为GUC=GLC=12.5K/km,低的地温梯度使得地壳的粘度过大(5aii),即使岩浆高超压也无法克服地壳的强流变实现岩浆聚集,此时岩浆主要在岩石圈底部聚集。图5bi中,温地壳的地温梯度为GUC=GLC=15K/km,在地壳具有较高的地温梯度时,岩浆在莫霍面之上聚集形成椭圆状岩浆房,并且岩浆可以上升至上、下地壳界面。图5ci中,热地壳的地温梯度为GUC=GLC=17.5K/km,高的地温梯度使得地壳的粘度降低,此时岩浆仅在莫霍面之上聚集形成席状岩浆房,产生以上现象是基于岩浆的正浮力与相应地壳地温梯度下的地壳强度共同作用的结果。

图5 地壳地温梯度对岩浆侵位深度的影响

我们进一步分析了岩浆侵位对地形带来的影响,图6a-c分别为冷地壳、温地壳、热地壳背景下地形随时间的变化,温地壳与热地壳的地形演化呈现相同的趋势:侵位均导致地表抬升,初始中心隆起后转变为中心凹陷两侧隆起形态(图6b, c),而冷地壳时,地表的地形无明显起伏(图6a)。速度骤增均对应多期岩浆脉冲的岩浆供应时刻(图6f)。图6e为 13.1kyr时的地形垂直剖面,地表最大隆起值发生在温地壳下,地表最高高程约为3km,热地壳次之,地表最高高程约为1.5km,冷地壳时,地表的相对高差小于140m。

图6 地壳地温梯度对地形的影响

综上,地壳地温梯度对岩浆的侵位深度具有重大影响,温地壳有利于岩浆上升到更浅的深度,地壳的地温梯度过高或过低都不利于岩浆的上升。并且,在岩浆多期次侵位的过程中,地形多期次抬升,壳内岩浆房的高度对地表的隆起高度贡献更大。

3.3 地壳强度的影响

为了探究不同的地壳强度对岩浆聚集的影响,我们引入了粘度因子方法调节地壳强度,在参考模型的基础上,系统测试了fUC=1、0.5、0.1、0.05、0.01,fLC=1、0.5、0.2、0.1时的地壳强度对岩浆上升过程的影响(model2～21;见表2)。在这里,我们展示model2、model17、model6三个端元模型的演化结果,其他模拟结果见电子版附图2。图7ai为fUC=1、fLC=1(model2)的模拟结果,在地壳强流变的情况下,岩浆主体在岩石圈底部聚集,仅有小体积的岩浆在下地壳底部侵位。当固定上地壳强度、仅变弱下地壳强度时(fUC=1,fLC=0.1,model17,图7bi),岩浆在莫霍面之上形成了扁平状岩浆房,外部为老期次岩浆内部为新期次岩浆,上、下地壳界面发生宽缓的隆升。当固定下地壳强度、仅变弱上地壳强度时(fUC=0.01,fLC=1,model6,图7ci),岩浆倾向于向上传播,破裂前的下地壳发育球状岩浆房,岩浆由岩墙穿过下地壳后(见图1),高超压驱动下的岩浆在上地壳自发上升至地表,上、下地壳界面处的韧性区仅有少量岩浆聚集,地表喷发主要为新期次的岩浆,包裹极少量的下地壳岩墙、上地壳岩石等物质。

图7 地壳强度对岩浆上升的影响

综上,我们基于不同的上、下地壳强度背景下模拟岩浆上涌的过程,表明上、下地壳的相对强度对岩浆的聚集样式具有重大影响,地壳强度对岩浆房的横向宽度影响更大,地壳越弱,岩浆房的宽度越宽。

4 讨论

4.1 岩浆聚集样式的控制因素

根据模型测试,我们发现有两种因素会影响壳内岩浆的聚集样式。首先是地壳的地温梯度,地壳越冷,岩浆无法在壳内聚集,地壳越热,岩浆越易在莫霍面之上聚集;其次是上、下地壳的相对强度,岩浆主体在流变较弱的深度聚集。在模拟过程中,我们设定来自地幔深部的岩浆热源提供六期岩浆脉冲,添加到岩石圈内部的岩浆体积是相同的,排除了由于岩浆体积不同带来的影响,因此岩浆的聚集样式主要取决于地壳的强度和热状态。

对于岩浆聚集样式的控制因素,Gerya and Burg (2007)通过模拟单期幔源岩浆在均质下地壳聚集的过程,强调围岩内聚力和内摩擦系数在控制侵入体的高度和形态上起着重要作用,韧性的下地壳有利于岩浆的横向扩张,塑性下地壳使得岩浆沿断层上涌至地表。Schubertetal.(2013)通过模拟镁铁质岩浆注入预设长英质岩浆房的过程,发现地壳的塑性强度决定了岩浆上升的高度,岩浆侵位低强度的地壳仅会停留在预设的下地壳岩浆房,岩浆侵位中等强度的地壳可以沿断层上升至上地壳,岩浆侵位高强度的地壳仅会上升至上、下地壳界面。Gorczyk and Vogt (2018)利用三维热-力学模拟的方法探究地壳温度和流变对岩浆侵位形态的影响,发现地壳的温度越高,岩浆更易在莫霍面聚集形成岩浆房,地壳的塑性强度越小,岩浆侵位的深度越浅。但他们的结果均没有形成分层聚集的岩浆房。我们的模拟结果与前人的研究结果基本一致,但我们通过引入粘度因子来调节地壳的粘性流变性质,发现地壳粘性流变分层性在控制岩浆多层级的聚集过程中具有重要影响。

以参考模型的地温梯度为基础,我们系统地测试了不同上地壳、下地壳强度时岩浆的上升过程,并将模拟结果总结在图8中。结果显示,当上地壳与下地壳均具有较强的流变时(fUC≥1、fLC≥1),幔源岩浆主体在岩石圈底部聚集;当下地壳的强度较弱(fUC<1、fLC<1)时,在莫霍界面处可以形成下地壳岩浆房;进一步,当上地壳的强度较弱时(fUC≤0.05、fLC≤0.5),在上、下地壳界面和莫霍面可以形成双层聚集的上地壳岩浆房和下地壳岩浆房。岩浆在岩石圈、莫霍面、上、下地壳界面等不同深度聚集与上、下地壳的强度呈现了一定的趋势,强度越弱岩浆房越扁平,上地壳岩浆房更难形成,并且上地壳岩浆房的厚度总是小于下地壳岩浆房的厚度。

图8 上、下地壳强度对岩浆聚集样式的影响在这些模式中,上地壳与下地壳的粘度因子是不同的,观测到岩浆聚集的三种模式:仅在岩石圈聚集、仅在下地壳聚集与双层聚集

4.2 对长白山火山活动的启示

长白山火山历史上曾发生多次剧烈喷发(Weietal., 2013),目前该火山仍存在着较高的潜在喷发危险(刘嘉麒等, 2015; Zhangetal., 2018)。前人运用了不同的地震学探测方法探究了长白山火山区下方的地壳结构,结果显示,该火山下方的地壳内存在多层聚集的岩浆房,在下地壳表现为球状岩浆房、在上地壳表现为扁平状岩浆房(张先康等, 2002; Songetal., 2007; Choietal., 2013; 陈棋福等,2019; Fanetal., 2022)。大地电磁结果显示该地区火山下的地壳存在显著的电导率不均一的特征,对应为岩石成分或熔流体含量的不均一(汤吉等, 2001; 仇根根等, 2014; Yangetal., 2021),预示着壳内存在多层级的流变性质分界面。来自地幔深部的低波速、高温物质为长白山地区的岩浆活动持续提供热源(Lei and Zhao, 2005; Tangetal., 2014; Lietal., 2020),维持了长白山火山现今较热的地温状态(平均热流值大于70mW/m2,Jiangetal., 2019)。我们的模拟结果表明,当地壳围岩流变分层性越显著、地温梯度越高时,更容易形成多层级聚集的岩浆房,长白山壳内结构特征与附图2biv所示的结果具有很好的一致性。因此,我们推测长白山地区壳内岩浆的聚集样式主要受围岩地温梯度和地壳内强烈的流变分层性的影响。

4.3 幔源岩浆侵位过程中地壳热状态对岩浆活动区地形演化的影响

在这里,我们通过设定不同的背景地壳地温梯度模拟了幔源岩浆的上升过程,进一步分析地表的地形演化结果发现:在不同的背景地壳地温梯度下,岩浆活动区的地形演化具有明显的分类特征,冷地壳背景下岩浆底侵对壳内及地表基本没有影响,而温地壳和热地壳背景下,由于壳内岩浆多期次聚集形成岩浆房,岩浆房之上的地表地形出现了多期次的抬升。这种地壳热状态与地表高程之间的对应关系与许多地质实例相吻合。例如,我国大同火山的背景热流值约为50～60mW/m2(Jiangetal., 2019),属于中低热流区,对应为我们的模拟中的冷地壳背景,该地区的地形相对高差较小(电子版附图3a),与我们针对冷地壳的地形模拟结果具有高度的一致性。长白山火山区的背景热流值约为70～80mW/m2(Lucazeau, 2019;附图3b),对应为我们的模拟中的温地壳背景,相对于周围的地形,该火山区地形的相对高差约2km(附图2);黄石超级火山是美国西部的一个地幔柱成因火山(Huangetal., 2015),该地区的平均热流值高于100mW/m2(Lucazeau, 2019),为典型的热地壳,黄石超级火山区作为岩浆底侵区域(Colónetal., 2018),相对于斯内克河平原,相对高差约1.5km(Peng and Humphreys, 1998, 附图3c)。这两个火山区均由于岩浆底侵导致了地表抬升,与我们针对温地壳和热地壳的地形模拟结果高度吻合。因此,在幔源岩浆侵位的过程中,地壳的热状态与地表的地形演化之间具有极大的相关关系,结合我们对不同地壳地温梯度背景下岩浆侵位过程中岩浆活动区地形的演化结果,可以得出:岩浆活动区的地形演化能够反映岩浆侵位时地壳的热状态,地壳越冷,地形的相对变化越小,地壳越热,壳内聚集的岩浆房会在地表形成相应的地形抬升,抬升幅度与岩浆房的高度正相关。

4.4 模型局限性

我们采用高分辨率的二维热-力学模型模拟了幔源岩浆上升的精细动力学过程。但是,我们的模型也存在一些局限性,首先是岩浆的粘度,由于岩浆熔流体含量的复杂性,我们考虑部分熔融岩石的流变时对粘度做了一定的简化,这与实际地质过程中熔融岩石的粘度范围存在差别(Dingwell, 2006),熔融岩石的粘度与温度和压力有关,模拟过程中综合考虑以上因素会对模拟结果带来影响,但是不会影响我们对壳内岩浆的聚集样式的探究。其次是岩墙的宽度,实际地质过程中岩墙的宽度为厘米或者米级尺度(Geshietal., 2010, 2020),受计算条件的限制,我们模拟的岩墙宽度超过了实际的分辨率。虽然我们的模型不能解释自然界所有的复杂性,但是我们可以在极短时间内模拟更为接近实际的过程,并给定了岩浆上升的一般性模型。未来我们将添加两相流的模拟(Kelleretal., 2013; Rässetal., 2019),考虑挥发分出溶带来的影响(Tian and Buck, 2022)。

5 结论

我们通过二维热-力学方法模拟了幔源岩浆上涌的动力学过程,探究了地壳地温梯度和地壳强度对岩浆上升及其聚集过程的影响,主要获得以下结论:

(1)地壳地温梯度对岩浆的侵位深度有重要影响。岩浆侵位冷地壳,岩浆主体在岩石圈深度聚集;岩浆侵位温地壳形成下地壳岩浆房,并且岩浆可以侵位到上、下地壳界面;岩浆侵位热地壳仅形成下地壳席状岩浆房。同时,伴随岩浆脉冲的多期次侵位,地表呈现多期次快速抬升的特征。在温地壳与热地壳条件下,岩浆房上方的地表地形呈现中心凹陷两侧隆起的形态;在冷地壳条件下,地表基本无起伏。

(2)地壳流变性质的分层性决定了岩浆侵位的样式。同一地壳地温梯度下,强的上地壳与强的下地壳(fUC≥1、fLC≥1)使得岩浆主体在岩石圈聚集,无法上升至地壳;当上地壳的强度强、下地壳的强度弱时(fUC≥1、fLC<1),岩浆在下地壳聚集形成岩浆房;当上地壳的强度弱、下地壳的强度强时(fUC≤0.05、fLC≥1),岩浆倾向于在地表喷发;在较弱的上地壳与较弱的下地壳的情况下(fUC≤0.05、fLC≤0.5),岩浆在壳内更容易形成双层聚集的岩浆房。此外,我们还发现在壳内岩浆多层聚集过程中,周期性岩浆脉冲的补给会导致岩浆房内部超压值呈现瞬时骤增随后下降的特征。

(3)根据我们的分析,长白山火山壳内岩浆的分层聚集样式主要受控于地壳的热状态和流变性质分层性;火山区地形的相对变化能够反映岩浆侵位时地壳的热状态,地壳越冷,地形的相对变化越小,地壳越热,壳内聚集的岩浆房会在地表形成显著的地形抬升。我们的模拟结果对理解岩浆在壳内的聚集样式和岩浆活动区的地形演化具有启示意义。

致谢本研究使用的MATLAB源程序来自www.cambridge.org/9781107143142,所有的图片均使用GMT绘制(https://www.soest.hawaii.edu/gmt/; Wesseletal., 2019)。感谢三位审稿专家和编辑提出的建设性意见,这些意见极大地帮助我们提高了文章的质量。感谢北京大学的李扬教授、中国科学院地质与地球物理研究所的刘博达副研究员对本文前期的研究工作提供的建议,感谢课题组颜智勇、唐嘉萱、谢仁先、司润港、向宵对本文的帮助。感谢国家超级计算天津中心提供“天河一号”计算平台,感谢北京超级云计算中心(BSCC)提供的高性能计算资源。