松潘—甘孜地块地壳和上地幔顶部现今变形模式数值模拟研究

万永魁, 刘峡, 刘瑞丰, 张扬, 沈小七, 郑智江

1 中国地震局地球物理研究所， 北京 100081 2 中国地震局第一监测中心， 天津 300180 3 中国地震局地质研究所， 北京 100029

0 引言

松潘—甘孜地块隆升机制主要有侧向挤出(Molnar and Tapponnier, 1975; Tapponnier et al., 2001; Hubbard and Shaw, 2009; Wang et al., 2011; Zhang, 2013)和中下地壳流(Royden et al., 1997; Clark and Royden, 2000; Clark et al., 2005a; Burchfiel et al., 2008)两种模式.汶川MS8.0地震在龙门山断裂带中南段和北段分别形成6.2±0.5 m、6.5±0.5 m的最大垂直位移，对应地壳缩短量约7.0 m，直接支持侧向挤出模型(徐锡伟等，2010).然而侧向挤出模式不能解释汶川地震余震多集中在5～20 km，25～40 km处的中下地壳仅有少量地震发生(黄媛等，2008；朱艾斓等，2008；陈九辉等，2009；易桂喜等，2012).龙门山前陆盆地新生代沉积分布范围有限，且发育不全，缺乏同期显著的构造挠曲，不是正常的盆—山耦合关系，这也与侧向挤出逆冲推覆模式的预测相矛盾(Burchfiel et al., 1995; Chen et al., 2000; Zhang et al., 2004; Li et al., 2012; Sun et al., 2018).

龙门山断裂两侧壳幔速度差异显著，相对四川盆地，松潘—甘孜地块普遍具有低速、高导的地球物理场特征，指示该区中下地壳和上地幔顶部可能存在壳幔流(Zhao et al., 2012; Liu et al., 2014；王绪本等，2017；王志等，2017；朱介寿等2017).松潘—甘孜地块深20～25 km处存在厚约5 km的低阻低速层(Li et al., 2009；朱介寿，2008；刘启元等，2009)，构成上地壳与中下地壳“解耦”运动的滑脱层或剪切带.龙门山后山、中央、前山和山前隐伏断裂向下延伸倾角逐渐变缓，最终收敛于该剪切带(Hubbard and Shaw, 2009; Yu et al., 2010; Wan et al., 2017；许志琴等，2007；张培震，2008).因此，目前对松潘—甘孜地块隆升之争，主要集中在仅受控于低阻低速层“解耦”还是同时受壳幔流共同作用这两种主流观点上.

姚琪等(2012)将低阻低速层构建为软弱薄层，与上地壳底面构成接触单元，利用与速度相关的非线性接触摩擦准则，模拟了低阻低速层对高原东缘隆升的影响，结果显示，低阻低速层可以控制松潘—甘孜地块快速隆升，但模型未考虑中下地壳和上地幔顶部变形的影响.尹力和罗纲(2018)、庞亚瑾等(2019)采用连续模型对青藏高原东缘垂向变形控制因素展开讨论，结果显示地形起伏、地壳结构和密度差异、岩石圈流变性质均是该区垂向变形的重要控制因素.汪昌亮(2012)沙箱模拟实验结果表明，脆性上地壳与中下地壳流先后对松潘—甘孜地块隆升起到作用.Xie等(2020)指出岩石圈均衡和岩石圈挠曲的静态支撑，以及下地壳流和地幔对流的动力作用是控制松潘—甘孜地块高海拔地形的主要机制.滕吉文等(2008，2014)指出松潘—甘孜地块中下地壳和上地幔顶部物质分别以埋深约20 km处的低阻低速层为第一滑移面，上地幔软流层顶面(深100±10 km)为第二滑移面，整体向SE方向运动，受到四川盆地深部“刚性”物质阻挡后，向上运移，物质汇聚、应力集中，从而引发汶川地震.胡幸平等(2012)采用弹性有限元模型对高原东缘深部构造变形模式展开讨论，结果表明，松潘—甘孜地块一侧，即模型北、西边界深100 km处水平运动是地表5倍的前提下，所得理论震源机制解与实际震源机制解最为吻合.此外，已有研究表明地幔对流拖曳力作为中国陆地岩石圈构造运动的重要驱动力，对青藏高原内部地壳运动方向有明显控制作用(黄建平等，2008；祝爱玉等，2019).综上所述，松潘—甘孜地块隆升可能与低阻低速层“解耦”，地形起伏、壳幔密度和结构差异、岩石圈流变等因素引起地壳至上地幔顶部水平运动速率逐渐增加以及地幔流对岩石圈底部的拖曳作用均有关联.

