引用本文
赵由佳, 张国宏#cod#x0002A;, 张迎峰, 单新建, 屈春燕. 基于连续-离散单元法的汶川地震动力学二维自发破裂全周期模拟研究[J]. 地震地质, 2018,40(1): 12-15.
ZHAO You-jia, ZHANG Guo-hong, ZHANG Ying-feng, SHAN Xin-jian, QU Chun-yan. TWO-DIMENSIONAL WHOLE CYCLE SIMULATION OF SPONTANE-OUS RUPTURE OF THE 2008 WENCHUAN EARTHQUAKE USING THE CONTINUOUS-DISCRETE ELEMENT METHOD[J]. SEISMOLOGY AND GEOLOGY, 2018,40(1): 12-15.  
DOI:10.3969/j.issn.0253-4967.2018.01.002
Permissions
Copyright©2018, 《地震地质》编辑部
《地震地质》编辑部 所有
基于连续-离散单元法的汶川地震动力学二维自发破裂全周期模拟研究
赵由佳, 张国宏*, 张迎峰, 单新建, 屈春燕
中国地震局地质研究所, 地震动力学国家重点实验室, 北京 100029
*通讯作者: 张国宏, 副研究员, E-mail: zhanggh@ies.ac.ccn。

〔作者简介〕 赵由佳, 女, 1993年生, 固体地球物理学专业在读硕士研究生, 主要研究方向为震源动力学, 电话: 010-62009024, E-mail: 983531711@qq.com。

收稿日期: 2017-11-08
基金项目: 地震动力学国家重点实验室自主课题(LED2014A01)与国家自然科学基金(41474013)共同资助
摘要

首先利用地震地质调查、 地震剖面探测结果对汶川地震的发震断层几何形态及周边块体介质参数进行约束, 建立了3个包含有不同主断层的龙门山断裂带二维黏弹性有限元模型; 然后基于连续-离散单元法对汶川地震的孕育与发生进行了全周期模拟, 包括震间应力积累阶段、 同震应力释放阶段及震后恢复调整阶段。模型对比研究结果表明, 龙门山断裂带中由于断层几何结构及断层面剪切应力与正应力的综合作用, 使得最先达到断层破裂准则阈值的灌县-江油断层率先破裂, 初始破裂在深度15~20km处, 破裂进一步传导至映秀-北川断层; 同时由于灌县-江油断层与映秀-北川断层的破裂导致系统内应力大幅度降低, 使得灌县-江油断层不具备产生同震破裂的应力条件。然而, 尽管汶川-茂汶断层上的应力水平与汶川地震前相比有所减弱, 但仍积累了较高的震间应力, 可能预示其在汶川地震后具有较高的地震危险性。模拟结果还表明, 铲型逆冲断层系统的深部断层破裂能够推动浅部高倾角区域发生被动破裂, 因而其发震断层的倾角上限很可能比传统认识的大。

关键词: 连续-离散单元法; 汶川地震; 自发破裂; 全周期模拟; 多断层相互作用
中图分类号:P315.2 文献标志码:A 文章编号:0253-4967(2018)01-0012-15
TWO-DIMENSIONAL WHOLE CYCLE SIMULATION OF SPONTANE-OUS RUPTURE OF THE 2008 WENCHUAN EARTHQUAKE USING THE CONTINUOUS-DISCRETE ELEMENT METHOD
ZHAO You-jia, ZHANG Guo-hong, ZHANG Ying-feng, SHAN Xin-jian, QU Chun-yan
State Key Laboratory of Earthquake Dynamics, Institute of Geology,China Earthquake Administration, Beijing 100029, China
Abstract

The May 12, 2008 MS7.9 Wenchuan earthquake is ranked as one of the most devastating natural disasters ever occurred in modern Chinese history. The Longmenshan Fault(LMSF)zone is the seismogenic source structure, which consists of three sub-parallel faults, i.e., the Guanxian-Jiangyou Fault(GJF)in the frontal, the Yingxiu-Beichuan Fault(YBF)in the central fault and the Wenchuan-Maowen Fault(WMF)in the back of the LMSF. In this study, geological survey and seismic profiles are used to constrain the faults geometry and medium parameters. Three visco-elastic finite element models of the LMSF with different main faults are established. From the phase of interseismic stress accumulation to coseismic stress release and postseismic adjustment, the Wenchuan earthquake is simulated using Continuous-Discrete Element Method(CDEM). Modeling results show that before the 2008 Wenchuan earthquake, the GJF becomes unstable due to the interaction between its unique fault geometry and the tectonic stress loading. In the fault geometry model, the GJF is the most gently dipped fault among the three faults, which in return makes it having the smallest normal stress and the greatest shear stress. The continuous shear stress loading finally meets the fault failure criteria and the Wenchuan earthquake starts to initiate on the GJF at the depth of 15~20km. The earthquake rupture then propagated to the YBF. At the same time, due to the GJF and YBF rupture, the interseismic stress accumulation has been greatly reduced, causing the WMF failed to rupture. Although the stress accumulation in the WMF has been reduced significantly after the earthquake, yet it has not been released completely, which means that the WMF likely has with high seismic risk after the 2008 Wenchuan earthquake. We also find that the stress perturbation caused by gently dipping segment of the fault can promote the passive rupture in the steeply dipping segment, making the upper limit of dip angles larger than traditional assumption.

