地球物理学报  2012, Vol. 55 Issue (3): 863-874   PDF    
青藏高原东缘变形机制的讨论: 来自数值模拟结果的限定
姚琪1 , 徐锡伟1 , 邢会林2 , 张微3 , 高翔1,4     
1. 中国地震局地质研究所,国家地震活动断层研究中心,北京 100029;
2. Earth Systems Science Computational Centre(ESSCC),School of Earth Sciences, The University of Queensland,St.Lucia,QLD 4072.Australia;
3. 国土资源厅航遥中心,北京 100083;
4. 浙江大学地球科学系,杭州 310027
摘要: 青藏高原东缘的地壳结构是两种主流青藏高原隆升模式争辩的焦点之一.中下地壳流曾经被认为是高原东缘隆升的主要构造驱动力,但是中上地壳之间低阻低速层的发现及其与2008 MS8.0汶川地震良好的对应关系表明,高原东缘具有向东刚性挤出的可能性.然而大部分关于龙门山断裂的数值模拟仍建立在下地壳流的基础上,仅将低阻低速层作为断裂的延续或是弱化地壳物性参数的软弱层,而非能够控制块体滑动的"解耦层",也没有考虑到刚性块体变形中的断裂相互作用.本文建立了包含相互平行的龙门山断裂与龙日坝断裂的刚性上地壳模型,用极薄的低阻低速层作为块体滑动的解耦带,采用速率相关的非线性摩擦接触有限元方法,基于R最小策略控制时间步长,计算了在仅有侧向挤压力作用下,低阻低速层对青藏高原东缘的刚性块体变形和断裂活动的作用.计算结果显示,低阻低速层控制了刚性块体的垂直变形和水平变形分布特征.在侧向挤压力的持续作用下,在低阻低速层控制下的巴颜喀拉块体能够快速隆升,而缺乏低阻低速层的四川盆地隆升速度和隆升量均极小,隆升差异集中在龙门山断裂附近,使其发生应力积累乃至破裂.龙日坝断裂被两侧的刚性次级块体挟持着一起向南东方向运动,但该断裂的走滑运动分解了绝大部分施加在块体边界上的走滑量,使得相邻的龙门山次级块体的走滑分量遽然减少,也使得龙门山断裂表现出以逆冲为主,兼有少量走滑的运动性质.本文所得的这些计算结果显示了在缺乏中下地壳流,仅在低阻低速层解耦下刚性块体隆升过程及相关断裂活动,提供了青藏高原东缘刚性块体挤出的可行性,为青藏高原东缘隆升机制的研究讨论提供了重要依据.
关键词: 青藏高原东缘      低阻低速层      有限元      龙门山断裂      龙日坝断裂     
Deformation mechanism of the eastern Tibetan plateau: Insights from numerical models
YAO Qi1, XU Xi-Wei1, XING Hui-Lin2, ZHANG Wei3, GAO Xiang1,4     
1. National Center for Active Fault Studies,Institute of Geology,China Earthquake Administration,Beijing 100029,China;
2. Earth Systems Science Computational Centre(ESSCC),School of Earth Sciences, The University of Queensland,St.Lucia,QLD 4072.Australia;
3. China Aero Geophysical Survey & Remote Sensing Center for Land and Resources, Beijing 100083, China;
4. Earth Science Department of Zhejiang University, Hangzhou 310027, China
Abstract: The crustal structure of the eastern Tibetan plateau plays a key role in studying the uplift mechanism of Tibetan. Ductile flow in middle-lower crust was considered to be the main driving force in the uplifting of eastern Tibetan plateau. However, both the discovery of the low resistivity and low velocity layer between the mid-upper crust and its correspondence with the 2008 Ms8.0 Wenchuan earthquake show the possibility of southeastward splaying of coherent lithospheric blocks in eastern Tibetan plateau. Unfortunately, popular FEM models of eastern Tibetan plateau and Longmenshan fault zone in recent years were established with viscous flow of the crust. And the low resistivity and low velocity layer was set as a fault plane or a factor weakening the crust in these models. Few researchers have addressed the decoupling problem of the low resistivity and low velocity layer. The interaction of faults was ignored in previous FEM models either, for continuous thickening and deformation. In this paper, we present a coherent block model with elasto-plastic material, containing Longmenshan fault and Longriba fault which parallel to each other, setting a thin soft layer under the upper crust as the decoupling layer. A non-linear rate-independent finite element method (including the frictional contact) was carried out using R-minimum strategy to limit the step size, to simulate the block deformation of upper crust and fault behavior of the eastern Tibetan plateau. The results show that the low resistivity and low velocity layer controls the distribution of horizontal and vertical deformation on surface. The continuous lateral extrusion pressure causes quick and huge uplift of the Bayan Har block which is under the control of decoupling layer, while little uplift of the Sichuan Basin which is not controlled by decoupling layer. The difference of uplifts makes stress accumulating on the Longmenshan fault, which is the border of the two blocks, and results in large earthquake on the fault, as Wenchuan earthquake. The Longriba fault moves southeastward with the two rigid subblocks. While its shear deformation absorbs most of the strike-slip displacement added on the model boundary, and strongly reduces the strike-slip displacement in the adjacent subblock and Longmenshan fault. Therefore, the Longmenshan fault shows dominantly dip-slip reverse faulting, with a little right-lateral oblique faulting. This paper introduces deformations of coherent blocks and relative fault activities only controlled by the decoupling of the low resistivity and low velocity layer, without any other influence from ductile flow in middle-lower crust. It should be possible, therefore, to oblique stepwise rise in the eastern Tibetan plateau. This could eventually lead to more discussion about the growth of the Tibetan..
Key words: The eastern Tibetan plateau      Low resistivity and low velocity layer      Finite-element method      Longmenshan thrust belt      Longriba fault     
1 引言

