地质力学学报  2021, Vol. 27 Issue (2): 241-253
引用本文
孙云强, 罗纲, 黄禄渊. 基于分裂节点法的地震同震和震后形变数值模拟及其在汶川大地震中的应用[J]. 地质力学学报, 2021, 27(2): 241-253.
SUN Yunqiang, LUO Gang, HUANG Luyuan. Numerical simulation of coseismic and postseismic deformation through a node-splitting algorithm: A case study of the Wenchuan earthquake[J]. Journal of Geomechanics, 2021, 27(2): 241-253.
基于分裂节点法的地震同震和震后形变数值模拟及其在汶川大地震中的应用
孙云强1, 罗纲2, 黄禄渊3    
1. 福建农林大学交通与土木工程学院, 福建 福州 350002;
2. 武汉大学测绘学院, 湖北 武汉 430079;
3. 应急管理部国家自然灾害防治研究院, 北京 100085
摘要:大地震导致的同震及震后效应,对于分析不同地震之间的相互影响及区域地震危险性等有着重要的作用。文中开发了模拟地震同震及震后效应的三维黏弹性有限元程序,通过计算走滑断层震例(概念性模型)引起的同震及震后效应,并与解析/半解析解进行对比,验证了程序的可靠性。同时基于概念性模型,分析了不同介质参数对同震及震后的地表变形的影响。研究表明,地球介质的横向不均匀性对地震同震位移有显著的影响,而中下地壳上地幔的黏度对震后效应起着主要控制作用。最后将该程序应用于青藏高原东缘,计算分析了2008年MW7.9汶川大地震导致的同震及震后库仑应力变化对2013年MW6.6芦山地震及2017年MW6.5九寨沟地震的影响。结果显示,汶川地震导致的库仑应力变化在芦山地震震源附近(0.013 MPa)及九寨沟地震震源附近(0.009 MPa)都为正值,说明汶川地震可能使得两次地震提前发生。
关键词地震位错    震后形变    三维黏弹性    库仑应力    汶川地震    
DOI10.12090/j.issn.1006-6616.2021.27.02.023     文章编号:1006-6616(2021)02-0241-13
Numerical simulation of coseismic and postseismic deformation through a node-splitting algorithm: A case study of the Wenchuan earthquake
SUN Yunqiang1, LUO Gang2, HUANG Luyuan3    
1. College of Transportation and Civil Engineering, Fujian Agriculture and Forestry University, Fuzhou 350002, Fujian, China;
2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, Hubei, China;
3. National Institute of Natural Hazards, MEMC, Beijing 100085, China
Abstract: The coseismic and postseismic effects are crucial elements in analyzing fault interactions and the regional seismic risk. In this paper, we developed a three-dimensional viscoelastic finite element code to simulate the coseismic and postseismic deformation. We calculated the coseismic and postseismic deformations caused by the strike-slip fault with a conceptual model, and compared the results with the analytic and semi-analytic solutions so as to verify the reliability of the code. Meanwhile, we also analyzed the influence of different parameters on the coseismic and postseismic deformation, uncovering a significant effect of the earth's lateral heterogeneity on the coseismic displacement. And the viscosity of the middle-lower crust and upper mantle plays a major role in controlling the postseismic displacement. At last, we used the three-dimensional viscoelastic model to calculate the changes of the coseismic and postseismic Coulomb stress caused by the 2008 MW7.9 Wenchuan earthquake, and analyzed the subsequent effect on the 2013 MW6.6 Lushan earthquake and the 2017 MW6.5 Jiuzhaigou earthquake. The calculation results show that the Coulomb stress changes caused by the Wenchuan earthquake is positive near the hypocenters of the Lushan earthquake (0.013 MPa) and the Jiuzhaigou earthquake (0.009 MPa), indicating the Wenchuan earthquake might have triggered both the Lushan earthquake and Jiuzhaigou earthquake.
Key words: earthquake dislocation    postseismic deformation    three-dimensional viscoelasticity    Coulomb stress    Wenchuan earthquake    
0 引言

自1910年弹性回跳理论提出以来(Reid, 1910),人们逐渐认识到地震与断层之间的密切相关性。地震位错理论得到迅速的发展。1958年Steketee最早将弹性位错理论引入地震学,并推导出弹性半无限空间中点源位错的地表位移格林函数(Steketee, 1958)。随后,该理论被广泛应用于地震同震、震后应力应变分析中(Maruyama, 1964; Mansinha and Smylie, 1971; 陈运泰等,1979)。Okada(1985, 1992)在已有研究基础上,推导出均匀弹性半无限空间中三种位错源(走滑位错、倾滑位错、引张位错)引起的同震位移场及应力场的解析表达式,由于求解速度快,Okada模型已成为计算弹性半空间同震效应的经典方法。但是均匀弹性半无限空间模型不能够很好地处理地球介质的横向与纵向不均匀,在实际应用中仍具有较大局限性(邵志刚等,2008; 张贝等,2015; 黄禄渊等,2019)。另外,学者也逐渐认识到震后形变,即大地震发生后,地球介质的黏弹性松弛效应。为了使模型更符合真实的地球情况,同时考虑大地震同震及震后效应,研究者相继又开发了一系列的模型(Pollitz,1992; Sun,1992; Sabadini and Vermeersen, 1997; Wang et al., 2003; 黄禄渊等,2019)。Pollitz(1992)基于自由震荡简正模方法,计算了无重力情况下,黏弹性地球模型的地震位移场与应变场。Sun(1992)提出了球对称分层介质的位错理论。Wang et al.(2003)提出利用矩阵传播算法计算地震应力场格林函数,开发了黏弹性的垂向分层模型下的地震同震和震后形变计算程序PSGRN/PSCMP等。上述解析/半解析的方法计算速度快,但是在处理地球介质的横向不均匀性、复杂的断层面几何形状时仍然遇到了困难。随着计算机技术的发展,有限元数值模拟得到快速发展。有限元方法能够根据已有的观测资料建立更符合实际情况的地质模型,在地震同震及震后形变数值模拟研究中有着更广泛的应用(邵志刚等,2008; 李玉江等,2009; Zhu and Zhang, 2013; Xu et al., 2020)。

地震的发生是断层上应力释放和调整的过程,会导致区域应力场发生变化,进而影响区域的地震活动性。由于地球介质的复杂性和地球深部的不可入性(陈运泰, 2009),获取地壳的绝对应力值仍较为困难。但是基于弹性位错理论,可以计算每次地震产生的应力张量变化及库仑应力变化。大地震导致的区域最优破裂面及断层面上的库仑应力变化(同震库仑应力变化及震后的黏弹性松弛过程导致的库仑应力变化)作为一个重要的指标,在分析不同地震之间的相互关系、解释余震序列分布及区域地震危险性评估中都发挥了重要的作用(King et al., 1994; Stein et al., 1997; Harris and Simpson, 1998; Freed, 2005; 万永革等,2007; 胡才博和蔡永恩, 2016; 汪建军和许才军, 2017; 朱晓杰和何建坤, 2019; Xu et al., 2020)。例如,King et al.(1994)通过计算1992年M 7.4 Landers地震造成的同震库仑应力变化,认为Landers地震触发了随后发生的M 6.5 Big Bear地震,且Landers地震的余震序列大部分发生在库仑应力变化为正的区域。Stein et al.(1997)通过计算库仑应力变化成功解释了土耳其1939—1992年的地震序列。万永革等(2007)计算了青藏高原东北部的库仑应力演化,发现后发大地震大部分都发生在先发大地震导致的库仑应力为正的区域,触发率达85%。Freed and Lin(2001)分析了1992年Landers地震的震后库仑应力变化,发现由于中下地壳上地幔的黏弹性松弛造成的库仑应力变化在1999年M 7.1 Hector Mine地震震中附近增加了约0.1~0.2 MPa,从而触发了该地震的发生。朱晓杰和何建坤(2019)通过计算1970年通海MS7.7地震的同震及震后库仑应力变化,认为该地震使得小江断裂与红河断裂上的库仑应力增加,潜在地震活动性增强。

