基于时间域半空间 Green 函数建立的边界积分方程及其在动态破裂模拟中的应用

刘煜杭, 钱峰, 冯禧, 张海明

北京大学地球与空间科学学院地球物理系, 北京 100871

0 引言

地震断层的自发破裂传播模拟是通过震源过程了解地球介质应力状态的重要手段.对地面设施破坏力较强的大地震理所当然受到更多关注,在研究这一类地震时,往往需要重点考虑其破裂直接扩展到地表以及断层规模较大的特征.当破裂直接扩展到地表时,其传播会明显受到地表的影响,因此不能用简单的全空间模型来考虑,必须引入自由界面来进行研究(Zhang and Chen, 2006a,b; Ando and Kaneko, 2018);而较大的断层规模则往往伴随复杂的断层几何结构,这对破裂的传播也会产生很大的影响(Aochi and Madariaga, 2003; 钱峰等, 2019).边界积分方程方法(Boundary Integral Equation Method,BIEM)可以将域内的问题转化到边界上,适于解决这类具有复杂边界的问题,而且这种做法将三维问题转化为二维问题,可以大大减少计算量.同时与有限差分法和有限元法等域方法相比,边界积分方程方法中自由表面和断层几何形态等介质信息可以单独在事先计算好的积分核中体现,不需要耦合在破裂过程中,因此在需要针对同一断层模型进行大量正演计算的情况下只需要计算一次积分核,避免了重复计算Green函数造成的计算资源损耗.目前边界积分方程方法已经在需要大量动力学正演计算的破裂相图相关研究中得到应用(Xu et al., 2015; 郑玲珑等,2021).

因为边界积分方程方法在解决复杂边界问题中具有优势,所以很多地震学家将其应用到复杂断层的震源自发破裂传播问题中.Das和Aki(1977)首先利用位移表示定理建立了位移边界积分方程,但位移边界积分方程只能应用于对平面断层的研究.为了更好地体现断层复杂几何结构的影响,地震学家渐渐将目光转向了牵引力边界积分方程(Fukuyama and Madariaga, 1995, 1998).利用牵引力边界积分方程不仅能处理平面断层的问题(Yamashita and Fukuyama, 1996; Kame and Yamashita, 1997),也能解决任意复杂断层的震源自发破裂传播问题(Aochi et al. 2000a, b; Tada et al., 2000).

因为牵引力边界积分方程基于全空间 Green 函数建立,所以只适用于全无限空间介质中的破裂传播问题.目前考察自由表面影响的方法主要有三种:第一种是利用几何对称性引入镜像虚拟震源(Lapusta et al., 2000; Aochi and Fukuyama, 2002),但是这种方法只适用于垂直断层;第二种是引入满足无牵引力条件的虚拟断层面来模拟自由表面(Ando and Okuyama, 2010; Hok et al., 2011),这种方法适用于任意复杂断层,但是虚拟断层面会导致模型空间变大,大大增加了计算量;最后一种则是将全空间 Green 函数替换为半空间 Green 函数来建立边界积分方程 (Zhang and Chen, 2006a,b),这种方法理论上适用于任意形状的断层,但相关的研究使用的是频率域半空间 Green 函数,在高频区域需要额外的采样点来保证精度,造成了非必要的计算浪费.

此前研究使用频率域半空间Green 函数计算积分核 (Zhang and Chen, 2006a,b),需要大量频率域采样点来保证变换回时间域后Green函数的精度.本研究直接从半无限空间Green 函数的时间域积分形式解(Johnson, 1974)出发,利用数值积分方法计算半空间Green 函数积分核,在保证精度的同时减少了频率域采样点较多导致的额外计算量,提高了计算效率.并在此基础上改进了钱峰等(2019)的方案,将复杂断层系统下的自发破裂传播模拟方法从全无限空间介质模型拓展到半无限空间介质模型中.在验证了方法的正确性之后,本文探究了自由表面对破裂传播过程的影响范围,并初步展示了自由表面影响下弯折断层的破裂过程.

1 方法介绍

本文从半空间 Green 函数(Johnson, 1974)出发,推导出离散化积分核的积分形式表达式,在用数值积分的方法预先计算得到积分核之后,将这些积分核代入边界积分方程中,联立适当的摩擦准则即可计算出自由界面影响下的地震断层自发破裂传播过程.

1.1 半空间Green函数积分核的理论公式