地壳内部复杂的结构及变形既是力学作用的结果,也能作为力学载体改变地壳变形和断裂活动特征.青藏高原地壳结构是两种主流高原隆升机制---中下地壳管流(Channel Flow)说[1-2]和侧向逃逸学说(Oblique Step wise Rise)[3-4]争辩的焦点,而青藏高原东缘作为地壳物质向东逃逸的必经之路,其地壳结构一直以来都受到广泛的关注.这两种隆升模式对高原内部地壳结构和相应的断层活动有不同的设想:Channel Flow 模式认为印度板块向北推挤导致韧性下地壳增厚,中下地壳流挤压上地壳,造就高原东缘巨大的高程差,上地壳发育大量缓慢滑动断层,发生分布式的连续变形[15];而侧向逃逸模式[2]更注重大型斜向走滑断裂在高原扩展过程中所起的作用,关注的是刚性中上地壳的变形及这种变形的传递扩展.

地球物理勘探和岩石学实验证明青藏高原地壳中部或是下部存在一个软弱带[6],成为Channel Flow 模式的重要依据.基于GPS 数据与龙门山地区的陡峻地形,Channel Flow 模式提出龙门山断裂上盘受到了下地壳流的向上挤压,龙门山及其邻近地区的隆升是分布式连续变形的结果[17].2008 年发生龙门山构造带上的MS8.0 汶川地震则表明青藏高原东缘的边界能够积累大量的应力并产生大地震,而非Channel Flow 所说的分布式连续变形.大量地球物理勘探资料提出在龙门山构造带深部约20km 处可能存在一个极薄的高导低阻低速层,被夹在相对强硬的层中间,成为大地震的解耦带[8-11].BaiD.H.等[12]所完成的4 条横跨青藏高原东缘、东南缘的电性剖面显示,所谓低电阻率的软弱下地壳,在羌塘块体和拉萨块体可能是存在的,但不存在于巴颜喀拉块体深部.在该块体仅上下地壳之间存在一薄层的低电阻层,其深度与分布与上述低阻低速层对应.

由此,青藏高原地壳结构之争集中到了低阻低速层与中下地壳流体的作用上.对青藏高原东缘相关计算模拟来说,两种隆升模式中地壳结构的差异反映出构造驱动力的差异,即该地区的变形是主要受控于低阻低速层的滑脱(解耦)作用还是主要受控于中下地壳流体的向上隆升作用.这种驱动力的争辩还直接关系到2008年汶川地震发震构造、龙门山断裂的几何形态、川西高原的隆升变形等重要问题.前人在关于青藏高原东缘相关计算模拟中,大部分把下地壳向上挤压的力作为汶川地震主要孕震机制[13],以满足地表破裂带调查所得的龙门山断裂的高倾角[14].即使有部分考虑到了低阻低速层的存在,也多将其考虑为受摩擦作用影响的断面[15],或是减弱上地壳杨氏模量的存在[16],或是存在于均质上地壳内的软弱层[17].尽管这些计算都取得了能够与同震GPS数据或者是地表露头对应的结果,但是这些计算模型均以中下地壳流向上挤压力作为主要驱动力,仅仅将低阻低速层的存在作为一个补充.因此,这些计算从本质上来说还是基于Channel Flow模式的讨论.虽然有了大量的地球物理勘探资料和地质资料来证实侧向逃逸模式在青藏高原东缘的可能性,但这些讨论仅制止步于地质学家们的猜想,尚未有基于侧向逃逸模式的有限元计算来讨论青藏高原东缘隆升机制和断裂活动.

Channel Flow 模式中,下地壳流变形驱动刚性上地壳,发生分布式连续变形,并不重视脆性断裂之间的相互作用.侧向逃逸模式中,刚性块体沿着边界断裂挤出的过程中,大型断裂的活动可以改变刚性块体中的应力,进而影响到块体周边或是块体内部其他断裂的活动.因此,探讨低阻低速层控制下的刚性块体变形时,仅关注一条断层是不够的.然而,目前对青藏高原东缘隆升机制的研究多集中于单条断裂,汶川地震相关数值计算也仅针对龙门山断裂的活动[1315-17],将龙门山断裂与其他断裂孤立,很少考虑到多条断裂的联动.

由此,本文建立一个西起阿坝,东至简阳,垂直于巴颜喀拉块体南东边界的青藏高原东缘上地壳的三维构造模型.模型包含了两条互相平行的断裂:以逆冲为主,含有少量走滑分量的龙门山断裂带中段,以及以走滑为主的龙日坝断裂北段.在深度上,模型仅限于上地壳,通过设置厚仅2km 的软弱薄层来模拟中上地壳之间的低阻低速层,以此来控制这两条目标断层的活动.利用澳大利亚昆士兰大学地球科学计算中心(ESSCC,Earth Systems Science Computational Centre)邢会林教授研发的基于R-minimum 的非线性摩擦接触有限元程序,计算在较短的时间内(如一个地震复发间隔)龙门山断裂和龙日坝断裂的接触行为,以及相关的刚性上地壳变形特征,论证刚性块体挤出在青藏高原东缘的可能性,并探讨块体边界断裂与块体内部断裂的活动差异,及其在块体挤出中所起的作用,为侧向逃逸模式提供了重要参考,对汶川地震孕震机制、发震构造研究也具有重要意义.

2 区域构造背景

本文研究区集中在地形最为陡峻的龙门山地区,横跨巴颜喀拉块体和华南块体(四川盆地),包含了块体边界断裂---龙门山断裂,以及块体内部走滑断裂---龙日坝断裂带,如图 1a所示.本文模型则位于该地区断层活动研究程度较高的龙门山断裂中段,即自造山带向盆地方向看,垂直于断层走向,两条目标断层近于平行,且边界GPS速度变化不大的区域(图 1b).

图 1 青藏高原东缘构造简图(a, b)(视图方向为自造山带向前陆盆地),深地震反射剖面(c)(剖面位置与图a中模型位置近似,图中数字为P波速度,单位km/S)「22,25],及研究区位置图(d) Fig. 1 Simplified geology map in the eastern margin of Tibetan Plate(a?b) (view direction is from the organ to foreland basin),deep seismic profile(c) (profile location is as the model location in Fig.1a) ,and the location of study area (d)