汶川地震后，针对龙门山断裂带及其周缘区域开展了一系列地球物理探测，新成像技术的广泛应用，使得获取更为精细的内部结构图像成为可能(Zhang et al., 2010；朱介寿，2008；朱介寿等，2017；王绪本等，2017；王志等，2017)，为构建二维有限元精细模型提供了条件.多期GPS观测、跨龙门山断裂水准观测以及最新的低温热年学剥蚀速率相关研究成果(Kirby et al., 2002; Godard et al., 2009; Li et al., 2012; Tian et al., 2013, 2015; Tan et al., 2017; Wang and Shen, 2020)，为模型边界加载提供了有效约束，同时也为模拟结果合理性检验提供了有力支持.为获得松潘—甘孜地块地壳和上地幔顶部现今变形模式，本文依据跨龙门山断裂带探测剖面，构建了长300 km、宽104 km的二维有限元接触模型(龙门山断裂带设置为接触单元)，在考虑岩石蠕变特性的前提下，以1991—2016年长期GPS观测结果为初始约束(图1a黑色箭头，Wang and Shen, 2020)，通过改变低阻低速层蠕变参数(即是否考虑低阻低速层的存在)以及模型西边界和底边界加载条件，共计开展了11项数值模拟实验.将模拟结果与综合多学科、不同时间尺度获得的隆升速率进行对比，各模型结果横向对比，讨论了松潘—甘孜地块地壳至上地幔顶部变形及物质运移模式，有助于进一步理解松潘—甘孜地块隆升以及汶川MS8.0地震孕育和发震机制.

1 模型

1.1 模型构建

青藏高原东缘由西向东依次为高海拔强烈变形变质作用的松潘—甘孜褶皱带、地势起伏剧烈的龙门山冲断带和低海拔弱变形的四川盆地(刘树根等，2019).晚三叠世以来，受印度板块与欧亚大陆强烈碰撞及中国陆地主体拼合的影响，松潘—甘孜地块和龙门山造山带发生强烈变形和冲断隆升(刘树根等，2009)，在东西宽约30～50 km范围内，西侧松潘—甘孜地块与东侧四川盆地形成高达约4 km地形差异，成为整个青藏高原地形梯度最大地区(Kirby et al., 2002; Li et al., 2012；李勇等，2005).龙门山断裂带呈北东—南西方向，全长约500 km，宽30～50 km，自北西向南东方向发育4条主干断裂，即汶川—茂县断裂(后山断裂)、映秀—北川断裂(中央断裂)、灌县—江油断裂(前山断裂)和大邑—广元断裂(山前隐伏断裂).断裂呈叠瓦状向四川盆地逆冲推覆，浅部断层倾角为60°～70°，向下延伸，倾角逐渐减小并收敛合并于埋深约20 km处的低阻低速层(Burchfiel et al., 2008；Yu et al., 2010；许志琴等，2007).低阻低速层厚约5 km，埋深于松潘—甘孜地块20～25 km处(Yu et al., 2010；刘启元等，2009).龙门山造山带南段、中段和北段分别以宝兴杂岩、彭灌杂岩和轿子顶杂岩为核心(刘树根等，2009)，具有强度大、能够积累高密度弹性应变能的特性(张培震等，2008).阿坝—双流人工地震探测剖面显示龙门山断裂带两侧地壳厚度变化显著，松潘—甘孜地块至四川盆地宽约80 km范围，地壳厚度由约60 km迅速递减至约45 km(朱介寿，2008).

依据阿坝—双流人工地震探测剖面(图1黑色粗实线)，参考前人数值模拟相关模型构架(Liu et al., 2015；朱守彪和张培震，2009；张竹琪等，2010；柳畅等，2014；陈棋福等，2015；祝爱玉等，2016；万永魁等，2017)，本文构建了长300 km、宽104 km的二维有限元接触模型.模型中龙门山4条主干断裂简化为1条，按照上陡下缓“铲式”结构，在20 km深度处与低阻低速层顶界相切.松潘—甘孜地块与四川盆地存在显著地形高差，中下地壳埋深差异更为显著，故模型上地壳设置了宽30 km过渡带，中下地壳设置了宽80 km过渡带.上地壳在过渡带内以强度较高的彭灌杂岩、宝兴杂岩及轿子顶杂岩为核心，其弹性模量为四川盆地的1.2倍.松潘—甘孜地块0～4 km范围代表高海拔地形，埋深20～25 km范围为低阻低速层(图2a黄色区域)，用于模拟上地壳与中下地壳的“解耦”.模型整体构架见图2a.

