2008年5月12日在青藏高原东缘龙门山造山带上发生了8.0级强烈破坏性大地震，已有众多研究者从不同角度对地震“孕育”、发生和发展的机理进行了研究和解释，取得了一系列重要的研究进展^[1-11].前人的多项研究表明：该研究区大震活动是青藏高原本体与邻域特异大地构造背景下，并在几个大块体的共同作用下，即印度板块北向强烈俯冲所导致的流变动力学过程所致^[12-18].

龙门山造山带恰置松潘-甘孜地块与四川盆地之间，从地形地貌上表现为山地向平原过渡的区域，地表高程变化剧烈，龙门山与山前地区的高差大于4000m，而地形陡变带的宽度则仅为15~20km.据1970-2006年多期水准观测结果，相对四川盆地，川西地区（包括龙门山断裂带）均表现为大面积的快速隆升，其上升速率一般在2.5~5.8 mm/a，相当于-0.8~-1.8 mGal/a的重力异常变化，反映出区域地壳垂直运动强烈，并具有短期垂向运动剧烈的波动特点^[19].重力场的分布表明，布格重力异常值变化剧烈，由负几十毫伽突降到负几百毫伽.根据Airy均衡模式，青藏高原东缘为弱负均衡重力异常区，几乎处于均衡状态；龙门山造山带则为正均衡重力异常区，均衡重力异常可达100×10^-5 m/s²；四川盆地为负均衡重力异常区^[20]，青藏高原和扬子克拉通重力已基本达到均衡，而龙门山造山带重力却尚未达均衡.由于龙门山地带处在明显的不均衡状态，故深部物质与能量存在强烈交换，具有特异的深层动力过程^{[1, 21-22]}，其动力演化过程复杂，不仅受到印度板块俯冲的近东西向水平拉张与大构造背景的影响，而且在深部垂向存在物质的重新分异和剧烈的重力均衡调整作用.因此，青藏高原东缘重力均衡作用在动力效应过程中的演化势必十分引人注目.

本文采用有限元的数值方法，对龙门山造山带及邻区在印度板块向北俯冲的作用下，依据位移场、应力场、应变场及应变能密度分布的数值模拟结果，研究重力均衡效应在横贯松潘-甘孜、龙门山和四川盆地剖面辖区构造运动的动力学效应，并据此探讨龙门山地区地壳的流变结构，进而为研究汶川大地震的孕震环境与深部流变结构提供参考依据.

本文取青藏高原东缘松潘-甘孜块体、龙门山造山带及其东侧扬子克拉通西北缘的四川盆地作为研究域.青藏高原东部已有人工源地震深部探测剖面8条，而穿越龙门山断裂系，并与松潘-甘孜块体和四川盆地相连的剖面主要有3条：阿尔金-都兰-龙门山剖面^[23]、巴塘（竹巴龙）-四川资中剖面^[24]和阿坝-理县-双流剖面^[25].这三条剖面从不同的方向穿过龙门山造山带并进入四川盆地，所得该区深部壳、幔结构有明显变化，其速度分布亦存在显著差异.结合天然地震观测^[26]以及电性结构特征等方面的工作，滕吉文等^{[1, 21]}通过对该区高精度人工源深部地震探测剖面的综合研究，提出了研究区壳、幔介质结构的参数分布表，如表 1所示.

表 1(Table 1)

表 1 研究区壳、幔介质几何结构参数 Table 1 The geometrical structure parameters of the crust and mantle in the model

表 1 研究区壳、幔介质几何结构参数 Table 1 The geometrical structure parameters of the crust and mantle in the model

现取研究区域垂直深度为150km的空间，并将其划分为地壳、岩石圈地幔和软流圈3层，松潘-甘孜块体地壳和岩石圈深度分别为60km和80km，四川盆地取为42km及105km，龙门山处的深度为在二者之间过渡.考虑到青藏高原东缘深约20km处存在一厚达5~10km的低速层，因此可将松潘-甘孜块体的地壳划分为上地壳厚15km、低速层厚5km及下地壳厚40km.在地表水平方向上，松潘-甘孜块体的长度取200km，龙门山造山带宽度取为60km，四川盆地取140km，模型水平向总长度为400km.在地势形态上，为便于单元划分，四川盆地东侧高程取1km，与龙门山交界的西缘过渡带为2.5km，龙门山块体以东青藏高原松潘-甘孜地带海拔高度取为5km.为简化起见，除龙门山地区的地壳及岩石圈区域的过渡型边界外，其它边界均考虑为直立或是水平状.