Keyword: continuous-discrete element method; Wenchuan earthquake; spontaneous rupture; complete period simulation; multi-fault interaction
0 引言

2008年5月12日四川西部龙门山构造带发生MS8.0汶川大地震。龙门山断裂带作为发震断层, 位于青藏高原的东缘与四川盆地的西缘, 属于松潘-甘孜造山带的前缘逆冲断裂带, 具有典型的逆冲推覆构造特征(图1; 许志琴等, 2007)。龙门山地区海拔变化极其剧烈, 从四川盆地到龙门山地区约10km的距离内, 海拔高度从约500m陡增至约5, 000m; 自西向东发育汶川-茂汶断层(后山断层)、 映秀-北川断层(中央断层)和灌县-江油断层(前山断层)3条主断层(邓起东等, 1994)。野外地质考察发现, 汶川地震造成了十分复杂的地表破裂带(徐锡伟等, 2008)。同时, 考察还发现, 此次地震仅沿映秀-北川断层和灌县-江油断层破裂了龙门山断裂带3条主断层中的2条, 距震中仅十几km的汶川-茂汶断层却未发生地表破裂。

图 1 龙门山地区地震构造图(a)与青藏高原东缘及龙门山地区剖面断层结构示意图(b)
a 黑线表示龙门山地区主要的活动断层, 红线表示汶川地震地表破裂带, 蓝色空心圆表示汶川地震的余震分布(陈九辉等, 2009); b 红线表示2008年汶川地震发震断层, 黑色虚线描绘深处的韧性剪切带(Feng et al., 2016); WMF 汶川-茂汶断层, YBF 映秀-北川断层, GJF 灌县-江油断层, XYDF 小鱼洞断层Fig. 1 Seismotectonic map of the Longmen Shan region(a) and tectonic profile around Longmen Shan(Feng et al., 2016)(b).

发震断层的几何结构对震源破裂传播具有控制性影响, 同时多断层分支存在破裂路径选择及扩展方向等问题。汶川地震发生已近10a, 已有的汶川地震震源动力学破裂过程模拟主要针对主破裂映秀-北川断层进行自发破裂研究(Duan, 2010; 陶玮等, 2011; 朱守彪等, 2016), 均没有考虑现存3条主断层之间的相互作用机制。因而, 对汶川地震震源破裂过程的发震机制、 多条断层之间的相互作用以及汶川地震中未发生地表破裂的汶川-茂汶断层的地震危险性等认识仍然十分有限。另外, 传统有限元方法一般采用负位错与正位错来描述断层面的相对错动(Duan, 2010), 由此带来的问题是无法模拟地震孕育至发震的全周期过程, 因而我们采用离散元与有限元耦合数值计算方法(Continuous-Discontinuous Element Method, CDEM; Wang et al., 2013; 王杰等, 2013, 2015), 即断层单元采用离散元进行模拟, 块体内部则采用有限元方式进行模拟, 从而实现了地震震间加载、 同震破裂及震后恢复的全周期自发破裂过程。

本研究首先利用地震地质调查、 地震剖面探测结果对发震断层几何形态及介质参数进行约束, 建立了龙门山断裂带的二维有限元块体模型。基于虚、 实质量切换计算方法, 实现了震间期准静态缓慢加载到地震能量快速释放期的自动转换。依据龙门山断层带的构造环境设置边界条件得到初始应力场, 基于CDEM方法模拟了汶川地震孕育及发震的全周期过程, 探讨了龙门山断裂带3条主断层之间的应力分配机制、 相互作用关系、 震前应力积累水平和震后应力释放程度等, 分析了汶川地震后汶川-茂汶主断层的地震危险性。

1 CDEM简介及二维构造模型建立

连续-离散单元法(Continuous-Discontinuous Element Method, CDEM)是1种有限元(FEM)和离散元(DEM)耦合的数值计算方法(Wang et al., 2013; 王杰等, 2013, 2015)。这种方法将断层带区域离散为块体单元, 在块体内部进行有限元计算, 在块体边界(断层面)进行离散元计算, 模拟断层面上从连续到离散的破坏演化过程, 因而既能分析连续体的弹塑性变形问题, 又能集成离散体的接触破坏和大变形过程, 具有两者的优势(王杰等, 2013)。

计算模型宽500km, 深度200km, 在深度20km左右设有滑脱层, 层厚约为5km; 莫霍面位于左侧深度60~70km过渡至右侧深度40km处。网格划分采用三角形单元, 单元最大边长为2km; 为获得更精细的应力和位移变化分布, 在断层带处对网格进行适当加密, 最小边长为200m左右。断层倾角参照地震剖面图及地表调查结果综合分析后进行设置(Xu et al., 2009; Zhang et al., 2010; Feng et al., 2016), 其中, 汶川-茂汶断层倾角从深部40° 延伸至地表75° ; 映秀-北川断层从深部倾角近乎水平连接滑脱面延伸至地表70° ; 灌县-江油断层倾角从深部20° 延伸至地表40° 。

2 介质参数、 边界加载条件及断层破裂准则