图1 研究区水准观测路线、GPS测站速度矢量(a)及垂直于龙门山断裂速度剖面(b)Fig.1 Leveling observation route and velocity vector of GPS in the study area(a) and velocity profile perpendicular to Longmenshan Fault (b)

1.2 参数设置

岩石在长时间载荷作用下应力、应变符合幂指数关系(Kirby,1983)，满足

(1)

B=Aexp(-E/RT),(2)

对于处于柔性变形阶段的岩石，其等效黏滞系数可利用公式(3)计算：

(3)

中国陆地地壳和上地幔顶部等效黏滞系数一般在1019～1024Pa·s，相对稳定，因此可以根据等效黏滞系数计算蠕变系数，即

(4)

本文在应变率设定为10-14s-1前提下(石耀霖和曹建玲，2008)，计算得到松潘—甘孜地块和四川盆地岩石圈各层蠕变系数.模型各层相关参数见表1和表2.

表1 松潘—甘孜地块岩石圈物质参数Table 1 Material parameters of rocks in the lithosphere in Songpan-Garzê block

表2 四川盆地岩石圈物质参数Table 2 Material parameters of rocks in the lithosphere in Sichuan basin

1.3 网格划分

模型由面单元plane182组成，龙门山断裂由二维接触单元contact171、targe169组成(图2a红色部位).采用三角形单元对模型进行网格划分，网格尺寸约1 km，网格划分结果见图2b，单元总数31749，节点总数32116个.

图2 模型构架及重力加载过程边界约束(a)和网格划分结果(b)Fig.2 Model structure and boundaries constraint of gravity loading process (a) and meshing results (b)

2 加载计算

2.1 加载

本文通过改变低阻低速层蠕变参数(是否考虑低阻低速层的存在)、模型西边界和底边界加载条件，共计开展了11项数值模拟实验，详情见表3.

表3 11项模拟实验低阻低速层和边界设置Table 3 Soft layer and boundary settings of 11 simulated experiments

加载分两步进行：(1)重力加载；(2)位移(构造)加载.重力加载时，模型东、西边界(x=300及x=0 km)水平向固定、垂向自由，底边界(y=-100 km)垂向固定、水平向自由，施加重力后维持1 Ma自由变形.位移加载是在重力加载的基础上完成的，一方面更接近实际，另一方面在构造变形过程中重力的影响不可忽略.依据1991—2016年GPS观测结果，模型300 km范围内，垂直于龙门山断裂带水平压缩速率约3 mm·a-1(图1b黄色透明区域)，故位移加载，首先考虑在模型东边界和底边界约束不变的前提下，西边界施加3 mm·a-1水平位移，作为基础模型(Model 1).基于基础模型，通过改变低阻低速层蠕变参数(参考已有研究由低阻低速层黏滞系数计算获得蠕变参数)，体现低阻低速层“解耦”对高原隆升的影响(Model 2).参考胡幸平等(2012)由地表至深100 km处，加载速率线性增加5倍，所得龙门山地区理论震源机制解与实际震源机制解最为吻合，刘峡等(2010)在华北地区现今地壳运动模拟研究中设定由地表至深100 km处，加载速率线性增加1.2倍，模拟结果与GPS观测结果最为一致，本文通过模型西边界加载速率随深度线性增加(1.2、1.5、1.8、2.0、2.5倍或3.0倍)，体现岩石圈水平向运动随深度增加对高原隆升的影响(Model 3—Model 8).最后，基于Model 5的模拟结果，通过改变底边界位移加载条件(类比地幔拖曳力)，体现地幔对流拖曳力对高原隆升的影响(Model 9—Model 10).此外，不考虑Model 5中低阻低速层的作用，模拟高原隆升速率(Model 11)，通过分析两模型隆升速率差异，进一步研究低阻低速层“解耦”对高原隆升的作用.

2.2 计算

蠕变和接触涉及材料非线性、几何非线性并与变形过程密切相关，接触状态事前通常未知，因此有限元处理蠕变和接触问题通常采用增量法、自动时间步长和N-R(Newton-Raphson)迭代联合求解.增量法首先将载荷分为Q0，Q1，Q2，…，若干步，对应位移分别为a0，a1，a2，…，假设已知第m步载荷Qm和相应的位移am，当载荷增加至Qm+1(Qm+1=Qm+ΔQm)，求解位移am+1(am+1=am+Δam)，理论上如果载荷增量ΔQm足够小，则可以计算出相应的位移增量.选择自动时间步长求解，ANSYS会依据载荷增量和时间步长自动将每一载荷步划分为若干子步进行求解，如果计算结果仍难以收敛，会通过增加子步数、减小子步时间，使子步载荷增量ΔQ足够小，从而实现收敛.N-R迭代的目的是为了改进计算精度，对于m+1次增量步的第n+1次迭代可表示为