在利用边界积分方程方法模拟断层的自发破裂过程时,需要先计算包含了介质信息的积分核,即应力 Green 函数.为了研究自由表面对破裂传播过程的影响,本文使用半空间 Green 函数(Johnson, 1974)计算积分核,再代入破裂过程计算(钱峰等, 2019)中,最终得到了自由表面影响下的断层自发破裂过程.

在此前的研究(Zhang and Chen, 2006a,b) 中,应力分量的表达式是在地表坐标系Ox1x2x3中给出的,其中自由表面在x3=0的平面内,x3轴指向地心(图 1).

但是由于其整个计算过程都是在频率域进行的,为了减少计算量的损失,首先反变换得到时间域内的应力表达式

其中

(2)

(3)

(4)

(5)

在无限空间介质下,Green函数的空间导数满足(Aki and Richards, 2002)

(6)

但将自由界面的影响纳入考虑范围后,式(6)在x3和ξ3方向不再成立.为表示方便,式(2)—(5)中

(m,q=1,2,3)

(7)

为更加方便地研究断层的破裂行为,选择震源坐标系O′x′1x′2x′3,使得断层面在x′3=0平面内,通过简单的坐标转换即可得到震源坐标系O′x′1x′2x′3下的应力分量.仅考虑相对普遍的断层剪切破裂的情况,自然有断层沿x′3方向的位错 Δu′3=0.

以走滑断层为例,在计算破裂过程时,需要使用断层上的正应力τ′33分量和剪应力τ′31分量,若考虑断层的弯折、分叉等复杂情况,还需要增加应力τ′11分量进行应力的合成.具体到震源坐标系O′x′1x′2x′3中,对于三个位错分量来说只有Δu′1非零.在此假设下,将各应力分量由地表坐标系Ox1x2x3中转换到震源坐标系O′x′1x′2x′3中,可以得到当前问题所关心的应力τ′11,τ′31,τ′33分量的表达式:

+2Δu1[η0(G11,21+G12,11)-(G11,31+G13,11)]}dτdS,

(8)

+η1η2(-G′11,23-G′12,13+G31,21+G32,11+G21,31+G23,11+G11,32+G13,21)]}dτdS,

(9)

(10)

在得到三个应力分量表达式后,为了通过数值方法进行地震断层破裂的自发传播模拟,需要将断层离散化.与之前研究使用的结构化网格(Zhang and Chen, 2006a, b)不同,本研究使用三角形离散化方案 (钱峰等, 2019),这种离散化方案在研究复杂断层结构时比较方便.

设三角形离散化后的应力分量可以分别表示为

(11a)

(11b)

(11c)

-η1(G′31,33+G′33,13-G11,13-G13,11-G21,23-G23,12)

+2η1(G11,13+G13,11)-η2(G11,12+G12,11)]}dS′,

(12)

(13)

(14)

特别地,当δ=90°即断层为垂直走滑断层时,(12)—(14)式退化为

(15)

(16)

(17)

至此,得到了离散化的半空间Green函数积分核的表达式.对于具体的断层离散化方案,利用数值积分方法计算得到积分核,便可以将其代入破裂过程的计算中.需要指出的是,在计算积分核的过程中可以将全空间Green函数积分核与自由表面部分的影响分离,即可以使用闭合解计算全空间Green函数积分核(Tada et al., 2000),再使用数值积分方法计算自由表面项,两者相加即为所需的半空间Green 函数积分核.

1.2 边界积分方程与摩擦准则

边界积分方程方法最早由Das和Aki(1977)引入到动力学破裂模拟中,本研究使用的边界积分方程为:

×V(ξ,τ)dτ(x∈Γ).

(18)

同时,本研究使用滑动弱化摩擦准则 (Ida, 1972),将其表示为如下形式:

(19)

式中,τ是当前时刻单元的剪应力,τu表示破裂强度,τf为残余应力,Du和Dc分别为当前时刻单元的总滑动量和临界滑动弱化位移.将断层离散化后,在给定的断层初始应力分布下,联立摩擦准则求解边界积分方程,可以得到各网格点上的位错变化历史,即断层破裂的自发传播过程.

2 正确性检验

因为边界积分方程方法可以将计算积分核和破裂模拟两个过程完全分离,所以可以采取先检验积分核,再检验破裂过程的逐层检验方法来验证本文所得结果的正确性.

2.1 积分核的检验

为检验积分核的正确性,选取张海明(2004)中使用结构化网格给出的一个半空间Green函数积分核计算实例,为与其统一,本文使用两个等腰直角三角形拼成的正方形单元上的积分核与其计算结果进行对比(图 2).

2.2 近地表平面断层的自发破裂过程

