2014, Vol. 57
2. Laboratoire de Géologie, CNRS-UMR8538, école Normale Supérieure, Paris 75231;
3. 中国科学院大学, 中国科学院计算地球动力学重点实验室, 北京 100049;
4. 陕西省地震局, 西安 710068
2. Laboratoire de Géologie, CNRS-UMR8538, école Normale Supérieure, Paris 75231, France;
3. Laboratory of Computational Geodynamics, University of the Chinese Academy of Sciences, Beijing 100049, China;
4. Earthquake Administration of Shaanxi Province, Xi'an 710068, China
2008年5月12日14时28分在青藏高原东缘与四川盆地交汇处的龙门山断裂带上发生了汶川Mw7.9级强烈地震.地震造成的西南向北东方向发展的破裂带长度约350 km，破裂持续时间长达90 s，整个断层面上的平均滑动量约2.4 m，最大滑动量达7.3 m(张勇等，2008；张培震等，2008)；并且沿北川—映秀断裂和彭县—灌县断裂分别形成了长达240 km和72 km的地表破裂带，最大垂直位移量约为6.5 m，右旋走滑位移量约为4.9 m(徐锡伟等，2013).震源机制解显示，汶川地震主震以逆冲为主兼有少量右旋走滑.汶川地震较大余震的“缺失”分析认为汶川地震发生在龙门山断裂带北东段从映秀到青川之间大约 350 km的断裂上, 而留下了西南段从汶川西南到芦定之间大约 120 km 的地震亏空区(陈运泰等，2013).5年后，2013年4月20日8时2分在汶川地震的亏空区发生了芦山地震.地震破裂过程研究表明，断层滑动面长宽约为30 km×30 km，最大滑动量为1.6 m；破裂起始点接近于地震滑动量集中区，破裂面从震源处向下延伸至20 km左右深的地方，破裂面并未到达地表，破裂过程没有明显的方向性；震源机制解显示，芦山地震主要以逆冲为主兼有非常少量的右旋走滑分量(陈运泰等，2013).结合地震勘探成果所揭示的龙门山断裂带地壳结构，余震震源定位工作表明龙门山区域地震空间分布特征为：纵向上，龙门山断裂带这两次地震主震均发生在上地壳的底部，汶川地震的震源深度为19 km(USGS)，芦山地震的震源深度为14 km(USGS)，全球其他各个研究机构对芦山地震的震源深度定位结果在12~19 km深度范围之间(杜方等，2013)；两次地震的绝大部分余震均发生在上地壳范围(5~25 km)(张瑞青等，2008；黄媛等，2008；刘巧霞等，2008；陈九辉等，2009；张广伟等，2013；吕坚等，2013；陈运泰等，2013)，而在中、下地壳深度范围内鲜见余震发生；横向上，地震(Mw>3)在龙门山断裂带青藏高原一侧密集分布，且有历史大 震发生，如：1630年M6.5级虎牙地震、1913年M7.0 级叠溪地震、1932年M7.5级叠溪地震、1960年M6.7级松潘地震与1976年3次震级为6.6< M<7.3级的松潘地震群；而四川盆地地震 (Mw>3) 稀少(李勇等，2006；滕吉文等，2008；Robert et al.，2010；雷兴林等，2013).
在低地形变速率(约3 mm/a)的龙门山断裂带相隔5年发生两次强震：为探讨龙门山断裂带地震孕震成因，不同学者就青藏高原东缘的构造运动模式，提出了两种不同的概念性地震地质孕震模型.Tapponnier等(2001)主张龙门山断裂带可能贯穿青藏高原东缘整个脆性地壳；Hubbard等(2009)认为汶川地震是龙门山断裂带在构造挤压环境下的再一次活动结果，并强调地壳缩短是地震发生的首要机制.更多的研究(张培震等，2008；滕吉文等，2008； Royden et al.，2008; Burchfiel et al.，2008；Stone et al.，2009；Robert et al.，2010；柳畅等，2012b)则支持如下的观点：认为印度板块对欧亚板块的推挤作用造成了青藏高原物质的东向运动，高原柔性的中、下地壳物质在龙门山断裂带处遭到相对坚硬的四川盆地的阻挡之后，部分中、下地壳物质在龙门山断裂带下堆积、产生应力集中而导致龙门山断裂带地震的发生.这一观点所描述的青藏高原脆性上地壳与韧性的中、下地壳组合的地壳分层流变结构，区别于Tapponnier等(2001)所认为的全脆性地壳流变结构.总结以上的观点，其共同点在于认同青藏高原的东向挤压作用，但又强调了完全不同的地壳流变结构.如何对这两种概念性地震地质孕震模型从数值模拟的角度作出验证，进一步对龙门山断裂带地震发生机理加以更深层次的认识，是一个尚未解决且需要探讨的问题.为此，本文将建立四种不同流变结构的龙门山断裂带三维岩石圈模型，数值模拟不同岩石圈模型中的应力积累过程与分布；并探讨青藏高原中、下地壳不同的黏滞系数，以及青藏高原下地壳流Channel flow(Royden et al.，1997；Clark et al，2000，2005；Beaumont et al，2001，2004；石耀霖等，2008)对应力积累过程的影响.
 | 图1 青藏高原东缘与四川盆地区域构造与地表GPS水平速度场 红色圆点代表2008汶川地震Mw5.0级以上余震，粉色圆点代表2013芦山地震Mw5.0级以上的余震.黄色箭头代表地表GPS水平速度(Shen et al.，2005).黑色小三角代表地名，分别为：1, 汶川; 2, 成都; 3, 叠溪; 4, 松潘; 5, 卢霍; 6, 康定; 7, 芦山.黑色虚线框表示数值模拟的模型区域.D-S fault代表大川—双石断裂 Fig.1 Map of the eastern Tibetan Plateau and the Sichuan Basin showing the Longmenshan fault and the horizontal GPS velocity The epicentres of the main shock of 2008 Wenchuan earthquake and 2013 Lushan earthquake, and aftershocks (Mw>5.0) of the Wenchuan earthquake (red dots) and Lushan earthquake (pink dots). The yellow arrow indicates the GPS velocity(Shen et al.，2005). The black triangles indicate the locations of cities. The numbers mark the city names: 1, Wenchuan; 2, Chengdu; 3, Diexi; 4, Songpan; 5, Luhuo; 6, Kangding; 7, Lushan. The black solid lines indicate the faults. D-S fault indicates the Dachuan-Shuangshi fault. The black dashed lines forming a rectangle mark the area where stress accumulation is simulated. |