文中首先基于Abaqus的二次开发平台,开发了三维Maxwell(马克斯威尔体)黏弹性有限元程序,并引入分裂节点法来模拟地震位错。其次,通过计算走滑断层的地震同震及震后效应,并与解析/半解析解对比,验证了程序的正确性与可靠性。最后,将程序应用于青藏高原东缘,计算了2008年MW7.9汶川地震导致的同震应力变化及震后中下地壳上地幔的黏弹性应力调整,分析了汶川地震对2013年MW6.6芦山地震及2017年MW6.5九寨沟地震的影响。

1 计算方法 1.1 分裂节点法

地震同震位移可以处理为断层两侧大小相等、方向相反的两部分位移(在断层两侧同时施加地震同震位移一半的位移约束,且方向相反)。在有限元数值模拟中,研究者采用了多种不同方式来处理断层,包括采用一定宽度的弱化带(Luo and Liu, 2010)、摩擦接触(李玉江等,2013; Zhu and Zhang, 2013)、分裂节点法(Melosh and Raefsky, 1981; 邵志刚等,2008; 黄禄渊等,2018)等。其中分裂节点法能够简单有效地将断层位移以初始位移的形式引入有限元数值模拟中,是一种简单有效且物理意义明确的方法(Melosh and Raefsky, 1981)。因此文中拟使用分裂节点法来模拟断层地震。首先以一个包含两个单元的一维有限元模型为例,介绍分裂节点法的具体处理过程(图 1)。

箭头表示地震断层错动方向;U2a=-ΔU1b,表示地震断层的同震位移;数字表示单元节点号 The arrows show the direction of the fault slip. ΔU2a=-ΔU1b represents coseismic slip. The numbers represent the nodes 图 1 分裂节点法的一维有限元示意图 Fig. 1 1-Dimensional finite element model with the node-splitting algorithm

假设模型遵循静力平衡方程。在有限元分析中,平衡方程可以表示为:

$ \mathit{\boldsymbol{KU}} = \mathit{\boldsymbol{F}} $ (1)

其中,K为单元刚度矩阵;U为位移;F为等效节点荷载。

用节点位移表示的单元a的节点平衡方程为:

$ \left[ {\begin{array}{*{20}{c}} {k_{11}^{\rm{a}}}&{k_{12}^{\rm{a}}}\\ {k_{21}^{\rm{a}}}&{k_{22}^{\rm{a}}} \end{array}} \right]\;\left[ {\begin{array}{*{20}{c}} {U_1^{\rm{a}}}\\ {U_2^{\rm{a}} + \Delta U_2^{\rm{a}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {F_1^{\rm{a}}}\\ {F_2^{\rm{a}}} \end{array}} \right] $ (2)

将地震引起的位移量移动到方程式的右侧:

$ \left[ {\begin{array}{*{20}{c}} {k_{11}^{\rm{a}}}&{k_{12}^{\rm{a}}}\\ {k_{21}^{\rm{a}}}&{k_{22}^{\rm{a}}} \end{array}} \right]\;\left[ {\begin{array}{*{20}{c}} {U_1^{\rm{a}}}\\ {U_2^{\rm{a}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {F_1^{\rm{a}} - k_{12}^{\rm{a}}\Delta U_2^{\rm{a}}}\\ {F_2^{\rm{a}} - k_{22}^{\rm{a}}\Delta U_2^{\rm{a}}} \end{array}} \right] $ (3)

同理,可以得到单元b的节点平衡方程:

$ \left[ {\begin{array}{*{20}{c}} {k_{11}^{\rm{b}}}&{k_{12}^{\rm{b}}}\\ {k_{21}^{\rm{b}}}&{k_{22}^{\rm{b}}} \end{array}} \right]\;\left[ {\begin{array}{*{20}{c}} {U_1^{\rm{b}}}\\ {U_2^{\rm{b}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {F_1^{\rm{b}} - k_{12}^{\rm{b}}\Delta U_1^{\rm{b}}}\\ {F_2^{\rm{b}} - k_{22}^{\rm{b}}\Delta U_1^{\rm{b}}} \end{array}} \right] $ (4)

将单元刚度矩阵组装到全局,形成总体刚度矩阵,得到整个模型的节点平衡方程组:

$ \left[ {\begin{array}{*{20}{c}} {k_{11}^{\rm{a}}}&{k_{12}^{\rm{a}}}&0\\ {k_{21}^{\rm{a}}}&{k_{22}^{\rm{a}} + k_{11}^{\rm{b}}}&{k_{12}^{\rm{b}}}\\ 0&{k_{21}^{\rm{b}}}&{k_{22}^{\rm{b}}} \end{array}} \right]\;\left[ {\begin{array}{*{20}{c}} {{U_1}}\\ {{U_2}}\\ {{U_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{F_1} - k_{21}^{\rm{a}}\Delta U_2^{\rm{a}}}\\ {{F_2} - k_{22}^1\Delta U_2^{\rm{a}} - k_{11}^2\Delta U_1^{\rm{b}}}\\ {{F_3} - k_{21}^{\rm{b}}\Delta U_1^{\rm{b}}} \end{array}} \right] $ (5)

其中U3=U2b; F3=F2b

通过上式,可以看到:在不改变整体刚度矩阵的情况下,分裂节点法将地震位移以初始位移边界的形式加载到右端项,可以有效处理地震不连续面。

1.2 三维黏弹性本构关系

地震同震的持续时间较短,同震变形可以视为地球介质的弹性响应。而地震发生后,形变随着时间还将继续积累。研究指出这是由于不同机制,如震后余滑、黏弹性松弛等共同作用的影响(邵志刚等,2008),其中黏弹性松弛现象被广泛应用于地震的震后形变研究(Pollitz, 1992; Wang et al., 2001; Hu et al., 2004; 邵志刚等,2008; Shao et al., 2011)。因此,作为简化,文中只考虑黏弹性松弛的影响,采用Maxwell黏弹性材料来模拟中下地壳上地幔物质。模型采用增量法求解每一时间步的应变增量。应变增量由黏性应变增量和弹性应变增量组成:

$ \left\{ {{\rm{d}}\varepsilon } \right\} = \left\{ {{\rm{d}}{\varepsilon ^{\rm{v}}}} \right\} + \left\{ {{\rm{d}}{\varepsilon ^{\rm{e}}}} \right\} $ (6)

其中{dεv}和{dεe}分别表示黏性和弹性的应变增量;{}为张量形式。

三维黏弹性本构关系描述为:

$ \begin{array}{l} \left\{ {{\rm{d}}{\varepsilon ^{\rm{v}}}} \right\} = {\mathit{\boldsymbol{Q}}^{ - 1}}\left\{ {{\sigma ^t}} \right\}{\rm{d}}t = {\mathit{\boldsymbol{Q}}^{ - 1}}(\left\{ {{\sigma ^{t - {\rm{d}}t}}} \right\} + \\ \;\;\;\;\;\;\;\;\left\{ {{\rm{d}}\sigma } \right\}){\rm{d}}t\left\{ {{\rm{d}}{\varepsilon ^{\rm{e}}}} \right\} = {\mathit{\boldsymbol{D}}^{ - 1}}\left\{ {{\rm{d}}\sigma } \right\} \end{array} $ (7)

其中,{σt}为t时刻的应力张量;dt为时间增量;{dσ}为应力张量增量;D为弹性材料矩阵;Q为与黏度相关的材料矩阵。

$ \begin{array}{l} \mathit{\boldsymbol{D}} = \frac{E}{{\left({1 + v} \right)\left({1 - 2v} \right)}}\\ \left[ {\begin{array}{*{20}{c}} {1 - v}&v&v&0&0&0\\ v&{1 - v}&v&0&0&0\\ v&v&{1 - v}&0&0&0\\ 0&0&0&{0.5 - v}&0&0\\ 0&0&0&0&{0.5 - v}&0\\ 0&0&0&0&0&{0.5 - v} \end{array}} \right] \end{array} $ (8)
$ {\mathit{\boldsymbol{Q}}^{ - 1}} = \frac{1}{\eta }\left[ {\begin{array}{*{20}{c}} {\frac{1}{3}}&{ - \frac{1}{6}}&{ - \frac{1}{6}}&0&0&0\\ { - \frac{1}{6}}&{\frac{1}{3}}&{ - \frac{1}{6}}&0&0&0\\ { - \frac{1}{6}}&{ - \frac{1}{6}}&{\frac{1}{3}}&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1 \end{array}} \right] $ (9)

其中,E为杨氏模量; υ为泊松比; η为黏度。

2 概念性模型

根据上述理论,文中开发了模拟地震同震及震后效应的三维黏弹性有限元程序。首先以一个概念性模型(包含一条直立的走滑断层)为例(图 2),对程序的可靠性进行验证。模型尺寸设置为100 km×100 km×100 km。在垂向方向上模型介质分为两层:上层20 km为弹性层,模拟上地壳物质;杨氏模量和泊松比分别设置为8.75×1010 Pa和0.25。往下80 km为黏弹性层,模拟中下地壳上地幔;杨氏模量和泊松比分别为1.1×1011 Pa和0.25,黏度设置为1×1020 Pa·s(算例1)。模型的断层长度为20 km,断层延伸深度为20 km;通过在断层两侧各施加0.5 m的位移来模拟同震位移为1 m的右旋走滑地震(其中断层左侧施加0.5 m的同震位移,断层右侧施加-0.5 m的同震位移)。在模型的四个侧面边界施加滑移边界条件(法向固定,切线方向可以自由移动),底部边界固定。

白色箭头表示模型的断层为右旋走滑;AA′为剖面位置 The white arrows indicate a right-lateral strike-slip fault. AA′ is the profile location in Fig. 4. 图 2 走滑断层的有限元模型 Fig. 2 Conceptual finite element model with single strike-slip fault

模型计算的地表同震位移分布结果显示:同震位移分布表现出明显的关于断层成对称(图 3a)或反对称(图 3b3c)分布的花瓣图,与已有计算得到的结果也基本一致(邵志刚等,2008; 黄禄渊等,2018)。为了更进一步地验证程序的正确性,便于与解析/半解析解进行比较,文中选取跨断层的一个剖面AA′(图 2),并画出剖面节点上与断层走向平行的同震位移Uy(图 4a蓝色线段)。同时文中利用Okada模型(Okada, 1985)计算了平行断层走向方向的同震位移Uy(图 4a橙色虚线)。从图中可以看到有限元模型计算的结果与解析解得到的结果是非常吻合的(图 4a)。这也验证了文中所编写程序的正确性。

图中白色线段表示断层位置;UxUyUz分别表示x方向、y方向和z方向的同震位移分布
a—x方向同震位移分布图;b—y方向同震位移分布图;c—z方向同震位移分布图
The white line shows the location of fault in the model. Ux, Uy, and Uz represent the coseismic slip along the x-direction, y-direction, and z-direction, respectively. 图 3 走滑断层的地表同震位移分布图 Fig. 3 Surface coseismic slip distribution of the strike-slip fault. (a) Coseismic slip along the x-direction. (b) Coseismic slip along the y-direction. (c) Coseismic slip along the z-direction

a—地震同震位移分布(实线表示有限元模型的结果,虚线为Okada模型计算得到的结果);b—不同算例的地震同震位移;c—震后200年剖面同震位移变化图;d—同震和震后200年的位移总变化量 (a) Coseismic slip. The solid lines show the results of the finite element model. The dashed lines show the calculation results using the Okada model. (b) Coseismic slip of difference cases. (c) Postseismic slip after 200 years. (d) Total slip (Coseismic and postseismic slip after 200 years) on the profile AA′. 图 4 地震同震及震后在剖面AA′上沿着断层走向(y方向)的位移分布图(剖面AA′的位置见图 2) Fig. 4 Coseismic and postseismic along-strike slip on the profile AA′. The location of the profile AA′ in the model is shown in Fig. 2

另外,为了探讨地球介质横向不均匀性对地震同震位移的影响,文中在算例1的基础上又设置了两组不同的算例(表 1):算例2将位于断层左侧的上地壳杨氏模量设置为1.75×1011 Pa;算例3将位于断层左侧的上地壳杨氏模量设置为4.375×1010 Pa;其他参数保持不变。同样取剖面AA′上平行断层走向的同震位移的结果,如图 4b所示。可以看到,断层两侧杨氏模量的不同,会改变沿着断层走向方向的同震位移关于断层成反对称的趋势。断层左侧上地壳的杨氏模量增大,位于断层左侧上地壳的同震位移减小,而断层右侧上地壳的同震位移增大(图 4b紫色线);断层左侧上地壳的杨氏模量减小,位于断层左侧上地壳的同震位移增大,而断层右侧上地壳的同震位移减小(图 4b黄色线)。可见地球介质的横向不均匀对地震同震位移有显著的影响,在实际的应用中不容忽视。

表 1 单断层有限元模型的部分材料参数 Table 1 Material parameters of the finite element model with single fault

同时,文中分析了中下地壳上地幔的黏弹性松弛对上地壳形变的影响,以算例1为计算模型,计算剖面AA′的同震及震后200年累积的形变Uy(图 4c4d蓝色线),可以看到地震后上地壳的位移仍继续积累(图 4c4d)。为了进一步分析不同中下地壳上地幔黏度对震后效应的影响,在算例1的基础上,文中还分别设置了中下地壳上地幔黏度为1×1018 Pa·s(算例4)、1×1019 Pa·s(算例5)、1×1021 Pa·s(算例6)的模型,并计算不同算例的同震及震后形变(图 4c4d表 1)。通过对比不同的算例结果,可以看到,中下地壳上地幔的黏度对震后形变有显著的影响。并且在震后相同的时间段里,中下地壳上地幔的黏度越小,上地壳的震后形变持续积累越大;而中下地壳上地幔的黏度越大,上地壳的震后形变持续积累越小。这个结果(图 4)与相关学者用PSGRN/PSCMP程序计算得到的结果(徐昊等,2018)也比较一致。

在验证了程序可靠性的基础上,接下来文中将该程序应用于青藏高原东缘实例中;计算了2008年MW7.9汶川地震的同震及震后库仑应力变化,分析汶川地震对芦山地震及九寨沟地震的影响。

3 三维黏弹性有限元程序在青藏高原东缘的应用 3.1 青藏高原东缘地震活动

青藏高原东缘的川西地区是中国重要的地震活动带,区域构造变形强烈,活动断层展布,控制着一系列历史强震的发生(徐锡伟等,2005; 张培震等,2008; 郑文俊等,2020)。2008年5月12日在这个区域的龙门山断裂上发生了汶川MW7.9大地震,造成了巨大的人员伤亡和经济损失。随后在2013年4月20日发生了MW6.6芦山地震;2017年8月8日在九寨沟地区又发生了MW6.5九寨沟地震等(图 5)。如此频繁的地震也引起了社会和学界极大的关注。研究者在这个区域开展了大量的研究,包括地震破裂过程、地震复发周期等,积累了丰富的资料(王卫民等,2008; 张培震等,2008; Wang et al., 2011)。基于这些观测的数据资料建立模型,并计算地震导致的库仑应力变化,在分析不同地震之间的相互关系、解释余震序列分布及区域未来地震危险性等方面发挥了重要的作用(解朝娣等,2010; 贾科和周仕勇, 2018; Luo and Liu, 2018; 董培育等,2019; 黄禄渊等,2019)。许多学者基于不同模型计算了汶川地震、芦山地震与九寨沟地震导致的应力变化及几次地震之间的相互关系(Parsons et al., 2008; Toda et al., 2008; 单斌等,2009; Luo and Liu, 2018; 黄禄渊等,2019; Xu et al., 2020)。Parsons et al.(2008)采用地震波反演得到的均匀弹性有限断层模型,计算了汶川地震导致的库仑应力变化,结果显示汶川地震引起的同震库仑应力变化在龙门山断裂带南段增加了0.1 MPa,而该区域在2013年发生了MW6.6芦山地震。单斌等(2009, 2017)利用弹性位错理论计算,结果显示汶川地震造成芦山地震、九寨沟地震震中库仑应力变化为正,促进了这两次地震的发生等。黄禄渊等(2019)基于三维黏弹性有限元模型的计算结果认为,汶川地震使得芦山地震与九寨沟地震提前发生。贾科和周仕勇(2018)基于库仑应力变化及背景地震活动变化的结果认为,汶川地震明显地触发了芦山地震的发生,但却在一定程度上延迟了九寨沟地震的发生。可见不同学者得到的研究结果都显示了汶川地震使得芦山地震的库仑应力变化增加,即汶川地震触发了芦山地震的发生。而关于汶川地震是否触发九寨沟地震,不同研究者得到的结果仍存在差异。

(震源机制解数据来自于GCMT;https://www.globalcmt.org/) (Focal mechanism solutions from GCMT; https://www.globalcmt.org/) 图 5 青藏高原东缘构造背景 Fig. 5 Tectonic background of the eastern margin of the Tibetan Plateau.

文中综合考虑活动构造、地球物理等多学科观测资料,建立青藏高原东缘的三维黏弹性有限元模型;使用Ji and Hayes(2008)基于远场体波资料和有限断层反演得到的汶川地震同震滑动位移模型,数值模拟汶川地震造成的同震及震后应力场,分析汶川地震对芦山地震及九寨沟地震的影响。

3.2 青藏高原东缘有限元模型设置

文中建立了青藏高原东缘的三维黏弹性有限元模型(图 6)。模型东西向为800 km;南北向为1000 km;深度为100 km。根据活动块体划分(邓起东等,2003; 张培震等,2003),在横向上将模型划分为:青藏高原东北缘、巴颜喀拉块体、川滇块体、华南块体这4个块体(图 6)。在纵向上,采用半空间分层模型,划分为4层,包括上地壳、中地壳、下地壳和上地幔层,各分层的深度见表 2。模型的上地壳为弹性层;中地壳、下地壳及上地幔为黏弹性层。根据研究区的三维速度结构(https://igppweb.ucsd.edu/~gabi/crust2.html;吴建平等,2006; 王椿镛等,2008; 朱介寿, 2008)计算得到研究区各块体区域的平均介质参数(杨氏模量、泊松比):

a—青藏高原东缘的三维有限元模型;b—汶川地震同震破裂模型 图 6 青藏高原东缘的三维有限元模型及汶川地震同震破裂模型 Fig. 6 Three-dimensional finite element model of the eastern margin of the Tibetan Plateau (a) and coseismic rupture model of the Wenchuan earthquake (b).

表 2 青藏高原东缘有限元模型的参数设置 Table 2 Material parameters set for the finite element model of the eastern margin of the Tibetan Plateau
$ E = \frac{{V_{\rm{S}}^2\rho (3V_{\rm{P}}^2 - 4V_{\rm{S}}^2)}}{{V_{\rm{P}}^2 - V_{\rm{S}}^2}} $ (10)
$ \upsilon = \frac{{V_{\rm{P}}^2 - 2V_{\rm{S}}^2}}{{2\left({V_{\rm{P}}^2 - V_{\rm{S}}^2} \right)}} $ (11)

其中,E为杨氏模量;υ为泊松比;ρ为平均密度;VS为S波速度;VP为P波速度。各参数的设置如表 2所示。

中下地壳上地幔的黏度结构对震后形变有着重要的作用。许多学者基于不同的研究方法对该区域的中下地壳上地幔的黏度结构做了定量分析(石耀霖和曹建玲, 2008; 张晁军等,2008; Shao et al., 2011)。如Royden(1996)认为青藏高原地区的中下地壳较为柔软,并基于解析解,得到其黏度约为1×1020 Pa·s。Shao et al.(2011)基于二维黏弹性有限元模型分析了汶川地震的震后形变,认为川西地区的中下地壳黏度约为4×1017 Pa·s,四川盆地的中下地壳黏度约为1×1020 Pa·s。石耀霖和曹建玲(2008)利用实验室流变结果估算了中国大陆岩石圈的等效流变结构,其中青藏高原下地壳等效黏度在1×1019~1×1020 Pa·s之间;四川盆地中地壳黏度在1×1021~1×1023 Pa·s之间,下地壳黏度1×1021~1×1022 Pa·s之间等。虽然不同研究者得到的结果存在差异,但是相关研究都显示了青藏高原的中下地壳较为柔软,而华南地块的中下地壳硬度较大。综合已有研究成果,文中设置华南地块的中下地壳上地幔黏度为1×1022 Pa·s;其他块体(包括青藏高原东北缘、巴颜喀拉块体、川滇地块)的中下地壳黏度为1×1019 Pa·s,上地幔黏度为1×1020 Pa·s(表 2)。

另外文中假设地震对远场的影响较小,在整个有限元模型的四个侧面边界施加法向固定,切向自由的位移边界条件。底部边界为水平方向可以自由移动,垂向固定。

3.3 汶川地震导致的同震库仑应力变化

汶川地震发生后,研究者根据GPS观测、地震波形记录等资料反演了地震破裂的滑动位移分布(Ji and Hayes, 2008; 王卫民等,2008; Wang et al., 2011)。由于选取的观测数据或反演过程中的参数选择等不同,不同研究得到的结果存在一定的差别。解朝娣等(2010)通过对比5个不同的汶川地震震源模型导致的同震库仑应力变化,发现其中有两个模型计算得到的应力变化空间分布与余震的空间分布对应较好,Ji and Hayes(2008)研究得到的模型为其中之一。因此,文中使用Ji and Hayes(2008)得到的汶川地震同震破裂模型,在模型的断层处施加相应的同震位移约束,计算汶川地震导致的同震及震后库仑应力变化。库仑应力变化ΔCSC定义为:

$ \Delta CSC = \Delta \tau + \mu {'}\Delta {\sigma _{\rm{n}}} $ (12)

其中,Δτ为断层面上的剪切应力变化,Δσn为断层面上的法向应力变化,μ′为有效摩擦系数。基于地震位错理论可以计算模型每个节点的应力张量变化;再对应力张量进行空间坐标变换即可得到断层面的法向应力变化和剪切应力变化。μ′为有效摩擦系数,其取值范围通常在0~0.8之间(King et al., 1994)。ΔCSC大于0,表示岩石更趋近于破裂,断层的地震危险性增加;ΔCSC小于0,表示岩石更加远离破裂,断层的地震危险性减小。

地震导致的应力变化是一个张量。计算库仑应力变化时的应力投影方向对结果有很大的影响。文中采用将应力投影到后续发震断层面上的方式来计算库仑应力变化。有效摩擦系数取0.4,计算深度为10 km。另外,由于文中重点关注汶川地震对芦山地震及九寨沟地震的库仑应力影响,因此选取芦山地震与九寨沟地震的震源附近位置,并绘制其在汶川地震后的库仑应力随时间的演化(图 7)。图中显示汶川地震导致的库仑应力变化在芦山地震的震源附近为正值。汶川地震导致的库仑应力变化在芦山地震震源附近增加了约0.013 MPa(其中同震库仑应力变化增加了0.011 MPa;恰在芦山地震发生之前的震后库仑应力增加了0.002 MPa)。可见汶川地震导致的库仑应力变化在芦山地震震源附近已经超过地震触发的典型值0.01 MPa(King et al., 1994; Freed, 2005)。即2008年汶川地震造成的同震及震后库仑应力变化触发了芦山地震的发生。汶川地震导致的库仑应力变化在九寨沟地震的震源附近也为正值(0.009 MPa;图 7):其中,同震库仑应力变化在九寨沟地震震源附近增加了0.006 MPa;震后库仑应力增加了0.003 MPa。可见汶川地震导致的同震及震后库仑应力变化使得九寨沟地震提前发生。

a—汶川地震发生后芦山地震震源附近的应力随时间的变化;(应力投影方向为芦山地震震源机制,走向212°、倾角42°、滑动角100°);b—汶川地震发生后芦山地震震源附近及九寨沟地震震源附近的应力随时间的变化(应力投影方向为九寨沟地震震源机制,走向150°、倾角78°、滑动角-13°) (a) Stress evolution near the hypocenter of the Lushan earthquake after the Wenchuan earthquake. The stress is projected to the plane determined by the focal mechanism of the Lushan earthquake (strike: 212°, dip: 42°, rake: 100°). (b) Stress evolution near the hypocenter of the Jiuzhaigou earthquake after the Wenchuan earthquake. The stress is projected to the plane determined by the focal mechanism of the Jiuzhaigou earthquake (strike: 150°, dip: 78°, rake: -13°). 图 7 汶川地震发生后芦山地震震源附近及九寨沟地震震源附近应力随时间的变化 Fig. 7 Stress evolution near the hypocenter of the Lushan earthquake and that of the Jiuzhangou earthquake after the Wenchuan earthquake.

另外,文中也计算了不同的有效摩擦系数对结果的影响。分别取有效摩擦系数为0.0,0.4和0.8,并计算相应的库仑应力变化(图 7)。有效摩擦系数增大,也就是正应力在库仑应力变化中的权重增大,库仑应力增加的量值减小。反之,有效摩擦系数减小,正应力在库仑应力变化中的权重减小,库仑应力增加的量值增大(图 7)。但是图中的结果显示,库仑应力变化的趋势及正负并没有随着有效摩擦系数的不同而发生改变。

4 讨论

随着地震位错理论的进一步发展,现代地震学的很多研究都基于精细的地震位错理论(Sun, 1992; Wang et al., 2003; 王启欣等,2015)。实际的地球模型是非常复杂的;地球介质除了有垂向上的分层特征,横向上也有较大的不均匀性。有限元数值模拟可以建立更接近实际情况的地质模型,在处理地球介质的横向及纵向不均匀性都具有天然的优势(邵志刚等,2008; 黄禄渊等,2018)。在有限元数值模拟中,对断层地震位错的处理,Melosh and Raefsky(1981)提出了分裂节点法。该方法将地震引起的位移看作是外力中的一部分;在有限元方程中其他部分不做任何改变的情况下,将断层引起的地震引入到数值计算中,是一种简单有效且物理意义明确的方法(Melosh and Raefsky, 1981)。

文中分析了走滑断层实例的同震及震后效应;研究结果显示,地球介质的横向不均匀性对地震同震位移有显著的影响,而中下地壳上地幔的黏弹性松弛对震后效应起着主要控制作用。在实际的震例计算中,为了分析大地震引起的同震及震后效应,应该同时考虑地球介质的横向不均匀性及中下地壳上地幔的黏弹性松弛效应。最后,根据上述理论,文中将所开发的程序应用到青藏高原东缘,计算了2008年MW7.9汶川大地震导致的同震应力变化及震后应力调整。研究结果显示了汶川地震导致的库仑应力变化在芦山地震及九寨沟地震的震源附近都为正值(图 7),说明了汶川地震可能使得芦山地震与九寨沟地震都提前发生。文中计算的结果与大部分研究得到的结果是比较一致的(汪建军和许才军, 2017; 黄禄渊等,2019)。而也有部分研究指出汶川地震延迟了九寨沟地震的发生(贾科和周仕勇, 2018),这可能与计算地震库仑应力变化时选取的模型参数、应力投影方向等有关。但是相关研究都显示了地震之间的这种相互作用的存在。地震之间的相互作用是地震时空迁移与丛集的重要影响因素(Luo and Liu, 2018; Xu et al., 2020; 赵根模等,2020),定量地研究这种相互作用有助于更进一步地分析区域的地震危险性。

文中重点讨论了地震引起的静态的同震效应及中下地壳上地幔的黏弹性松弛导致的震后效应,在计算中没有考虑地震波的动态传播、震后余滑等因素的影响,更进一步的分析还需要结合多个方面共同作用的影响。另外,断层上的构造加载是地震孕育发生的主要因素。在下一步的工作中拟结合区域的构造加载来模拟历史地震序列的时空演化,为区域地震危险性分析提供参考。

5 结论

文中开发了模拟地震同震及震后效应的三维黏弹性有限元程序,并分析了地壳横向不均匀性、中下地壳上地幔的黏度对地震同震及震后效应的影响;通过建立青藏高原东缘的三维有限元模型,计算了汶川地震导致的同震及震后库仑应力变化,探讨了汶川地震对芦山地震和九寨沟地震的影响,得到以下结论。

(1) 地球介质的横向不均匀性对地震同震位移有显著的影响。如断层两侧的杨氏模量不一致,地震同震形变关于断层反对称的性质会发生改变;且杨氏模量大的一侧同震位移较小;杨氏模量小的一侧同震位移较大。

(2) 中下地壳上地幔的黏度对大地震震后形变起着主要控制作用。黏度越小,震后形变积累量越大;而黏度越大,震后形变积累量则越小。

(3) 汶川地震引起的库仑应力变化在芦山地震震源附近为0.013 MPa,在九寨沟地震震源附近为0.009 MPa。即汶川地震的发生有助于芦山地震和九寨沟地震提前发生。

参考文献/References
CHEN Y T, LIN B H, WANG X H, et al., 1979. A dislocation model of the Tangshan earthquake of 1976 from the inversion of geodetic data[J]. Acta Geophysica Sinica, 22(3): 201-217. (in Chinese with English abstract)
CHEN Y T, 2009. Earthquake prediction: Retrospect and prospect[J]. Science China Series D: Earth Sciences, 39(12): 1633-1658. (in Chinese with English abstract)
DENG Q D, ZHANG P Z, RAN Y K, et al., 2003. Active tectonics and earthquake activities in China[J]. Earth Science Frontiers, 10(S1): 66-73. (in Chinese with English abstract)
DONG P Y, LIU C, SHI Y L, 2019. Relationship between the 2008 MS8.0 Wenchuan earthquake and the 2017 MS7.0 Jiuzhaigou earthquake[J]. Journal of Geodesy and Geodynamics, 39(8): 777-782. (in Chinese with English abstract)
FREED A M, LIN J, 2001. Delayed triggering of the 1999 Hector Mine earthquake by viscoelastic stress transfer[J]. Nature, 411(6834): 180-183. DOI:10.1038/35075548
FREED A M, 2005. Earthquake triggering by static, dynamic, and postseismic stress transfer[J]. Annual Review of Earth and Planetary Sciences, 33: 335-367. DOI:10.1146/annurev.earth.33.092203.122505
HARRIS R A, SIMPSON R W, 1998. Suppression of large earthquakes by stress shadows: A comparison of Coulomb and rate-and-state failure[J]. Journal of Geophysical Research: Solid Earth, 103(B10): 24439-24451. DOI:10.1029/98JB00793
HU C B, CAI Y E, 2016. A possible mechanism of Omori-Utsu's law through an example of the great Tangshan earthquake[J]. Acta Seismologica Sinica, 38(4): 580-589. (in Chinese with English abstract)
HU Y, WANG K L, HE J, et al., 2004. Three-dimensional viscoelastic finite element model for postseismic deformation of the great 1960 Chile earthquake[J]. Journal of Geophysical Research: Solid Earth, 109(B12): B12403. DOI:10.1029/2004JB003163
HUANG L Y, WANG C H, YANG S X, 2018. Numerical modeling of earthquake dislocation using MSC. Marc[J]. North China Earthquake Sciences, 36(1): 1-9. (in Chinese with English abstract)
HUANG L Y, CHENG H H, ZHANG H, et al., 2019. Coseismic and postseismic stress evolution caused by the 2008 Wenchuan earthquake and its effects on the 2017 MS7.0 Jiuzhaigou earthquake[J]. Chinese Journal of Geophysics, 62(4): 1268-1281. (in Chinese with English abstract)
JI C, HAYES G. 2008. Preliminary result of the May 12, 2008 MW7.9 Eastern Sichuan, China earthquake[J/OL]. http://ji.faculty.geol.ucsb.edu/big_earthquakes/2008/05/12/Shichuan.html.
JIA K, ZHOU S Y, 2018. Triggering relationship in strong earthquake sequence around the Bayan Har block and its tectonic significance based on Coulomb stress changes and seismicity[J]. Acta Seismologica Sinica, 40(3): 291-303. (in Chinese with English abstract)
KING G C P, STEIN R S, LIN J, 1994. Static Stress Changes and the Triggering of Earthquakes[J]. Bulletin of the Seismological Society of America, 84(3): 935-953.
LI Y J, CHEN L W, YE J Y, 2009. Application and development of numerical simulation in stress field evolvement and seismology[J]. Progress in Geophysics, 24(2): 418-431. (in Chinese with English abstract)
LI Y J, CHEN L W, LU Y Z, et al., 2013. Numerical simulation on influences of Wenchuan earthquake on the stability of faults in the neighborhood[J]. Earth Science-Journal of China University of Geosciences, 38(2): 398-410. (in Chinese with English abstract) DOI:10.3799/dqkx.2013.039
LUO G, LIU M, 2010. Stress evolution and fault interactions before and after the 2008 Great Wenchuan earthquake[J]. Tectonophysics, 491(1-4): 127-140. DOI:10.1016/j.tecto.2009.12.019
LUO G, LIU M, 2018. Stressing rates and seismicity on the major faults in eastern Tibetan Plateau[J]. Journal of Geophysical Research: Solid Earth, 123(12): 10968-10986. DOI:10.1029/2018JB015532
MANSINHA L, SMYLIE D E, 1971. The displacement fields of inclined faults[J]. Bulletin of the Seismological Society of America, 61(5): 1433-1440.
MARUYAMA T, 1964. Statical elastic dislocations in an infinite and semi-infinite medium[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 42(2): 289-368.
MELOSH H J, RAEFSKY A, 1981. A simple and efficient method for introducing faults into finite element computation[J]. Bulletin of the Seismological Society of America, 71(5): 1391-1400.
OKADA Y, 1985. Surface deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 75(4): 1135-1154.
OKADA Y, 1992. Internal deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 82(2): 1018-1040.
PARSONS T, JI C, KIRBY E, 2008. Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan basin[J]. Nature, 454(7203): 509-510. DOI:10.1038/nature07177
POLLITZ F F, 1992. Postseismic relaxation theory on the spherical earth[J]. Bulletin of the Seismological Society of America, 82(1): 422-453.
REID H F, 1910. The Mechanics of the Earthquake. Vol. 2 of the California Earthquake of April 18, 1906[R]. Report of the State Earthquake Investigation Commission. Washington D C: Carnegie Institution of Washington, 2.
ROYDEN L, 1996. Coupling and decoupling of crust and mantle in convergent orogens: Implications for strain partitioning in the crust[J]. Journal of Geophysical Research: Solid Earth, 101(B8): 17679-17705. DOI:10.1029/96JB00951
SABADINI R, VERMEERSEN L L A, 1997. Influence of lithospheric and mantle stratification on global post-seismic deformation[J]. Geophysical Research Letters, 24(16): 2075-2078. DOI:10.1029/97GL01979
SHAN B, XIONG X, ZHENG Y, et al., 2009. Stress changes on major faults caused by MW7.9 Wenchuan earthquake, May 12, 2008[J]. Science in China Series D: Earth Sciences, 52(5): 593-601. DOI:10.1007/s11430-009-0060-9
ZHENG Y, LIU C L, et al., 2017. Coseismic Coulomb failure stress changes caused by the 2017 M7.0 Jiuzhaigou earthquake, and its relationship with the 2008 Wenchuan earthquake[J]. Science China Earth Sciences, 60(12): 2181-2189. DOI:10.1007/s11430-017-9125-2
SHAO Z G, FU R S, XUE T X, et al., 2008. The numerical simulation and discussion on mechanism of postseismic deformation after Kunlun MS8. 1 earthquake[J]. Chinese Journal of Geophysics, 51(3): 805-816. (in Chinese with English abstract)
SHAO Z G, WANG R J, WU Y Q, et al., 2011. Rapid afterslip and short-term viscoelastic relaxation following the 2008 MW7.9 Wenchuan earthquake[J]. Earthquake Science, 24(2): 163-175. DOI:10.1007/s11589-010-0781-z
SHI Y L, CAO J L, 2008. Effective viscosity of China continental lithosphere[J]. Earth Science Frontiers, 15(3): 82-95. (in Chinese with English abstract) DOI:10.1016/S1872-5791(08)60064-0
STEIN R S, BARKA A A, DIETERICH J H, 1997. Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering[J]. Geophysical Journal International, 128(3): 594-604. DOI:10.1111/j.1365-246X.1997.tb05321.x
STEKETEE J A, 1958. On volterra's dislocations in a semi-infinite elastic medium[J]. Canadian Journal of Physics, 36(2): 192-205. DOI:10.1139/p58-024
SUN W K, 1992. Potential and gravity changes caused by dislocations in spherically symmetric earth models[J]. Bulletin of the Earthquake Research Institute, University of Tokyo, 67(2): 89-238.
TODA S J, LIN J, MEGHRAOUI M, et al., 2008. 12 May 2008 M=7.9 Wenchuan, China, earthquake calculated to increase failure stress and seismicity rate on three major fault systems[J]. Geophysical Research Letters, 35(17): L17305. DOI:10.1029/2008GL034903
WAN Y G, SHEN Z K, ZENG Y H, et al., 2007. Evolution of cumulative Coulomb failure stress in northeastern Qinghai-Xizang (Tibetan) Plateau and its effect on large earthquake occurrence[J]. Acta Seismologica Sinica, 29(2): 115-129. (in Chinese with English abstract) DOI:10.1007/s11589-007-0117-9
WANG C Y, LOU H, LV Z Y, et al., 2008. S-wave crustal and upper mantle's velocity structure in the eastern Tibetan Plateau-Deep environment of lower crustal flow[J]. Science in China Series D: Earth Sciences, 51(2): 263-274. DOI:10.1007/s11430-008-0008-5
WANG J J, XU C J, 2017. Coseismic Coulomb stress changes associated with the 2017 MW6.5 Jiuzhaigou earthquake (China) and its impacts on surrounding major faults[J]. Chinese Journal of Geophysics, 60(11): 4398-4420. (in Chinese with English abstract)
WANG K L, HE J H, DRAGERT H, et al., 2001. Three-dimensional viscoelastic interseismic deformation model for the Cascadia subduction zone[J]. Earth, Planets and Space, 53(4): 295-306. DOI:10.1186/BF03352386
WANG Q, QIAO X J, LAN Q G, et al., 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan[J]. Nature Geoscience, 4(9): 634-640. DOI:10.1038/ngeo1210
WANG Q X, JIANG Z S, WU Y Q, et al., 2015. A review on comparison and progress in applications of earthquake dislocation theories based on different models[J]. Acta Seismologica Sinica, 37(4): 690-704. (in Chinese with English abstract)
WANG R J, MARTÍN F L, ROTH F, 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP[J]. Computers & Geosciences, 29(2): 195-207.
WANG W M, ZHAO L F, LI J, et al., 2008. Rupture process of the MS 8.0 Wenchuan earthquake of Sichuan, China[J]. Chinese Journal of Geophysics, 51(5): 1403-1410. (in Chinese with English abstract)
WU J P, MING Y H, WANG C Y, 2006. Regional waveform inversion for crustal and upper mantle velocity structure below CHuandian region[J]. Chinese Journal of Geophysics, 49(5): 1369-1376. (in Chinese with English abstract)
XIE Z D, ZHU Y Q, LEI X L, et al., 2010. Pattern of stress change and its effect on seismicity rate caused by MS8.0 Wenchuan earthquake[J]. Science China Earth Sciences, 53(9): 1260-1270. DOI:10.1007/s11430-010-4025-9
XU D Y, XIAO J, HE J K, et al., 2020. Strong earthquake clustering around the eastern Tibetan Plateau after the 2008 MW7.9 Wenchuan earthquake[J]. Science China Earth Sciences, 63(7): 999-1012. DOI:10.1007/s11430-019-9581-x
XU H, SUN Y J, WU Z H, 2018. The effect of lithospheric structure on the seismic deformation-taking the 1976 Tangshan earthquake and 2001 Kunlunshan earthquake as an example[J]. Chinese Journal of Geophysics, 61(8): 3170-3184. (in Chinese with English abstract)
XU X W, ZHANG P Z, WEN X Z, et al., 2005. Features of active tectonics and recurrence behaviors of strong earthquakes in the western Sichuan Province and its adjacent regions[J]. Seismology and Geology, 27(3): 446-461. (in Chinese with English abstract)
ZHANG B, ZHANG H, SHI Y L, 2015. Equivalent-bodyforce approach on modeling elastic dislocation problem using finite element method[J]. Chinese Journal of Geophysics, 58(5): 1666-1674. (in Chinese with English abstract)
ZHANG C J, CAO J L, SHI Y L, 2008. Studying the viscosity of lower crust of Qinghai-Tibet Plateau according to post-seismic deformation[J]. Science in China Series D: Earth Sciences, 52(3): 411-429.
ZHANG P Z, DENG Q D, ZHANG G M, et al., 2003. Active tectonic blocks and strong earthquakes in the continent of China[J]. Science in China Series D: Earth Sciences, 46(2): 13-24.
ZHANG P Z, XU X W, WEN X Z, et al., 2008. Slip rates and recurrence intervals of the Longmen shan active fault zone, and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China[J]. Chinese Journal of Geophysics, 51(4): 1066-1073. (in Chinese with English abstract)
ZHAO G M, WU Z H, LIU J, 2020. The types, characteristics and mechanism of seismic migration[J]. Journal of Geomechanics, 26(1): 13-32. (in Chinese with English abstract)
ZHENG W J, ZHANG P Z, YUAN D Y, et al., 2020. Basic characteristics of active tectonics and associated geodynamic processes in continental China[J]. Journal of Geomechanics, 25(5): 699-721. (in Chinese with English abstract)
ZHU J S, 2008. The Wenchuan earthquake occurrence background in deep structure and dynamics of lithosphere[J]. Journal of Chengdu University of Technology (Science and Technology Edition), 35(4): 348-356. (in Chinese with English abstract)
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. DOI:10.1016/j.tecto.2012.06.024
ZHU X J, HE J K, 2019. Coulomb stress change on the Xiaojiang and the Red River faults, southeastern Tibetan Plateau, from the 1970 Tonghai MS7.7 earthquake[J]. Journal of Geodesy and Geodynamics, 39(12): 1223-1227. (in Chinese with English abstract)
陈运泰, 林邦慧, 王新华, 等, 1979. 用大地测量资料反演的1976年唐山地震的位错模式[J]. 地球物理学报, 22(3): 201-217. DOI:10.3321/j.issn:0001-5733.1979.03.001
陈运泰, 2009. 地震预测: 回顾与展望[J]. 中国科学D辑: 地球科学, 39(12): 1633-1658.
邓起东, 张培震, 冉勇康, 等, 2003. 中国活动构造与地震活动[J]. 地学前缘, 10(S1): 66-73.
董培育, 柳畅, 石耀霖, 2019. 2008年汶川MS8.0地震与2017年九寨沟MS7.0地震成因关系探讨[J]. 大地测量与地球动力学, 39(8): 777-782.
胡才博, 蔡永恩, 2016. 大森-宇津定律的一种可能机制: 以唐山大地震为例[J]. 地震学报, 38(4): 580-589.
黄禄渊, 王成虎, 杨树新, 2018. 地震位错模拟在有限元软件MSC. Marc中的实现[J]. 华北地震科学, 36(1): 1-9. DOI:10.3969/j.issn.1003-1375.2018.01.001
黄禄渊, 程惠红, 张怀, 等, 2019. 2008年汶川地震同震-震后应力演化及其对2017年九寨沟MS7.0地震的影响[J]. 地球物理学报, 62(4): 1268-1281.
贾科, 周仕勇, 2018. 基于库仑应力改变和地震活动性研究巴颜喀拉块体周缘强震序列的触发关系及其构造意义[J]. 地震学报, 40(3): 291-303.
李玉江, 陈连旺, 叶际阳, 2009. 数值模拟方法在应力场演化及地震科学中的研究进展[J]. 地球物理学进展, 24(2): 418-431. DOI:10.3969/j.issn.1004-2903.2009.02.007
李玉江, 陈连旺, 陆远忠, 等, 2013. 汶川地震的发生对周围断层稳定性影响的数值模拟[J]. 地球科学-中国地质大学学报, 38(2): 398-410.
单斌, 熊熊, 郑勇, 等, 2009. 2008年5月12日MW7.9汶川地震导致的周边断层应力变化[J]. 中国科学D辑: 地球科学, 39(5): 537-545.
单斌, 郑勇, 刘成利, 等, 2017. 2017年M7.0级九寨沟地震同震库仑应力变化及其与2008年汶川地震的关系[J]. 中国科学: 地球科学, 47(11): 1329-1338.
邵志刚, 傅容珊, 薛霆虓, 等, 2008. 昆仑山MS8.1级地震震后变形场数值模拟与成因机理探讨[J]. 地球物理学报, 51(3): 805-816. DOI:10.3321/j.issn:0001-5733.2008.03.021
石耀霖, 曹建玲, 2008. 中国大陆岩石圈等效粘滞系数的计算和讨论[J]. 地学前缘, 15(3): 82-95. DOI:10.3321/j.issn:1005-2321.2008.03.006
万永革, 沈正康, 曾跃华, 等, 2007. 青藏高原东北部的库仑应力积累演化对大地震发生的影响[J]. 地震学报, 29(2): 115-129. DOI:10.3321/j.issn:0253-3782.2007.02.001
王椿镛, 楼海, 吕智勇, 等, 2008. 青藏高原东部地壳上地幔S波速度结构: 下地壳流的深部环境[J]. 中国科学D辑: 地球科学, 38(1): 22-32. DOI:10.3321/j.issn:1006-9267.2008.01.003
汪建军, 许才军, 2017. 2017年MW6.5九寨沟地震激发的同震库仑应力变化及其对周边断层的影响[J]. 地球物理学报, 60(11): 4398-4420. DOI:10.6038/cjg20171127
王启欣, 江在森, 武艳强, 等, 2015. 不同模型下地震位错理论的对比及其应用进展综述[J]. 地震学报, 37(4): 690-704.
王卫民, 赵连锋, 李娟, 等, 2008. 四川汶川8.0级地震震源过程[J]. 地球物理学报, 51(5): 1403-1410. DOI:10.3321/j.issn:0001-5733.2008.05.013
吴建平, 明跃红, 王椿镛, 2006. 川滇地区速度结构的区域地震波形反演研究[J]. 地球物理学报, 49(5): 1369-1376. DOI:10.3321/j.issn:0001-5733.2006.05.016
解朝娣, 朱元清, LEI X L, 等, 2010. MS8.0汶川地震产生的应力变化空间分布及其对地震活动性的影响[J]. 中国科学: 地球科学, 40(6): 688-698.
徐锡伟, 张培震, 闻学泽, 等, 2005. 川西及其邻近地区活动构造基本特征与强震复发模型[J]. 地震地质, 27(3): 446-461. DOI:10.3969/j.issn.0253-4967.2005.03.010
徐昊, 孙玉军, 吴中海, 2018. 岩石圈结构对大地震震后形变的影响: 以1976年唐山大地震和2001年昆仑山大地震为例[J]. 地球物理学报, 61(8): 3170-3184.
张贝, 张怀, 石耀霖, 2015. 有限元模拟弹性位错的等效体力方法[J]. 地球物理学报, 58(5): 1666-1674.
张培震, 邓启东, 张国民, 等, 2003. 中国大陆的强震活动与活动地块[J]. 中国科学(D辑), 33(S1): 12-20.
张培震, 徐锡伟, 闻学泽, 等, 2008. 2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因[J]. 地球物理学报, 51(4): 1066-1073. DOI:10.3321/j.issn:0001-5733.2008.04.015
张晁军, 曹建玲, 石耀霖, 2008. 从震后形变探讨青藏高原下地壳黏滞系数[J]. 中国科学D辑: 地球科学, 38(10): 1250-1257. DOI:10.3321/j.issn:1006-9267.2008.10.008
赵根模, 吴中海, 刘杰, 2020. 地震迁移的类型、特征及机制讨论[J]. 地质力学学报, 26(1): 13-32.
郑文俊, 张培震, 袁道阳, 等, 2020. 中国大陆活动构造基本特征及其对区域动力过程的控制[J]. 地质力学学报, 25(5): 699-721.
朱介寿, 2008. 汶川地震的岩石圈深部结构与动力学背景[J]. 成都理工大学学报(自然科学版), 35(4): 348-356. DOI:10.3969/j.issn.1671-9727.2008.04.002
朱晓杰, 何建坤, 2019. 1970年通海MS7.7地震后青藏高原东南部小江断裂带和红河断裂带的库仑应力变化研究[J]. 大地测量与地球动力学, 39(12): 1223-1227.