3 模拟结果

3.1 重力和摩擦系数对模拟结果的影响

施加重力后模型会产生一定程度塌陷，将Model 5接触单元(龙门山断裂)摩擦系数设定为0.6，分别提取仅重力作用下0.5 Ma、1 Ma、2 Ma、4 Ma和6 Ma时模型各节点累计位移，计算0.5～1 Ma、>1～2 Ma、>2～4 Ma和>4～6 Ma四个时段因施加重力节点运动速率(图3)，结果显示，0.5～1 Ma和>1～2 Ma模型中、上地壳有一定程度的下沉，量值分别<0.06 mm·a-1和<0.03 mm·a-1，下地壳和上地幔存在由松潘—甘孜地块向四川盆地流动的趋势，量值分别<0.03 mm·a-1和<0.02 mm·a-1.>2～4 Ma和>4～6 Ma中、上地壳下沉和下地壳和上地幔东向流动趋势进一步减小，中、上地壳下沉速率分别<0.015 mm·a-1和<0.01 mm·a-1，反映出重力加载后随时间增加模型逐渐趋于稳定，这一结果与前人关于重力势作用可能是控制青藏高原边缘动力学变形的主要因素相矛盾，主要是因为模型边界采用了强制约束，这与实际边界条件存有差异.提取松潘—甘孜地块和四川盆地地表平均累计塌陷量，结果显示，0.5 Ma松潘—甘孜地块平均累计塌陷量为1197.6 m，四川盆地为699.7 m，1.0 Ma松潘—甘孜地块为1286.1 m，四川盆地为741.7 m(图4a)，由此得到0.5～1.0 Ma内松潘—甘孜地块和四川盆地因重力加载导致的垂向变形速率分别为0.057 mm·a-1和0.020 mm·a-1(图4b)，远小于两区域隆升速率，即加载重力后经过1 Ma的变形，因重力作用导致的垂向变形基本可以忽略，所以本文在重力加载1 Ma模型相对稳定后进行位移(构造)加载，但此时松潘—甘孜地块与四川盆地地形高差减小至约3500 m，与4000 m实际地形高差相差约500 m，因此本文模拟结果低估了重力势能的影响.针对龙门山断裂摩擦系数的选取，同样利用Model 5，尝试将接触单元摩擦系数分别设置为0.2、0.4、0.6、0.8和1.0，模拟结果显示，位移加载后节点174(x≈150 km)和节点21913(x≈250 km)在不同摩擦系数下隆升趋势几乎一致(图4c、4d).地表各节点垂直于龙门山断裂水平向压缩速率(图4e)和垂向隆升速率(图4f)也较为一致，反映出摩擦系数对模拟结果的影响并不显著，据此本文11项模拟实验摩擦系数统一设定为0.6.

图3 重力加载后不同时段节点平均运动速率(a) 0.5～1 Ma； (b) >1～2 Ma； (c) >2～4 Ma；(d) >4～6 Ma.Fig.3 Average velocity of nodes in different periods after gravity loading

图4 重力加载后不同时间地表相对高程(a)和相应时间段地表垂向速率(b)；构造加载后Model 5在不同摩擦系数下节点174(c)、节点21913(d)相对高程变化；构造加载20万年Model 5在不同摩擦系数下地表节点水平压缩速率(e)和垂向隆升速率(f)Fig.4 Relative elevation of surface nodes at different times after gravity loading (a) and vertical velocity in corresponding time period; relative elevation change of node 174 (c) and node 21913 (d) in Model 5 after displacement loading under different friction coefficients; surface nodes in Model 5 horizontal compression rate (e) and vertical rate (f) under different friction coefficients at loading 0.2 Ma.

图5 Model 5地表节点隆升位移随位移加载时间变化曲线Fig.5 The graph of uplift displacement of surface nodes with tectonic loading time of Model 5

图6 Model1-Model11构造加载20万年时地表水平向速率及GPS拟合结果(a)、(c)和隆升速率及水准结果(b)、(d)Fig.6 Simulation results of 11 models at tectonic loading of 0.2 Ma, surface horizontal rates and GPS (a)、(c) and uplift rates and leveling data (b)、(d)

