-
岩石圈地幔的破坏、减薄是一个复杂的深部地幔动力学过程，其破坏原因可能是由地幔柱[1-9]或俯冲作用导致岩石圈地幔的重力不稳定，从而发生拆沉、减薄、侵蚀。其中俯冲作用对岩石圈地幔的破坏可能有2种方式：第1种方式由俯冲引起的地幔对流可导致岩石圈结构破坏[10-11]，该破坏减薄过程往往需要岩石圈地幔的加水弱化，即俯冲过程中洋壳—洋幔会发生变质脱水和含水过程[12-17]，由此引起上覆岩石圈地幔的弱化和相变。除此之外，前人的数值模拟结果[3,18-20]发现，岩石圈地幔破坏减薄还可能存在另外一种作用机制，即板块深俯冲过程中形成的俯冲板片熔融柱对上覆岩石圈地幔的破坏作用。在大洋板块俯冲过程中，因温度压力的升高，俯冲板片可能会发生部分熔融形成俯冲板片熔融柱。该熔融柱能够在浮力作用下快速上涌，并与上覆板块岩石圈地幔底面发生相互作用，从而导致岩石圈的破坏减薄。目前，对于俯冲板片熔融柱的形成及其对岩石圈破坏机制、程度及过程并不清楚。
类似于地幔柱作用，俯冲板片熔融柱对岩石圈的破坏在地表也一定存在一种响应关系，这种响应关系对约束地球内部的地幔动力学过程以及岩浆活动至关重要[21]。前人曾通过地震层析成像和三维数值模拟手段揭示出上升的地幔柱与岩石圈地幔相互作用过程，并试图与地表变化联系起来，但得到的结果十分不清晰，从而引发动力地形变化与深部地幔动力学过程间耦合关系的争论[21-22]。这种争论当涉及到与俯冲过程相关的地形变化时，就变得更加复杂不清。这是因为俯冲带地形变化既包括由区域构造变形导致的地形变化，又包含了由深部地幔动力学过程导致的动力地形变化[21]，这里特别涉及到地表动力地形变化与深部破坏作用间的响应关系方面。
为了解决上述问题，本文主要采用I2VIS二维地球动力学数值模拟方法，在给定边界条件下，在物理方程和流变方程控制下，模拟与俯冲相关形成的俯冲板片熔融柱破坏岩石圈地幔结构的动力学演化过程，提取地形数据，结合地质证据深入分析其影响地表地形的变化过程。
-
基于I2VIS[16]二维地球动力学数值算法结合marker-in-cell[23-27]技术建立数值模型。主要是对3组数值模拟方程进行求解：斯托克斯流体动量守恒方程、物质守恒方程以及热量守恒方程[28]：
（1）斯托克斯方程：
$\frac{{\partial \sigma' _{{{ij}}}}}{{\partial {x_j}}} - \frac{{\partial P}}{{\partial {x_i}}}={g_i}\rho (C,M,P,T)$
(1) 式中，
$ \sigma' _{{{ij}}}$ 表示偏应力张量，g表示重力加速度，其中密度ρ与岩性成分（C）、部分熔融体积比例参数（M）、压强（P）和温度（T）相关。（2）质量守恒方程可以近似看做是不可压缩性流体方程：
$\frac{{\partial {v_i}}}{{\partial {x_i}}}=0$
(2) 式中，v是速度分量。
（3）热守恒方程
$\begin{aligned} \rho {C_{\rm{p}}}\frac{{DT}}{{Dt}}=- \frac{{\partial {q_i}}}{{\partial {x_i}}} + H \\ {q_i}={\rm{ - }}k(T,P,C)\frac{{\partial T}}{{\partial {x_i}}} \end{aligned} $
(3) 式中，重复出现的索引代表每个坐标系（
$ { x},{ y},{ \textit z}$ ）计算热通量的总和；ρ表示密度（kg/m3）；Cp表示恒定压力下的热容量（J·kg−1·K−1）；H表示产生的体积热（W/m3），包含产生的放射热、剪切热、绝热以及潜热。$\dfrac{{DT}}{{D{\rm{t}}}}$ 是对应标准拉格朗日-欧拉方程下温度对时间的导数。k是热导率，与温度（T）、压力（P）、物质组成以及物质结构（C）相关[8,29]。与先前的数值模型进行了数值对比，部分类型岩石在一定的温度和压力条件下会发生部分熔融，因此，本模型中考虑到了部分岩石部分熔融的计算[30-31]。基于实验岩石学条件的限制，对温度和岩石的部分熔融体积比例做了近似处理，假设两者之间是一种线性关系。因此，地幔减压熔融也可以用数值代码进行模拟[32]，这一过程与密度和岩石组成的变化有关。岩石的部分熔融可以用下列方程进行计算：
$ \begin{aligned} &M=0,T \text{≤} {T_{{{{\text{固态}}}}}}\\ &M=\frac{{(T - {T_{{\text{固态}}}})}}{{({T_{{\text{液态}}}} - {T_{{\text{固态}}}})}},\quad{T_{{\text{固态}}}} \text{＜}T \text{＜} {T_{{\text{液态}}}}\\ &M=1,T \text{≥} {T_{{\text{液态}}}} \end{aligned} $
(4) 式中，M表示岩石部分熔融体积比例。T固态、T液态分别表示根据实验获得的湿固相线和干固相线温度。
如上文所述，物质组成、熔融体积分数、温度、压力与岩石的密度紧密相关。对于部分熔融的岩石，有效密度遵循下述状态方程：
$ \begin{aligned} {\rho _{\text{有效}}}={\rho _{{\text{固体}}}} - M({\rho _{{\text{固体}}}} - {\rho _{{\text{熔融}}}})\\ {\rho _{{\text{固体}}}}=\rho \left[ {1 - \alpha (T - {T_0})} \right]\left[ {1 + \beta (P - {P_0})} \right] \end{aligned} $
(5) 式中，ρ熔融和ρ固体表示岩石的熔融密度和固体密度。ρ0代表岩石在标况下的标准密度（表1）。热膨胀系数α为3×10−5 K−1，热压缩系数β为1×10−5 K−1。
表 1 数值模拟采用的黏滞性流变参数（据文献[35-38]）
Table 1. Parameters of viscous flow in the numerical experiments (after references[35-38])
标号 流变性质 E/(K·J·mol−1) V/(J·MPa−1·mol−1) n AD/(MPa−n·s−1) η0a/(Pa·s) A* 空气/水 0 0 1.0 1.0×10−12 1×1018 B* 湿石英（强） 154 0 2.3 3.2×10−6 1.97×1019 C* 斜长An75（强） 238 0 3.2 3.3×10−6 4.80×1024 D* 斜长石 An75 238 0 3.2 3.3×10−4 4.80×1022 E* 无水橄榄岩 532 8 3.5 2.5×104 3.98×1016 F*b 湿橄榄岩 470 8 4.0 2.0×103 5.01×1020 G*b 长英质熔体 0 0 1.0 2.0×10−9 5.00×1014 H 铁镁质熔体 0 0 1.0 1.0×10−7 1.00×1013 注：:a η0表示为有效黏滞系数，计算公式为：η0=(1/AD)×106n;
b 熔融的长英质熔体表示的是熔融的沉积物和地壳。当岩石处于塑性变形阶段时，模型将遵循莫尔-库伦屈服准则[33]：
$\begin{aligned} \sigma _{\text{屈服}}={C_0} + P\sin ({\phi _{\rm{eff}}}) \\ {\eta _{\text{摩尔-库伦}}}=\frac{{{\sigma_{\text{屈服}}}}}{{2{\varepsilon'_{{\scriptsize{\text Π}}} } }} \\ \end{aligned} $
(6) 式中，
$ {{\rm{\sigma }}_{\text{屈服}}}$ 表示屈服应力，$ {\varepsilon '_{{\scriptsize{\text Π}}}}$ 表示第二不变偏应变率张量，P表示压力，C0表示标况下的凝聚力，η表示黏度系数，${\phi _{\rm{eff}}}$ 表示有效内摩擦角。在模型中还会涉及到地形的变化，地形变化会受到侵蚀和沉积速率的影响并遵循下述状态方程[27]：
$ \begin{aligned} \frac{{\partial {{\rm{y}}_{\rm{es}}}}}{{\partial t}}=&{v_{\rm y}} - {v_{\rm x}}\frac{{\partial {y_{\rm {es}}}}}{{\partial y}} - {v_{\rm s}} + {v_{\rm e}}\\ {v_{\rm s}}=& 0\;{\rm{mm/a}},\quad{\nu _{\rm{e}}}={v_{{\rm e}0}},\;{\text{当}}\;y_{\rm{es}}\text{＜}{y_{{\text{水}}}}\\ {v_{\rm s}}=&{v_{{\rm s}0}}\;{\rm{mm/a}},\quad{v_{\rm e}}=0,\;{\text{当}}\;y_{\rm{es}}\text{＞} {y_{{\text{水}}}} \end{aligned} $
(7) 式中，
${y_{\rm{es}}}$ 表示相对水平方向距离地表垂直方向的函数，${\nu _{\rm x}}$ 和${\nu _{\rm y}}$ 分别表示在地表物质运动速度的水平分量和垂直分量。${\nu _{{\rm e}0}}$ 和${\nu _{{\rm s}0}}$ 分别是侵蚀速率常数和沉积速率常数，${\nu _{\rm e}}$ 和${\nu _{\rm s}}$ 分别表示侵蚀速率和沉积速率。 -
本文中设置的模型主要是在一定的边界条件下模拟因板片俯冲导致上覆大陆岩石圈地幔破坏以及对浅表的影响过程。初始模型的长度为4 000 km，宽度为670 km（图1）。使用699×134非均一但规则的网格对模型进行离散化。在模型设计中俯冲带区域分辨率为2 km×2 km，模型边界区域分辨率设置为30 km×30 km。在模型构建过程中约有700万个拉格朗日自由移动点，表示各个岩石圈结构及属性，模型中岩体的物性参数属性以及温度等可以通过这些拉格朗日点进行传递。根据板块俯冲的动力学特征以及数值模拟[34]，构建初始模型主要包含两大部分（图1）：大陆板块和大洋板块。模型中设置大陆岩石圈厚度为130 km，地壳厚度为34 km，包含上层沉积物6 km，上地壳厚度为14 km，下地壳厚度为15 km。大洋板块岩石圈厚度是100 km，洋壳的厚度为8 km。不同岩石类型的物性参数参考表1和表2。
图 1 初始模型设置
Figure 1. Initial model configuration
表 2 数值模型中的主要材料参数
Table 2. Parameters of the materials in the numerical models
物质 状态 ρ0
/(kg·m−3)ρe
/(kg·m−3)Cp
/(J·kg−1·K−1)Ka
/(W·m−1·K−1)$T_{\simfont\text{固相}}^{\rm b} $/K $T_{\simfont\text{液相}}^{\rm b} $/K Hr
/(μW·m−3)α
/K−1β
/MPa黏滞性流变参数 塑性性质
Sin (φeff)空气 — 1 — 100 20 — — 0 0 0 A* 0 水 — 1 000 — 3 330 20 — — 0 0 0 A* 0 沉积物
(6 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 上地壳
(14 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 下地壳
(15 km)固态 3 000 — 1 000 K2 TS2 TL2 0.5 3×10−5 1×10−5 C* 0.15 熔融 2 500 G* 0.06 洋壳(8 km) 固态 3 000 3 800 1 000 K2 TS2 TL2 0.25 3×10−5 1×10−5 D* 0.15 熔融 2 900 H* 0.06 岩石圈—软流圈地幔 固态 3 300 — 1 000 K3 — — 0.022 3×10−5 1×10−5 E* 0.6 熔融 2 700 0.06 水化地幔 固态 3 200 — 1 000 K3 — — 0.022 3×10−5 1×10−5 F* 0.6 熔融 2 700 0.06 注：a. K1=[0.64+807/(TK+77)]exp(0.000 04P); K2=[1.18+474/(TK+77)]exp(0.000 04P); K3=[0.73+1 293/(TK+77)]exp(0.000 04P);
b. 当P＜1 200 MPa, TS1=889+17 900/(P+54)+20 200/(P+54)2; 当P＞1 200 MPa, TS1=831+0.06P. TL1=1 262+0.09P; 当P＜1 600 MPa, TS2=973–70 400/(P+354)+778×105/(P+354)2; 当P＞1 600 MPa, TS2=935+0.003 5P+0.000 006 2P2. TL2=1 423+0.105P模型设置过程中需要根据研究的问题加载多种边界条件。首先，模型的初始位移条件，本文模型计算域内顶部、左边界和右边界均为自由滑移边界[39-41]，底部边界为渗透性边界[42]以保证满足物质守恒。在模型内部边界中，对大洋岩石圈施加俯冲速率，大陆岩石圈施加朝向内陆的仰冲速率，并且2个速度在整个计算过程中保持不变；其次，模型的温度条件设置为：垂向上，模型的顶部边界温度为0 ℃，莫霍面温度为450 ℃，岩石圈地幔底部边界温度为1 300 ℃，软流圈地幔深度温度梯度是0.5 ℃/km。横向上，模型左右两边界的温度梯度为0 ℃/km，即保证模型在横向上零热流的散失[42]；最后，在地壳表面和自由滑动模型的顶部界面设计一层黏性较低的空气层（大陆地壳以上的黏性空气层厚度设计为10 km，洋壳上方的黏性空气层设置为12 km），该层能够用来分析地壳地表的地形变化[33]，也用来限定地表侵蚀和沉积[42-45]过程对地形的影响。
-
从模拟结果可以看出大洋板块俯冲初始阶段，大洋板块在低角度俯冲过程中可以从上覆大陆板片的底部刮掉约20 km厚的岩石圈地幔。因密度变大，刮下的岩石圈地幔随板片继续俯冲下沉进入到软流圈地幔并触发小尺度地幔对流。同时，大陆前缘因受到强烈挤压而不断增厚（图2a）。随着俯冲的进行，俯冲角度逐渐由低角度变为高角度俯冲，绝大部分的沉积物只能俯冲至200 km深（地幔楔与软流圈的交界处），但有少部分的沉积物伴随洋壳能够俯冲至500 km深。随着温度-压力的逐渐升高，这部分沉积物和洋壳能够发生部分熔融，并混合在一起，形成一个与周围地幔物质存在密度差的低密度熔融柱。该熔融柱在重力不稳定（浮力）作用的驱动下能够产生垂直向上的运移速度，不断上涌（图2b）。当模拟演化至44.5 Ma（图2c），熔融柱到达上覆岩石圈地幔底部。由于压力的快速降低，熔融柱能够发生减压熔融事件，并以柱顶部为中心成蘑菇状横向侵蚀岩石圈底部，并进一步导致岩石圈地幔的升温熔融。在模拟演化后期（46.9～52.9 Ma），持续俯冲的沉积物和洋壳不断混合熔融，并沿着原熔融柱运移通道不断上涌，大范围侵蚀上覆大陆岩石圈地幔底部（图2d—e）。在55～88.9 Ma期间，滞留在200 km以内浅的沉积物开始熔融形成岩浆，并沿着俯冲通道折返至增生楔中。这部分壳源性质的岩浆能够迅速积累形成岩浆底辟，并上侵至地表形成穹隆结构（图2f）。值得注意的是，在整个俯冲过程中，岩石圈地幔底部破坏减薄作用可以不断拓展，最终可达300 km。
图 2 与大洋板块俯冲相关形成的熔融柱对岩石圈地幔破坏的物质场演化过程
Figure 2. Material field evolution of mantle plume damage to lithospheric mantle associated with subduction of an oceanic plate
大洋板块俯冲过程以及熔融柱对岩石圈地幔破坏作用的深部过程直接导致了浅部的地形变化（图3a）。板块俯冲初始阶段，上覆大陆壳受到强烈的挤压，变形范围拓展宽达300 km，地形抬升最大高程为8 km（18.5 Ma）。从黏度场和第二不变偏应变率场（图4）看出因板块的持续俯冲（22.4～88.9 Ma），低黏度区和高应变区逐渐向内陆拓展，与之相对应的褶皱变形带向内陆拓展宽达700 km，地表隆升最大高度为8 km。但俯冲至88.9 Ma时，上覆大陆板块前缘的地形变化为隆升―沉降―隆升，该变化与沉积物岩浆底辟上侵过程相关（图2f）。更重要的是，深部熔融柱与大陆岩石圈作用的范围内，地形逐渐抬升，且变形范围也局限在300 km宽，与深部作用范围一致。
图 3 不同的时间动力地形瞬时变化（a）及地形随时间演化过程（b）
Figure 3. a. Different temporal dynamic terrain instantaneous changes (in meters) corresponding to Fig. 2, b. Represents the evolution of terrain over time
图 4 黏度场和第二不变偏应变率
Figure 4. Viscosity field and second constant deviation strain rate
-
对于大陆岩石圈地幔结构遭到强烈的破坏的现象，一个共同的认识是，由热浮力驱动的地幔柱以及俯冲作用可造成岩石圈的拆沉与减薄[10-11,46-48]。但本文中的数值模拟结果与前人的研究成果存在一定程度的差异。主要表现在热熔融柱形成来源不同。在上述数值模型中，由于大洋板块的俯冲，被带入软流圈的沉积物、少量陆壳以及洋壳物质发生部分熔融。该熔融物质与周围软流圈地幔存在温度和密度上的差异，从而产生扰动[3,28,49]，并在热浮力作用下向上侵蚀大陆岩石圈地幔底部。除此之外，俯冲过程中形成的俯冲板片熔融柱物质能够以脉动的形式快速上涌[5]，从图2可观察到，在约22 Ma时洋壳俯冲至深地幔（约600 km），由于物质的流变性导致洋壳部分熔融，从而与俯冲板块发生拆离[28,49]，形成第一股熔融柱。随着俯冲的继续，时隔约22 Ma之后，俯冲的部分沉积物和下地壳能够再度熔融混合，可形成第二股熔融柱（44.5 Ma），其可垂直向上运动，与大陆岩石圈地幔相互作用，并加速熔融侵蚀岩石圈地幔。熔融柱脉动式上涌的原因，可能与板块俯冲速率及温度结构有关。
从本文模拟结果来看，俯冲板片熔融柱的上涌及其对大陆岩石圈地幔破坏作用过程主要分为2个方向：垂向上，俯冲板片熔融柱在地幔中持续上升，到达岩石圈底部，热熔柱以柱顶为中心引发地幔对流，从而与岩石圈地幔发生热交换[50-53]，向上热熔侵蚀，最多可减薄20 km的岩石圈地幔。值得注意的是，熔融柱头部分到达岩石圈底部能够促使岩石圈地幔熔融，从而加速了破坏效率。其次，在横向运动上，熔融柱与岩石圈地幔接触时，开始以柱顶为中心不断向两侧水平方向大规模侵蚀岩石圈地幔底部，且破坏范围不断增大。但是当熔融柱横向扩展至俯冲带地幔楔时，侵蚀作用明显受阻，并停止破坏。这是因为在板块俯冲过程中，冷的洋壳以及俯冲脱水作用能够明显降低地幔楔的温度[54]，从而间接提高了地幔楔的强度，这一现象在本文数值模拟结果中表现为地幔楔的高黏度、低应变特征（图4）。因此，受到地幔楔的约束作用，熔融柱向俯冲带方向对岩石圈地幔的破坏最大只能拓展至地幔楔处，向陆内方向由于无明显约束可大规模侵蚀，且最终的破坏范围局限在约300 km。
-
地表地形的演化过程是构造变形和动力地形联合作用的结果[5]。其中，构造变形是板块间的相互作用导致岩石圈（或地壳）的挤压或拉伸，从而引起地形上的变化[55]，其可控制地表的山川、盆地的分布，波及范围可达数百千米。动力地形源于深部地幔的运动过程，是俯冲带（或地幔柱）地幔对流导致的地表一阶尺度上的地形变化[21]，其波及范围根据地幔对流情况可达数百至上千千米。对于动力地形变化对地幔对流的响应关系，前人利用地球物理数据，并结合Gplates和Citcoms地球动力学模拟软件进行了深入研究[56-57]，揭示了板块俯冲过程中深部地幔对流活动对动力地形的影响，但这些研究成果由于受到模拟算法的限制，往往忽略了构造变形对地形的控制。然后，本文的模拟结果却能够很好地揭示出俯冲带地区动力地形和构造地形间的叠合关系，即俯冲带岩石圈构造变形与深部地幔动力学过程（由熔融柱引起）间的叠合关系。
在本文模拟结果中的构造变形区，地表地形变化明显受板块俯冲的影响。在18～89 Ma期间（图3a），俯冲带的前陆地区因持续受到强烈的横向挤压作用发生变形，变形范围可达300 km，地形抬升最高处达8 km，主要位于陆缘弧。随着板块持续的俯冲，约89 Ma时在俯冲渠道内的沉积物能够发生部分熔融，并沿着俯冲通道发生折返，侵入至俯冲带增生楔中，形成隆起构造，地形表现为双峰状（图3a的黑色线）。值得注意的是，大洋板块俯冲的整个过程中，构造变形的影响范围不断缩小，且由陆域逐渐向海沟方向迁移，最后主要集中于俯冲带陆缘弧地区（约200 km宽）（图3b）。紧临陆缘弧，在原构造变形区存在一个相对稳定的高地形隆起带（图3a，2 800～2 900 km），地形高程稳定在4 km以内，该区域很好地限定了后期的构造变形范围，可称为构造变形和动力地形的叠合区。这是因为，在叠合区内，构造变形的影响较弱，而受俯冲板片熔融柱控制的动力地形又无法穿过低温地幔楔（高黏度低应变区）去影响构造变形区。因此，叠合区内的地形变化相对稳定。
动力地形的演化与板块俯冲引起的熔融柱活动密切相关，在俯冲带后方形成广泛的动力隆起带，并随着时间而不断发生变化（图5）：自22 Ma开始，俯冲板片熔融柱开始上升作用于岩石圈地幔，对应的地形开始抬升。随着熔融柱持续上升，熔融柱对岩石圈地幔底部的侵蚀作用逐渐增强，动力地形变化幅度增大，持续动力抬升，抬升至最高4 km（图3），且动力地形变化上水平范围也局限在300 km，这与岩石圈地幔的破坏范围（300 km）保持一致。在俯冲板片熔融柱与大陆岩石圈地幔相互作用的整个过程中，由于受到地幔楔（高黏度、低应变、低温度）的约束作用，熔融柱向陆内方向大规模破坏岩石圈地幔，动力地形抬升中心也向内陆方向迁移，这与构造地形变化方向相反。其次，在动力地形变化范围内还分布着紧密间隔的平行偏移断层（图3b），说明在这种情况下断裂的主要驱动力来自于活跃的熔融柱与岩石圈的相互作用[58]。因此，熔融柱活动的动力学过程与地表地形变化存在很强的耦合性[57]。基于以上分析，地表地形变化大致分为3个阶段：板块俯冲挤压上覆大陆壳，地形抬升（阶段1）；然后缓慢剥蚀，地形沉降并趋于稳定（阶段2）。随着板块持续俯冲，在深部形成的俯冲板片熔融柱在热浮力作用下向上侵蚀大陆岩石圈地幔，动力抬升，沉降中心逐渐向内陆迁移（阶段3），并且在地表发育堑垒构造[56]。由此可知，地表地形的变化可反映出大洋板块俯冲形成的熔融柱活动参与大陆板块构造演化的过程，而不是单纯的板块俯冲运动控制地表变化。
图 5 深部俯冲动力学过程与浅部地表变化的响应
Figure 5. Responses of shallow topographic changes to deep subduction dynamics
-
（1）大洋板块俯冲形成的熔融柱在热浮力作用下，减压熔融纵向侵蚀大陆岩石圈地幔，岩石圈地幔随之受热熔融减薄；在水平方向上，岩石圈地幔随着热熔柱横向扩展，其熔融侵蚀减薄范围增加。
（2）熔融柱与大陆岩石圈地幔的相互作用可引起地表地形的变化。在构造变形活跃区，陆缘弧地区因持续受到强烈的挤压发生变形，变形范围达300 km，地形抬升高达8 km。俯冲板片熔融柱对岩石圈地幔底部的侵蚀作用逐渐增强，动力地形变化幅度增大，持续动力抬升，抬升至最高4 km，且动力地形变化上水平范围也局限在300 km，这与岩石圈地幔的破坏范围保持一致。
Numerical simulation of subduction-induced molten plume: Destruction of overriding plate and its dynamic topographic responses
More Information-
摘要: 洋壳俯冲过程中温度、压力升高和密度差异，可导致俯冲板片熔融柱的快速上涌，并作用在上覆板块岩石圈地幔底部，从而导致岩石圈的破坏减薄以及地表形态的剧烈变化，该过程类似于地幔柱对岩石圈的破坏作用。目前，对于俯冲板片熔融柱的形成及其对岩石圈破坏程度的研究相对较少，特别是地表动力地形变化与深部岩石圈破坏作用之间的响应关系依然不清楚。为此，本文将利用I2VIS有限差分方法，基于质量守恒方程、动量守恒方程以及能量守恒方程，通过给定材料参数和一定边界条件，计算揭示俯冲洋壳在不同时间和不同深度下发生部分熔融并形成俯冲板片熔融柱的过程，从而模拟再现该熔融柱对上覆板块岩石圈的破坏过程，并进一步分析其导致的浅部地表地形变化响应。数值模拟结果显示，在大洋板片俯冲过程中，由俯冲的陆源沉积物以及洋壳形成的混合熔融柱垂向侵蚀岩石圈底部，造成岩石圈减薄。在熔融柱的横向侵蚀过程中，岩石圈地幔熔融范围增加，可达300 km。在地形变化方面，板块俯冲造成大陆前缘受挤压变形，引起构造变形，构造变形范围可达300 km。同时，与俯冲相关形成的熔融柱对岩石圈地幔底部的侵蚀作用逐渐增强，动力地形变化幅度增大，并持续抬升，最终可垂向抬升至4 km。动力地形的变化范围局限在300 km以内，这与岩石圈地幔的破坏范围保持一致。Abstract: In the process of oceanic crust subduction, with the increase in temperature and pressure and the difference in density, the subduction-induced molten plume will rise rapidly and act on the lithospheric mantle bottom of the overriding plate, which may lead to the decrease in the lithospheric damage and the drastic change of surface morphology. This process is similar to the destruction of the lithosphere caused by mantle plume. So far, there have been little studies on the formation of subduction-induced molten plumes and their damage to the lithosphere, especially on the responses of surface dynamic topographic changes to the deep destruction. Based on the conservation equations of matter, momentum and energy, the I2VIS finite difference method is adopted by the authors to calculate and reveal the partial melting of the subducted oceanic crust at different times and depths under given material parameters and boundary conditions. The process of forming a subduction-induced molten plume is obtained, and then the process of molten plume-lithosphere interactions is further simulated, and the response of shallow topographic changes are analyzed. The numerical simulation results show that in the process of oceanic plate subduction, the composite molten plumes formed by subducted terrigenous sediment and oceanic crust eroded the bottom of the lithosphere longitudinally, and resulted in lithospheric thinning. During the transverse erosion of the molten plumes, the melting range of the lithospheric mantle increases up to 300 km. In terms of geomorphic change, plate subduction results in compression deformation of the continental front, which may reach 300 km. At the same time, the erosion of the molten plumes associated with subduction to the bottom of the lithospheric mantle is gradually strengthened, and the dynamic topographic changes increased, while uplifting continued, and ultimately reached a figure of 4 km. The variation range of dynamic topography is limited to 300 km, which is consistent with the damage range of lithospheric mantle.
-
图 1 初始模型设置
a.模型主要区域（4 000 km×670 km）的初始物质场和温度场以及边界条件，图中白色的线条是以400 ℃为梯度绘制的等温线，黄色箭头是大洋俯冲速率和大陆仰冲速率；b.模型中不同颜色代表不同的岩石组成，分别是：1.空气层，2.海水，3.沉积物，4.洋壳，5.上地壳，6.下地壳，7.岩石圈地幔，8.软流圈，9.薄弱带，10、11.熔融的沉积物，12.熔融的洋壳，13.熔融的上地壳，14.熔融的下地壳
Figure 1. Initial model configuration
a.Enlargement (4 000 km×670 km) of the numerical model shows composition field and boundary conditions. The isotherms (white lines) are plotted for each 400 °C increment. Yellow arrow represents the subduction rates of oceanic plate and obduction rates of continental plate，b.The colored grid for different rock types: 1—air; 2—water; 3 —sedimentary cover; 4—oceanic crust; 5—upper continental crust; 6—lower continental crust; 7—lithospheric mantle; 8—athenospheric mantle; 9—weak zone mantle; 10 and 11—partially molten sediment; 12—partially molten oceanic crust;13 and 14—partially molten continental crust(5and 6)
图 2 与大洋板块俯冲相关形成的熔融柱对岩石圈地幔破坏的物质场演化过程
黑色数字表示实验积累时间，单位是Ma；白色的实线是温度线，数字表示温度，单位 ℃
Figure 2. Material field evolution of mantle plume damage to lithospheric mantle associated with subduction of an oceanic plate
The black number indicates the accumulated time of the experiment in millions of years. The solid white line is the temperature line, the number indicates the temperature, unit ℃
图 3 不同的时间动力地形瞬时变化（a）及地形随时间演化过程（b）
时间与图2相对应。红色虚线表示在不同时间的地形剖面
Figure 3. a. Different temporal dynamic terrain instantaneous changes (in meters) corresponding to Fig. 2, b. Represents the evolution of terrain over time
The dotted red line represents the topographic profile at different times, which is consistent with a
图 5 深部俯冲动力学过程与浅部地表变化的响应
（红色虚线表示动力地形扩展范围，黑色虚线表示构造变形范围，虚线重合区带表示地形变化稳定区）
Figure 5. Responses of shallow topographic changes to deep subduction dynamics
(The red dashed line indicates the extension range of dynamic topography, the black dashed line indicates the range of deformation, the dotted line superposition area indicates the stable area of topographic change)
表 1 数值模拟采用的黏滞性流变参数（据文献[35-38]）
Table 1. Parameters of viscous flow in the numerical experiments (after references[35-38])
标号 流变性质 E/(K·J·mol−1) V/(J·MPa−1·mol−1) n AD/(MPa−n·s−1) η0a/(Pa·s) A* 空气/水 0 0 1.0 1.0×10−12 1×1018 B* 湿石英（强） 154 0 2.3 3.2×10−6 1.97×1019 C* 斜长An75（强） 238 0 3.2 3.3×10−6 4.80×1024 D* 斜长石 An75 238 0 3.2 3.3×10−4 4.80×1022 E* 无水橄榄岩 532 8 3.5 2.5×104 3.98×1016 F*b 湿橄榄岩 470 8 4.0 2.0×103 5.01×1020 G*b 长英质熔体 0 0 1.0 2.0×10−9 5.00×1014 H 铁镁质熔体 0 0 1.0 1.0×10−7 1.00×1013 注：:a η0表示为有效黏滞系数，计算公式为：η0=(1/AD)×106n;
b 熔融的长英质熔体表示的是熔融的沉积物和地壳。
下载: 导出CSV
表 2 数值模型中的主要材料参数
Table 2. Parameters of the materials in the numerical models
物质 状态 ρ0
/(kg·m−3)ρe
/(kg·m−3)Cp
/(J·kg−1·K−1)Ka
/(W·m−1·K−1)$T_{\simfont\text{固相}}^{\rm b} $/K $T_{\simfont\text{液相}}^{\rm b} $/K Hr
/(μW·m−3)α
/K−1β
/MPa黏滞性流变参数 塑性性质
Sin (φeff)空气 — 1 — 100 20 — — 0 0 0 A* 0 水 — 1 000 — 3 330 20 — — 0 0 0 A* 0 沉积物
(6 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 上地壳
(14 km)固态 2 700 — 1 000 K1 TS1 TL1 2 3×10−5 1×10−5 B* 0.15 熔融 2 500 G* 0.06 下地壳
(15 km)固态 3 000 — 1 000 K2 TS2 TL2 0.5 3×10−5 1×10−5 C* 0.15 熔融 2 500 G* 0.06 洋壳(8 km) 固态 3 000 3 800 1 000 K2 TS2 TL2 0.25 3×10−5 1×10−5 D* 0.15 熔融 2 900 H* 0.06 岩石圈—软流圈地幔 固态 3 300 — 1 000 K3 — — 0.022 3×10−5 1×10−5 E* 0.6 熔融 2 700 0.06 水化地幔 固态 3 200 — 1 000 K3 — — 0.022 3×10−5 1×10−5 F* 0.6 熔融 2 700 0.06 注：a. K1=[0.64+807/(TK+77)]exp(0.000 04P); K2=[1.18+474/(TK+77)]exp(0.000 04P); K3=[0.73+1 293/(TK+77)]exp(0.000 04P);
b. 当P＜1 200 MPa, TS1=889+17 900/(P+54)+20 200/(P+54)2; 当P＞1 200 MPa, TS1=831+0.06P. TL1=1 262+0.09P; 当P＜1 600 MPa, TS2=973–70 400/(P+354)+778×105/(P+354)2; 当P＞1 600 MPa, TS2=935+0.003 5P+0.000 006 2P2. TL2=1 423+0.105P
下载: 导出CSV
-
[1] Morgan W J. Convection plumes in the lower mantle [J]. Nature, 1971, 230(5288): 42-43. doi: 10.1038/230042a0 [2] Wilson J T. A possible origin of the Hawaiian islands [J]. Canadian Journal of Physics, 1963, 41(6): 863-870. doi: 10.1139/p63-094 [3] Burov E, Guillou-Frottier L, d'Acremont E, et al. Plume head-lithosphere interactions near intra-continental plate boundaries [J]. Tectonophysics, 2007, 434(1-4): 15-38. doi: 10.1016/j.tecto.2007.01.002 [4] Christensen U R, Harder H. 3-D convection with variable viscosity [J]. Geophysical Journal International, 1991, 104(1): 213-220. doi: 10.1111/j.1365-246X.1991.tb02505.x [5] 李建康, 王登红. 地幔柱数值模拟研究进展[J]. 地质科技情报, 2005, 24(4):13-20. [LI Jiankang, WANG Denghong. Advances in the numerical simulation of the mantle plume [J]. Geological Science and Technology Information, 2005, 24(4): 13-20. doi: 10.3969/j.issn.1000-7849.2005.04.003 [6] 卢记仁. 峨眉地幔柱的动力学特征[J]. 地球学报, 1996, 17(4):424-438. [LU Jiren. Dynamic characteristics of EMEI mantle plume [J]. Acta Geoscientia Sinica, 1996, 17(4): 424-438. [7] 徐义刚. 地幔柱构造、大火成岩省及其地质效应[J]. 地学前缘, 2002, 9(4):341-353. [XU Yigang. Mantle plumes, large igneous provinces and their geologic consequences [J]. Earth Science Frontiers, 2002, 9(4): 341-353. doi: 10.3321/j.issn:1005-2321.2002.04.014 [8] Gorczyk W, Hobbs B, Gessner K, et al. Intracratonic geodynamics [J]. Gondwana Research, 2013, 24(3-4): 838-848. doi: 10.1016/j.gr.2013.01.006 [9] King S D, Ritsema J. African hot spot volcanism: small-scale convection in the upper mantle beneath cratons [J]. Science, 2000, 290(5494): 1137-1140. doi: 10.1126/science.290.5494.1137 [10] Liao J, Gerya T. Influence of lithospheric mantle stratification on craton extension: Insight from two-dimensional thermo-mechanical modeling [J]. Tectonophysics, 2014, 631: 50-64. doi: 10.1016/j.tecto.2014.01.020 [11] Yang S H, Li Z H, Gerya T, et al. Dynamics of terrane accretion during seaward continental drifting and oceanic subduction: Numerical modeling and implications for the Jurassic crustal growth of the Lhasa Terrane, Tibet [J]. Tectonophysics, 2018, 746: 212-228. doi: 10.1016/j.tecto.2017.07.018 [12] van Avendonk H J A, Holbrook W S, Lizarralde D, et al. Structure and serpentinization of the subducting Cocos plate offshore Nicaragua and Costa Rica [J]. Geochemistry, Geophysics, Geosystems, 2011, 12(6): Q06009. [13] Contreras-Reyes E, Grevemeyer I, Watts A B, et al. Deep seismic structure of the Tonga subduction zone: Implications for mantle hydration, tectonic erosion, and arc magmatism [J]. Journal of Geophysical Research: Solid Earth, 2011, 116(B10): B10103. doi: 10.1029/2011JB008434 [14] Fumagalli P, Stixrude L, Poli S, et al. The 10Å phase: A high-pressure expandable sheet silicate stable during subduction of hydrated lithosphere [J]. Earth & Planetary Science Letters, 2001, 186(2): 125-141. [15] Irifune T, Kubo N, Isshiki M, et al. Phase transformations in serpentine and transportation of water into the lower mantle [J]. Geophysical Research Letters, 1998, 25(2): 203-206. doi: 10.1029/97GL03572 [16] Ranero C R, Weinrebe W, Grevemeyer I, et al. Tectonic structure of the Middle America Pacific margin and incoming Cocos Plate from Costa Rica to Guatemala[C]//American Geophysical Union, Fall Meeting 2003. AGU, 2003. [17] Sano A, Ohtani E, Kubo T, et al. Effect of water on garnet-perovskite phase transformation in MORB system[C]//American Geophysical Union, Fall Meeting 2004. AGU, 2004. [18] Li Z H, Xu Z Q, Gerya T V. Flat versus steep subduction: Contrasting modes for the formation and exhumation of high- to ultrahigh-pressure rocks in continental collision zones [J]. Earth & Planetary Science Letters, 2011, 301(1-2): 65-77. [19] Campbell I H, Griffiths R W. Implications of mantle plume structure for the evolution of flood basalts [J]. Earth & Planetary Science Letters, 1990, 99(1-2): 79-93. [20] d’Acremont E, Leroy S, Burov E B. Numerical modelling of a mantle plume: the plume head-lithosphere interaction in the formation of an oceanic large igneous province [J]. Earth and Planetary Science Letters, 2003, 206(3-4): 379-396. doi: 10.1016/S0012-821X(02)01058-0 [21] Pekeris C L. Thermal convection in the interior of the earth [J]. Geophysical Journal, 1935, 3(8): 343-367. [22] Flament N, Gurnis M, Müller R D. A review of observations and models of dynamic topography [J]. Lithosphere, 2013, 5(2): 189-210. doi: 10.1130/L245.1 [23] Duretz T, Gerya T V. Slab detachment during continental collision: Influence of crustal rheology and interaction with lithospheric delamination [J]. Tectonophysics, 2013, 602: 124-140. doi: 10.1016/j.tecto.2012.12.024 [24] Gerya T V, Yuen D A. Rayleigh-Taylor instabilities from hydration and melting propel ‘cold plumes’ at subduction zones [J]. Earth & Planetary Science Letters, 2003, 212(1-2): 47-62. [25] Gerya T V, Yuen D A. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems [J]. Physics of the Earth and Planetary Interiors, 2007, 163(1-4): 83-105. doi: 10.1016/j.pepi.2007.04.015 [26] Vogt K, Gerya T V, Castro A. Crustal growth at active continental margins: Numerical modeling [J]. Physics of the Earth and Planetary Interiors, 2012, 192-193: 1-20. doi: 10.1016/j.pepi.2011.12.003 [27] Gerya T. Future directions in subduction modeling [J]. Journal of Geodynamics, 2011, 52(5): 344-378. doi: 10.1016/j.jog.2011.06.005 [28] Toussaint G, Burov E, Jolivet L. Continental plate collision: Unstable vs. stable slab dynamics [J]. Geology, 2004, 32(1): 33-36. doi: 10.1130/G19883.1 [29] Gorczyk W, Willner A P, Gerya T V, et al. Physical controls of magmatic productivity at Pacific-type convergent margins: Numerical modelling [J]. Physics of the Earth and Planetary Interiors, 2007, 163(1-4): 209-232. doi: 10.1016/j.pepi.2007.05.010 [30] Gerya T V, Meilick F I. Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts [J]. Journal of Metamorphic Geology, 2011, 29(1): 7-31. doi: 10.1111/j.1525-1314.2010.00904.x [31] Li Z H, Gerya TV. Polyphase formation and exhumation of high- to ultrahigh-pressure rocks in continental subduction zone: Numerical modeling and application to the Sulu ultrahigh-pressure terrane in eastern China [J]. Journal of Geophysical Research: Solid Earth, 2009, 114(B9): B09406. [32] Liao J, Wang Q, Gerya T, et al. Modeling craton destruction by hydration-induced weakening of the upper mantle [J]. Journal of Geophysical Research: Solid Earth, 2017, 122(9): 7449-7466. doi: 10.1002/2017JB014157 [33] Manglik A, Singh R N. Rheological stratification of the Indian continental lithosphere: Role of diffusion creep [J]. Proceedings of the Indian Academy of Sciences - Earth and Planetary Sciences, 1999, 108(1): 15-21. [34] Gerya T V, Stöckhert B, Perchuk A L. Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation [J]. Tectonics, 2002, 21(6): 6-1-6-19. [35] Kirby S H, Kronenberg A K. Rheology of the lithosphere: selected topics [J]. Reviews of Geophysics, 1987, 25(6): 1219-1244. doi: 10.1029/RG025i006p01219 [36] Kirby S H. Rheology of the lithosphere [J]. Reviews of Geophysics, 1983, 21(6): 1458-1487. doi: 10.1029/RG021i006p01458 [37] Ranalli G, Murphy D C. Rheological stratification of the lithosphere [J]. Tectonophysics, 1987, 132(4): 281-295. doi: 10.1016/0040-1951(87)90348-9 [38] Ranalli G. Rheology of the Earth[M]. 2nd ed. Netherlands: Springer, 1995. [39] Burg J P, Gerya T V. The role of viscous heating in Barrovian metamorphism of collisional orogens: thermomechanical models and application to the Lepontine Dome in the Central Alps [J]. Journal of Metamorphic Geology, 2005, 23(2): 75-95. doi: 10.1111/j.1525-1314.2005.00563.x [40] Li Z H, Gerya T V, Burg J P. Influence of tectonic overpressure on P-T paths of HP-UHP rocks in continental collision zones: thermomechanical modelling [J]. Journal of Metamorphic Geology, 2010, 28(3): 227-247. doi: 10.1111/j.1525-1314.2009.00864.x [41] Li Z H, Xu Z Q, Gerya T, et al. Collision of continental corner from 3-D numerical modeling [J]. Earth and Planetary Science Letters, 2013, 380: 98-111. doi: 10.1016/j.jpgl.2013.08.034 [42] Turcotte B, Schubert J. Geodynamics [J]. Geological curtain: English version, 2002, 450(2): 136-136. [43] Gerya T V, Yuen D A. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties [J]. Physics of the Earth and Planetary Interiors, 2003, 140(4): 293-318. doi: 10.1016/j.pepi.2003.09.006 [44] Huangfu P P, Wang Y J, Cawood P A, et al. Thermo-mechanical controls of flat subduction: insights from numerical modeling [J]. Gondwana Research, 2016, 40: 170-183. doi: 10.1016/j.gr.2016.08.012 [45] Schmeling H, Babeyko A Y, Ennsa A, et al. A benchmark comparison of spontaneous subduction models-Towards a free surface [J]. Physics of the Earth and Planetary Interiors, 2008, 171(1-4): 198-223. doi: 10.1016/j.pepi.2008.06.028 [46] Larsen T B, Yeun D A. Fast plumeheads: Temperature-dependent versus non-Newtonian rheology [J]. Geophysical Research Letters, 1997, 24(16): 1995-1998. doi: 10.1029/97GL01886 [47] Manga M, Stone H A, O'Connell R J. The interaction of plume heads with compositional discontinuities in the Earth's mantle [J]. Journal of Geophysical Research, 1993, 98(B11): 19979-19990. doi: 10.1029/93JB00441 [48] Van Keken P. Evolution of starting mantle plumes: A comparison between numerical and laboratory models [J]. Earth & Planetary Science Letters, 1997, 148(1-2): 1-11. [49] Burov E, Guillou-Frottier L. The plume head-continental lithosphere interaction using a tectonically realistic formulation for the lithosphere [J]. Geophysical Journal International, 2005, 161(2): 469-490. doi: 10.1111/j.1365-246X.2005.02588.x [50] 郭慧丽, 徐佩芬, 张福勤. 华北克拉通及东邻西太平洋活动大陆边缘地区的P波速度结构: 对岩石圈减薄动力学过程的探讨[J]. 地球物理学报, 2014, 57(7):2352-2361. [GUO Huili, XU Peifen, ZHANG Fuqin. P wave velocity structure of the North China Craton and West Pacific active continental margin: exploration for dynamic processes of lithospheric thinning [J]. Chinese Journal of Geophysics, 2014, 57(7): 2352-2361. doi: 10.6038/cjg20140729 [51] 李三忠, 索艳慧, 李玺瑶, 等. 西太平洋中生代板块俯冲过程与东亚洋陆过渡带构造-岩浆响应[J]. 科学通报, 2018, 63(16):1550-1593. [LI Sanzhong, SUO Yanhui, LI Xiyao, et al. Mesozoic Plate Subduction in West Pacific and Tectono-magmatic Response in the East Asian Ocean-Continent Connection Zone [J]. Chinese Science Bulletin, 2018, 63(16): 1550-1593. [52] Yu S Y, Li S Z, Zhang J X, et al. Multistage anatexis during tectonic evolution from oceanic subduction to continental collision: A review of the North Qaidam UHP Belt, NW China [J]. Earth-Science Reviews, 2019, 191: 190-211. doi: 10.1016/j.earscirev.2019.02.016 [53] Yu S Y, Zhang J X, Li S Z, et al. TTG-Adakitic-like (Tonalitic-Trondhjemitic) magmas resulting from partial melting of Metagabbro under high-pressure condition during continental collision in the North Qaidam UHP Terrane, Western China [J]. Tectonics, 2019, 38(3): 791-822. doi: 10.1029/2018TC005259 [54] 郑永飞, 陈仁旭, 徐峥, 等. 俯冲带中的水迁移[J]. 中国科学: 地球科学, 2016, 59(4):651-682. [ZHENG Yongfei, CHEN Renxu, XU Zheng, et al. The transport of water in subduction zones [J]. Science China Earth Sciences, 2016, 59(4): 651-682. [55] Liu S F, Nummedal D, Gurnis M. Dynamic versus flexural controls of Late Cretaceous Western Interior Basin, USA [J]. Earth & Planetary Science Letters, 2014, 389: 221-229. [56] Liu L J, Spasojević S, Gurnis M. Reconstructing Farallon plate subduction Beneath North America back to the Late Cretaceous [J]. Science, 2008, 322(5903): 934-938. doi: 10.1126/science.1162921 [57] 刘少峰, 王成善. 构造古地理重建与动力地形[J]. 地学前缘, 2016, 23(6):61-79. [LIU Shaofeng, WANG Chengshan. Reconstruction of tectono-paleogeography and dynamic topography [J]. Earth Science Frontiers, 2016, 23(6): 61-79. [58] Burov E, Gerya T. Asymmetric three-dimensional topography over mantle plumes [J]. Nature, 2014, 513(7516): 85-89. doi: 10.1038/nature13703 -
点击查看大图
图(5) / 表(2)
计量
- 文章访问数: 3085
- HTML全文浏览量: 554
- PDF下载量: 87
- 被引次数: 0
