强震孕育、发生及其复发循环过程的有限单元法模拟

袁杰， 崔泽飞， 朱守彪*， 王进廷

1 应急管理部国家自然灾害防治研究院， 北京 100085 2 清华大学水利水电工程系， 北京 100084

0 引言

定量地模拟地震、孕育发生及其复发循环过程对于深入认识地震及地震预测预报等具有非常重要的科学意义.但是，地震从孕育到发生是个十分复杂的力学过程.通常，强震孕育需要数十年、上百年甚至几千年的时间(朱守彪和张培震，2009；Zhu and Zhang, 2010, 2013)，而地震的发生则只需数十秒；此外，由于地震孕育过程缓慢、变形速率小，该过程可视为准静态过程，可用静态力学方程来求解；但地震发生时，断层产生错动，发射地震波，地面产生强烈震动，只能用动力学过程来描述.在地震、孕育发生的模拟过程中，不仅涉及力学状态的改变，同时，如何来选择计算中的时间步长，也面临着更大的困难：若按照漫长的地震孕育过程来选择时间步长(如：步长为半年或一个月)，则无法模拟短暂的地震发生过程(数十秒的范围)；若以地震发生过程来选择时间步长(如：0.0001 s)，则无法模拟漫长的地震孕育过程(数千年).因此，模拟过程中，还必须改变计算的时间步长.

为克服上述模拟中的困难，不少研究人员进行了许多探索：利用静态计算方法模拟断层的孕震过程，利用动态计算方法继续模拟断层的同震破裂过程.在断层处于临界状态时，将静态过程的计算结果作为初始条件，并改变计算方法(力学状态)从而继续计算，进而连续模拟孕震-同震过程 (Duan and Oglesby 2005a,b).然而，人为改变计算过程中的力学状态会将应力、位移、摩擦本构等之间的复杂耦合关系解耦，进而会干扰断层失稳演化的自然力学过程；同时，此方法也无法精确的判断断层的失稳临界点.针对这种情况，Lapusta等(2000)在前人的研究之上(Tse and Rice, 1986; Rice and Ben-Zion, 1996; Ben-Zion and Rice, 1997),基于边界积分方程方法，开发了一种可连续模拟断层的孕震-同震及交替循环的计算方法.随后又将此方法的应用范围由二维扩展到三维 (Lapusta and Liu, 2009).特别值得指出的是，利用该方法， Barbot等(2012)模拟了San Andreas 断层的Parkfield段的强震孕育发生过程，模拟结果再现了Parkfield断层段上1966年至2004年间发生的MW6.0以上的所有地震事件，并且模拟给出的每次地震事件的发生时刻、成核位置、震级大小都与实际记录的观测资料很好的吻合.可以推测，若利用该方法继续计算，则可以预测出San Andreas断层的Parkfield段下一次大地震发生的时、空、强三要素；从而实现数值地震预测预报.

然而，边界积分方程方法使问题的维数降低一维，例如：三维问题降为二维问题，二维降为一维问题，使得求解的自由度下降.同时，边界元仅适应规则区域及边界条件的模型，适应于求解线性、均匀材料属性的计算问题.这样，在处理实际的地质问题时(如：复杂的几何边界、复杂的物性结构等)会遇到很大的困难.尽管对边界积分方程方法进行了不断地优化(如：Sathiakumar et al., 2020)，但目前进展不大.因此，对于复杂情况的实际问题，模拟断层的孕震-同震及其循环过程，最好还是选择能够适应复杂几何形状、复杂边界条件、复杂的材料属性的有限单元方法.虽然很多研究人员(如：Duan and Oglesby, 2005a, b; 朱守彪等，2008；朱守彪和张培震, 2009; Zhu and Zhang, 2010, 2013; Zhu, 2013; 马林飞等，2018)利用有限单元法模拟了断层的发震周期及断层破裂过程，但没有很好地连续模拟断层孕震-同震(准静态-动态)及其循环的全过程.

1 计算原理及有限元模型

1.1 动力学方程

利用有限元方法研究断层自发破裂动力学过程问题，即整个区域内求解二维波动方程(王勖成，2003)：

σij,j+fi-ρui,tt-cui,t=0.

(1)

(1)式又称为平衡方程.其中σij是应力张量，fi是体力，ρ是弹性介质的密度，c是阻尼系数，ui表示位移，ui,t和ui,tt分别是ui对t的一次导数和二次导数，即分别是速度和加速度；-ρui,tt和-cui,t分别代表惯性力和阻尼力.具体的动力学基本方程和有限元方法见我们之前的文章(袁杰和朱守彪，2014)，以及相关有限元理论文献(王勖成，2003)，此处不再赘述.