由于不同区域的介质性质存在差异, 本文将模型分为17个块体。横向以龙日坝断层和龙门山断裂带为界划分成松潘地块、 龙门山以及四川盆地3部分; 纵向划分依次为上地壳、 中地壳、 下地壳、 上地幔。其中, 上地壳采用线弹性模型, 中、 下地壳与上地幔采用蠕变Bugers与Mohr-Coulomb耦合黏弹塑性蠕变模型, 在原有蠕变Bugers体基础上添加M-C塑性元件用以判断是否进入塑性变形阶段, 模型的黏弹特性由Burgers体实现, 塑性特性由Mohr-Coulomb屈服准则实现(袁海平等, 2006)。

表1 模型材料块体参数 Table1 Model parameters used in simuation

根据龙门山断裂带的构造环境设置位移边界条件来模拟Ma尺度的垂向加厚和横向缩短以及ka尺度的推挤作用。在模型西侧施加水平位移场, 从而模拟青藏高原的E向推挤作用; 模型右侧水平向位移固定, 从而模拟四川盆地的阻挡作用(陶玮等, 2011)。水平位移场量级的选取保证模拟得到的地壳中初始应力量值与实际地壳内的应力量级基本一致, 大约在100~1, 000MPa量级之间(黄福明, 2013)。边界条件设置如图2。

图 2 二维模型块体划分及边界条件Fig. 2 Block partition and boundary conditions of two-dimensional model.

上述分层设置将模型分为17个块体, 得到多个接触面, 在可能出现断层的地方设置节理面模拟断层面。 τ和σn为作用在断层面上的剪切应力和正应力。利用Mohr-Coulomb断裂准则结合摩擦本构方程判断初始破裂点, 初始应力场通过重力均衡后的平衡沉降位移获取(朱守彪等, 2009), 结合岩石的重力和水平位移场, 当应力增量超过断裂准则的阈值时断层实现自发破裂。

τ≤c+μ·σn　μ=μs=0.6(黏结状态)τ> c+μ·σn　μ=μd=0.1(滑动状态)1

式(1)中, c、μ分别为黏聚力和内摩擦系数。

3 模拟结果

为了研究汶川地震多断层破裂问题, 我们建立了3大类对比模型, 分别是包含映秀-北川断层(以下简称YBF)的模型I; 包含汶川-茂汶(以下简称WMF)与映秀-北川2条断层的模型Ⅱ ; 以及包含汶川-茂汶、 映秀-北川与灌县-江油(以下简称GJF)3条断层的模型Ⅲ 。在摩擦系数、 介质参数设置等条件都相同的情况下, 模拟3类模型从震间期断层面应力积累情况(孕震过程)到同震期自发破裂过程, 再到震后应力残余情况的全周期过程。我们将全周期地震问题分为5个主要的阶段: 重力均衡阶段、 水平加载平衡沉降阶段、 震间应力加载阶段(孕震期)、 地震破裂阶段(发震期)及震后恢复期。同时, 将重力均衡阶段过渡至水平加载平衡沉降阶段称为过渡阶段1; 水平加载平衡沉降阶段过渡至震间应力加载阶段称为过渡阶段2; 震间应力加载阶段过渡至地震破裂阶段称为过渡阶段3, 即震前阶段, 如图4所示。我们重点关注地震周期过程中的3个典型阶段即震间、 同震及震后阶段。通过模型对比和应力场分析, 探讨汶川地震在多断层系统选择破裂扩展路径时的可能性, 进一步认识多断层逆冲地震的力学演化与破裂特征。

图 3 地震模拟阶段说明Fig. 3 Schematic earthquake phases in simulation.

3.1 震间期断层应力积累

震间期应力积累是地震在闭锁期时的能量积累, 可以用于判定地震危险性及未来可能的发震位置等。当震间应力加载阶段应力积累阈值超过断裂准则时, 产生地震进入发震阶段(同震期)。模拟孕育期断层面上的正应力与剪切应力累积情况, 能帮助更好地理解同震期多断层破裂模式。

图 4 震前阶段断层面正应力与剪应力累积图Fig. 4 The plot of accumulation of normal and shear stress from loading to release.

图4显示了3类模型震前阶段断层累积的剪切应力与正应力结果。可以发现, 模型I中剪切应力与正应力高值区位于断层浅部高倾角区域(B— C段), 最大剪切应力位于深部低倾角至浅部高倾角过渡处, 即图中O点, 最大正应力位于近地表; 而剪应力与正应力的低值区都处于深部低倾角处(平缓区域, A— B段所示)。模型Ⅱ 中YBF(右)模拟结果与模型I类似, WMF(左)正应力高值区处于浅部高倾角区域, 断层向下扩展正应力逐渐减小; 而剪切应力高值区处于断层深部, 最大剪切应力位于图中M点( 图4)。模型Ⅲ 中YBF(中)与WMF(左)模拟应力结果与模型Ⅱ 基本相同; 而GJF(右)正应力与剪切应力高值区都处于断层浅部。模拟结果证实了同1条断层不同深度的应力积累量确实是不同的(陶玮等, 2011)。

