基于分裂节点法的地震同震和震后形变数值模拟及其在汶川大地震中的应用

2021-05-19孙云强黄禄渊

地质力学学报 2021年2期

孙云强，罗纲，黄禄渊

1.福建农林大学交通与土木工程学院，福建福州 350002；2.武汉大学测绘学院，湖北武汉 430079；3.应急管理部国家自然灾害防治研究院，北京 100085

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年M7.4 Landers地震造成的同震库仑应力变化，认为Landers地震触发了随后发生的M6.5 Big Bear地震，且Landers地震的余震序列大部分发生在库仑应力变化为正的区域。Stein et al.(1997)通过计算库仑应力变化成功解释了土耳其1939—1992年的地震序列。万永革等(2007)计算了青藏高原东北部的库仑应力演化，发现后发大地震大部分都发生在先发大地震导致的库仑应力为正的区域，触发率达85%。Freed and Lin(2001)分析了1992年Landers地震的震后库仑应力变化，发现由于中下地壳上地幔的黏弹性松弛造成的库仑应力变化在1999年M7.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)。

箭头表示地震断层错动方向；表示地震断层的同震位移；数字表示单元节点号

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

KU=F

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

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

1.2 三维黏弹性本构关系

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

{dε}={dεv}+{dεe}

(6)

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

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

{dεv}=Q-1{σt}dt=Q-1({σt-dt}+

{dσ})dt{dεe}=D-1{dσ}

(7)

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

(8)

(9)

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

2 概念性模型

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

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

白色箭头表示模型的断层为右旋走滑；AA′为剖面位置

图中白色线段表示断层位置；Ux，Uy，Uz分别表示x方向、y方向和z方向的同震位移分布

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

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

a—地震同震位移分布(实线表示有限元模型的结果，虚线为Okada模型计算得到的结果)；b—不同算例的地震同震位移；c—震后200年剖面同震位移变化图；d—同震和震后200年的位移总变化量

表1 单断层有限元模型的部分材料参数

在验证了程序可靠性的基础上，接下来文中将该程序应用于青藏高原东缘实例中；计算了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)。许多学者基于不同模型计算了汶川地震、芦山地震与九寨沟地震导致的应力变化及几次地震之间的相互关系(Parsonset 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)基于库仑应力变化及背景地震活动变化的结果认为，汶川地震明显地触发了芦山地震的发生，但却在一定程度上延迟了九寨沟地震的发生。可见不同学者得到的研究结果都显示了汶川地震使得芦山地震的库仑应力变化增加，即汶川地震触发了芦山地震的发生。而关于汶川地震是否触发九寨沟地震，不同研究者得到的结果仍存在差异。

图5 青藏高原东缘构造背景(震源机制解数据来自于GCMT；https://www.globalcmt.org/)

文中综合考虑活动构造、地球物理等多学科观测资料，建立青藏高原东缘的三维黏弹性有限元模型；使用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—汶川地震同震破裂模型

表2 青藏高原东缘有限元模型的参数设置

其他块体包括青藏高原东北缘、川滇块体和巴颜喀拉块体；上地壳的平均密度设置为2800 kg/m3，中下地壳上地幔的平均密度设置为3200 kg/m3

(10)

(11)

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

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

ΔCSC=Δτ+μ′Δσ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°)

另外，文中也计算了不同的有效摩擦系数对结果的影响。分别取有效摩擦系数为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)，定量地研究这种相互作用有助于更进一步地分析区域的地震危险性。

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