1.2 摩擦本构关系

断层上的摩擦关系非常复杂，研究人员通过拟合实验室实验数据，参照断层破裂运动学反演结果，提出了数十种断层面上的本构关系假设(Bizzarri,2011).其中，运用较为广泛的是速率-状态相关摩擦的关系和滑移弱化摩擦的关系.虽然利用这些摩擦本构关系模拟断层的破裂传播过程，其模拟结果也能很好解释一些地震现象，但使用不同的摩擦本构关系，其模拟结果也不尽相同(Ryan, 2012; Ryan and Oglesby, 2014; Yuan et al., 2020)，到底哪一种摩擦本构关系最贴合断层的实际状况，目前还没有定论.本文的主要研究目的是模拟断层孕震-同震准周期循环过程，因此本文中所有模型拟使用形式简单且易于实现、得到非常广泛的应用的滑移弱化的摩擦本构关系.

在经典的滑移弱化摩擦本构关系中，当断层两侧的相对滑移距离达到特征滑移距离后，断层上的摩擦系数将一直不变，维持动摩擦系数状态，导致一次地震事件之后，断层无法继续积累能量，也就无法产生下一次地震(Duan and Oglesby, 2005a).为了能连续模拟出断层孕震-同震周期循环过程，我们采用改进后的滑移弱化摩擦本构关系(Olsen-Kettle et al., 2008; 袁杰和朱守彪,2014; Yuan et al., 2020):

(2)

利用改进后的滑移弱化摩擦本构关系模拟断层破裂过程时，在一次地震事件完成之后，断层的摩擦系数将会再次强化，由动摩擦系数变为静摩擦系数.具体在计算中实现的做法是：当断层面上对应两节点之间相对滑移速度降为零，即停止滑移时，将摩擦系数变为静摩擦系数.为保证计算过程的稳定、自洽，本文中的摩擦系数再次强化是一个缓慢平滑的过程，而不是直接由动摩擦系数重置为静摩擦系数.

1.3 动力学方程求解方法

显式算法和隐式算法是利用有限单元法求解动力学问题的两种主要算法，其主要对应的中心差分法和Newmark法，相关的理论介绍见有限元理论文献(王勖成，2003).表1简略介绍了两种算法的优缺点.

模型中，i和t分别表示不同省份和年份；C1、C2、C4分别表示投资、技术创新和人口对产业结构升级的作用效应；C3表示非老龄人口比例对产业结构升级的影响，当C3＞0时，表示老龄化程度越高，越不利于产业结构升级；C5表示分配给劳动人口的资源比例对产业结构升级的影响，当C5＜0时，表示分配给老龄人口的资源越多，越能促进产业结构升级。

表1 显式算法和隐式算法的优缺点Table 1 Advantages and disadvantages of explicit and implicit algorithms

两种算法都可以用来模拟断层的自发破裂过程(同震过程)，但是目前研究中主要采用的是显式算法.这是由于断层的瞬时破裂过程是高度非线性的问题，采用隐式算法是需要很小的时间步长能收敛，这时就体现不出隐式算法的优点(时间步长Δt不影响解的稳定性).实际对比两种算法模拟同震过程的结果时，我们发现：只有当隐式算法的时间步长小于0.001 s时，其计算精度才能达到显式算法模拟结果.然而时间步长同样采用0.001 s时，隐式算法计算同震过程的求解时间是显式算法的数十倍，所以单纯模拟同震过程时，主要还是采用显式算法.

但是在模拟断层孕震过程时，由于显式算法需要满足稳定条件，即时间步长Δt必须小于Δtcr，基本上不可能模拟几十年甚至上千年的孕震过程.因此本文基于Newmark隐式时间积分法，发展出一种时间步长可以随着断层状态自动平滑缩放的计算方法，进而模拟同时包含孕震-同震两种状态的地震复发模型.

1.4 模型几何及有限元网格剖分

为简单起见，本文将实际逆冲断层的三维情况抽象为二维模型，并将断层简化为直线.图1为研究所用的模型几何及边界条件：模型空间尺度为100 km×40 km的矩形，断层的倾向长度为40 km，倾角为30°.模型底部边界条件是水平方向自由运动，垂直方向固定；右侧边界条件是：垂直方向上自由运动，水平方向上固定；左侧边界条件为施加15 mm·a-1的水平构造载荷；于自由地表至地下2 km处施加重力加速度.此外，为防止地震波通过边界反射而影响结果，模型在四周设置了吸收边界.