震间期断层面的剪应力及正应力共同决定了未来地震的发震位置及发震能力。可以通过震间期断层面上应力的分布模式初步判断未来地震的初始破裂点。由Mohr-Coulomb断裂准则(式(1))可以看出, 破裂启动及扩展同时与断层面的剪切应力和法向力都有关系。因此, 模型Ⅱ 与Ⅰ 的最可能的破裂起始点可能位于YBF断层(右)深部低倾角至浅部高倾角过渡处, 倾角30° 左右处, 为图中O点。而模型Ⅲ , 最可能成为初始破裂点的位置有3处, 分别位于3条断层各自的最大剪切应力范围内, 即图中M, O, X点所示, 再结合正应力低值区的范围可以大致判断最可能的破裂起始点位于GJF(右)深部低倾角至浅部高倾角过渡处, 即图中X点。总之, 破裂起始点既不是最大正应力也不是最大剪应力处, 而是由断层面上的剪应力高值区与正应力低值区共同决定的。

3.2 同震期应力与位移动态变化

图5显示3种模型在同震阶段剪应力随破裂扩展的动态变化图。 可以看出, 初始发震点处首先聚集能量, 在断层两侧形成正、 负2组方向相反的最大剪应力区域; 随破裂扩展, 最大剪应力区域也跟着扩展, 显示出断层破裂时的扩展路径。同时, 图5中破裂起始点的位置与震间阶段应力积累模式非常吻合。

模型I显示, 在YBF约18km深度首先出现正、 负2组方向相反的最大剪应力区域, 即为初始发震点; 随后破裂沿断层同时向深处与浅表传播, 正负2区域的最大剪应力区域也不断扩展直至传播到整条断层, 最终断层完全破裂。模型Ⅱ 的初始发震点与模型I类似, 破裂传播跨断层传至WMF(左), 最终2条断层发生完全破裂。模型Ⅲ 显示初始发震点位于GJF(右)约15km深度, 即平缓断层到浅处陡峭断层的倾角过渡处, 随后破裂沿断层同时向深处与浅处传递, 在沿原断层扩展的同时, 破裂传播跨断层传至YBF(中), 最终2条断层破裂至地表, 逐步形成大地震, 但是同震期间WMF(左)始终未参与破裂。

图 5 3种模型同震剪应力分布的动态变化
模型Ⅰ 显示地震在YBF位于深处平缓断层到浅处陡峭断层过渡处约深度18km左右发震; 模型Ⅱ 破裂起始点同模型Ⅰ , 再传至WMF; 模型Ⅲ 起始破裂点位于GJF倾角过渡处, 接着传至YBFFig. 5 Snapshots of coseismic shear stress distribution in three models.

图 6 3种模型同震位移矢量分布的动态变化Fig. 6 Snapshots of coseismic displacement distribution in three models.

图6显示3种模型同震阶段位移矢量随破裂扩展的动态变化图, 用以描绘逆冲系统发震位移场的演化过程。箭头方向指示滑移方向, 箭头长短表示滑移距离大小; 上盘箭头方向向上, 下盘箭头方向向下, 指示断层逆冲运动。震间期断层处于闭锁状态, 当剪切应力超过断裂准则的阈值时, 首先在初始发震点滑移, 断层两侧运动趋势相反; 随破裂扩展, 形成断层上盘上升、 下盘相对下降的逆冲运动。在相同深度上, 断层面上盘的同震位移变形大于下盘, 与野外调查资料相符(徐锡伟等, 2008)。模型I中最大同震位移位于断层近地表附近。模型Ⅱ 显示YBF(右)与WMF(左)最大值都处于近地表附近。模型Ⅲ 显示YBF(中)与GJF(右)2条断层破裂完全, WMF(左)未发生破裂; 2条发震断层最大同震位移也出现在近地表。

由于断层面采用离散元法模拟, 破裂过程不像大部分数值模拟只降低一部分应力, 而是将模拟积累的Ma尺度的水平位移场1次性释放, 相当于数次大地震释放的总应力降, 得到1次大震数倍的滑移量。因此3种模型的剪应力数值大约为108~109Pa, 位移数值大约在100m量级, 根据剪应力与位移矢量大小估计发震造成的震级为模型Ⅱ > Ⅲ > I。虽然模拟结果量级偏大, 但是可以看成多次特征地震的叠加, 并不影响运动场性质的体现。同时也可以解释为什么龙门山地区海拔从约500m陡增至最高峰约5, 000m, 其中一部分隆升可能是由数ka来多次地震滑移叠加造成的。

3.3 震后期断层面上的应力残余

震后应力释放程度有助于判定断层的地震危险性。图7为3种模型震后期断层面上的剪应力与正应力残余量。从图中可以看出模型I、 Ⅱ 应力基本清零; 模型Ⅲ 2条发震断层的应力基本清零, 而WMF虽与震前阶段(图4)相比应力有所减小, 但其残余应力数值仍然是应力的高储能区; 这部分残余能量极可能以震后黏弹性形变或震后余滑方式释放。尽管如此, 由于WMF的震后残余应力明显高于其震前应力积累水平, 似乎预示其未来可能的地震危险性。

图 7 震后阶段断层面正应力与剪应力累积图Fig. 7 The plot of accumulation of normal and shear stress in post-earthquake stage.