在计算出积分核之后,可以直接将其应用到破裂过程的模拟中.为考察半空间 Green 函数积分核在破裂过程中的应用效果,本文选择The SCEC/USGS Code Verification Web Server 项目中的TPV5问题作为参考算例进行比较,该项目的网站地址为:https:∥strike.scec.org/cvws/.

考虑一个30 km×15 km的与自由表面相连的垂直走滑断层,破裂成核区是位于7.5 km深度处的边长为3 km的正方形区域,区域内的初始应力为81.60 MPa,成核区两侧分别有同样大小的高初始应力区和低初始应力区,这两个区域的初始应力分别为78.00 MPa和62.00 MPa,其余区域的初始应力设为70.00 MPa(图4).断层的破裂强度为81.24 MPa,残余应力τf=63.00 MPa,P波速度和S波速度分别取cL=6000 m·s-1和cT=3464 m·s-1,介质密度ρ=2670 kg·m-3.在模拟破裂过程时,滑动弱化距离Dc设为0.40 m.如图5所示,断层使用三角形网格离散.

本文得到的破裂前锋演变图与USGS 提供的破裂前锋演变图在整体行为上一致(图 6).首先,破裂前锋从断层中心的正方形成核区开始传播,并逐渐向外扩展;当破裂前锋传播到左侧的高初始应力区,破裂速率有明显加快,表现为破裂前锋外凸,相应地,在右侧的低初始应力区,破裂速率减慢,对应到破裂前锋内凹;破裂前锋抵达自由表面时,由于直达震相与反射震相的耦合作用,自由表面处破裂传播速率会明显加快,表现为破裂前锋向外弯曲,这个效应随着破裂向两侧传播会越来越明显,尤其是在断层右侧由于低初始应力区的存在,导致直达震相速率减慢,使得自由表面震相更加明显.需要注意的是,本文得到的结果与USGS提供的结果相比,破裂的整体传播速率会相对慢一些,这可能是由破裂模拟的网格相对较大导致的.

3 分析与讨论

与全空间Green函数积分核相比,半空间Green函数积分核中增加了自由表面的影响.为探讨当前方案下自由表面对破裂过程模拟的影响程度,本节将对自由表面项的作用做出进一步讨论.

3.1 积分核的行为

理论上,自由表面的作用可以直接在积分核的行为中体现出来,具体表现为反射震相导致的积分核数值波动.为更加清晰地显示自由表面对积分核行为的影响,本文在图7显示模型下探究积分核行为随观测点深度的变化.断层面假设为一个x2=0平面内的等腰直角三角形单元ABC,三角形单元的顶点分别为A(0,2)km,B(1,2)km和C(0,3)km,观测点在x1=10 km的垂直于自由表面的直线上移动,P波速度和S波速度分别取cL=6000 m·s-1和cT=3464 m·s-1,比较观测点在不同深度时的积分核行为.

如图8所示,全空间Green函数积分核仅有直达P波和直达S波导致的积分核数值波动,而半空间Green函数积分核的自由表面项则包含了反射 PP,PS,SP,SS四种震相导致的积分核数值波动.在观测点深度逐渐增大的过程中,各个震相导致的积分核数值波动由于反射波的传播距离增大有明显的时间延迟,而且波动的幅度有明显的减小,这也印证了自由表面的影响会随着与自由表面距离的增大而逐渐削弱.

3.2 破裂过程随断层深度的变化

自由表面会导致积分核行为上出现反射震相导致的数值波动,同时随着震源或者观测点随自由表面距离的增加,自由表面的效应会明显减弱,可以预见,将半空间 Green 函数积分核应用到破裂过程模拟中时,也会出现类似的特征.为验证自由表面对破裂过程的影响,选取图9所示的简单平面断层模型进行破裂过程模拟.假设28 km×15 km的断层上沿到自由表面的距离为h,h不断变化,破裂成核区为断层左侧一个半径为1.5 km的圆形区域,成核区到断层左端的距离为8 km,到断层上沿的距离为7.5 km.不同深度的断层模型根据三角形离散化方案离散化后,首先计算出半空间Green函数积分核,再将积分核代入破裂过程模拟程序中,可以绘制出不同断层深度下的破裂过程快照.断层上的正应力和剪应力分别为120 MPa和70 MPa,滑动弱化距离设为3 m,成核区的初始剪应力和断层的破裂强度与2.2节相同.