依据上述几何模型，本项工作采用Ansys多物理场有限元分析软件构建了如图 1所示的单元模型，不同块体的单元用不同的颜色来表示，各构造单元表征着不同的空间状态和物理-力学属性.模型采用了8节点Plane183的四边形平面单元，网格密度在软流圈内均设置为5×5km，其余层位均设置为3×3km，但对龙门山地带的地壳和岩石圈地幔以及松潘-甘孜低速层的单元网格进行加密，以得到更为精确的解答.整个有限元模型单元数9387，节点数28585.选用软件自动选择的稀疏矩阵求解器，它建立在与迭代法相对应的直接消元法基础上，即使存在病态矩阵也不会引起求解困难，但该求解器对内存及硬盘的配置要求较高.

图 1

Fig. 1

图 1 计算单元模型图 Fig. 1 Sketch of calculation elements of the numerical model

在计算过程中选择平面应变模式，二维模型边界条件的施加方式对于计算结果有重要影响，即不合适的边界条件易带来不合实际乃至错误的结果解释.因此，二维模型不仅需要在几何尺度上尽量减小因区域过小而导致的边界条件影响，同时要根据实际地质构造来合理约束其边界条件.因而，模型东部边界是稳定的四川盆地，基于区域地表构造展布，东部不是变形的主要区域，故此处边界条件简化为在模型的右边界施加全约束，而模型底部的上地幔底界面只施加竖直向约束，使其水平向能自由移动.

在单元本构关系的构建中，青藏高原东北缘松潘-甘孜地区上地壳及四川盆地地壳设置为弹性，其余部分选择为Maxwell体的黏弹性本构关系.为了表征不同地体在构造演化中的不同影响，本文对各地层采用了不同的物理-力学参数.其中，四川盆地作为力学性质稳定且坚硬的构造体，地壳与岩石圈地幔介质的弹性模量及黏滞系数均高于青藏高原东缘及龙门山地域.考虑到青藏高原东缘下地壳层处于温度较高的状态，速度低且发育有高导层，介质力学属性偏软，体现了较强的流变性，因而松潘-甘孜块体下地壳部分选取的黏滞系数小于龙门山地壳.

为了探讨研究区构造环境下的重力均衡效应，本文采用了不同的初始力源，分别为只考虑重力效应、只考虑水平力源效应，以及综合考虑水平力源和重力效应双重影响三种模型.第一种模型是在零初始位移和零初始应力条件下，仅仅考虑由重力作用产生的应力场作为初始场，这种模型常见于巨型高原重力坍塌对周边区域应力场作用的数值模拟工作^[27]；第二种模型则不考虑重力的影响，只选取研究区的运动学边界条件，即速度约束.模型左端设定的速度应满足青藏高原东缘GPS资料获得的速度场^[28]，约2.5cm·a^-1，这种模型常见于多数青藏高原隆升的数值模拟工作中^[29]；第三种模型则综合考虑了重力效应和速度边界条件的影响.为了进行对比，三种模型中对应单元所用的物理-力学参数相同，所选取的初始力学参数示于表 2.各模型总的求解时限均为0.1 Ma，即远大于岩石圈和软流圈的黏弹性松弛时间（约数十至数千年），以达表征出各圈层重力均衡调整的效应.各求解所选用的时间步长依据计算收敛情况在5000~15000步之间不等，运算时间步长均小于相应黏度条件下的应力松弛时间.

表 2(Table 2)

表 2 计算模型中的力学参数 Table 2 The mechanical parameters of materials in the model

表 2 计算模型中的力学参数 Table 2 The mechanical parameters of materials in the model

在三种模型所分别计算得到的竖直方向的应变场和位移矢量场分布中（图 2），虚线表征变形前的模型几何形状（后同）.由垂向应变结果可见，四川盆地地壳与软流圈的应变场大小趋近于零，保持了稳定的状态，体现了坚硬的力学性质.

图 2

Fig. 2