图1 模型几何及边界条件Fig.1 Model geometry and boundary conditions

模型全部采用三角形单元来剖分.断层附近为研究的重点区域，为了保证精度，在该区域对网格进行细化，单元边长为100 m，此外，离断层越远地方的单元尺度越大，模型最外围部分单元的边长为500 m.模型中单元的节点数和总数分别为99714和196955.模型中的材料参数、摩擦系数等见表2.

表2 模型参数Table 2 Model parameters

2 模拟结果

模拟孕震-同震准周期循环过程最大的难点是精准地找到孕震-同震之间的临界点，即断层的临界失稳状态.此外，模拟程序中时间步长如何随着断层的滑动状态做出相应的自动调整也是此研究中的难点之一.本文利用新的有限单元计算方法可以很好的解决上述难点，进而连续稳定的模拟断层孕震-同震及其循环过程，模拟结果如下.

图2展示了断层上盘上的典型点(中间点)的位移随时间的变化曲线.随着构造载荷的不断挤压，孕震期间典型点的位移随时间呈线性增加，在整个计算的2000年时间里，该点的位移发生了5次跳跃式的突然增加.位移的每次跳跃式增加，表示发生了一次地震事件，两次地震事件之间的时间间隔即为孕震时间.图3是断层典型点位置上下盘之间的位错随时间的变化曲线.可以看出：(1)在时间较长的孕震期间断层两盘对应点之间位错变化缓慢，而在时间很短同震期间位错突然增加.(2)孕震时间与同震位错近似呈正比关系，当孕震时间较长时，紧接着到来的地震也会相对较大.这与前人(马林飞等,2018)利用有限单元法模拟断层准周期行为的结果类似.但模拟方法不同，本文采用方法中的时间步长可以自动随着断层状态做出相应的调整，还可以模拟断层的瞬时破裂过程.

图2 断层上盘上的典型点(中间点)的位移随时间的变化曲线Fig.2 Displacements change with time at the typical point on the hanging wall of the fault

图3 断层典型点位置上下盘之间的位错随时间的变化曲线Fig.3 Sliding distance at the typical point on the fault vary with time

实现时间步自动调整的思路是：(1)用较大的时间步长计算断层孕震过程.当在某一计算时间步内(以时间步长107s为例)，断层两盘的相对滑移速率发生突变时，即认为在这一时间步内发生了地震事件.但由于时间步长较大，无法捕捉同震过程，因此放弃这一步的计算结果，并将时间步长减小(以时间步长106s为例)，再进行计算.此时，若断层两盘的相对滑移速率未发生突变，则认为是孕震过程，继续计算；若发生了突变，继续减小时间步长.以此类推，直至时间步长降到10-3s以下，进而计算断层的同震过程.(2)用较小的时间步长计算断层同震过程.当在较长一段时间内(以时间步长10-3s为例)，断层两盘的相对滑移速率约为零，即认为同震过程结束，此时逐渐放大时间步长，直至时间步长增至107s，进而计算断层的孕震过程.图4展示的是地震事件一和地震事件二之间，孕震过程和同震过程的时间步长变化.

图4 孕震期间和同震期间的自动增量步时间的变化Fig.4 The change of the automatic time increments during interseismic and coseismic periods

为了更好的展示断层两侧相对滑移距离随时间的变化，图5a中忽略孕震期间断层两盘位错变化，只展示2000年时间内五次地震事件的断层两盘位错随时间的变化，图中相邻曲线的时间间隔是0.5 s，五次地震事件共造成的最大累计位错约为7 m.此外，图5b—f分别展示了每次地震事件断层两侧相对滑移距离随时间的变化，每次地震的成核位置、破裂持续时间、最大位错等都是不同的，这说明这些地震事件不是简单的重复.结合图3和图5可以看出，同一条断层，当孕震时间越长，其下一次地震事件造成的最大位错就会越大，即震级越大.图6给出地震事件五中断层在破裂过程中产生的地震波在介质中传播时，质点振动速度云图在不同时刻的快照.破裂在断层上成核后，向地表方向传播，约6 s后传至地表.图5和图6都表明本文所使用的方法可以模拟断层详细的同震破裂瞬时过程.

从上述模拟结果看，本文所使用的方法可以连续模拟断层孕震-同震准周期循环过程，不仅可以实现强震复发周期，还可以模拟断层的同震破裂过程.此外，相较于其他数值计算方法，有限单元法可以很容易处理复杂几何边界、复杂物性结构、流固耦合等等高度非线性问题，因此今后利用此方法可以将更加简单、方便的研究一些复杂条件下的地震序列问题.