当前离散化方案使用三角形的网格平均边长为 1 km.破裂过程快照显示,自由表面对破裂过程的影响随断层上沿到自由表面距离的增加有明显的减弱.在该离散化方案下,当断层与自由表面直接相交,可以明显看到自由表面产生的反射震相,以及自由表面处的滑动速率和破裂传播速率增加;当断层上沿距自由表面的距离为 1 km 时自由表面的影响已经比较微弱,几乎只能看到自由表面导致的微弱反射震相;而在断层上沿距自由表面的距离为 5 km 时则几乎看不到自由表面的影响(图 10).

自由表面对断层自发破裂过程传播的影响主要体现在与地表直接相交的断层上,这与预期结果吻合.在给定的离散化方案下,当断层逐渐深入地下,自由表面的影响会逐渐削弱而且削弱速率非常快,因而对于断层距离地表较近尤其是与地表直接相交的情况,自由表面对破裂过程传播的影响是必须考虑在内的,而对于断层较深的情况,自由表面的影响显著减弱.

3.3 近地表弯折断层的自发破裂过程

在无限半空间介质中,可以清晰地看到自由表面对简单平面断层上自发破裂过程的影响,包括反射震相和地表处的破裂传播速率增加.为了进一步探究半空间Green 函数积分核在复杂断层系统破裂传播过程中的应用,考虑将图9所示平面断层划分为 16 km×15 km和12 km×15 km两个断层组成的弯折断层系统,两断层之间的弯折角为30°,断层上沿与自由表面相交.动力学参数设置与 3.2 节相同(图 11).

为了更直观地显示加入自由表面影响后破裂传播过程的变化,除直接使用半空间 Green 函数积分核模拟破裂外,还将半空间 Green 函数积分核中的自由表面项去除,直接使用全空间 Green 函数积分核计算破裂过程进行对比.

破裂过程快照结果显示,应力首先在成核区积累,然后从成核区向外传播(图12).在开始阶段,使用两种积分核模拟出来的破裂过程是完全相同的.当破裂传播到自由表面时,两种介质模型下的破裂传播过程开始发生变化,使用全空间积分核模拟出的破裂在上下界面处的行为是完全一致的;而使用半空间积分核的破裂过程模拟结果表明,破裂在自由表面处有明显的加速,导致破裂先接触断层上沿即自由表面.破裂继续传播,相较全空间介质模型下的结果,在半空间介质模型下,破裂传播出现非常明显的反射震相.同时由于反射震相和直达震相的叠加,自由表面破裂前锋处的滑动速率有明显增大,破裂速率也随之加快.当破裂从主平面经过弯折部位向弯折面方向传播时,破裂会产生短暂停留,之后继续沿弯折面方向传播,在这个过程中,以上所述自由表面的影响仍然明显存在.无论在主平面还是弯折面,反射震相都会随着距离自由表面的距离增大而逐渐削弱,当传播到断层下沿时已经几乎没有显现.

图1 地表坐标系Ox1x2x3和震源坐标系O′x′1x′2x′3示意图浅灰色区域为自由表面,深灰色区域表示断层面,断层上沿到自由表面的距离为h,倾角为δ,Ox1x2x3和O′x′1x′2x′3分别为地表坐标系和震源坐标系.

图2 用于积分核正确性检验的断层模型断层为x2=0平面内与自由表面相交的由两个等腰直角三角形ABD和BCD组成的正方形,观测点R在图中星号所示位置.

图3 积分核对比结果上方为本文结果,下方为张海明(2004)结果,其中k-n表示时间步.

图4 TPV5问题的断层模型问题采用如图所示30 km×15 km的垂直走滑断层,成核区为断层中心的3 km×3 km的正方形区域,在成核区的左右两侧分别有一个同样形状大小的高初始应力区和低初始应力区.

图5 离散化后的TPV5问题断层模型采用三角形离散化方案离散的TPV5 问题的断层模型,网格颜色越深表示该区域的初始应力越大.

图6 破裂前锋随时间的演变图像上方为本文结果(网格尺度为 750 m),下方为USGS(United States Geological Survey,美国地质调查局)提供的结果(有限元法,网格尺度为 100 m).

图7 用于探究积分核行为的断层单元模型断层为x2=0平面内的等腰直角三角形单元ABC,星号所示为观测点R,观测点R在x1=10 km 射线上移动.

图8 不同观测点深度下的积分核随时间变化图像(a) 黑线为全空间积分核,红线为半空间积分核的自由表面项; (b) 放大5倍后的半空间积分核自由表面项.