巴颜喀拉地块是青藏高原中东部被NWW 向东昆仑断裂带和甘孜-玉树-鲜水河断裂带、NE向龙门山断裂带中-南段和近南北向岷江断裂带围限的长条状构造域(图 1a).该块体是中国主要地震活动块体之一,1997 年西藏玛尼MS7.5级地震、2001年青海昆仑山MS8.1 级地震、2008 年3 月新疆于田MS7.3级地震和5月四川汶川MS8.0级地震,及2010 年青海玉树MS7.1 级地震均相继发生于巴颜喀拉块体四周边界活动断裂带上[18].

龙门山断裂带为巴颜喀拉块体向南东滑动与华南块体发生碰撞而形成,代表了青藏亚板块与华南板块的碰撞边界,以其陡峻的地形而著名,高地形的川西高原与低地形的四川盆地之间地表高差可达4.5~5.0km[19].斜滑逆冲型的汶川地震表明龙门山断裂带现今和第四纪时期以地壳缩短为主[20],且具有低活动速率、长复发间隔的特征[21-23].

龙日坝断裂带将巴颜喀拉块东段体划分为西部阿坝和东部龙门山两个次级块体.龙日坝断裂北东段由北支龙日曲断裂和南支毛尔盖断裂组成,两支均为与龙门山断裂带中段近于平行的NE 向断裂,相距约30km.其中龙日曲断裂为带有逆冲分量的右旋走滑断裂,毛尔盖断裂为一条纯右旋走滑断裂.GPS测量显示龙日坝断裂变形速率达4~6mm/a[22],断错地貌研究显示该断裂带北东段晚更新世以来平均右旋滑动速率为(5.4±2.0)mm/a, 垂直滑动速率约0.7 mm/a, 地壳缩短率约0.55 mm/a[24],表明龙日坝断裂带的活动以剪切变形为主,垂直滑动仅占很小的一部分.

横跨龙门山地区的深地震反射剖面(图 1c)显示,在巴颜喀拉块体东段,深度约20km 处存在一低阻低速层,为一个相对柔软的层,且该层分布于龙门山断裂带和龙日坝断裂带之下,止步于龙门山断裂带[25],即所谓的滑脱层(解耦面).除了该二维剖面外,研究区域的P波三维速度结构[26]同样显示了龙门山断裂中段地区北西侧的上地壳具有低速特征,S波速度结构也显示该地区在14~50km 深度范围内存在S波低速异常[11].而从断裂的地表分布来看(图 1a),新生代强烈隆生的岷江断裂成为了分隔龙门山断裂带中段与北段的边界,以岷江断裂为界(即高川阶区),地表破裂带显示其两侧的断裂活动分别以逆冲为主和以走滑为主[14].地形上来看,岷江断裂与龙门山断裂带中段以西均为高地形,成为了川西高原的北东边界,显示该两组断裂共同构成了隆升边界.由此我们可以推测,若低阻低速层为控制该地区上地壳内断层逆冲活动和地形隆升的关键因素,则至少在龙门山断裂中段以西,岷江断裂以南的区域有所分布.因此,本文模型(图 1b中紫色区域)正处于该低阻低速层所覆盖的范围.

3 摩擦接触有限单元数值模拟方法

大部分的地震是变形块体之间的断层上断裂摩擦失稳的结果,地震的孕育、发生过程也就是断层闭锁、滑动过程.Xing和Makinouchi[27-28]应用统一的数学公式表述了速率相关的摩擦接触中黏着(sticking)和滑移(sliding)这两种不同的运动状态,有限元计算中采用“node-to-point"接触算法和静力显示的时间积分方法,基于R最小策略控制时间步长以捕捉上述黏着和滑移运动状态的变化(即地震)并保持每一时间步内力学状态变化稳定,从而保证有限元计算过程平稳、收敛.该程序能够计算在外边界载荷下,可变形模型内部的接触面上摩擦状态黏着和滑移两种过程的变化,即断层上闭锁和解锁两种不同运动状态的变化特征,也就是地震的孕育、发生的过程,已经在变形块体间的断层非线性摩擦行为、单断层以及多断层破裂过程、地震复发周期的模拟中取得了良好的结果,并已实际应用在澳大利亚南部和南加州断裂行为的模拟[29-30].这里我们用该程序来模拟在现今的GPS速度场下、一个地震复发周期内、上地壳在低阻低速层的控制下的块体变形速度,以及龙日坝断裂与龙门山断裂在低阻低速层控制下的接触行为和地震活动,用以分析这两条断层的相互作用,讨论模拟所得的块体变形特征和青藏高原东南缘隆升之间的联系.

4 模型设置

龙门山断裂带距其西北侧NE 走向的龙日坝断裂带约180km, 距其东南侧NE 走向的龙泉山背斜约100km.鉴于龙泉山背斜是川西前陆盆地沉积边缘,代表了迄今造山带对前陆的最大影响范围.由此,本文建立横跨龙日坝断裂带与龙门山断裂带中段的三维模型,包括了阿坝次级块体、龙门山次级块体和四川盆地的一部分,模型总长度为350km, 其中龙日坝断裂上盘地表长度60km, 两断裂之间地表距离180km, 龙门山断裂下盘地表长度110km.鉴于汶川地震震源深度为19±5km, 深地震反射剖面显示的低阻低速层位置约为20±5km, 模型高度设置为26km, 其中断块高度24km, 如图 2所示.