龙门山断裂带处于年轻活跃的青藏高原与古老稳定的扬子克拉通的交汇地带，为NE-SW走向上长约470 km、宽约50 km的活动断裂带.该断裂带在空间上呈NE-NNE向展布，且以北NW-SE方向逆冲为主而兼具少量右旋走滑分量；是由汶川—茂汶逆断裂、映秀—北川逆断裂、彭县—灌县断裂、龙门山山前断裂、山前隐伏断裂段和相应的推覆体组成的一组断裂系.2008年汶川地震发生在映秀—北川断裂带上，2013年芦山地震发生在山前断裂南段大川—双石断裂上.
为揭示龙门山断裂带的岩石圈结构，该地域进行过大量的地震探测工作：如层析成像工作(Huang et al.，2007；郭飚等，2009)，接收函数工作(刘启元等，2009；Robert et al.，2010；郑勇等，2013)，人工地震勘探工作(Wang et al.，2005；Wang et al.，2010)和一些地质勘查工作(Burchfiel，1995，2008).以上研究结果表明从青藏高原东缘到四川盆地地壳结构变化强烈，青藏高原一侧地壳厚度在63 km左右，而四川盆地一侧地壳厚度在45 km左右.在龙门山区域横向50 km范围内地壳厚度的垂直变化幅度可达18 km左右，地形高差在龙门山陡降近5 km；因此，结合Moho面形状与地表地形看，龙门山区域地壳结构呈瓶颈状(图2 ).
 | 图2 龙门山断裂带三维岩石圈模型 模型尺寸为500 km×100 km×100 km(X×Y×Z). Y轴方向平行于龙门山断裂带走向.不同颜色代表不同岩石圈层位的不同物质，依据对应的黏滞系数,青藏高原中、下地壳物质较软而四川盆地中、下地壳物质相对较硬.黑色星星代表2013芦山Mw6.6级地震的震中位置 Fig.2 Lithospheric model of the Longmen Shan The size of this model is 500 km×100 km×100 km (X×Y×Z). The Y axis is parallel to the strike of the Longmen Shan fault. The different colours refer to the different material parameters of layers in the lithosphere. The materials of the middle and lower crusts of the Tibetan Plateau are softer than those of the Sichuan Basin, corresponding to the level of viscosity. Black star indicates the epicentre of the 2013 Lushan earthquake. |
大陆岩石圈一般存在脆性的上地壳、柔性的下地壳和较强上地幔这种三明治式的分层流变结构，尤其在青藏高原这种地壳厚、地热流量密度高的地区会更为显著(Royden et al.，1997).岩石流变性质与温度密切相关，安美建利用地震波速估算了上地幔50~200 km的温度(安美建等，2007)；石耀霖据此温度模型给出了中国大陆岩石圈不同层位的等效黏滞系数(石耀霖等，2008)，得到了相应区域的三维流变结构.胡圣标研究了川滇地区相对偏高的平均大地热流值(胡圣标等，2001)，认为川滇地区的中、下地壳较热，介质强度可能相对较软.GPS观测结果也显示, 川滇地区的现今运动模式支持较硬的脆性上地壳和软弱的柔性中、下地壳的分层流变结构(Shen et al.，2005).
P波层析成像工作所揭示的四川盆地下方地壳速度结构，认为扬子克拉通的稳定速度结构可能延伸至250 km深度(Li et al，2006，2008).并且，中生代与新生代的构造活动对四川盆地周边构造单元产生了显著的变形作用，而四川盆地并未受到明显的影响(Li et al，2006，2008)，从而表明四川盆地下方岩石圈的力学强度较周边构造单元较强.
本文建立如图1 中黑色虚线方框所示区域的三维黏弹性Maxwell体岩石圈模型(见图2 ).坐标系X轴(SE41°)垂直于龙门山断裂带走向，长500 km；Y轴(NE49°)平行于龙门山断裂带走向，长100 km，包含了龙门山断裂带在2008汶川地震未发生破裂的西南段部分；Z轴向上(岩石圈的深度在Z轴负向)，模型中深度为从地面至地下100 km；坐标原点位于模型上表面的西南端.模型具体考虑龙门山断裂带青藏高原与四川盆地的地形高差、地壳厚度在龙门山断裂带下方的突变与中、下地壳的黏滞系数在高原与盆地下方的差异.
本文计算所涉及的黏弹性介质的3个物质参数分别为：杨氏模量E，泊松比ν和黏滞系数η.岩石杨氏模量E与泊松比ν依据波速计算.一般随深度增加，地震波速增加，岩石的杨氏模量亦随深度增加.依据所选择物质参数，龙门山断裂带岩石圈100 km 深度范围内介质杨氏模量随深度的变化曲线见图3 .
 | 图3 龙门山断裂带岩石圈物质杨氏模量随深度分布图(Wang et al.，2005；Huang et al.，2007；石耀霖等，2008； Robert et al.，2009； Wang et al.，2010；Zhang et al.，2011) Fig.3 Relationship of Young′s modulus with depth in the lithosphere of the Longmen Shan area |
岩石圈模型中青藏高原与四川盆地的上地壳及岩石圈上地幔部分具有相同的流变性质与参数，而中、下地壳流变性质有非常明显的差异.对于黏滞系数而言，青藏高原中、下地壳较软；而四川盆地下方中、下地壳相对坚硬.参照石耀霖等(2008)的中国大陆岩石圈等效黏滞系数的研究结果，相应的层位黏滞系数取为: 高原下方中地壳为3.0×1020 ~5.0×1020 Pa·s, 下地壳为3.0×1019 ~6.0×1019Pa·s，四川盆地中地壳为1.0×1022 ~1.2×1022Pa·s，下地壳为7.0×1021 ~8.0×1021Pa·s.相应岩石圈各层位黏弹性物质的松驰时间(介质黏滞系数与杨氏模量之比的二分之一)分别为：青藏高原与四川盆地上地壳为30000a左右；青藏高原中、下地壳分别为200a、50a左右，四川盆地中、下地壳分别为1000a、500a左右；青藏高原与四川盆地岩石圈上地幔为300a左右.模型岩石圈物质黏滞系数随深度分布的变化曲线见图4 .
 | 图4 龙门山断裂带岩石圈各层位岩石介质黏滞系数在青藏高原一侧与四川盆地一侧随深度分布图(石耀霖等，2008) Fig.4 Relationship of viscosity with depth in the lithosphere of the Longmen Shan area |