3.2 模拟结果稳定性分析

为观察模型表面隆升位移随时间变化，在模型表面自西向东每隔约20 km取一个节点(node01-node 15).图5为模型Model 5位移加载1 Ma内各节点隆升位移随时间变化曲线(图5a黑色虚线框放大部分见图5c，图5b黑色虚线框放大部分见图5d).结果显示，位移加载开始至约8万年，各节点隆升速率(曲线斜率)逐渐增大，加载20万年后，隆升速率接近线性变化(图5c、5d)，反映出模拟结果基本稳定.因此，本文选取11项模型位移加载20万年时的模拟结果展开分析.

图6为11项模拟计算给出的地表各节点水平向和垂向运动速率，垂直于龙门山断裂GPS速度剖面的拟合曲线(图1b红色虚线，扣除整体运动)，以及基于1960—2010年水准观测数据采用动态平差方法获得的观测点垂向速率也显示在图6中(中国陆地现代垂直形变图集的编制与资料整编项目).模拟得到的水平速率自西而东逐步减小，在松潘—甘孜地块内减小缓慢，跨过龙门山断裂后，四川盆地内减小速度加快.从曲线形态看，Model 1-Model 8由折线逐渐向曲线过渡(图6a)，Model 9-Model 11与Model 5基本一致(图6c).GPS拟合曲线在松潘—甘孜地块内与Model 6最为一致，四川盆地内与Model 5最为一致.模拟给出的松潘—甘孜地块垂向隆升速率明显大于四川盆地，Model 2-Model 8垂向隆升速率逐渐增大，说明模型西边界加载位移越大，越能够促使地表隆升加速.水准观测结果与Model 5结果最为一致，而Model 9-Model 11结果与Model 5略有差异，反映出底边界加载和低阻低速层对松潘—甘孜地块隆升均有影响.考虑到GPS计算结果存在0.1～1.2 mm·a-1的误差，且大部分测站误差集中在0.2～0.6 mm·a-1，而Model 5水平向模拟结果与GPS拟合曲线最大误差小于0.2 mm，我们认为Model 5模拟结果与实际的水平向和垂向运动最为接近.

4 讨论

4.1 松潘—甘孜地块长期稳定隆升速率

Hao等(2014)、杜方等(2009)依据1975—1997年跨龙门山断裂带水准观测数据(图1白色圆点)，计算得到汶川地震前松潘—甘孜地块相对于四川盆地隆升速率高达4～6 mm·a-1.青藏高原东缘经历了由中生代—早新生代平缓至晚新生代突然加速的剥蚀过程(Arne et al., 1997).基于磷灰石裂变径迹(AFT)、锆石裂变径迹(ZFT)和(U-Th)/He等热年代学技术，前人研究了青藏高原东缘快速剥蚀的时间和速率(表4)，最新研究成果显示快速剥蚀始于约10 Ma，最大剥蚀速率约1.0 mm·a-1.如果隆升速率大于剥蚀速率，则地表隆升为正值，反之为负值.假定青藏高原东缘快速隆升始于距今10 Ma，在不考虑地震破裂引起地表垂向位移的前提下，取最大剥蚀速率1.0 mm·a-1，松潘—甘孜地块与四川盆地现今地形高差约为4000 m，计算得到最大隆升速率为1.4 mm·a-1(地表实际隆升速率=最大隆升速率-最大剥蚀速率，为0.4 mm·a-1)，与汶川地震前跨龙门山断裂带水准观测结果差异巨大.

基于中国陆地现代垂直形变图集的编制与资料整编项目，采用动态平差方法，中国地震局第一监测中心完成了中国陆地多期垂直形变图，其中1960—2010年跨龙门山断裂带长时段垂直形变结果显示松潘—甘孜地块相对四川盆地最大隆升速率约为1.5 mm·a-1(图6b、6d黑色虚线)，与低温热年代学结合现今地貌计算得到的约1.4 mm·a-1最大隆升速率较为接近.因此，Hao等(2014)、杜方等(2009)给出的汶川地震前松潘—甘孜地块相对四川盆地4～6 mm·a-1的快速隆升，作为长期动力地质学数值模拟的检验标准有待商榷.综合多学科、不同时间尺度研究成果，松潘—甘孜地块长期稳定最大隆升速率为1.4～1.5 mm·a-1可能更为客观.

4.2 低阻低速层对高原隆升的影响