图 2 几何模型及边界条件设置 为所有方向均固定,(7n为1与Z方向固定,为Y与Z方向固定,为自由边界,Y方向固定,而在X与Z方向分别以恒定的速率6X10-4Rd施加位移边界条件,Ref为相对速度. Fig. 2 Models and boundary conditions Here Uxyz - fixed at all the directions; Uxz - fixed along the X and Z directions; Uyz - fixed along the Y and Z directions; Ur - free along the all direction; DUxz - displacement applied in the positive X direction and Z direction at the constant velocity of6X 10-4Vref but fixed in the Y direction, Vref is the reference velocity .for the entire process.

北川-映秀断裂的破裂面几何形态一直存在争议.汶川震源机制解则显示最佳双力偶解的节面I倾角可能为33°、59°、32°、39°,与之对应的节面II倾角为70°、47°、63°、57°[31-32].USGS地震记录显示,汶川地震起始破裂处距北川-映秀地表破裂带垂向距离约17km, 而地震震源深度约19km, 若简化断面为平直的面,则断裂倾角为48°,重新定位后汶川地震余震线性密集带包络线推测的逆断层倾角约为47°[20].考虑到在静水孔隙压力条件下,45°是高倾角逆冲断裂的倾角下限,且尚处于汶川地震震源机制解的倾角范围,本文用45°倾角的平直断面做为北川-映秀断裂的主要几何形态,并用与之相切的曲面将断面与深部的滑脱构造连接起来(图 2),连接位置为19km 深处,即汶川地震震源深度.野外地震探测和人工地震勘探资料均显示在同样的应力作用下,较小倾角的断层具有倾滑的性质,高倾角的断层则走滑分量较大.为了表现出龙日坝断裂带的纯走滑特征,我们将其简化为倾角90°的平直断层,使其无法产生倾向上的滑动分解.GPS数据显示龙日坝断裂有效降低了龙门山次级块体的水平位移速率,表明该断裂并非一条浅薄的次级断裂,然而龙日坝断裂作为板块内部断裂,其深度不足以切入中下地壳,因此本文模型中将龙日坝断裂的深度设置为24km, 即延伸至低阻低速层.

鉴于龙日坝断裂带长度约为80km, 汶川地震北川-映秀破裂带龙门山镇-清平段长度约62km, 而汉旺-白鹿地表破裂带长度约72km[20],本文模型宽度设置为74km, 其中断裂宽度则为64km.有限元建模及计算结果后处理利用MSC.Patran商业软件进行,采用Paver网格对模型进行划分,网格长度为2km, 模型中有8节点六面体单元92093个,节点数为103094个.

虽然青藏高原内部具有弱负均衡重力异常,而龙门山地区为正均衡重力异常,四川盆地则为负均衡重力异常区[33-34],这种重力不均衡的分布表现了龙门山地区的高地形,以及对应的地壳厚度自龙门山地区向四川盆地减薄,但本文主要模拟对象为深度仅20余公里的上地壳的弹性变形,主要模拟目标是验证短时间周期内侧向挤压条件下快速隆升的可能性,重力不均衡对一个地震周期内块体隆升的影响是有限的.因此本文计算中初始地形为水平面,不考虑重力作用.上地壳材料假定为线弹性,密度取为2.60g/cm3,杨氏模量为44.80GPa, 泊松比为0.22(印支期花岗闪长岩,青海共和龙羊峡水电站岩石力学测试)[35].

接触面上的摩擦依速度而定:μ=0.60+(0.01-0.025)ln(V/Vref),μ 为摩擦系数,VVref分别为滑动速度及相对速度.GPS速度场[22]显示青藏高原东缘巴颜喀拉块体的水平变形速率达到4.2mm/a, 且与龙日坝断裂呈几乎45°的夹角,由此本文将此斜向的挤压分解为X方向和Z方向以相同的恒定速率挤压,挤压速度均为6×10-4Vref(≈6mm/a).深地震反射剖面显示低阻低速层在巴颜喀拉块体深处是广泛存在的,是一个相对软弱的薄层,因此我们在巴颜喀拉块体深部设置一个2km 厚材料软弱的薄层,而在四川盆地深部则设置2km 厚的强硬薄层,两个薄层的边界条件均设置为Y方向固定,X方向自由.这样不仅可以使得块体沿着该软弱薄层发生滑动,并止步于强硬薄层,避免造成畸变或是块体扭曲.其他模型边界设置如图 2所示.

5 模拟结果与讨论

在上述边界条件的作用下,本模型中的龙门山断裂和龙日坝断裂均发生了破裂,其中龙日坝断裂表现为纯走滑变形,而龙门山断裂带则以逆冲变形为主,兼有极其少量的走滑变形,且破裂时间晚于龙日坝断裂的破裂时间,两者之间的比值约为7∶11,表明这两条断裂的复发间隔不会相差一个数量级.

在整个应力累积直至断层破裂的过程中,X轴方向(垂直于断裂走向,即北西-南东方向)的变形速度在模型中的分布基本不变,且表现为自北西向南东减小:在龙日坝断裂西北侧,由于距离模型边界最近,变形速度为快速减小,自龙日坝断裂开始缓慢减速,至龙门山断裂东南侧(240km),X轴方向的变形速度极小,仅相当于龙日坝断裂东南侧的变形速度的7%左右,且缓慢趋近于0.地表节点(剖面I)的X轴方向变形速度曲线(图 3)显示,当龙日坝断裂发生破裂时(图 3b),X轴方向的变形速度在断裂附近约8km 的范围内缓慢减小,而当龙门山断裂发生破裂时(图 3d),该断裂成为了变形速度场的间断面,虽然此时地表节点变形速度曲线自北西向南东减小的趋势没有发生改变,但在龙门山断裂两盘出现了间断,断裂上盘相对变形速度为(2.90×10-5)/6×Vrefmm/a, 下盘为(1.00×10-5)/6×Vrefmm/a, 而该处也是断块下部软弱层的终点.这表明X轴方向的变形受到断层破裂的一些影响.龙门山断裂作为块体下部软弱层的终点,X轴方向的变形在传递到龙门山断裂时受阻,转化为应力积累于断层面上并使之发生破裂.而龙日坝断裂作为一条位于软弱层之上的走滑断层,对X轴方向的变形不仅没有起到削弱作用,反而传递了X轴方向的变形.

图 3 X轴方向变形速度场(曲线剖面位置如图 2中剖面I所示) 加载时间T=3×106 R/Vref, 其中(a)加载时间内步长的累积R=1.00×10-7 ;b)R=6.97×10-6 ;c)R=9.90×10-6 ;d)R=1.16×10-5. Fig. 3 Deformation velocity along X coordinate (the curve of deformation velocity is the surface nodes in the profile I as shown in Fig.2) Loading time=3 × 106R/Vref refers to the accumulated time step in the loading time, in which (a)R=1.00×10-7 ;(b)R=6.97×10-6 ;(c)R=9.90×10-6 ;(d)R=1.16×10-5.