针对以上黏弹性Maxwell体模型，在数千年的时间尺度内考察岩石圈的应力积累过程中，不论是青藏高原还是四川盆地，在地壳的缩短过程中，上地壳将表现为脆性，而中、下地壳介质的流变性质则可以充分得到体现(这取决于黏弹性介质的松驰时间).由于Maxwell体具有这种可以自然处理脆-韧性转变的优点，本文有限元模型实际可以有效地描述弹性上地壳覆盖在黏弹性中、下地壳和岩石圈上地幔的应力积累过程.
我们将该岩石圈模型命名为参考模型0，以方便与下文中另外的几种岩石圈模型加以区别.模型0从上到下一共分为13层，其中1—7层为上地层，8—9层为中地壳，10—11层为下地壳，12—13层为 岩石圈上地幔，图2 中不颜色代表不同分层的物质参数，其具体介质参数见表1 与表2 .用六面体单元对模型进行网格划分，单元总数为182400个, 节点总数为279600个.
|
表1 青藏高原东缘岩石圈物质参数 Table 1 ESRMaterial parameters of rocks in the lithosphere of the eastern margin of the Tibetan Plateau |
|
|
表2 四川盆地岩石圈物质参数 Table 2 Material parameters of rocks in the lithosphere of the Yangtze craton |
大量的研究对青藏高原的流变结构进行过探讨(Royden et al，1997，2008；Clark et al，2000，2005；Beaumont et al，2001，2004)，并对青藏高原下地壳的运动模式持有不同的见解(Royden et al.，2008).Clark等(2000，2005)认为青藏高原重力驱动作用下的低黏度(黏滞系数为2.0×1018Pa·s)的下地壳流动速度大约在80 mm/a左右(为地表运动速度的8倍左右).Cao等(2009)通过建立青藏高原的三维岩石圈模型拟合了青藏高原地表的水平运动速度，其结果表明青藏高原下地壳的流动速度仅比地表运动速度快8 mm/a左右，并且所估计的下地壳黏滞系数也比Royden与Clark的估计结果高出一个量级，为3.0×1019Pa·s.而Bendick等(2007)和Wang等(2008)则支持下地壳与上、中地壳和岩石圈上地幔之间无差异运动的观点.因此，我们在有限元模型中对青藏高原一侧的物质运动采用了三种不同的位移边界条件，以满足以上几种不同的观点，见图5 .
边界条件1 (BC1)：地壳运动观测网络给出了研究地区1998—2004年的地表GPS速度场结果，见图1 .将GPS 在研究区域边界附近实测速度值插值到四个侧边界，作为水平速度约束条件，且假定从地表到100 km深度保持一致，以拟合藏青高原下地壳的无差异运动的情况(Bendick et al.，2007；Wang et al.，2008)；垂直方向位移可以保持自由.上表面为自由边界，即法向应力和剪应力均为零.对于底部边界，鉴于目前对该区域岩石圈上地幔运动状态的未知，暂且将底面垂直方向速度约束为0，而水平方向自由.在黏弹性问题中，边界条件随时间的变化也是重要的问题，在目前的模拟中，我们初步假定边界位移速率不随时间变化，见图5 .
 | 图5 模型边界条件 剖面1垂直于龙门山断裂带走向，剖面2平行于龙门山断裂带走向.BC1，BC2和BC3分别为边界条件1，2和3的缩写 Fig.5 Boundary conditions Transection 1 is perpendicular to the strike of the Longmen Shan fault, while transection 2 is parallel to the strike of the Longmen Shan fault. BC1, BC2 and BC3 are abbreviations of boundary conditions 1, 2 and 3, respectively |
边界条件2 (BC2)：我们采用Cao等(2009)的观点，认为青藏高原下地壳以Channel flow形式的运动速度V1比地表运动速度V快8 mm/a.其他边界条件同边界条件1，见图5 .
边界条件3 (BC3)：青藏高原下地壳以Channel flow形式的运动速度V2为地表运动速度V的3倍.需要指出的是，V2的取值并没有理论依据，目的是为了测试下地壳不同的流动速度对模型构造应力积累的影响.其他边界条件同边界条件1，见图5 .
初始条件是构造应力场模拟中最困难的问题，尽管现今在龙门山断裂带进行了不同程度的应力测量结果，但基本都限于沉积层深度的钻孔应力测量，而对深部三维应力分布和应力演变历史几乎仍然没有定量的资料.在这种情况下，我们只能先假定初始应力为0，然后计算在定长的边界位移速率条件下应力的演变.虽然我们不可能模拟现今真实的应力状态，但我们可以了解在定长边界位移速率下的应力增长率.应力增长率高的地方，未必一定是目前应力最高的地方，但如果初始应力类似，则较高应力增长率的地方则更有可能是现今应力绝对值较大的区域.本文主要就应力增长率与地震活动性的关系进行讨论.
计算所使用的程序是利用“飞箭有限元程序自动生成系统(FEPG)”生成的Maxwell体三维有限元计算程序，程序计算的可靠性已经在大量事例中通过验证.在保证程序可靠性的前提下，对模型进行计算是合理的.本文采用1a为一个时间步长，根据时间步长逐步加载边界位移约束.本文计算主要的目标是计算龙门山断裂带岩石圈各层位在数千年时间尺度以上的应力积累过程，在计算过程中只考虑模型在边界条件作用下产生的构造应力的积累变化，并不考虑重力的因素.
黏弹性Maxwell体在外部载荷下的变形，不仅与边界条件及其随时间的变化有关，而且与初始条件和以前的应力演变历史有关.鉴于初始条件缺乏测量资料，只能假定为零应力、初始应变速率为0.边界条件也假定了位移速率为常量、不随时间变化.在这种条件下，开始的数百年内，下地壳、岩石圈上地幔等黏滞系数较小、弛豫时间较短(数十到数百年)的柔性可以占主导的层位，在压缩位移下的应力增长与柔性介质内的应力松弛将达到平衡，应力维持在一个较低水平而不再增加.相反，上地壳黏滞系数高(弛豫时间达数万年)的层位，在几百到几万年的期间内，弹性仍然占主导，应力随压缩几乎可以接近线性的速率增长.该结论在我们的另一研究中已经得到证明(柳畅等，2012b).因此，在模型中各层位应力积累达到稳定速率增长后，我们分别取经过芦山地震震源的纵剖面(图5 中剖面1)和芦山地震震源深度的横切面上的应力积累速率分布，以分析龙门山断裂带数千年以上的应力积累及其与区域地震活动性之间的关系.
在边界条件1(下地壳无差异运动)作用下参考模型0，垂直于龙门断裂带方向的压应力积累分布情况如下.
 | 图6 (a) 边界条件1(BC1)作用下模型0中经过芦山地震震中且垂直于龙门山断裂带的剖面上的应力积累速率分布状况；(b)边界条件2(BC2)作用下该剖面上的应力积累速率分布状况.应力在龙门山断裂带下集中积累，在龙门山断裂带上地壳底部有最大增长率.黑色星 星表示汶川地震和芦山地震主震震源在该剖面上的投影 TP表示青藏高原，LMS表示龙门山，SB表示四川盆地 Fig.6 (a) Normal stress accumulation rate distribution in the cross section of model 0 in the boundary condition 1 (BC1)；(b) Normal stress accumulation rate distribution in the boundary condition 2 (BC2)；The cross section is perpendicular to the strike of the Longmen Shan fault. Stress concentrates at the bottom of the upper crust of the Longmen Shan fault. The largest stress accumulation rate is located at the bottom of the upper crust of the Longmen Shan fault. The black stars indicate the hypocenters of the 2008 Wenchuan earthquake and the 2013 Lushan earthquake. TP indicates the Tibetan Plateau, LMS indicates the Longmen Shan, and SB indicates the Sichuan Basin |