图 2 不同模型计算得到的竖直向应变分布（a1-c1）及位移矢量场分布（a2-c2） （A）模型一；（B）模型二；（C）模型三. Fig. 2 Distribution of strain field in the vertical direction (a1-c1) and displacement vector field (a2-c2) (A)model1;(B)model2;(c)model3.

从竖直向应变场的计算结果看，模型一中（图 2a1），由于只受到重力作用，四川盆地整体变形较小，在龙门山地域及松潘-甘孜地块深部应变则呈明显的层状分布，竖直向应变均基本属于压缩状态，数值大小随深度递增；模型二（图 2b1）受到水平向的压缩力源，计算结果呈现出明显的块状分布，显现了张应力的状态，数值大小明显地体现出三个不同块体之间的差别.松潘-甘孜块体弹性的上地壳竖直向应变值大小趋近于零，自低速层以下应变值显著增加，但基本上属于一致变形的趋势.相比应变场相对简单的松潘-甘孜块体和四川盆地，龙门山地域竖直向应变总体上呈现出从青藏高原东缘到四川盆地之间过渡型的复杂而又混乱的形态.地壳上部，尤其是与青藏高原东缘低速层抬升尖灭区域相接触的部位，为整个模型中竖直向应变值最大的区域，属于应力最为集中的地方，也正是汶川8.0级大地震发生的边界区；与模型一和模型二中比较单一的压缩或拉张的垂向应变场不同，模型三（图 2c1）垂向应变场则不但有压缩向，也存在拉张效应的区域，而且形态分布更为复杂，各块体内部存在着比模型二明显的层状分布，呈现出压缩力源作用下的重力分层调整机制.同样，龙门山区域在与松潘-甘孜块体尖灭的低速层相连的地方依然是应力最为集中、应变值最大的区域.与模型二不同的是龙门山下方岩石圈与软流圈分界的区域存在明显的应力、应变集中效应，这个区域的集中突变性在模型一中亦有所体现.

从位移场的矢量计算结果看，模型一（图 2a2）中位移量最大的区域发生在松潘-甘孜上地壳，其下方呈现层状分布，位移变形量逐渐减小，体现了深部岩石的静水压力效应.由于右端固支，而且四川盆地壳、幔介质力学性质相对坚硬，呈现出向左端位移量逐渐增大的带状趋势，而在此之间的龙门山区域却具有明显的从层状过渡到带状的效应.龙门山下方自岩石圈地幔深度开始，位移矢量方向便由重力作用下的向下方向逐渐发生偏转，在龙门山与四川盆地岩石圈地幔相接的地方出现明显向右的位移方向，呈现出高原东缘及龙门山软流圈在四川盆地下方东向运移并下沉的趋势；模型二（图 2b2）中，松潘-甘孜块体在龙门山区域呈现出明显的推覆构造，受到四川盆地的阻碍，而四川盆地位移变化量很小；模型三（图 2c2）在地壳深度范围内呈现出与模型二一样的趋势，但在岩石圈地幔及软流圈深度处却表现出与模型一类似的效应，即体现在由于受到边界压缩效应，青藏高原东缘在向东扩展的过程中，其深部物质在遇到坚硬的四川盆地时，部分物质形成分叉，向下绕过四川盆地较厚的岩石圈后逐渐下沉，而这也正是模型一和模型二两者之间明显的差异之一.

在地球物理深部探测所得结果的基础上由于简化而得的物理-力学模型，与实际地质构造之间尚存在一定的差异，二维模型中各块体初始几何形态同样会给计算带来一定的影响.若四川盆地岩石圈厚度明显增加，不仅四川盆地深部块体的初始几何形态会改变，龙门山下方岩石圈形态亦会随之改变.模拟计算结果显示，各块体的位移场、应力和应变场其计算结果在数值大小上虽然有一定程度的改变，但其演化形态的相互关系基本不变.