Model 1与Model 2唯一差别是前者没有考虑低阻低速层，本研究通过对比这两项模拟结果，讨论低阻低速层的作用.从隆升速率看，Model 1和Model 2在松潘—甘孜东缘均出现了局部“鼓包”，最大隆升速率约为1.0 mm·a-1，与1.4～1.5 mm·a-1长期稳定最大隆升速率存有较大差异.而横坐标0～75 km模型范围，Model 1的地表垂向隆升速率比Model 2小约0.2 mm·a-1(图6b).

表4 青藏高原东缘快速剥蚀的时间和速率Table 4 Time and rate of rapid denudation in the eastern Margin of the Tibet Plateau

根据图7a—7d显示的水平向和垂向速率随深度变化，两模型运动速率的整体趋势基本一致.如松潘—甘孜内部隆升速率为0.4～0.7 mm·a-1，松潘—甘孜东缘至龙门山断裂带隆升速率最强为0.7～1.0 mm·a-1，四川盆地相对微弱为0.1～0.4 mm·a-1.所不同的是Model 2垂向隆升的范围大于Model 1，如隆升大于0.6 mm·a-1的区域，Model 1仅分布在地表横坐标50～200 km范围内，Model 2则覆盖了整个松潘—甘孜地块；又如隆升速率大于0.8 mm·a-1的区域，Model 1分布于横坐标100～175 km范围内，而Model 2扩展至模型100 km以西区域.上述结果揭示低阻低速层容易变形，有利于松潘—甘孜内部整体隆升，但不足以形成松潘—甘孜地块如此规模的强烈的地表抬升.

4.3 松潘—甘孜地块地壳和上地幔顶部变形模式

Model 2与Model 5唯一差别是后者西边界加载位移随深度线性增加1.8倍，即由地表3 mm·a-1至深100 km处增至5.4 mm·a-1.根据图7e—7f显示的Model 5水平向和垂向速率随深度变化，松潘—甘孜内部，水平向运动在低阻低速层发生了“解耦”，上地壳运动速率小于中下地壳.垂向速率与Model 2相比存有较大差异：(1)Model 2地表最大隆升速率约为1.0 mm·a-1，而Model 5约为1.44 mm·a-1，隆升速率显著增加；(2)Model 2仅在上地壳150 km附近隆升速率接近1.0 mm·a-1，Model 5隆升速率大于1.0 mm·a-1的区域基本扩展至整个松潘—甘孜地块，强烈隆升区域范围显著扩展；(3)Model 2隆升中心位于模型150 km附近，Model 5向西转移至模型100 km附近，与长时段水准观测结果更为一致(图6b).

图7 Model 1、Model 2、Model 5、Model 9-Model 11位移加载20万年水平向和垂向速率随深度变化(a) Model 1水平向速率； (b) Model 1垂向速率； (c) Model 2水平向速率； (d) Model 2垂向速率； (e) Model 5水平向速率； (f) Model 5垂向速率； (g) Model 9水平向速率； (h) Model 9垂向速率； (i) Model 10水平向速率； (j) Model 10垂向速率； (k) Model 11水平向速率；(l) Model 11垂向速率.Fig.7 Model 1, Model 2, Model 5 and Model 9-Model 11 simulation results of horizontal and uplift rate vary with depth at tectonic loading of 0.2 Ma(a) Model 1 horizontal rate; (b) Model 1 uplift rate; (c) Model 2 horizontal rate; (d) Model 2 uplift rate; (e) Model 5 horizontal rate; (f) Model 5 uplift rate; (g) Model 9 horizontal rate; (h) Model 9 uplift rate; (i) Model 10 horizontal rate; (j) Model 10 uplift rate; (k) Model 11 horizontal rate; (l) Model 11 uplift rate.

Model 1、Model 2和Model 5速度矢量见图8.Model 1和Model 2速率减小区域主要集中在松潘—甘孜东缘及龙门山断裂带中下地壳和上地幔顶部(图8a、8b)，故松潘—甘孜东缘地表隆升速率出现局部“鼓包”.与Model 1、Model 2相比，Model 5最显著的特征是松潘—甘孜内速度矢量旋转幅度显著增强，中下地壳和上地幔水平向运动速率快速减小区域集中在模型50～150 km处，松潘—甘孜东缘和龙门山断裂带(模型150～200 km)，速率减小已不显著，这与水准观测隆升中心并非在龙门山断裂带，而是在其西侧约100 km处的川西高原，即模型50～150 km范围相吻合.中下地壳和上地幔物质水平向运动速率快速减小与四川盆地的阻挡固然有关，但与地形起伏、地壳结构和密度差异、岩石圈流变及地幔对流拖曳力可能均有关联，至于每一项因素的影响程度还需开展进一步研究.