4 讨论
4.1 多断层相互作用机制及汶川-茂汶断层未来地震危险性初步讨论

现有汶川地震的动力学研究仅仅涉及主破裂映秀-北川断层的研究模拟, 未涉及汶川-茂汶断层、 映秀-北川断层、 灌县-江油断层的空间组合关系与应力相互作用(Duan, 2010; 陶玮等, 2011; 朱守彪等, 2016)。我们通过动力学进行全周期演化模拟, 来初步探讨多断层逆冲系统的发震破裂模式。

地表调查结果表明, 汶川地震仅仅在龙门山断裂带的中央断裂和前山断裂即映秀-北川断层及灌县-江油断层产生大规模地表破裂, 而后山断裂即汶川-茂汶断层未发生明显破裂(徐锡伟等, 2008)。然而, 通过对比模拟结果, 发现汶川-茂汶断层并非不能发生破裂。模型Ⅱ 与Ⅲ 中的汶川-茂汶断层几何结构相同, 但是只有模型Ⅱ 中汶川-茂汶断层发生了破裂。模型Ⅱ 发震点位于映秀-北川断层(右)深部平缓与浅部陡峭过渡处, 破裂至地表倾角70° 左右区域, 接着破裂跨断层传播至汶川-茂汶断层(左), 能量释放集中在龙门山中央断裂和后山断裂; 而模型Ⅲ 初始发震点位于灌县-江油断层(右)深处平缓断层到浅处陡峭断层过渡处, 接着破裂传播跨断层传至映秀-北川断层(中), 能量释放集中在龙门山中央断裂和前山断裂。虽然模型Ⅱ 、 Ⅲ 最终都造成了2条断层破裂, 但破裂模式不同; 这2种不同的破裂模式显然与模型Ⅲ 中独有的灌县-江油断层的几何结构相关。可以认为破裂总是选择断层系统中最脆弱的区域开始破裂, 直至剩余未破裂部分的应力量级不足以使断层发生破裂为止。由于模型Ⅲ 中灌县-江油断层的几何结构更容易累积构造应力、 并且具有更弱的正应力, 使得它在汶川地震过程中率先破裂, 并快速传导至映秀-北川断层; 而汶川-茂汶断层虽然具有较高的剪切应力水平, 但断层面上超高的正应力阻止了破裂的进一步传播; 汶川地震的发生过程中, 汶川-茂汶断层震间期累积的构造应力在同震阶段来不及释放, 系统的整体应力水平就已经显著降低, 从而阻止了地震进一步扩展至汶川-茂汶断层。

从多断层模拟结果对比分析来看, 断层数量对于应力积累量及地震发生的早晚有很大影响, 但不能简单认为断层数量多就更容易触发地震或逆冲系统更不稳定。比如, 从震间期积累的应力来看, 模型Ⅱ 最大, 发震时间却最晚; 模型I次之, 模型Ⅲ 最小。可以认为, 模型Ⅲ 中的灌县-江油断层加剧了整个区域的不稳定性; 而模型Ⅱ 中的映秀-北川断层则使系统更加稳定; 这似乎再次印证了灌县-江油断层在汶川地震的孕震与破裂过程中所具备的至关重要的作用。可以推测, 模型断层数量与断层几何结构可能共同影响孕震及发震过程。

另外, 关于汶川-茂汶断裂未破裂的原因和其地震危险性是汶川地震后需要深入研究的问题。因为真实初始应力场无从知晓, 且汶川-茂汶断层的发震周期不一定就与龙门山其余断层一致。本研究假设3条断层的复发周期一致且同时开始加载积聚能量。动力学模型Ⅲ 重现了汶川地震中3条主断层的破裂情况, 破裂依次发生在灌县-江油断层和映秀-北川断层上, 汶川-茂汶断层未破裂, 原因可能在于汶川-茂汶断层整体倾角较前2条断层更陡, 且越靠近浅部倾角越大, 即整条断层都处于高倾角状态。震后阶段模型Ⅲ 断层面上的应力分布表明(图7), 较震前而言(图4), 汶川地震后汶川-茂汶断层可能仍然存在较大的地震危险性。

4.2 断层破裂顺序的讨论

现有汶川地震运动学反演结果描述的破裂模式大致可分为2种: 1)汶川地震在映秀-北川断层地下深度10~20km左右开始发震, 然后破裂至灌县-江油断层(Ji et al., 2008; 王卫民等, 2008; 张勇等, 2008; Nakamura et al., 2010; Zhang et al., 2011; Fielding et al., 2013; Zhang et al., 2014); 2)汶川地震在灌县-江油断层先发震, 再破裂至映秀-北川断层上(Shao et al., 2010; Hartzell et al., 2013), 即认为灌县-江油断层优先产生破裂, 然后传递至映秀-北川断层。