Y轴方向(即垂直方向)的变形速度场(图 4)显示了块体的隆升速率在持续的侧向挤压过程中的变化:在应力开始加载时,初始速度场显示龙日坝断裂西北侧是隆升速度最快的地方,龙日坝断裂与龙门山断裂之间的块体隆升速度有所减慢,龙门山断裂下盘隆升速度极小.在龙日坝断裂发生走滑破裂之后(图 4b),剖面I显示龙日坝断裂北西侧约4km的范围内出现隆升速度的低值,而断裂南东侧则出现了隆升速度的较高值,该低值是断层右旋走滑造成的隆升中心的转移,即断裂西北侧的隆升中心向北东迁移,断裂南东侧的隆升中心则向南西迁移.在龙门山断裂发生破裂之前夕(图 4c),龙日坝断裂附近的隆升速度保持恒定而断裂造成的隆升中心偏移更大,龙门山断裂上盘的隆升速度迅速增大,几乎与龙日坝断裂南东侧近似,而龙门山断裂下盘的隆升速度保持不变.龙门山断裂发生破裂时(图 4d),龙日坝断裂附近的隆升速度基本不变,而龙门山断裂带上盘的隆升速度远远高于龙日坝断裂东南侧的变形速度.龙门山断裂的破裂是从断裂的中部开始(图 4c),向两端扩展,且垂直方向(Y轴)的变形速度自南西向北东逐渐变小(图 4d).与龙门山断裂上盘成鲜明对比的是,断裂下盘则仍然保持了相当低的隆升速度.此外,在整个过程中,由于倾角为90°,龙日坝断裂始终没有发生具有逆冲性质的破裂,仅有走滑性质的破裂,这与模型建立时的设想一致.

图 4 Y轴方向变形速度场(曲线剖面位置如图 2中剖面I所示) 加载时间T=3×106R/Vref, 其中(a)R=1.00×10-7 ;(b)R=6.97×10-6 ;(c)R=9.90×10-6 ;(d)R=1.16×10-5 Fig. 4 Deformation velocity along X coordinate (the curve of deformation velocity is the surface nodes in the protile I as shown in Fig. 2) Loading time T=3 × 106R/Vref, in which: (a)R=1.00×10-7;(b)R=6.97×10-6;(c)R=9.90×10-6;(d)R=1.16×10-5.

在龙门山断裂带发生破裂时Z轴方向的变形(如图 5所示),即走滑变形量,集中在龙日坝断裂附近,并经过该断裂之后迅速减小,但龙门山断裂带仍有少量的右旋走滑变形量.剖面I和剖面II显示了模型边缘部分和模型中部的地表节点Z轴的变形速度:剖面I中由于边界条件的影响,龙日坝断裂两侧虽然有较快的变形速度差异,但尚未破裂,龙门山断裂上盘近断层沿Z轴的变形速度分别为(-1.65×10-5)/6×Vrefmm/a, 下盘为(0.00×10-5)/6×Vref;剖面II显示大部分的走滑变形发生在龙日坝断裂上(距离60km),断层两侧相对变形速度分别为(7.2 ×10-5)/6×Vrefmm/a和(-9.60×10-5)/6×Vref, 而龙门山断裂带上(距离240km)两侧相对变形速度分别是(-1.70×10-5)/6×Vref和(0.00×10-5)/6×Vref, 仅约为龙日坝断裂上变形量的10%.

图 5 龙门山断裂破裂时Z轴方向变形速度场 加载时间T =3×106R/Vref, 其中R=1.16×10-5. Fig. 5 Deformation velocity along Z coordinate when the Longmenshan thrust fault slipping Loading time T=3 × 106R/Vref, in which R=1.16×10-5.
6 结论与讨论

本文选择了以低阻低速层作为主要控制因素,以侧向挤压力为驱动力,模拟了青藏高原东缘块体的隆升,以及相应的龙日坝断裂与龙门山断裂的活动,取得了以下结果:

(1) 在侧向挤压力的作用下,青藏高原东缘中上地壳之间的低阻低速层作为一个相对软弱的层,能够成为块体在侧向挤压力作用下发生相对运动的滑脱层(解耦带),造成青藏高原东缘的快速隆升.(2)在侧向挤压力的持续作用下,低阻低速层边缘产生了垂直变形(隆升)速度间断和水平变形速度间断,表明低阻低速层同时控制了水平形变和垂直形变的分布范围,而这种变形速度的间断使得低阻低速层边界上的龙门山断裂成为应力累积的场所.(3)位于低阻低速层之上的龙日坝断裂产生了剪切破裂,并以此分解了绝大部分施加于边界上沿Z轴方向的滑移量,导致其南侧的龙门山次级块体上仅保留了较小走滑分量,而龙门山断裂的活动也相应地表现为逆冲为主,兼有少量的走滑分量.(4)在龙日坝断裂的剪切变形和边界的侧向挤压力双重作用下,龙门山断裂由断裂的中部开始破裂,同时向两端扩展.龙门山断裂上盘近断层的垂直变形速度远远大于走滑变形速度,且垂直变形速度向具有沿走向向北东方向减小的趋势,这与长周期地震台垂直向记录的P波波形资料反演得出汶川地震的时空破裂过程一致[32].

结合地质背景和断裂特性,我们取得了一些青藏高原东缘隆升相关的断裂活动、块体变形特征的认识:

(1) 断裂几何形态和相对于低阻低速层的位置变化造成了不同的断裂活动特征.沿X轴方向上的变形速度(水平变形)显示,龙日坝断裂的纯剪变形并没有阻止水平变形向南东的扩展,相反,龙日坝断裂被两个刚性的次级块体挟持着一起向南东运动.地表的水平变形在龙门山断裂上盘呈自北西向南东缓慢下降,在龙门山断裂下盘则变形极小.龙门山断裂位于低阻低速层边缘,受该层影响,断裂下盘的水平变形受到下盘强硬层不易变形的制约,而断裂上盘变形则受益于块体底部柔软的低阻低速层的较大较快的变形.这两种变形的差异导致了块体变形的差异,使得应力积累在龙门山断裂带上,产生逆冲破裂和大地震,并使其两侧的水平变形产生更大的差异.