由此可见，尽管龙门山地区地表地形变化剧烈，深部地壳厚度及地壳密度差异明显，不仅导致其与邻近的四川盆地存在巨大的横向重力偏应力梯度差，同时它又远离印度板块呈北向挤压的区域，其主应力方向已经由喜马拉雅造山带的近南北向转化为近东西向.本文数值模拟结果表明，就最主要的动力学效应而言，印度板块向北的运动推力依然是龙门山地区深部动力学响应的主导因素，它决定了该研究区的主应力方向，并且促使龙门山地带的地壳上部与青藏高原东缘低速层抬升、尖灭部位相接处成为应力、应变最为集中的区域.然而，深部重力均衡效应导致了复杂的垂向动力学效应，引起了深部物质极为复杂的运移过程.

基于模型三的研究还发现，研究区岩石圈地幔及软流圈流变性的改变，尤其是黏度的改变，对整体动力学过程影响并不明显，可以说是并不敏感.但龙门山地带地壳黏度的改变，却会引起研究区均衡重力场的动力学效应发生明显的变化.可见，基于研究区均衡重力场动力学过程的数值模拟，乃是探索研究区龙门山造山带深部地壳流变结构及动力学响应的重要参考手段.

因此，在进一步的研究中，仍选用上述重力场效应及速度边界条件共存的第三种模型，在其所采用的初始力学参数中，各块体的本构关系不变，弹性参数的选取亦不变，黏度系数的选取则分别改变为松潘-甘孜地壳低速层为1.0×10²⁰Pa·s，下地壳为5.0×10¹⁹Pa·s，龙门山地壳为1.0×10²²Pa·s，松潘-甘孜、龙门山及四川盆地的岩石圈地幔黏度为5.0×10²¹Pa·s，软流圈黏度为1.0×10²⁰Pa·s.左边界速度场依然为2.5 cm·a^-1，模型总的求解时间为1 Ma，所选取的时间步数为35000步，其余求解设置不变.

图 3分别是0.25Ma（图 3a）、0.50Ma（图 3b）、0.75 Ma（图 3c）及1.0 Ma（图 3d）时计算所得的位移场分布图，而图 4则是各个阶段所对应的位移矢量场分布图.与虚线所代表的块体初始几何形态相对照，由图见，在整个构造变形中，松潘-甘孜地块及四川盆地位移场总的变形趋势变化不大，而龙门山地带变化明显.图中显示，松潘-甘孜块体及四川盆地所夹持的龙门山地带表现出明显的地壳缩短作用，推覆构造增长发育明显，反映出该区承载了青藏高原东缘现今地表的构造变形特征.同样，从位移矢量场的分布特征可以看出，龙门山区域地壳在缩短的过程中，受到坚硬四川盆地块体的阻挡，青藏高原东缘地壳低速层在压缩过程中不仅向上抬升，而且下地壳物质受到低速层影响，该低速层为滑脱面，向上运移，形成地表的推覆构造体系.而更深部位物质则形成了以软流圈顶部为滑脱面，在向东运移的过程中向下潜没于四川盆地岩石圈深处.

图 3

Fig. 3

图 3 位移场分布图 (a)0.25 Ma; (b)0.50 Ma; (c)0.75 Ma; (d)1.0 Ma. Fig. 3 Distribution of displacement field, (a) after 0.25 Ma, (b)after 0.50 Ma, (c)after 0.75 Ma, (d)after 1.0 Ma.

图 4

Fig. 4

图 4 位移矢量场分布图 (a)0.25 Ma; (b)0.50 Ma; (c)0.75 Ma; (d)1.0 Ma. Fig. 4 Distribution of displacement vector field, (a) after 0.25 Ma, (b) after 0.50 Ma, (c) after 0.75 Ma, (d) after 1.0 Ma.

图 5和图 6分别展示了四个阶段应变能及弹性应变能密度的分布.通过计算可得该研究区地壳内能量最为聚集的地方出现在龙门山地区深度约20km处，即汶川地震发震深度15±5km处^[1].在龙门山地带的地壳流变性变化时，则不论是增大还是降低其黏度达到一阶，也不论是在龙门山，还是在整个研究区的地壳深度范围内，不同地质时期，均无明显的能量聚集区，即无法形成孕震能量储集的优势层位.

图 5

Fig. 5

图 5 应变能分布图 (a) 0.25 Ma; (b) 0.50 Ma; (c) 0.75 Ma; (d) 1.0 Ma. Fig. 5 Distribution of strain energy field, (a) after 0.25 Ma, (b) after 0.50 Ma, (c) after 0.75 Ma, (d) after 1.0 Ma.