本研究模型Ⅲ 应力模拟结果与第2种更符合。但是, 由于我们无法全面了解地下的真实情况, 如是否存在障碍体或是凹凸体、 紫坪铺水库蓄放水是否导致岩石孔隙度和断层带处含水量的变化、 3条断层的地震复发周期是否一致以及对于地下断层结构和龙门山地区真实初始应力场的不完整认识等问题, 都会影响应力模拟结果及其发震破裂机制解释。本研究并未考虑以上这些问题, 对龙门山复杂的实际构造背景进行了简化, 且认为3条断层的复发周期一致并同时开始加载积聚能量。然而, 模型Ⅲ 模拟的是一定假设条件下的龙门山地区Ma尺度的地壳垂向增厚, 水平向缩短, 并最终导致汶川地震孕震及发震的情形, 基于模拟结果显示的是在岩石自重、 水平位移场与模型倾角相互作用下, 结合断裂准则得出的断层应力积累、 演化及释放过程; 显示了在一定深度低倾角断层浅部比深部更易破裂的力学机制与原因; 这一点与震源运动学反演获得的断层破裂顺序有本质差别。

4.3 断层倾角对于发震的影响

已有研究表明, 逆冲地震发震断层的倾角上限为50° ~60° (Sibson et al., 1998); 高倾角断层理论上难以发震, 然而映秀-北川断层上部倾角接近70° , 颠覆了我们对高倾角断层地震危险性的认识。关于映秀-北川断层整条断层全部破裂的原因, 我们认为主要是由于平缓的断层区域破裂驱动陡峭断层区域的扩展破裂。与Zhu等(2010, 2013)和Hartzell等(2013)等的观点类似, 模型I与Ⅱ 中初始发震的区域是映秀-北川断层深部缓倾角区域, 接着传至映秀-北川断层浅部高倾角处; 模型Ⅲ 中首先产生破裂扩展的区域是灌县-江油断层缓倾角区域, 接着传至映秀-北川深部缓倾角区域, 再传至高倾角处导致整条断层产生破裂。也就是说, 无论龙门山断裂带包含几条高倾角断层, 地震总是从倾角较为平缓的部分开始破裂, 随之向陡峭的断层上部传播。

虽然浅部的高倾角几何结构导致断层更加难以破裂, 但在震前积累的能量相对于缓倾角断层而言更高(Feng et al., 2016; 图4)。同时, 这也并不意味着高倾角断层破裂时释放的能量就更大, 2008年和2009年2次MW6.3大柴旦地震就证实了这一点(Eliiot et al., 2011; 温扬茂等, 2014; 温少妍等, 2016)。当然, 这种高倾角断层的破裂除了受到断层几何形状的影响外, 还强烈依赖于断层的介质性质和应力积累情况, 当浅部断层介质强度较弱时, 难以积累足够的应力, 即使存在深部断层破裂的推动, 浅部断层也难以破裂, 这也解释了俯冲带地震难以破裂至地表的原因(Hu et al., 2011)。

5 结论

本研究利用地震地质调查、 地震剖面探测结果对发震断层的几何形态及介质参数进行约束, 从而建立了龙门山断裂带的二维黏弹性有限元模型, 并基于连续-离散单元法对汶川地震的震间应力积累、 同震应力释放及震后应力恢复的全周期进行了模拟重现, 探讨了龙门山断裂带3条主断层之间的相互作用关系, 分析了汶川地震后汶川-茂汶断层的地震危险性, 主要结论如下:

(1)破裂总是选择断层系统中最脆弱的点开始, 直至剩余未破裂部分的应力不足以使断层发生破裂为止。由于模型中灌县-江油断层的几何结构更容易累积构造剪应力、 并且具有更弱的正应力, 使得它在汶川地震过程中率先开始破裂, 并快速传导至映秀-北川断层; 而汶川-茂汶断层虽然具有较高的剪切应力水平, 但断层面上超高的正应力阻止了破裂的进一步传播。

(2)由于汶川-茂汶断层的高倾角几何结构十分不利于同震破裂, 但却有利于震间能量的累积, 因而在汶川地震发生后仍然残留较高能量, 预示着汶川-茂汶断层可能存在较大的地震危险性。

(3)逆冲断层的几何形状呈铲型使得深部断层破裂造成的扰动推动原来处于闭锁状态的浅部高倾角区域发生被动破裂, 因而逆冲地震发震断层的倾角上限可能突破了传统认识的50° ~60° , 达到70° 。

致谢 感谢中国科学院非连续介质力学及工程灾害联合实验室提供的CDEM软件。

The authors have declared that no competing interests exist.