(2) 在持续的侧向挤压力作用下,阿坝次级块体和龙门山次级块体产生快速隆升,隆升速率自北西向南东减小,至低阻低速层边缘,即龙门山断裂,产生了隆升速度间断.龙门山断裂发生破裂时,缓慢的大尺度地壳构造运动转化为急剧的小尺度断层失稳,断裂上盘的隆升速度远远大于模型其他位置,而龙门山断裂下盘仍然保持了极低的隆升速度,使得该隆升间断进一步发育.可见在龙门山断裂大地震反复的孕育、发生过程中,这种隆升间断不断的积累,能够形成造山带与前陆盆地之间巨大的高程差.

龙日坝断裂两侧在低阻低速层的控制下也产生了快速隆升(其数值大小受边界条件影响较大),而龙日坝断裂的走滑使其两侧的隆升中心发生偏移:在龙日坝断裂西北侧,隆升分布较为均匀,且由于走滑作用向模型北边界聚集,这与以逆冲活动为主的龙日曲断裂和毛儿盖断裂[24]相对应;龙日坝断裂南东侧的隆升速度集中于断裂中部,且远小于大地震发生时龙门山断裂上盘隆升量,这与龙门山地区的局部高程差分布[35]类似,显示造山带和块体内部变形的差异.

(3) 位于低阻低速层之上的龙日坝断裂分解了绝大部分的剪切量,使得龙门山断裂的活动仅兼有少量走滑分量.多期GPS复测显示,经龙日坝断裂,巴颜喀拉地块的SE 向运动速率减少约3mm/a[22],而经龙门山断裂,SE 向运动速率减少得极其有限(图 1).巴颜喀拉地块东段相邻两个次级块体向南东滑移速率的减少量应大部分转化为龙日坝断裂带的右旋走滑和逆冲运动,剩余部分可能被龙门山次级块体的构造变形吸收[24].这与我们的计算结果一致.

在低阻低速层的控制下,次级块体边界的走滑断裂对块体边界断裂的滑动分解作用值得我们深思.GPS观测数据显示巴颜喀拉块体东段存在着左旋走滑速率,而龙门山断裂的活动表现为逆冲为主,兼有少量走滑,因此必然需要通过块体内部的变形来分解绝大部分的走滑量.从本文的计算结果来看,巴颜喀拉块体东段的走滑变形基本集中在龙日坝断裂带上,并不存在宽阔广泛的左旋剪切带[7].也就是说,在低阻低速层的控制下,通过刚性的阿坝次级块体与刚性的龙门山次级块体沿着龙日坝断裂的滑移,可以分解绝大部分施加于巴颜喀拉块体之上的左旋走滑量.

侧向逃逸模式强调大型断裂而忽略小断裂上可能发生的变形,但是该模式提出的块体边界走滑速率远大于现今GPS观测所得的断裂滑动速率[736],也高于晚第四纪以来的断裂活动速率[37].这里面固然有长期构造活动中断裂活动速率变化的因素,也可能有较大型的次级断裂分担了一部分走滑量的因素.至少对青藏高原东缘来说,估算块体挤出量时龙日坝断裂的活动不可忽视.此外,在侧向挤压力方向与断裂走向呈一定角度的情况下,较大规模的走滑断层能起到调整方向的作用,这些调节作用有没有分担一部分的挤出量也是值得考虑的.但是,何种规模的走滑断裂才能分担挤出量,该类走滑断裂在青藏高原的分布及其几何特征、运动特征尚需要进一步的研究,这对于进一步论证侧向逃逸模式具有极其重要的意义.

(4) 龙门山断裂的破裂时间晚于龙日坝断裂的破裂时间,两者之间的比值约为7∶11,表明这两条断裂的大地震复发间隔不会相差一个数量级.探槽古地震研究显示,龙门山断裂带的汶川地震震级的地震复发间隔为2300~3300a[23],由此推断龙日坝断裂的大地震复发间隔约为1463~2100a, 但尚需要古地震探测工作进行验证.

本文三维线弹性有限元模型仍存在以下的局限性:(1)没有考虑重力作用.重力在长时间尺度内的横向作用对于断层深部力学平衡具有很大影响,但是对短时间尺度的计算来说,在没有引入黏弹性介质性质的情况下,重力作用在弹性介质系统内部只是表现为垂向作用及侧压强.在弹性模型阶段,重力作用会对速度场和位移场量值有一定影响,但是并不影响断面接触点滑动的相对时间,以及速度场时空变化的性质.(2)本模型没有包含黏弹性的下地壳或者下地壳流的作用.为了突出上地壳的脆性断层活动,本文采用了弹性介质,将低阻低速层设置为一极薄层,并将其下边界设置为X方向上(NW-SE 方向)自由而Y方向上(垂向)固定,以保证低阻低速层能够传递X方向上的变形,而本身不会产生较大的垂向变形.这种模型设置仅计算了低阻低速层在侧向挤压力作用下对块体变形和断层活动的控制作用,验证了在缺乏下地壳流的条件下青藏高原东缘脆性上地壳隆升的可能性.如果下地壳流能产生如Channel Flow 假设中的作用,那么下地壳流至少能够对龙门山次级块体中上地壳产生向上的挤压力.在一个地震周期的短时间尺度内,下地壳流对上地壳的作用方式、挤压力的大小、与初始应力场的关系等都是未知的.低阻低速层与Channel Flow 是否是非此即彼,或是两者可以共存,以及侧向挤压力和下地壳流的向上两者在青藏高原隆升中的作用孰轻孰重,尚需要进一步工作进行研究.

尽管存在以上局限性,但在地震周期时间尺度内,对于现今青藏高原南缘上地壳变形来说,本文的模型基本能够反映低阻低速层控制下的龙日坝断裂和龙门山断裂中段的断裂活动和刚性块体变形特征.

致谢

本研究是在澳大利亚昆士兰大学地球系统科学计算中心的超级计算机上利用该中心邢会林研发的ESyS_Crustal软件进行的,计算中得到工作人员的大力协助,在此表示感谢.感谢各位评审老师,对文章的建议和意见.