图8 Model 1、Model 2和Model 5构造加载20万年速度矢量模拟结果(a) Model 1速度矢量； (b) Model 2速度矢量； (c) Model 5速度矢量Fig.8 Model 1, Model 2 and Model 5 simulation results of velocity vector at tectonic loading of 0.2 Ma(a) Model 1 velocity vector； (b) Model 2 velocity vector； (c) Model 5 velocity vector.

与Model 5相比，Model 9和Model 10对模型底边界进行了位移加载，用于类比地幔拖曳力，Model 9底边界加载速率采用线性递减方式，由西边界5.4 mm·a-1至东边界递减为0，Model 10采用非线性递减方式，速率减小主要集中在松潘—甘孜东缘和龙门山断裂带(模型100～200 km范围)，两模型底边界加载速率变化曲线见图9.模拟结果显示，Model 9松潘—甘孜内地表隆升速率约1.0 mm·a-1，四川盆地内约0.50 mm·a-1，Model 10松潘—甘孜内地表隆升速率由0逐渐增至约1.51 mm·a-1，而后快速减小，四川盆地内约0.35 mm·a-1，但最大隆升速率位于松潘—甘孜东缘(模型约150 km处).整体而言，Model 9和Model 10地表隆升速率与水准观测结果均存有了一定差异(图6d)，垂向速率随深度变化结果与Model 5模拟结果差异也较为显著(图7f、7h、7j).上述结果反映出隆升速率和隆升中心的位置分布对底边界加载条件较为敏感.

图9 Model 9、Model 10底边界加载速率变化曲线Fig.9 Bottom boundary loading rate curve of Model 9 and Model 10

与Model 5相比，Model 11不考虑低阻低速层的作用，模拟的水平向和垂向运动速率见图7k、7l，与Model 5相比，Model 11中上地壳与中下地壳水平方向未发生“解耦”，隆升中心向龙门山断裂带一侧偏移，隆升速率大于1.2 mm·a-1的区域向深部延伸.模型西侧(约0～75 km)垂向隆升速率显著减小，地表隆升速率与水准观测结果产生了一定差异(图6d)，这一差异与Model 1和Model 2在模型西侧地表隆升速率存有的差异是类似的(图6b)，再次揭示低阻低速层易于变形，可以作为中下地壳和上地幔物质运移的一个滑移面，构成上地壳与中下地壳的“解耦”带，有利于松潘—甘孜地块整体隆升.

4.4 青藏高原东缘地壳和上地幔顶部物质运移模式

Model 5在构造加载50万年和100万年时水平向和垂向速率随深度变化见图10.与构造加载20万年结果相比，构造加载50万年和100万年时松潘—甘孜地块整体运动趋势基本一致，即水平方向上地壳与中下地壳“解耦”，垂向方向整体隆升，隆升中心位于模型约100 km处.但值得注意的是最大隆升速率却不断增加，构造加载20万年最大隆升速率约1.44 mm·a-1，50万年增至约1.8 mm·a-1，100万年增至约2.0 mm·a-1.Clark and Royden(2000)指出青藏高原东流物质受坚硬四川盆地阻挡后向两方向流出，一部分转至北东向，另一部分转至南东向.邹镇宇等(2015)计算的多期GPS结果显示，相对于稳定华南块体GPS速度场在青藏高原东缘逐渐衰减，并发生明显转向，在龙门山断裂带中北段附近向北东方向偏转，在鲜水河断裂带附近向南东方向偏转.快剪切波偏振优势方向在龙门山断裂北段为NE向，南段为NW向(石玉涛等，2009).汶川MS8.0强余震震源机制解P轴分布也表现出类似的特征(胡幸平等，2008；郑勇等，2009).上述研究说明青藏高原东缘物质汇聚到一定程度后，可以向两侧运移，这与本文数值模拟获取的隆升速率随时间逐渐增加的结果有所不同，考虑到本文2D模型具有一定封闭性，汇聚物质并不能向周缘扩散，导致应力不断增加，隆升不断加快，因此本文选定稳定初期(20万年)模拟结果展开分析.另外，模拟过程也未考虑地壳增厚导致的拆沉和地幔物质的底侵作用，同样会对模拟结果造成影响.实际松潘—甘孜地块地壳和上地幔顶部物质，受四川盆地阻挡，当物质汇聚到一定程度后会被迫向两侧运动，同时也会发生地壳拆沉和地幔物质底侵作用，实现物质运移的动态平衡(图11).本文模拟结果随构造加载不断进行，隆升速率逐渐增大，也是对青藏高原东缘物质汇聚、扩展，存在稳定开放物质流动系统的侧面印证.此外，构造加载50万年和100万年时，模型西边界中下地壳和上地幔出现了不同程度的塌陷(图10b、10d黑色区域)，也是因为物质东移，在本文模型中缺少相应的物质补充导致.