3 问题讨论

3.1 初始场的设定

与前人(Duan and Oglesby, 2005a, b; Lapusta and Liu, 2009; Noda and Lapusta, 2010)模拟地震序列问题直接给定初始应力场不同，本文中设定的应力场是通过施加重力和构造载荷共同作用的结果.具体的做法是：先在自由地表至地下2 km处施加重力加速度，通过地应力平衡的方法消除模型在受到重力作用后产生的大变形；然后再施加构造载荷，随着构造载荷的不断累积，直至发生地震事件.

重力作用在地学数值模拟中是个非常重要的问题，对比本文研究的地震序列问题，以往的研究者(Duan and Oglesby, 2005a,b; Lapusta and Liu, 2009; Noda and Lapusta, 2010)都不予考虑.究其原因，如图7a所示，当对整个模型施加重力后，随着深度的增加，其应力大小也在不断的增加.然而地球内部的应力场绝不仅仅是一种简单的、线性的力学过程，而是多圈层相互作用、受多种物理和化学因素控制的、非线性作用过程.通常认为由于孔隙压、构造应力的作用，达一定深度后，断层面上的有效正应力就不再随着深度的增加而增加.然而若不考虑重力作用，整个模型特别是自由地表附近没有了重力加速度的约束，显然也是不合理的.因此参照前人在地下2 km以下给定均匀初始应力场的做法(Lapusta and Liu, 2009)，本文在自由地表至地下2 km处施加重力加速度，如图7b所示，兼顾深部应力场大小和重力加速度的作用.

3.2 摩擦系数对地震准周期的影响

为了考察摩擦系数对地震准周期的影响，本文通过改变静摩擦系数和动摩擦系数计算多组模型进行对比分析，图8给出了其中4种典型情况下断层典型点位置上下盘之间的位错随时间的变化.根据图8，比较4种情况下的计算结果：

图5 每次地震事件断层上各点滑移随时间的变化Fig.5 Snapshots of the slip profiles at every node along the fault vary with time in each event

图6 地震事件五中不同时刻断层周边介质的振动速度云图Fig.6 Contour distributions of particle velocities at different times around the fault in event 5

图7 重力作用下，垂直于地表方向应力云图(a) 对整个模型施加重力; (b) 自由地表至地下2 km处施加重力.Fig.7 The vertical stresses of the models under applying gravity(a) Gravity is applied to the whole model; (b) Gravity is applied to the region between free surface and 2 km underground.

(1) 比较图8a、图8b和图8d，在动摩擦系数相同时，当断层上的静摩擦系数越大，断层的强震复发周期越长.

(2) 比较图8b、图8c和图3，在静摩擦系数相同时，当断层上的动摩擦系数越大，断层的强震复发周期越短.

(3) 比较图8c和图8d，可以看出，当断层上的静、动摩擦系数差值相同时，断层的强震复发周期也相当.

由上述分析可见：断层上的静、动摩擦系数差值直接影响着断层的强震复发周期，差值越小，复发周期越短；差值越大，复发周期越长.

此外，为了提高效率，研究本质问题，本文采用二维有限元模型进行计算，但本文所采用的方法是完全可以用来计算三维模型的.在今后分析复杂的实际地震问题时，会将此方法扩展应用到三维模型.

4 结论

本文基于有限单元动态隐式算法发展出一种新的计算方法，连续稳定的模拟了断层孕震-同震准周期循环过程.新的计算方法既可以模拟断层的强震复发周期，同时也可以模拟断层同震的详细破裂过程.此外，我们研究了摩擦系数对地震准周期的影响，发现断层上的静、动摩擦系数差值直接影响着断层的强震复发周期，差值越小，复发周期越短；差值越大，复发周期越长.

由于有限单元方法能够适应复杂的断层几何、复杂的边界条件、复杂的材料属性，因此相对于其他数值模拟方法在计算真实的地震、地质问题时更具有优势.在今后定量分析复杂的实际地震问题时，可以运用本文所发展的方法开展深入研究.因此本研究对于探索强震孕育发生的动力学机制及地震预测预报等有重要的科学意义及实际应用价值.

图8 不同摩擦系数下，断层典型点位置上下盘之间的位错随时间的变化Fig.8 Slip between hanging wall and footwall at the typical point on the fault vary with time in different models in which friction coefficients are different