图 6

Fig. 6

图 6 弹性应变能密度分布图 (a)0.25 Ma; (b)0.50 Ma; (c)0.75 Ma; (d)1.0 Ma. Fig. 6 Distribution of elastic strain energy density field, (a) after 0.25 Ma, (b) after 0.50 Ma, (c) after 0.75 Ma, (d) after 1.0 Ma.

可见，基于黏弹性的流变结构模式，在流变性的对比中，龙门山地带地壳的黏度从力学性质上并不同于松潘-甘孜地块下地壳那样“热”、“软”，至少在其地壳上部并不容易发生塑性变形，而体现出较为“硬”的变形特征.

图 7a显示1.0Ma之后模型垂向位移的等值线分布，沿龙门山地壳中间区域选取不同深度的5个点（点1：地表；点2：8.8km；点3：21.4km；点4：36.2km；点5：51.4km），其随时间演化的垂向位移如图 7b所示.在三个块体地壳部分中，龙门山地壳变形最为明显，它与松潘-甘孜壳内低速层相连接的部分变形尤为突出，低速层和下地壳前沿有向上抬升趋势，并出现尖灭的现象，尖灭处亦是整个模型垂向位移，即隆升变形最为剧烈的区域.而图 7b显示，在计算模型中，1 Ma的时间尺度上，龙门山地壳内除选取的地表节点隆升高度约12km，深部四个节点的隆升高度均在3~6km，其隆升速度约为3~6mm·a^-1.由此，模型计算所得的竖向隆升速度与王庆良等^[30]采用以相邻水准点间高差变化速率为观测值的垂直形变网平差方法获得的川西地区近30年来的长期地壳垂直运动速度场2.5~5.0mm·a^-1的值基本一致.然而，不论是空间的水平向和深度方向，还是时间尺度上，龙门山地壳隆升情况均表现为差异性隆升，并非一致隆升.

图 7

Fig. 7

图 7 (a)1.0 Ma时垂向位移等值线分布; (b)龙门山地壳内不同深度节点垂向位移与时间关系; (c)龙门山地壳内不同深度节点垂向速度与时间关系 Fig. 7 (a)Contour map of the vertical displacements after 1.0 Ma; (b)Vertical displacements and time relation graph of the nodes in different depths beneath the Longmenshan crust; (c)Vertical velocity and time relation graph of the nodes in different depths beneath the Longmenshan crust

均衡重力异常场是深部壳、幔物质运动产生的静力平衡在重力场中的表征，变化幅度常随着新构造运动强度的增加而增加，研究均衡破坏可以了解区域内部发生的作用并获得有关新构造运动信息.因此，青藏高原东缘地貌与重力均衡之间的关系，对该区地球动力学研究具有十分重要的意义，同时，由于研究区当前动力学意义的解释常根据水平方向的观测结果（如GPS）来构建，因而垂向运动特征亟待进一步加强（王谦身，2010会议交流）.

（1）本文基于青藏高原东缘横穿松潘-甘孜、龙门山地带和四川盆地剖面的结构中重力均衡的动力学意义探讨，王谦身等^[31]的研究表明：根据Airy大陆均衡理论，理论均衡地壳厚度D=T+t，不同地形高程时该处地壳具有的理论均衡厚度t=4.45× H（H为地形高程），本研究区标准均衡地壳厚度T可近似地确定为39km.带入各块体的地形高程，可得研究区地壳理论均衡厚度分别为：松潘-甘孜块体约为61km，四川盆地约为43km，龙门山造山带在50~60km之间.显见相比之下，本文模型中松潘-甘孜块体和四川盆地的地壳基本达到重力均衡，而龙门山地带急剧变化的地形导致该区重力场处于明显的不均衡状态之下^[31].在印度板块的水平向挤压与本地区重力均衡垂向上调节的共同作用下，深部物质的重新分异、调整和运动当必会导致其深部复杂的物质与能量的强烈交换^[32]，而垂向运动则会导致水平向上更为复杂的构造运动.这对龙门山地区构造运动的深层动力学效应具有明显的影响，形成了特异的深层动力过程.