参考文献
[1] Clark M K, Royden L H. Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. Geology , 2000, 28(8): 703-706. DOI:10.1130/0091-7613(2000)28<703:TOBTEM>2.0.CO;2
[2] Beaumont C, Jamieson R A, Nguyen M H, et al. Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation. Nature , 2001, 414(6865): 738-742. DOI:10.1038/414738a
[3] Tapponnier P, Peltzer G, Armijo R. Continent-Continent Collision: Himalayan-Alpine Belt on the mechanics of the collision between India and Asia. Geological Society, London, Special Publications , 1986, 19: 113-157. DOI:10.1144/GSL.SP.1986.019.01.07
[4] Tapponnier P, Xu Z Q, Roger F, et al. Oblique stepwise rise and growth of the Tibet Plateau. Science , 2001, 23(294): 1671-1677.
[5] Shapiro N M, Ritzwoller M H, Molnar P, et al. Thinning and flow of Tibetan crust constrained by seismic anisotropy. Science , 2004, 305(5681): 233-236. DOI:10.1126/science.1098276
[6] Yao H J, Beghein C, van der Hilst R D. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—II. Crustal and upper-mantle structure. Geophysical Journal International , 2008, 173(1): 205-219. DOI:10.1111/gji.2008.173.issue-1
[7] 张培震, 王敏, 甘卫军, 等. GPS观测的活动断裂滑动速率及其对现今大陆动力作用的制约. 地学前缘 , 2003, 10(Suppl.): 81–92. Zhang P Z, Wang M, Gan W J, et al. Slip rates along major active faults from GPS measurements and constrains on contemporary continental tectonics. Earth Science Frontiers (in Chinese) , 2003, 10(Suppl.): 81-92.
[8] 王椿镛, MooneyW D, 王溪莉, 等. 川滇地区地壳上地幔三维速度结构研究. 地震学报 , 2002, 24(1): 1–16. Wang C Y, Mooney W D, Wang X L, et al. Study on 3-D velocity structure of crust and upper mantle in Sichuan-Yunnan region, China. Acta Seismologica Sinica (in Chinese) , 2002, 24(1): 1-16.
[9] 王椿镛, 常利军, 吕智勇, 等. 青藏高原东部上地幔各向异性及相关的壳幔耦合型式. 中国科学D辑 , 2007, 37(4): 495–503. Wang C Y, Chang L J, Lü Z Y, et al. Seismic anisotropy of upper mantle in eastern Tibetan Plateau and related crust-mantle coupling pattern. Science in China (Series D) (in Chinese) , 2007, 37(4): 495-503.
[10] 蔡学林, 曹家敏, 朱介寿, 等. 龙门山岩石圈地壳三维结构及汶川大地震成因浅析. 成都理工大学学报: 自然科学版 , 2008, 35(4): 357–365. Cai X L, Cao J M, Zhu J S, et al. A preliminary study on the 3-D crust structure for the Longmen lithosphere and the genesis of the huge Wenchuan earthquake, Sichuan, China. Journal of Chengdu University of Technology (Science & Technology Edition) (in Chinese) , 2008, 35(4): 357-365.
[11] 刘启元, 李昱, 陈九辉, 等. 汶川Ms8.0地震: 地壳上地幔S波速度结构的初步研究. 地球物理学报 , 2009, 52(2): 309–319. Liu Q Y, Li Y, Chen J H, et al. Wenchuan Ms8.0 earthquake: preliminary study of the S-wave velocity structure of the crust and upper mantle. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 309-319.
[12] Bai D H, Unsworth M J, Meju M A, et al. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging. Nature Geoscience , 2010, 3(5): 358-362. DOI:10.1038/ngeo830
[13] 朱守彪, 张培震. 2008年汶川MS8.0地震发生过程的动力学机制研究. 地球物理学报 , 2009, 52(2): 418–427. Zhu S B, Zhang P Z. A study on the dynamical mechanisms of the Wenchuan MS8.0 earthquake, 2008. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 418-427.
[14] Xu X W, Wen X Z, Yu G H, et al. Coseismic reverse- and oblique-slip surface faulting generated by the 2008 Mw7.9 Wenchuan earthquake, China. Geology , 2009, 37(6): 515-518. DOI:10.1130/G25462A.1
[15] 张竹琪, 张培震, 王庆良. 龙门山高倾角逆断层结构与孕震机制. 地球物理学报 , 2010, 53(9): 2068–2082. Zhang Z Q, Zhang P Z, Wang Q L. The structure and seismogenic mechanism of Longmenshan high dip-angle reverse fault. Chinese J. Geophys. (in Chinese) , 2010, 53(9): 2068-2082.
[16] 陈祖安, 林邦慧, 白武明, 等. 2008年汶川8.0级地震孕震机理研究. 地球物理学报 , 2009, 52(2): 408–417. Chen Z A, Lin B H, Bai W M, et al. The mechanism of generation of May 12, 2008 MS8.0 Wenchuan earthquake. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 408-417.
[17] 陶玮, 胡才博, 万永革, 等. 铲形逆冲断层地震破裂动力学模型及其在汶川地震研究中的启示. 地球物理学报 , 2011, 54(5): 1260–1269. Tao W, Hu C B, Wan Y G, et al. Dynamic modeling of thrust earthquake on listric fault and its inference to study of Wenchuan earthquake. Chinese J. Geophys. (in Chinese) , 2011, 54(5): 1260-1269.
[18] 邓起东, 高翔, 陈桂华, 等. 青藏高原昆仑—汶川地震系列与巴颜喀喇断块的最新活动. 地学前缘 , 2010, 17(5): 163–178. Deng Q D, Gao X, Chen G H, et al. Recent tectonic activity of Bayankala fault-block and the Kunlun-Wenchuan earthquake series of the Tibetan Plateau. Earth Science Frontiers (in Chinese) , 2010, 17(5): 163-178.
[19] 贾秋鹏, 贾东, 朱艾斓, 等. 青藏高原东缘龙门山冲断带与四川盆地的现今构造表现: 数字地形和地震活动证据. 地质科学 , 2007, 42(1): 31–44. Jia Q P, Jia D, Zhu A L, et al. Active tectonics in the Longmen thrust belt to the eastern Qinghai-Tibetan Plateau and Sichuan basin: evidence from topography and seismicity. Chinese Journal of Geology (in Chinese) , 2007, 42(1): 31-44.
[20] 徐锡伟, 闻学泽, 叶建青, 等. 汶川MS8.0地震地表破裂带及其发震构造. 地震地质 , 2008, 30(3): 597–629. Xu X W, Wen X Z, Ye J Q, et al. The MS8.0 Wenchuan earthquake surface ruptures and its seismogenic structure. Seismology and Geology (in Chinese) , 2008, 30(3): 597-629.
[21] Zhang P Z, Shen Z K, Wang M, et al. Continuous deformation of the Tibetan Plateau from global positioning system data. Geology , 2004, 32(9): 809-812. DOI:10.1130/G20554.1
[22] Shen Z K, Lü J N, Wang M, et al. ,Contemporary crustal deformation around the southeast borderland of the Tibetan Plateau. J. Geophys. Res. , 2005, 110: B11409.
[23] Ran Y K, Chen L C, Chen J, et al. Paleoseismic evidence and repeat time of large earthquakes at three sites along the Longmenshan fault zone. Tectonophysics , 2010, 491(1-4): 141-153. DOI:10.1016/j.tecto.2010.01.009
[24] 徐锡伟, 闻学泽, 陈桂华, 等. 巴颜喀拉地块东部龙日坝断裂带的发现及其大地构造意义. 中国科学 (D辑: 地球科学) , 2008, 38(5): 529–542. Xu X W, Wen X Z, Chen G H, et al. Discovery of the Longriba faults, eastern part of the Bayankela tectonic block and its geodynamic implications. Sciences in China (Series D) (in Chinese) , 2008, 38(5): 529-542.
[25] Li Q S, Gao R, Wang H Y, et al. Deep background of Wenchuan Earthquake and the upper crust structure beneath the Longmen Shan and adjacent areas. Acta Geologica Sinica , 2009, 83(4): 733-739. DOI:10.1111/j.1755-6724.2009.00096.x
[26] 吴建平, 黄媛, 张天中, 等. 汶川Ms8.0级地震余震分布及周边区域P波三维速度结构研究. 地球物理学报 , 2009, 52(2): 320–328. Wu J P, Huang Y, Zhang T Z, et al. Aftershock distribution of the Ms8.0 Wenchuan earthquake and three dimensional P-wave velocity structure in and around source region. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 320-328.
[27] Xing H L, Makinouchi A. Finite-element modeling of multibody contact and its application to active faults. Concurrency and Computation: Practice and Experience , 2002, 14(6-7): 431-450. DOI:10.1002/(ISSN)1532-0634
[28] Xing H L, Makinouchi A. Finite element modelling of frictional instability between deformable rocks. Int. J. Numer. Anal. Meth. Geomech. , 2003, 27(12): 1005-1025. DOI:10.1002/(ISSN)1096-9853
[29] Xing H L, Makinouchi A, Mora P. Finite element modeling of interacting fault systems. Physics of the Earth and Planetary Interiors , 2007, 163(1-4): 106-121. DOI:10.1016/j.pepi.2007.05.006
[30] Xing H L, Mora P. Construction of an intraplate fault system model of South Australia, and simulation tool for the iSERVO institute seed project. Pure and Applied Geophysics , 2006, 163(11-12): 2297-2316. DOI:10.1007/s00024-006-0127-x
[31] 刘超, 张勇, 许力生, 等. 一种矩张量反演新方法及其对2008 年汶川MS8.0地震序列的应用. 地震学报 , 2008, 30(4): 329–339. Liu C, Zhang Y, Xu L S, et al. A new technique for moment tensor inversion with applications to the 2008 Wenchuan MS8.0 earthquake sequence. Acta Seismologica Sinica (in Chinese) , 2008, 30(4): 329-339.
[32] 张勇, 许力生, 陈运泰. 2008年汶川大地震震源机制的时空变化. 地球物理学报 , 2009, 52(2): 379–389. Zhang Y, Xu L S, Chen Y T. Spatio-temporal variation of the source mechanism of the 2008 great Wenchuan earthquake. Chinese J. Gephys. (in Chinese) , 2009, 52(2): 379-389. DOI:10.1002/cjg2.v52.2
[33] 李勇, 徐公达, 周荣军, 等. 龙门山均衡重力异常及其对青藏高原东缘山脉地壳隆升的约束. 地质通报 , 2005, 24(12): 1162–1168. Li Y, Xu G D, Zhou R J, et al. Isostatic gravity anomalies in the Longmen Mountains and their constraints on the crustal uplift below the mountains on the eastern margin of the Qinghai—Tibet Plateau. Geological Bulletin of China (in Chinese) , 2005, 24(12): 1162-1168.
[34] 王谦身, 滕吉文, 张永谦, 等. 四川中西部地区地壳结构与重力均衡. 地球物理学报 , 2009, 52(2): 579–583. Wang Q S, Teng J W, Zhang Y Q, et al. The crustal structure and gravity isostasy in the middle western Sichuan area. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 579-583.
[35] 叶金汉. 岩石力学参数手册. 北京: 水利水电出版社, 1991 . Ye J H. Manual of Rock Mass Parameters (in Chinese). Beijing: China Water Power Press, 1991 .
[36] 汪一鹏, 沈军, 王琪, 等. 川滇块体的侧向挤出问题. 地学前缘 , 2003, 10(Suppl.): 188–192. Wang Y P, Shen J, Wang Q, et al. On the lateral extrusion of Sichuan-Yunnan block (Chuandian block). Earth Science Frontiers (in Chinese) , 2003, 10(Suppl.): 188-192.
[37] 李传友. 青藏高原东北部几条主要断裂带的定量研究. 北京: 中国地震局地质研究所, 2005 . Li C Y. Quantitative studies on major active fault zones in Northeastern Qinghai-Tibet Plateau (in Chinese). Beijing: Institute of Geology, China Earthquake Administration, 2005 .