图10 Model 5构造加载50万年和100万年水平向和垂向隆升速率结果(a) 50万年水平向速率； (b) 50万年垂向速率； (c) 100万年水平向速率； (d) 100万年垂向速率.Fig.10 Model 5 simulation results of horizontal and uplift rate at tectonic loading of 0.5 Ma and 1 Ma(a) Horizontal rate at 0.5 Ma; (b) Uplift rate at 0.5 Ma; (c) Horizontal rate at 1 Ma; (d) Uplift rate at 1 Ma.

5 结论

本文基于二维有限元接触模型，在考虑地形起伏、地壳结构和密度差异、重力及岩石长期载荷蠕变作用的前提下，依据1991—2016年GPS观测结果，通过改变低阻低速层蠕变参数(是否考虑低阻低速层的存在)、模型西边界和底边界加载条件，共计开展了11项数值模拟实验，将模拟结果与1960—2010年长时段水准观测结果进行对比，结合低温热年代学隆升及剥蚀速率相关研究成果，讨论了松潘—甘孜地块地壳至上地幔顶部现今变形及物质运移模式.主要得到以下结论：

(1) 不考虑剥蚀作用的前提下，松潘—甘孜地块长期稳定最大隆升速率为1.4～1.5 mm·a-1可能更为客观.

(2) 松潘—甘孜块体内低阻低速层可以作为中下地壳和上地幔物质运移的一个滑移面，构成上地壳与中下地壳的“解耦”带，促进松潘—甘孜块体整体隆升，但低阻低速层对隆升影响有限，不足以形成如此规模的强烈的地表抬升.

图11 青藏高原东缘物质运移模式示意图Ⅰ 松潘—甘孜地块； Ⅱ 四川盆地； Ⅲ 西秦岭块体； Ⅳ 川滇“菱形块体”. ① 龙门山断裂；② 龙日坝断裂； ③ 鲜水河断裂；④ 东昆仑断裂； ⑤ 岷江断裂； ⑥ 虎牙断裂.Fig.11 The pattern of material transport mode in the eastern margin of the Tibet Plateau Ⅰ Songpan-Garzê block; Ⅱ Sichuan basin; Ⅲ West Qinling block; Ⅳ Sichuan-Yunnan rhombic block； ① Longmenshan fault; ② Longriba fault; ③ Xianshuihe fault; ④ East Kunlun fault; ⑤Minjiang fault; ⑥Huya fault.

(3) 考虑低阻低速层作用，模型西边界加载速率由3 mm·a-1随深度线性增加1.8倍的前提下，地表隆升速率与长时段水准观测结果及现今地貌结合低温热年代学获得的隆升速率最为吻合.在这一动力条件下，中下地壳和上地幔物质受四川盆地强烈阻挡，速率快速减小区域主要分布于松潘—甘孜地块内部(模型50～150 km范围)，松潘—甘孜地块东缘和龙门山断裂带速率减小已不显著，故地表隆升中心位于川西高原内而非龙门山断裂带.

(4) 青藏高原东缘物质汇聚到一定程度后，会被迫向两侧运动，从而实现物质运移的动态平衡.

本文模拟过程存在以下不足：(1)重力加载模型稳定后，松潘—甘孜地块相对四川盆地地形高差减至约3500 m，与4000 m实际地形高差相差约500 m，低估了重力势的影响；(2)模拟过程未考虑地壳增厚导致的拆沉和地幔物质的底侵作用.基于获得的松潘—甘孜地块至四川盆地地壳和上地幔顶部现今变形及物质运移模式，结合龙门山断裂带周缘精细速度图像结果，构建三维有限元模型，研究汶川MS8.0地震NW走滑型余震带(米亚罗断裂)动力学机制，揭示青藏高原东缘挤压动力学背景下，走滑型强震成因，可以作为今后工作的一个研究方向.

致谢图件由GMT绘制而成，GPS数据采用了中国地震局地质研究所王敏研究员计算结果，水准观测结果采用了中国陆地现代垂直形变图集的编制与资料整编项目组计算结果，在此一并表示感谢！