我们取经过芦山地震震中、且垂直于龙门山断裂带走向的纵剖面(图5 中剖面1)上的应力积累速率分布，见图6 a，以揭示应力积累与龙门山区域地震纵向空间分布特征之间的关系.结果显示在脆性的上地壳应力以近乎定值的速率持续增长，在龙门山断裂带上地壳底部出现应力集中现象，最大应力积累速率为-3.86 kPa/a.而中、下地壳及岩石圈上地幔的应力积累速率因物质的黏性应力松驰而为0.这种应力积累状态表明，在构造挤压作用下龙门山断裂带的上地壳内应力可以持续增长；在应力积累速率最大的龙门山断裂带上地壳底部，应力可以优先增长至该部分岩石破裂强度而导致主震的发生(汶川地震主震震源深度19 km，芦山地震主震震源深度14 km)，进而触发上地壳内部的余震(5~25 km)；而在中、下地壳内，当构造挤压作用下的应力增长与介质黏性应力松驰效应达到动态平衡后，其应力趋于稳定，不再增长，因此该层位应力积累很难达到岩石的破裂强度.汶川地震和芦山地震的余 震震源定位结果均表明，在龙门山断裂带中、下地壳深度范围内鲜见余震发生.
我们取上地壳深度14 km处(芦山地震源深度)横切面上的应力积累速率分布，见图7 .结果显示青藏高原上地壳应力增长率(约-2.5 kPa/a)远远高于四川盆地的(约-1.0 kPa/a)，而最大应力增长率(-3.4 kPa/a)则位于龙门山断裂带.这种应力积累分布现象有利于解释龙门山区域地震横向空间分布特征：龙门山断裂带青藏高原一侧地震(Mw> 3)密集分布(李勇等，2006；滕吉文等，2008；雷兴林等，2013)，且有历史大震发生，如：1630年M6.5级虎牙地震、1913年M7.0级叠溪地震、1932年M7.5级叠溪地震、1960年M6.7级松潘地震与1976年3次震级为6.6<M<7.3级的松潘地震群；最大级别地震——汶川Mw7.9级地震，发生在龙门山断裂带；而四川盆地地震稀少(Mw>3) (李勇等，2006；滕吉文等，2008；雷兴林等，2013).
 | 图7 边界条件1作用下模型0中深度14 km处横切面上的应力积累速率分布状况；青藏高原应力增长速率远远高于四川盆地，在龙门山断裂带集中积累有最大值.黑色星星表示Mw6.6级芦山地震主震.黑色三角标志城市位置.白色线表示断裂 Fig.7 Normal stress accumulation rate distribution in the map at the depth of 14 km in model 0 in the boundary condition 1 (BC1). The stress accumulation rate in the Tibetan Plateau is much higher than that of the Sichuan Basin. The stress concentrates at the Longmen Shan fault. The black star indicates the hypocentre of the 2013 Lushan earthquake. Black triangles indicate the location of cities |
我们计算了边界条件2(下地壳流动速度快于地表8 mm/a)与边界条件3(下地壳流速度为地表速度的3倍)作用下模型0中的应力积累速率分布状态，以探讨青藏高原下地壳流对龙门山断裂带应力积累过程的影响.计算结果表明，除了在下地壳的上、下底界面处(相应于图6 b中L层位的上、下界面) 因下地壳的快速流动而引起的拖曳力(张应力)外，边界条件2和3作用下模型的应力积累速率分布与边界条件1作用下的结果并无太大区别.因此，我们只展示边界条件2作用下，模型0中的应力积累速率分布状况(图6 b).结果显示龙门山断裂带上地壳底部同样出现应力集中现象，且该部位应力 增长速率有所增大.取芦山地震震源处的应力增长率分析，边界条件1，边界条件2和边界条件3作用下该处的应力增长率分别为-3.60 kPa/a，-3.75 kPa/a和-3.9 kPa/a， 见表5 .可见，下地壳流动有助于龙门山断裂带的应力积累，并且随流动速度增大而增大，但增大量并不显著.
为探讨青藏高原地壳流变结构对龙门山断裂带应力积累的影响，我们在参考模型0的基础上另外引入了四个不同流变结构的岩石圈模型(模型1、模型2、模型3和模型4).在这四个模型中，我们仅在参考模型0的基础上改变了岩石圈不同层位的黏滞系数，见表3 和表4 ，而其他物质参数保持不变，且均加载边界条件1.我们复述模型0的参数取值以方便与其他四个模型做比较.各个模型的具体黏滞系数随深度分布见图8 .
|
|
表3 五个不同模型中青藏高原部分不同岩石圈层位的黏滞系数(Pa·s) Table 3 Viscosity of rock layers of the eastern Tibetan Plateau in the five models |
|
|
表4 五个不同模型中四川盆地部分不同岩石圈层位的黏滞系数(Pa·s) Table 4 Viscosity of rock layers of the Sichuan Basin in the five models |
模型0 ：上地壳表现为脆性；中、下地壳表现为韧性，且横向上青藏高原与四川盆地的黏滞系数有差异；岩石圈上地幔表现为韧性.黏滞系数的取值在 上地壳为1.0×1023Pa·s；高原下方中地壳为3.0× 1020 ~5.0×1020Pa ·s, 下地壳为3.0×1019 ~6.0×1019Pa·s，四川盆地中地壳为1.0×1022~1.2×1022Pa·s，下地壳为7.0×1021~8.0×1021Pa·s； 岩石圈上地幔为4.0×1020~6.0×1020Pa·s.
模型1 ：整个岩石圈表现为脆性，岩石圈各层位有统一的较高的黏滞系数,为1.0×1023Pa·s，且横向上从青藏高原到四川盆地黏性无差异，具体黏滞系数随深度分布见图8e.
模型2 ：整个地壳表现为脆性，地壳各层位有统一的较高的黏滞系数,为1.0×1023Pa·s；岩石圈上地幔表现为韧性且黏滞系数与模型0中相同.具体黏滞系数随深度分布见图8f.该模型中所描述的岩石圈流变结构符合Tapponnier等(2001)支持的青藏高原全脆性地壳模型特征.
模型3 ：上地壳表现为脆性；中、下地壳和岩石圈上地幔表现在为韧性，且横向上青藏高原与四川盆地的黏滞系数无差异(这是该模型与模型0的唯一区别)，中、下地壳黏滞系数分别为3.0×1020 ~5.0×1020Pa·s和3.0×1019~6.0×1019Pa·s；其他层位黏滞系数均与模型0相同.具体黏滞系数随深度分布见图8g.
 | 图8 不同模型中经过芦山地震震中且垂直于龙门山断裂带的纵剖面上应力积累速率分布，与相应模型的黏滞系数随深度分布图 (a)、(b)、(c)和(d)分别为模型1、模型2、模型3和模型0中的应力分布图；(e)、(f)、(g)和(h)分别为模型1、模型2、模型3和模型0中黏 滞系数随深度分布图.模型0的龙门山断裂带上地壳底部出现应力集中现象.U表示上地壳，M表示中地壳，L表示下地壳.TP表示青藏高原，LMS表示龙门山，SB表示四川盆地.白色虚线表示岩石圈不同层位的分界面 Fig.8 Normal stress accumulation rate distribution in the transaction 1 (in Fig.5 ) in different models by boundary condition 1 (a) Model 1； (b) Model 2；(c) Model 3 and (d) Model 0. The viscosity distributes with depth in different models: (e) Model 1, (f) Model 2, (g) Model 3 and (h) Model 0. Stress concentrates heavily at the bottom of the upper crust of the Longmen Shan fault in Model 0. Black stars indicate the hypocentres of the 2008 Wenchuan earthquake and the 2013 Lushan earthquake. The black dot indicates the transferring point (marked as A) on the Moho surface. TP indicates the Tibetan Plateau. LMS indicates the Longmen Shan. SB indicates the Shichuan Basin. U indicates upper crust. M indicates middle crust. L indicates lower crust. White dash lines indicate the contact surface in the crust layers. |