（2）四川盆地对青藏高原东部深部物质向东流动所起到的阻滞作用不仅仅局限于地壳尺度，而是深达上地幔软流圈的顶部，即整个岩石圈（包括地壳和上地幔盖层）^{[1, 21, 33]}.

青藏高原壳、幔物质在东向运移的过程中，受到四川盆地刚性壳、幔物质阻滞后，一部分绕四川盆地向东南和东北方向运移，体现为深部物质呈水平方向运动，这也得到大地电磁测深等相关工作支持；而另外一部分的岩石圈深部物质以低速层（圈）为滑脱面，沿龙门山断裂系向上运移，并在龙门山下方软流圈深处东向运移过程中向下潜入地幔，这可表现为该研究区印度-欧亚大陆连续汇聚时壳、幔物质运移的另外一种补偿方式.

（3）相对于四川盆地，该研究区的川西高原及松潘-甘孜地块部分，即贡嘎山相对四川盆地的现今隆起速率至少为5.8 mm·a^-1，其隆升速率之大，与喜马拉雅山5~10 mm·a^-1的隆起速率相比亦不逞多让.位于雅江、理塘之间的剪子弯山-卡子拉山地区现今隆起速率亦达4.0~4.8 mm·a^-1，可见研究区除四川盆地之外，现今仍处在差异性的快速隆升阶段^[30].

基于青藏高原东缘的快速隆升伴随着大规模的侵蚀作用^[34]，自中新世中晚期（12 Ma）以来，估计有8~10km的地壳物质被剥蚀^[35]，其剥蚀量之大，所引起的地壳均衡反弹，对研究区动力学效应有着巨大的影响.当然稳定的四川盆地亦不平静，Richardson等^[36]的研究表明，四川盆地在过去40 Ma年里，其上覆的沉积盖层中有厚达1~4km的层位被剥蚀掉.这种邻近区域不平衡的剥蚀量，将必使构造演化的动力学效应趋于愈发复杂化.

由此可见，在基于该区地壳与上地幔深、浅部物质相互作用研究，并结合龙门山地区壳、幔介质属性和结构变异以及地表剥蚀等信息，由挤压缩短与地表剥蚀共同作用下的重力均衡调整模型，可以合理地解释龙门山-青藏高原周缘最陡地形梯级带地区的动力学响应并可提供重要的深层过程参考依据^[37-38].

（4）松潘-甘孜在与龙门山交界处低速层发生尖灭的地方容易形成能量储集的优势层位^{[1, 39]}，表明龙门山造山带东西两侧显著的地形高差、大幅度的地壳厚度变化、介质属性和分层差异变化对汶川地震的孕育、发生和发展在能量储集方面均起到了关键性的作用.可见，松潘-甘孜块体、龙门山造山带和四川盆地所组成的推覆构造系统不仅承载着青藏高原东缘现今的构造变形，同时其能量储集方式还控制着这里的强震活动和发生与发展的趋势.

此外，应当清晰地认识到，应变能量首先开始在龙门山造山带地带的地壳中与松潘-甘孜宽体交界的西段一侧储集，随着地表推覆构造的发育，其能量储集区逐渐向中部及东部扩展.本文前述有关龙门山造山带地带地壳的流变性特征研究表明，其上部地壳具备较强的力学刚度，变形破裂将以脆性破裂的方式为主.因此龙门山地壳中所发育的复杂深部断裂同样可以使得龙门山造山带中段及东段的上、中地壳也均能成为能量储集危险度很高的地段.地震容易在应力集中达临界时突发，而并不一定稳定在某一区段.同样，这种能量储集差异的动力学过程所揭示的现象和发震断层西侧上盘的逆冲规模大于东侧四川盆地俯冲的事实相一致，故显示了汶川大震逆冲型发震断层地带独特的异常响应.

也有学者^[3]的研究认为，龙门山西侧的应变率强度乃是东侧的近4倍，而且衰减频度东侧亦比西侧快.

（5）由于资料的限制，本文采用的力学模型在物性参数和几何关系上均做了一些必要的简化处理，这对计算结果的量值会有一定影响，但对探索其深层过程和动力学响应的轨迹是十分有益的.显然随着工作不断深入和数据的不断采集，认识将必会深化，并将会不断完善，取得更为逼近实际地球动力学过程的结果^[40-44].