图9 平面断层模型问题采用如图所示28 km×15 km的垂直走滑断层,成核区为如图所示半径为1.5 km 的圆形区域,成核区到断层左侧边缘的距离为8 km,到断层上沿的距离为7.5 km,断层上沿到自由表面的距离h不断变化.

图10 不同断层面深度下的不同时刻的滑动速率快照(a) 断层与自由表面相交; (b) 断层上沿距自由表面 1 km; (c) 断层上沿距自由表面 5 km.

图11 弯折断层模型问题采用断层由16 km×15 km的主平面和12 km×15 km的弯折面组成,成核区为如图所示半径为1.5 km的圆形区域,成核区到断层左侧边缘的距离为8 km,到断层上沿的距离为7.5 km,断层上沿与自由表面相交.

图12 不同介质模型下的不同时刻的滑动速率快照(a) 半空间模型; (b) 全空间模型.

在破裂传播过程快照中,可以清晰地看到破裂的整体传播路径以及变化趋势.但为了更深入了解传播过程中破裂前锋处破裂速率及滑动速率的变化情况,本文分别在断层的自由表面处及深度为 7.5 km 处取了两条网格单元线,绘制了断层破裂前锋传播图.

如图 13 所示,破裂前锋在约 3 s 时抵达自由表面并向两侧延伸.在主平面上,破裂前锋传播图的斜率绝对值很小,说明破裂传播速率很快,这与自由表面加速破裂传播的结论相吻合,同时滑动速率也比较大.在弯折面上,随着传播时间逐渐增加,反射震相由于传播速率较快,逐渐与直达震相分离,形成两条破裂前锋,其中反射震相导致的破裂前锋速率略快于直达震相导致的破裂前锋.此外,两条破裂前锋上的滑动速率也会由于解耦而减小.而深度为 7.5 km 的网格单元线刚好穿过成核区,因此可以直接看到初始时刻破裂在成核区产生并向两侧扩展的过程.同时随着破裂的传播,单元的滑动速率逐渐加快,直达震相经过弯折部位的短暂停顿,在弯折面继续传播.最早的反射震相在约 6 s 的时刻抵达成核区即x1=8 km处,但是由于离自由表面较远,反射震相明显较弱.

4 结论

本文在边界积分方程方法的框架下,从半无限空间Green函数的积分形式解出发,利用数值积分的方法推导出了三角形离散化的半空间Green函数积分核表达式,使其替代全空间Green函数积分核,模拟出自由表面影响下的断层自发破裂传播过程.并通过与前人结果的对比,验证了本文半空间积分核结果的正确性.为了探究半空间积分核中包含的自由表面效应对积分核以及破裂过程的影响,本文分别给出了积分核与破裂过程随观测点/断层深度变化的图像,同时拓展了半空间积分核在弯折断层破裂过程中的应用.

相比全空间积分核,半空间积分核中包含了自由表面的效应,这种效应会导致破裂传播过程中地表反射震相的产生以及自由表面处的破裂传播加速.但是在给定的断层离散化方案下,这种效应仅在断层直接与自由表面相交或者二者相距很近的情况下才能明显体现,在震源较深时自由表面的影响可以忽略.弯折断层的弯折部位往往会限制破裂往弯折面的传播,但是由于自由表面会导致反射震相的产生,反射震相可以与直达震相耦合,加快耦合处的断层滑动速率,从而促进破裂的传播,因此在二者的耦合作用下,破裂在弯折面继续传播的可能性也会增大.而在计算成本上,本文得到的积分核为时间域内的结果,相比频率域内的计算,可以大大减少高频区域采样点增加造成的浪费,很大程度上提高了计算效率.

使用半空间Green函数积分核模拟的破裂过程可以很好地体现自由表面的影响,非常有助于对浅源地震的研究,边界积分方程方法中采用的三角形离散化方法则有助于讨论各种复杂断层的情况,而浅源与断层模型复杂也恰好是破坏性大地震的共有特征.本文将自由表面的效应与复杂断层的效应成功结合起来,建立了一套较为完善的破裂过程正演体系,同时,随着Feng和Zhang(2021)给出Lamb问题的广义闭合解,计算积分核的效率也预期会在当前使用积分解的基础上有很大提升.在后续的研究中,将会进一步进行动力学反演的尝试,这有助于深层次了解断层破裂产生并传播的力学机制,对于了解破坏性大地震具有重要意义.

致谢本研究工作得到北京大学高性能计算校级公共平台支持,这使得计算积分核的效率大大提高;同时感谢北京超级云计算中心提供的计算资源支持.此外,对责任编委和匿名审稿人提出的宝贵意见表示衷心感谢.