模型4 ：仅将模型0中青藏高原中、下地壳的黏滞系数均减小为原来的20倍，则其中、下地壳黏滞系数分别为1.5×1019~2.5×1019Pa·s 和1.5×1018~3.0×1018Pa·s.这一取值仍在Clark等 (2005)、石耀霖等(2008)和Godard等(2009)对青藏高 原中、下壳的黏滞系数的估计范围之内，为合理取值.
引入模型3、模型0和模型4其目的是为了研究青藏高原与四川盆地中、下地壳黏滞系数差异从无到有、到增大差异这三种情况下，相应岩石圈模型应力积累分布的区别.
由于结果显示模型4与模型3的应力积累结果分布十分相似，因此我们只给出边界条件1作用下模型1、模型2、模型3和模型0中经过芦山地震震中且垂直于龙门山断裂带的纵剖面上的应力积累速率分布状况，见图8 .结果具体如下：
模型1中，如图8 a所示，整个地壳内部横向上应力积累速率均匀分布，没有应力集中的现象.纵向上应力积累速率随深度递增大，不同层位的应力增长率随层位加深、杨氏模量的增大而增加.芦山地震的震源处应力速率为-2.1 kPa/a.
模型2中，如图8 b所示，在龙门山断裂带下地壳的底部出现应力集中现象，且模型1中Moho面拐点处A应力积累速率从-2.8 kPa/a增大到模型2中的-4.3 kPa/a.可见，在Moho面起伏的岩石圈结构中，构造应力积累在很大程度上受控于岩石圈的流变结构，这一点在我们的另外一个研究中也同样得到验证(柳畅等，2012a).地壳其他部分应力积累分布与模型1中相似.韧性的岩石圈上地幔部分应力积累速率为0.
模型3中，如图8 c所示，上地壳内龙门山断裂带下方并无应力集中现象出现，并且在上地壳内青藏高原的应力积累速率小于四川盆地的.韧性的中、下地壳和岩石圈上地幔部分应力积累速率为0.
以上结果表明，模型1、模型2与模型3的应力积累结果均不能与龙门山地区的地震空间分布特征相对应(包括纵向与横向).从应力积累的角度而言，这三种模型既无法解释汶川地震与芦山地震的发生，也无法解释如上文所述的龙门山断裂带地震空间分布特征(李勇等，2006；滕吉文等，2008；雷兴林等，2013).
只有模型0中的应力积累分布特征与龙门山地区的地震空间分布特征有较好的吻合，这一点在上文中已得到详细分析，这里不再赘述.模型0与其他几个模型相比较，其最重要的特征就在于模型中青藏高原与四川盆地的中、下壳的黏性差异.
模型4中各层位的应力积累分布与模型0的结 果相似.但是由于模型4相对于模型0增大了青藏高原与四川盆地中、下地壳的黏性差异(降低青藏高原中、下地壳黏滞系数为原来的1/20)，使得龙门山断裂带的应力集中程度被加剧，相同部位(芦山地震 震源)的应力增长率从模型0中的-3.60 kPa/a增大到模型4中的-4.46 kPa/a.由此可见，青藏高原与四川盆地韧性层中更大的黏性差异将有利于龙门山断裂带的应力加速积累.5个不同模型中芦山地震震源处的应力增长率见表5 .
|
|
表5 三种边界条件作用下五个不同模型中芦山地震震源处的应力增长率(kPa/a) Table 5 Normal stress accumulation rate at the hypocentre of the Lushan earthquake in the 5 models by 3 kinds of boundary conditions |
在2008年汶川地震之后，大量的研究探讨了汶川地震是否为紫坪铺水库所触发的问题.这些研究大多是围绕水库水体的渗流作用对龙门山断裂带的库伦应力的影响，而进行的不同程度的探讨(Parson et al.，2008；Ge et al.，2009；Stone et al.，2009；Deng et al.，2010；Zhu et al.，2010)，且这些研究所支持的观点并不尽相同.这种因素的作用可能在某种程度上对龙门山断裂上的应力改变产生了一定的影响，但是导致汶川地震发生的主要原因还在于，长期构造挤压环境下龙门山断裂带应力的集中积累、释放的结果.
本研究中我们建立了四种不同流变结构的龙门山断裂岩石圈模型，计算了青藏高原下地壳在无差异运动和channel flow的情况下，该区域数千年尺度下的应力积累过程与分布状况.通过比较不同模型的应力积累结果表明，仅在青藏高原与四川盆地的中、下地壳存在黏性差异的模型0中，龙门山断裂带上地壳底部才出现应力集中现象，这一应力集中现象可以解释汶川地震与芦山地震的发生；同时，该模型中的应力积累空间分布状态能与地震空间分布特征有较好的吻合，可以解释龙门山区域地震空间分布特征：纵向上，龙门山断裂带这两次地震主震均发生在龙门山断裂带上地壳的底部(14~19 km)，绝大部分余震均发生在上地壳范围(5~25 km)，而在其中、下地壳深度范围内鲜见余震发生；横向上，地震(Mw>3)在龙门山断裂带青藏高原一侧密集分布且曾有大震发生，而四川盆地地震(Mw>3)稀少(李勇等，2006；滕吉文等，2008；雷兴林等，2013).除了模型0外，其他3个模型中龙门山断裂带均未见有应力集中的现象出现.由此推断，在构造挤压的条件下，青藏高原与四川盆地中、下地壳的黏滞系数的差异是龙门山断裂带应力集中积累、优先达到断层强度而导致地震发生的主要原因.
通过研究2008年汶川地震余震系列的震源定位结果，张瑞青等(2008)认为汶川地震的余震系列主要分布在8~22 km的上地壳范围内，进而指出龙门山断裂带只存在于上地壳22 km深度范围内，并没有穿透65 km左右的整个地壳.刘巧霞等(2010)的汶川地震余震震源定位结果也支持张瑞青等的这一推断.震源定位所指向的龙门山断裂带的发震层仅为上地壳深度范围，这一地震观测事实在很大程度上表明，青藏高原的中、下地壳物质表现为较低力学强度的韧性特征.以上对龙门山区域地壳的流变特征这一推断与Tapponnier等(2001)所认为的青藏高原东缘全脆性地壳模型并不一致.我们在模型2中模拟了整个脆性地壳上伏在韧性岩石圈上地幔的岩石圈结构中的应力积累分布，结果并没有显示龙门山断裂上地壳底部的应力集中现象；而是出现了龙门山断裂带下地壳底部的应力集中现象.这一应力积累结果与龙门山断裂带的两次地震主震与绝大部分余震的震源深度并不相符合.所以，从应力积累的角度而言，我们认为全脆性地壳的地震地质模型(Tapponnier et al.，2001；Hubbard et al.，2009)不能解释龙门山断裂带地震的发生.大量的地质学证据(Bird et al.，1991)、数值模拟青藏高 原地壳流变结构的研究(Clark 2000，2005；Beaumont et al，2001，2004)、 以及地球物理观测包括GPS地表运动监测(Shen et al.，2005)与深部大地电磁测量数据(Bai et al.，2010)所反映的龙门山断裂带岩石圈的流变结构，均支持青藏高原东缘有较弱的柔性中、下地壳这一观点.
Godard等(2009)通过龙门山断裂带岩石圈的温度计算推断青藏高原东缘中、下地壳的黏滞系数在1.0×1019~1.0×1021Pa·s之间.该估计与Hilley等(2005)对于青藏高原北部中、下地壳的黏滞系数的估计一致；但是比Clark等(2000，2005)通过模拟龙门山的地形而推断的2.0×1018Pa·s黏滞系数要高出一到两个量级.尽管关于青藏高原韧性中、下地壳的黏滞系数在各家研究中不尽相同，比较模型0与模型4的计算结果表明，更低的青藏高原中、下地壳黏滞系数则有利于加速龙门山断裂带的应力积累.
青藏高原下地壳的运动模式(无差异运动或下地壳流)在不同的研究中仍然存在较大的争议，但是我们的计算结果表明，青藏高原下地壳不论做何种运动模式，应力集中现象均会在模型0的相同部位出现.在假设存在下地壳流存在的情况下，我们的结果表明下地壳的快速流动将有利于应力在龙门山断裂带的加速积累.其原因可解释为青藏高原下地壳流的东向快速流动加剧了物质在龙门山断裂带下方的堆积程度，而有利于应力的加速增长.
需要指出的是，在本文的计算中并没有将龙门山断裂带的断层视为力学性质的薄弱带考虑到计算的模型当中，其原因有二：其一，本文计算讨论的主要目的是，发现龙门山断裂带在构造挤压的环境下，地壳的流变结构对应力积累的控制因素；其二，从单纯的数值模拟本身来看，如果在模型中加入龙门山地震带的断层薄弱带，那么薄弱带势必为周边非断层介质所包围.一个完整的介质内部夹入一个薄弱带的模型，在挤压的边界条件情况下，将会产生人为的应力集中现象，其结果未必与实际力学情况接近或者相符.所以并未将断层当作力学性质的薄弱带考虑到模型中来.
我们建立了四种不同流变结构的龙门山断裂岩石圈模型，计算并比较了这几种模型在数千年以上、长期匀速构造挤压作用下的应力积累特征，分析了该区域地震空间分布与构造应力积累速率的关系，得到如下结论：
(1)龙门山断裂带在数千年的应力积累过程中，脆性上地壳中应力表现近于恒定值的线性增长趋势，龙门山断裂带上地壳底部出现应力集中积累现象，这一应力集中现象可以解释龙门山断裂带汶川地震与芦山地震的主震的发生，及其大部分余震在脆性上地壳中的触发.而在柔性的中、下地壳内，应力积累在稳定状态之后其增长速率近于零，构造应力积累很难达到岩石破裂强度，而鲜见地震发生.
(2) 青藏高原一侧上地壳应力积累速率远远高于四川盆地的应力积累速率，这一应力积累分布现象可以解释龙门山区域青藏高原一侧地震密集而四川盆地地震稀少的地震空间分布特征.
(3)通过四种不同流变结构模型的应力计算结果比较，认为导致龙门山断裂带以上应力积累空间分布状态的重要控制因素在于青藏高原中、下地壳较低的黏滞系数与四川盆地中、下地壳较高的黏滞系数的差异.地壳各层位的应力增长率差异与地震成层分布的现象共同揭示了龙门山区域岩石圈分层流变结构：脆性上地壳、韧性中、下地壳(青藏高原一侧较弱，四川盆地一侧较强)、韧性岩石圈上地幔.
致谢 感谢苏黎世联邦理工学院的朱桂芝博士参与结果讨论.| [1] | An M J, Shi Y L. 2007. 3D temperature of crust and lithosphere mantle in Chinese continent.Science China (D) (in Chinese), 37(6): 736-745. |
| [2] | Bai D, Unsworth M J, Meju M A, et al. 2010. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging. Nature Geoscience, 3:359-362. |
| [3] | Beaumont C, Jamieson R A, Nguyen M H, et al. 2001. Himalayan tectonics explained by extrusion of a low viscosity crustal channel coupled to focused surface denudation. Nature, 414: 738-742. |
| [4] | Beaumont C, Jamieson R A, Nguyen M H, et al. 2004. Crustal channel flows: 1. Numerical models with applications to the tectonics of the Himalayan-Tibetan orogen. J. Geophys. R.,109(B6): B06406, doi: 10.1029/2003JB002809. |
| [5] | Bendick R, Flesch L. 2007. Reconciling lithospheric deformation and lower crustal flow beneath central Tibet. Geology, 35(10): 895-898. |
| [6] | Bird P. 1991. Lateral extrusion of lower crust from under high topography in the isostatic limit. J. Geophys. R., 96(B6): 10 275-10 286. |
| [7] | Burchfiel B C, Chen Z, Liu Y, et al. 1995. Tectonics of the Longmen Shan and adjacent regions, Central China. International Geology Review, 37(8): 661-735. |
| [8] | Burchfiel B C, Royden L H, van der Hilst R D, et al. 2008. A geological and geophysical context for the Wenchuan earthquake of 12 May 2008, Sichuan, People's Republic of China. GSA Today, 18(7): 4-11. doi: 10.1130/GSATG18A.1. |
| [9] | Cao J L, Shi Y L, Zhang H, et al. 2009. Numerical Simulation of GPS observed clockwise rotation around the eastern Himalayan syntax in the Tibetan Plateau. Chinese Science Bulletin, 54(8): 1398-1410. |
| [10] | Chen J H, Liu Q Y, Li S C, et al. 2009. Seismotectonic study by relocation of the Wenchuan MS8.0 earthquake sequence. Chinese J.Geophys. (in Chinese), 52(2): 390-397. |
| [11] | Chen Y T, Yang Z X, Zhang Y, et al. 2013. From Wenchuan earthquake to Lushan earthquake. Science China (D) (in Chinese), 43(6): 1064-1072. |
| [12] | Clark M K, Royden L H. 2000. Building the eastern margin of Tibet by lower crustal flow. Geology, 28(8): 703-706. |
| [13] | Clark M K, Bush J M, Royden L H. 2005. Dynamic topography produced by lower crustal flow against rheological strength heterogeneities bordering the Tibetan Plateau. Geophysical Journal International, 162(2): 575-590. |
| [14] | Deng K, Zhou S, Wang R, et al. 2010. Evidence that the 2008 Mw7.9 Wenchuan earthquake could not have been induced by the Zipingpu Reservoir. Bulletin of the Seismological Society of America, 100(5B): 2805-2814. |
| [15] | Du F, Long F, Ruan X, et al. 2013. The M7.0 Lushan earthquake and the relationship with the M8.0 Wenchuan earthquake in Sichuan, China. Chinese J. Geophys. (in Chinese), 56(5): 1772-1783. |
| [16] | Ge S M, Liu M, Lu N, et al. 2009. Did the Zipingpu Reservoir trigger the 2008 Wenchuan earthquake? Geophysical Research Letters, 36(20): L20315, doi: 10.1029/2009GL040349. |
| [17] | Godard V, Cattin R, Lave J, et al. 2009. Erosional control on the dynamics of low convergence rate continental plateau margins. Geophysical Journal International, 179(2): 763-777. |
| [18] | Guo B, Liu Q Y, Chen J H, et al. 2009. Teleseismic P-wave tomography of the crust and upper mantle in Longmenshan area, west Sichuan. Chinese J. Geophys. (in Chinese), 52(2): 346-355. |
| [19] | Hilley G E, Burgmann R, Zhang P Z, et al. 2005. Bayesian inference of plastosphere viscosities near the Kunlun Fault, northern Tibet. Geophysical Research Letters, 32(1): L01302.doi: 10.1029/2004GL021658. |
| [20] | Huang Y, Wu J P, Zhang T Z, et al. 2008. Relocation of the M8.0 Wenchuan earthquake and its aftershock sequence. Sciecne China (D) (in Chinese), 38(10): 1242-1249. |
| [21] | Hubbard J, Shaw J H. 2009. Uplift of the Longmen Shan and Tibetan Plateau, and the 2008 Wenchuan (M=7.9) earthquake. Nature, 458(7235): 194-197. |
| [22] | Huang J L, Zhao D P, Zheng S H. 2007. Lithospheric structure and it's relationship to seismic and volcanic activity in southwest China. J. Geophys. R., 107(B10): ESE 13-1-ESE 13-14. |
| [23] | Hu S B, He L J, Wang J Y. 2001. Compilation of heat flow data in the China continental area (the 3rd edition). Chinese J. Geophys. (in Chinese), 44(5): 611-626. |
| [24] | Lei X L, Ma S L, Su J R, et al. 2013. Inelastic trigering of the 2013 Mw6.6 Lushan earthquake by the 2008 Mw7.9 Wenchuan earthquake. Seismology and Geology (in Chinese), 35(2): 411-422. |
| [25] | Li C, van der Hilst R D,Toksoz N M. 2006. Constraining spatial variations in P-wave velocity in the upper mantle beneath SE Asia. Physics of the Earth and Planetary Interiors, 154(2): 180-195. |
| [26] | Li C, van der Hilst R D, Engdahl E R, et al. 2008. A new global model for P wave speed variations in Earth's mantle. Geochem. Geophys. Geosyst., 9(5): Q05018. doi: 10.1029/2007GC001806. |
| [27] | Li Y, Zhou R J, Densemore A L, et al. 2006. The geodynamics process and geology response at the eastern margin of the Tibetan Plateau. Beijing: Geology Publication House (in Chinese). |
| [28] | Liu C, Shi Y L, Zheng L, et al. 2012a. Relation between earthquake spatial distribution and tectonic stress accumulation in the North China Basin based on 3D visco-elastic modelling. Chinese J. Geophys. (in Chinese), 55(12): 3942-3957. |
| [29] | Liu C, Zhu B J, Shi Y L. 2012b. Numerical modelling stress accumulation on Longmen Shan fault and the recurrence interval of Wenchuan earthquake. Acta Geologica Sinica (in Chinese), 86(1): 157-169. |
| [30] | Liu Q X, Zhu J S, Cao J X, et al. 2010. Relocation of the Wenchuan MS8.0 earthquake sequence and its aftershocks and their sparial distribution characters. Quaternary Science (in Chinese), 30(4): 736-744. |
| [31] | Liu Q Y, Li L, Chen J H. 2009. Wenchuan MS8.0 earthquake: preliminary study of the S-wave velocity structure of the crust and upper mantle. Chinese J. Geophys. (in Chinese), 52(2): 309-319. |
| [32] | Lü J, Wang X S, Su J R, et al. 2013. Hypocentral location and source mechanism of the MS7.0 Lushan earthquake sequence. Chinese J. Geophys. (in Chinese), 56(5): 1753-1763. |
| [33] | Parsons T, Ji C, Kirby E, et al. 2008. Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan Basin. Nature, 454(7203): 509-510. |
| [34] | Robert A, Zhu J, Vergne J, et al. 2010. Crustal structures in the area of the 2008 Sichuan earthquake from seismologic and gravimetric data. Tectonophysics, 491(1-4): 205-210. |
| [35] | Royden L H, Burchfiel B C,King R W. 1997. Surface deformation and lower crustal flow in eastern Tibet. Science, 276(2): 788-790. |
| [36] | Royden L H, Burchfiel B C,van der Hilsta R D. 2008. The geological evolution of the Tibetan Plateau. Science, 321: 1054-1058. |
| [37] | Shen Z K, Lu J N, Wang M, et al. 2005. Contemporary crustal deformation around the southeast borderland of the Tibetan Plateau. J. Geophys. R., 110(B11409): doi:10.1029/2004JB 003421 |
| [38] | Shi Y L, Cao J L. 2008. Effective viscosity of Chinese continent lithosphere. Earth Sciences Frontier (in Chinese), 15(3): 83-95. |
| [39] | Stone R. 2009. Some unwelcome questions about big dams. Science, 324(5928): 714-714. |
| [40] | Tapponnier P, Xu Z Q, Roger F, et al. 2001. Oblique stepwise rise and growth of the Tibet Plateau. Science, 294(5547): 1671-1677. |
| [41] | Teng J W, Bai D H, Yang H, et al. The mechanism and dynamics of the generation and occurrence for Wenchuan Mw7.9 earthquake. Chinese J. Geophys. (in Chinese), 51(5): 1385-1402. |
| [42] | Wang Y X, Wang W D, Han H G. 2008. Crustal P-wave velocity structure from Altyn Tagh to Longmen mountains along the Taiwan-Altay geoscience transect. Chinese J. Geophys. (in Chinese), 2005, 48(1): 98-105. |
| [43] | Wang C Y, Flesch L M, Silver P G, et al. 2008. Evidence for mechanically coupled lithosphere in central Asia and resulting implications. Geology, 36(5): 363-366. |
| [44] | Wang C-Y, Lou H, Silver P G, et al. 2010. Crustal structure variation along 30°N in the eastern Tibetan Plateau and its tectonic implications. Earth and Planetary Science Letters, 289(3): 367-376. |
| [45] | Xu X W, Chen G H, Yu G H, et al. 2013. Seismogenic structure of Lushan earthquake and its relationship with Wenchuan earthquake. Earth Science Frontiers (in Chinese), 20(3): 11-20. |
| [46] | Zhang G W, Lei J S. 2013. Relocations of Lushan, Sichuan strong earthquake (MS7.0) and its aftershocks. Chinese J. Geophys. (in Chinese), 56(5): 1764-1771. |
| [47] | Zhang P Z, Xu X W, Wen X Z, et al. 2008. Slip rate and recurrence intervals of the Longmen Shan fault zone, and tectonic implication of the Wenchuan earthquake 2008. Chinese J. Geophys. (in Chinese), 51(4):1066-1073. |
| [48] | Zhang R Q, Wu Q J, Li Y H, et al. 2008. Location of the hypocenters of the Wenchuan earthquake after shockes and its application. Science China (in Chinese), 38(10): 1234-1241. |
| [49] | Zhang Y, Feng W P, Xu L S, et al. 2008. The temporal and spacial rapture process of the 2008 Wenchuan earthquake. Science China (D) (in Chinese) , 38(10): 1186-1194. |
| [50] | Zhang Z, Deng Y, Teng J, et al. 2011. An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. Journal of Asian Earth Sciences, 40(4): 977-989. |
| [51] | Zheng Y, Ge C, Xie Z J, et al. 2013. the structure and enviroment of earthquake generation of the crust and mantle in Wenchuan earthauke and Lushan earthquake zone. Science China (D) (in Chinese), 43(6): 1027-1037. |
| [52] | Zhu B J, Liu C, Shi Y L, et al. 2011. Application of flow driven pore-network crack model to Zipingpu reservoir and Longmenshan slip. Science China (E), Mechanics and Astronomy, 54(8): 1532-1540. |
| [53] | 安美建, 石耀霖. 2007. 中国大陆地壳和上地幔三维温度场. 中国科学(D), 37(6): 736-745. |
| [54] | 陈九辉, 刘启元, 李顺成等. 2009. 汶川MS8.0地震余震序列重新定位及其地震构造研究. 地球物理学报, 52(2): 390-397. |
| [55] | 陈运泰, 杨智娴, 张勇等. 2013. 从汶川地震到芦山地震. 中国科学(D), 43(6): 1064-1072. |
| [56] | 杜方, 龙锋, 阮祥等. 2013. 四川芦山7.0级地震及其与汶川8.0级地震的关系. 地球物理学报, 56(5): 1772-1783. |
| [57] | 郭飚, 刘启元, 陈九辉等. 2009. 川西龙门山及邻区地壳上地幔远震P波层析成像. 地球物理学报, 52(2): 346-355. |
| [58] | 黄媛, 吴建平, 张天中等. 2008. 汶川8.0级大地震及其余震序列重定位研究. 中国科学(D), 38(10): 1242-1249. |
| [58] | 胡圣标, 何丽娟, 汪集旸. 2001. 中国大陆地区大地热流数据汇编(第三版). 地球物理学报, 44(5): 611-626. |
| [59] | 雷兴林, 马胜利, 苏金蓉等. 2013. 汶川地震后中下地壳及上地幔的黏弹性效应引起的应力变化与芦山地震的发生机制. 地震地质, 35(2):411-422. |
| [60] | 柳畅, 石耀霖, 郑亮等. 2012a. 三维黏弹性数值模拟华北盆地地震空间分布与构造应力积累关系.地球物理学报, 55(12): 3942-3957. |
| [61] | 柳畅, 朱伯靖, 石耀霖. 2012b. 粘弹性数值模拟龙门山断裂带应力积累及大震复发周期. 地质学报, 86(1): 157-169. |
| [62] | 李勇, 周荣军, Densemore A L等. 2006. 青藏高原东缘大陆动力学过程与地质响应. 北京: 地质出版社. |
| [63] | 刘巧霞, 朱介寿, 曹俊兴等. 2010. 汶川MS8.0 级地震余震重新定位及其空间分布特征研究. 第四纪研究, 30(4): 736-744. |
| [64] | 刘启元, 李昱, 陈九辉. 2009. 汶川MS8.0地震:地壳上地幔S波速度结构的初步研究. 地球物理学报, 52(2): 309-319. |
| [65] | 吕坚,王晓山,苏金蓉等.2013.芦山7.0级地震序列的震源位置与震源机制解特征. 地球物理学报, 56(5): 1753-1763. |
| [66] | 石耀霖,曹建玲. 2008. 中国大陆岩石圈等效粘滞系数的计算和讨论. 地学前缘, 15(3): 83-95. |
| [67] | 滕吉文, 白登海, 杨辉等. 2008. 2008汶川MS8.0地震发生的深层过程和动力学响应. 地球物理学报, 51(5): 1385-1402. |
| [68] | 徐锡伟, 陈桂华, 于贵华等. 2013. 芦山地震发震构造及其与汶川地震关系讨论. 地学前缘, 20(3):11-20. |
| [69] | 张广伟, 雷建设. 2013. 四川芦山 7.0 级强震及其余震序列重定位. 地球物理学报, 56(5): 1764-1771. |
| [70] | 张培震, 徐锡伟, 闻学泽等. 2008. 2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因. 地球物理学报, 51(4): 1066-1073. |
| [71] | 张瑞青, 吴庆举, 李永华等. 2008. 汶川中强余震震源深度的确定及其意义. 中国科学(D), 38(10): 1234-1241. |
| [72] | 张勇, 冯万鹏, 许力生等. 2008. 2008年汶川大地震的时空破裂过程. 中国科学(D), 38(10):1186-1194. |
| [73] | 郑勇, 葛粲, 谢祖军等. 2013. 芦山与汶川地震震区地壳上地幔结构及深部孕震环境. 中国科学(D), 43(6): 1027-1037. |