参考文献
[1] 陈九辉, 刘启元, 李顺成, 等. 2009. 汶川 MS8. 0地震余震序列重新定位及其地震构造研究[J]. 地球物理学报, 52(2): 390—397.
CHEN Jiu-hui, LIU Qi-yuan, LI Shun-cheng, et al. 2009. Seismotectonic study by relocation of the Wenchuan MS8. 0 earthquake sequence[J]. Chinese Journal of Geophysics, 52(2): 390—397(in Chinese). [本文引用:1]
[2] 邓起东, 陈社发, 赵小麟, 等. 1994. 龙门山及其邻区的构造和地震活动及动力学[J]. 地震地质, 16(4): 389—402.
DENG Qi-dong, CHEN She-fa, ZHAO Xiao-lin, et al. 1994. Tectonics, seismicity and dynamics of Longmenshan mountains and its adjacent regions[J]. Seismology and Geology, 16(4): 389—402(in Chinese). [本文引用:1]
[3] 黄福明. 2013. 断层力学概论 [M]. 北京: 地震出版社.
HUANG Fu-ming. 2013. Introduction to Fault Mechanics [M]. Seismological Press, Beijing(in Chinese). [本文引用:1]
[4] 陶玮, 胡才博, 万永革, 等. 2011. 铲形逆冲断层地震破裂动力学模型及其在汶川地震研究中的启示[J]. 地球物理学报, 54(5): 1260—1269.
TAO Wei, HU Cai-bo, WAN Yong-ge, et al. 2011. Dynamic modeling of thrust earthquake on listric fault and its inference to study of Wenchuan earthquake[J]. Chinese Journal of Geophysics, 54(5): 1260—1269(in Chinese). [本文引用:1]
[5] 王杰, 李世海, 张青波. 2015. 基于单元破裂的岩石裂纹扩展模拟方法[J]. 力学学报, 47(1): 105—118.
WANG Jie, LI Shi-hai, ZHANG Qing-bo. 2015. Simulation of crack propagation of rock based on splitting elements[J]. Chinese Journal of Theoretical and Applied Mechanics, 47(1): 105—118(in Chinese). [本文引用:1]
[6] 王杰, 李世海, 周东, 等. 2013. 模拟岩石破裂过程的块体单元离散弹簧模型[J]. 岩土力学, 34(8): 2355—2362.
WANG Jie, LI Shi-hai, ZHOU Dong, et al. 2013. A block-discrete-spring model to simulate failure process of rock[J]. Rock and Soil Mechanics, 34(8): 2355—2362(in Chinese). [本文引用:1]
[7] 王卫民, 赵连锋, 李娟, 等. 2008. 四川汶川8. 0级地震震源过程[J]. 地球物理学报, 51(5): 1403—1410.
WANG Wei-min, ZHAO Lian-feng, LI Juan, et al. 2008. Rupture process of the MS8. 0 Wenchuan earthquake of Sichuan, China[J]. Chinese Journal of Geophysics, 51(5): 1403—1410(in Chinese). [本文引用:1]
[8] 温少妍, 单新建, 张迎峰, 等. 2016. 基于InSAR的青海大柴旦地震三维同震形变场获取与震源特征分析[J]. 地球物理学报, 59(3): 912—921.
WEN Shao-yan, SHAN Xin-jian, ZHANG Ying-feng, et al. 2016. Three-dimensional co-seismic deformation of the Da Qaidam, Qinghai Earthquakes derived from D-InSAR data and their source features[J]. Chinese Journal of Geophysics, 59(3): 912—921(in Chinese). [本文引用:1]
[9] 温扬茂, 许才军, 李振洪, 等. 2014. InSAR约束下的2008年汶川地震同震和震后形变分析[J]. 地球物理学报, 57(6): 1814—1824.
WEN Yang-mao, XU Cai-jun, LI Zhen-hong, et al. 2014. Coseismic and postseismic deformation of the 2008 Wenchuan earthquake from InSAR[J]. Chinese Journal of Geophysics, 57(6): 1814—1824(in Chinese). [本文引用:1]
[10] 徐锡伟, 闻学泽, 叶建青, 等. 2008. 汶川 MS8. 0地震地表破裂带及其发震构造[J]. 地震地质, 30(3): 597—629. doi: DOI: 103969/j. issn. 0253-4967. 2008. 03. 003.
XU Xi-wei, WEN Xue-ze, YE Jian-qing, et al. 2008. The MS8. 0 Wenchuan earthquake surface ruptures and its seismogenic structure[J]. Seismology and Geology, 30(3): 597—629(in Chinese). [本文引用:1]
[11] 许志琴, 李化启, 侯立炜, 等. 2007. 青藏高原东缘龙门-锦屏造山带的崛起: 大型拆离断层和挤出机制[J]. 地质通报, 26(10): 1262—1276.
XU Zhi-qin, LI Hua-qi, HOU Li-wei, et al. 2007. Uplift of the Longmen-Jinping orogenic belt along the eastern margin of the Qinghai-Tibet Plateau: large-scale detachment faulting and extrusion mechanism[J]. Geological Bulletin of China, 26(10): 1262—1276(in Chinese). [本文引用:1]
[12] 袁海平, 曹平, 许万忠, 等. 2006. 岩石粘弹塑性本构关系及改进的Burgers蠕变模型[J]. 岩土工程学报, 28(6): 796—799.
YUAN Hai-ping, CAO Ping, XU Wan-zhong, et al. 2006. Visco-elastop-lastic constitutive relationship of rock and modified Burgers creep model[J]. Chinese Journal of Geotechnical Engineering, 28(6): 796—799(in Chinese). [本文引用:1]
[13] 张勇, 冯万鹏, 许力生, 等. 2008. 2008年汶川大地震的时空破裂过程[J]. 中国科学(D辑), 38(10): 1186—1194.
ZHANG Yong, FENG Wan-peng, XU Li-sheng, et al. 2009. Spatio-temporal rupture process of the 2008 great Wenchuan earthquake[J]. Science in China(Ser D), 52(2): 145—154. [本文引用:1]
[14] 朱守彪, 袁杰. 2016. 2008年汶川大地震单侧破裂过程的动力学机制研究[J]. 地球物理学报, 59(11): 4063—4074.
ZHU Shou-biao, YUAN Jie. 2016. Mechanisms for the fault rupture of the 2008 Wenchuan earthquake( MS=8. 0)with predominately unilateral propagation[J]. Chinese Journal of Geophysics, 59(11): 4063—4074(in Chinese). [本文引用:1]
[15] 朱守彪, 张培震. 2009. 2008年汶川 MS8. 0地震发生过程的动力学机制研究[J]. 地球物理学报, 52(2): 418—427.
ZHU Shou-biao, ZHANG Pei-zhen. 2009. A study on the dynamical mechanisms of the Wenchuan MS8. 0 earthquake, 2008[J]. Chinese Journal of Geophysics, 52(2): 418—427(in Chinese). [本文引用:1]
[16] Duan B C. 2010. Role of initial stress rotations in rupture dynamics and ground motion: A case study with implications for the Wenchuan earthquake[J]. Journal of Geophysical Research: Solid Earth, 115(B5): B05301. [本文引用:1]
[17] Elliott J R, Parsons B, Jackson J A, et al. 2011. Depth segmentation of the seismogenic continental crust: The 2008 and 2009 Qaidam earthquakes[J]. Geophysical Research Letters, 38(6): L06305. [本文引用:1]
[18] Feng S Y, Zhang P Z, Liu B J, et al. 2016. Deep crustal deformation of the Longmen Shan, eastern margin of the Tibetan plateau, from seismic reflection and Finite Element modeling[J]. Journal of Geophysical Research: Solid Earth, 121(2): 767—787. [本文引用:1]
[19] Fielding E J, Sladen A, Li Z H, et al. 2013. Kinematic Fault slip evolution source models of the 2008 M7. 9 Wenchuan earthquake in China from SAR interferometry, GPS and teleseismic analysis and implications for Longmen Shan tectonics[J]. Geophysical Journal International, 194(2): 1138—1166. [本文引用:1]
[20] Hartzell S, Mendoza C, Ramirez-Guzman L, et al. 2013. Rupture history of the 2008 MW7. 9 Wenchuan, China, earthquake: evaluation of separate and joint inversions of geodetic, teleseismic, and strong-motion data[J]. Bulletin of the Seismological Society of America, 103(1): 353—370. [本文引用:1]
[21] Hu Y. 2011. Deformation processes in great subduction zone earthquake cycles [D]. University of Victoria, Victoria. [本文引用:1]
[22] Nakamura T, Tsuboi S, Kaneda Y, et al. 2010. Rupture process of the 2008 Wenchuan, China earthquake inferred from teleseismic waveform inversion and forward modeling of broadband seismic waves[J]. Tectonophysics, 491(1-4): 72—84. [本文引用:1]
[23] Shao G, Ji C, Lu Z, et al. 2010. Slip history of the 2008 MW7. 9 Wenchuan earthquake constrained by jointly inverting Seismic and geodetic observations. 2010 Fall Meeting, AGU[C]. Abstract S52B-04. San Frand sco. [本文引用:1]
[24] Sibson R H, Xie G Y. 1998. Dip range for intracontinental reverse fault ruptures: truth not stranger than friction?[J]. Bulletin of the Seismological Society of America, 88(4): 1014—1022. [本文引用:1]
[25] Wang L X, Li S H, Zhang G X, et al. 2013. A GPU-based parallel procedure for nonlinear analysis of complex structures using a coupled FEM/DEM approach[J]. Mathematical Problems in Engineering, 618980: 15. [本文引用:1]
[26] Xu X W, Wen X Z, Yu G H, et al. 2009. Coseismic reverse- and oblique-slip surface faulting generated by the 2008 MW7. 9 Wenchuan earthquake, China[J]. Geology, 37(6): 515—518. [本文引用:1]
[27] Zhang G H, Qu C Y, Shan X J, et al. 2011. Slip distribution of the 2008 Wenchuan MS7. 9 earthquake by joint inversion from GPS and InSAR measurements: A resolution test study[J]. Geophysical Journal International, 186(1): 207—220. [本文引用:1]
[28] Zhang P Z, Wen X Z, Shen Z K, et al. 2010. Oblique, high-angle, listric-reverse faulting and associated development of strain: The Wenchuan earthquake of May 12, 2008, Sichuan, China[J]. Annual Review of Earth and Planetary Sciences, 38: 353—382. [本文引用:1]
[29] Zhang Y, Wang R J, Zschau J, et al. 2014. Automatic imaging of earthquake rupture processes by iterative deconvolution and stacking of high-rate GPS and strong motion seismograms[J]. Journal of Geophysical Research: Solid Earth, 119(7): 5633—5650. [本文引用:1]
[30] Zhu S B, Zhang P Z. 2010. Numeric modeling of the strain accumulation and release of the 2008 Wenchuan, Sichuan, China, Earthquake[J]. Bulletin of the Seismological Society of America, 100(5B): 2825—2839. [本文引用:1]
[31] Zhu S B, Zhang P Z. 2013. FEM simulation of interseismic and coseismic deformation associated with the 2008 Wenchuan earthquake[J]. Tectonophysics, 584: 64—80. [本文引用:1]