• 首页关于本刊投稿须知期刊订阅编委会期刊合作查询检索English
层状介质对破裂过程的影响
层状介质对破裂过程的影响
齐梦雪*, 周红*
(中国地震局地球物理研究所,北京 100081)
 [收稿日期]: 2016-11-10
摘要

震源动力学中破裂产生的地震动在层状介质中的传播模拟,是地震学以及地震工程学研究的前沿课题之一。本文通过建立精确的三维模型,选取具备灵活网格、高精度高效率计算性能的谱元法,利用有效抑制伪震荡的时间域离散方法——加权速度Newmark方法以及多次透射人工边界条件,进行了SCEC/USGS基准项目中TPV5模型的地震破裂过程模拟,得到基于层状介质模型和均匀介质模型(后者采用相同破裂模型)的埋深2km的震源参数结果。将二者进行对比,并具体分析破裂面位错、地震矩、破裂传播时间、上升时间和地表位移,发现层状介质对破裂过程的传播影响较为明显:①层状介质的存在整体增加了破裂面上的位错,在层状介质模型下计算得到的地震矩约是均匀介质模型结果的1.3倍,因此认为层状介质增强了地震破裂过程中的能量释放;②层状介质的存在使得破裂传播至地表的速度减慢,并缩短了地表各点的上升时间,增强了地表的地震动响应;③层状介质对于地表位移有着明显的增加作用,同时协同破裂面上的初始应力异常区域对位移峰值中心的改变有显著影响。④介质分异面附近地震动强烈。对结果进行整理后发现,在具有地下层状介质的地区要充分考虑层状介质产生的场地效应,否则可能会低估该地区的地震危险性。



引言

动态破裂问题,即地震的孕育、发生、地震波的传播等问题,一直是地震学研究的重要课题之一。在破裂模型的研究中,许多学者运用二维模型展开了大量的模拟工作(Andrews等,1997Duan等,2006)。但基于地震破裂动力学的概念,断层的几何特性、介质的不均匀性等因素在破裂传播过程中起重要作用(Harris,2004),有必要合理地考虑地质构造、局部场地条件、介质各向异性以及地震波的传播特性等多方面因素。为了保证地震破裂模拟的真实性和可靠性,需要进一步考虑断层几何特性、传播介质的不均匀性等因素。目前,基于非均匀介质模型的破裂过程模拟还很不完善。例如,Andrews等(1997)对位于加利福尼亚地区的部分圣安德烈亚斯断层的二维模拟结果显示,该区域的地震波应该率先按一个方向传播,但调查显示该区域存在横向介质不均匀性,同时通过一系列真实地震的观测数据发现,地震波的传播在各个方向上均有呈现(Fletcher等,1998),因此需要更多的地质结构模型和大震观测研究来确定介质对于破裂动力学问题的影响。我国地域辽阔,地壳介质构造复杂,具有浅层多沉积、地下异常体构造复杂和地表起伏剧烈等特征,通过研究我们猜测这种复杂介质会对地震动的场地效应产生一定的影响,而本文的实验也验证了介质的不均匀性对地表的地震动有影响。因此,针对非均匀介质模型的地震破裂传播过程的研究不仅在地震学领域有重要影响,对于工程学领域更有着参考意义。

在地震应急响应的相关研究中,考察地震的动态破裂过程十分重要(刘成利等,2013),通过其结果不仅可以了解强地面震动(PGV和PGA)的情况、评估地震的影响范围,还能够调查地表破裂特征(白玉柱等,2009陈立春等,2010)。通过动态破裂传播过程中包含的破裂细节,如断层面的位错、破裂速度、断层破裂引起的地表运动以及地表破裂的尺度等,探究近断层地震动的基本特征和分布规律,从而预测未来发震断层附近的强地震动分布区域和地面运动时程,为结构抗震、城市规划及震害预防奠定基础(刘启方,2005杜晓晨等,2008)。近年来,动态破裂过程的研究结果常常用于考察近断层地震灾害的致灾机理(张勇等,2008刘成利等,2013),为抗震救灾及设防工作提供了重要参考。

由于计算手段的制约,20世纪的断层破裂模拟工作主要针对简单介质中简单断层的破裂开展研究(Burridge等,1974Richard,1976)。随着计算技术和设备的发展,复杂的断层破裂可以通过有限差分法、边界积分法等数值方法求解(Madariaga等,2000Zhang等,2004, 2005, 2006Oglesby等,2000a, 2000b)。在数值模拟方法丰富的情况下,学者们越来越关注破裂所在的介质模型。最初介质模型为简单的全空间均匀介质(Burridge等,1974Richard,1976),仅从物理学破裂传播方面建立模型模拟地震动,随后发展到在均匀半空间介质条件下开展研究(Oglesby等,2000a, 2000b, 2001)。同时,水平层状介质中的断层破裂研究也日益丰富(Andrews等,1997Harris等,1991)。Andrews等(1997)曾指出,在非均匀介质中展布的断层,其破裂位错在较软介质中更大,由此也可以看出介质性质对破裂传播有明显影响。目前,随着计算技术的发展,针对在更符合实际地壳情况的非均匀介质中的断层破裂行为开展的研究开始受到越来越多的关注,Zhang等(2014)利用Zhang等(2006)提出的曲面网格有限差分法模拟了地震波在任意不规则形态介质中的传播情况,并灵活运用曲线网格构建各种复杂断层模型,使得各种断层模型的研究在非均匀介质中可以系统地展开。但目前基于非均匀介质模型所做的工作还远远不够,例如对非均匀介质对断层破裂在地表产生的影响等问题的研究还较少,因此未来需要通过三维模型和大震观测开展更多的研究,以确定介质对于破裂动力学问题的重要性(Harris,2004)。

谱元法(SEM)由Patera(1984)引入流体问题研究中,用于求解Navier-Stokes方程的数值解。随后Seriani等(1992)首次将其应用于地震波传播模拟运算。谱元法的基本思想是将求解区域划分为由若干不均匀节点组成的有限个子单元,随后采用正交多项式展开式,形成谱单元,最后采用Galerkin方法求得全局近似解,即得到求解区域内动力学波传播过程的模拟。Ronquist等(1988)最早将Lagrange插值基函数引入谱元法中,与Gauss-Lobatto-Legendre(GLL)积分相结合,形成了Legendre谱元法,其基本思想是对以Legendre正交多项式为基的Lagrange插值进行Gauss-Lobatto-Legendre(GLL)积分,生成对角质量矩阵。其数值计算为简单的时间域内的显式时间积分,因此具有计算过程简单、效率高等特点。Legendre谱元法也适用于非均匀、各向异性以及复杂地质构造区域的地震动场模拟(Komatitsch等,2000, 2002, 2004王童奎等,2007李孝波等,2014)。基于如上所述,Legendre谱元法的计算结果精度较高,且具有网格介质的灵活性,因此本文考虑了断层信息以及分层介质信息,以简单层状介质为研究对象,探讨了不均匀介质对地震破裂在传播过程中的具体影响。

对于动态破裂问题,谱元法的数值模拟存在伪震荡的弊端,Zhou等(2015)提出了加权速度Newmark方法(The Weighted Velocity Newmark Scheme),该法抑制了谱元法的伪震荡效应,同时不影响破裂传播的模拟。本文利用引入了加权速度Newmark方法的谱元法程序代码,围绕不均匀介质模型与地震动特性展开研究工作。将破裂面位错、破裂传播时间、上升时间等断层震源参数作为研究因子来描述层状介质对破裂过程的影响,量化断层面信息,以便于归纳与分析。同时,选取简单层状介质为研究模型,降低了建立介质模型的难度,也使得对于各种震源参数的分析更易结合介质差异,更加有针对性。通过对上述断层震源参数的研究发现层状介质对破裂过程的传播影响较为明显。

1 理论及方法
1.1 谱元法理论

若将地震波近似为弹性波,则地下介质地震波的传播可以表示成Ω空间内带有初始条件和边界条件的线弹性波动方程:

$\left\{ \begin{matrix} \rho \bf{\ddot{u}}=\nabla \ \centerdot \ \bf{ \pmb{\mathsf{ σ}} }+\bf{f} \\ \bf{ \pmb{\mathsf{ σ}} }=\bf{C}:\mathit{\boldsymbol{ }}\!\!\varepsilon\!\!\rm{ } \\ \bf{ \pmb{\mathsf{ ε}} }=\frac{\bf{1}}{\bf{2}}[\nabla \bf{u}+{{(\nabla \bf{u})}^{\rm{T}}}] \\ \end{matrix} \right.$ (1)

式(1)中第一个等式称为波动方程的强化形式,式中ρ代表密度,u为位移矢量,σ为应力张量,f为震源项,C表示刚度矩阵,ε是应变张量。

为求得波动方程的数值解,需要在公式(1)的两端点乘一个任意测试向量w,在整个计算域Ω内和边界进行分部积分得到波动方程的弱化形式:

$\underset{\Omega }{\mathop{\int }}\,\rho \mathbf{w}\cdot \ \mathbf{\ddot{u}}\text{d }\!\!\Omega\!\!\text{ }+\underset{\Omega }{\mathop{\int }}\,\nabla \mathbf{w}:\mathbf{C}:\nabla \mathbf{u}\text{d }\!\!\Omega\!\!\text{ }=\underset{\Omega }{\mathop{\int }}\,\mathbf{w}\ \cdot \ \mathbf{f}\text{d }\!\!\Omega\!\!\text{ }+\underset{\Gamma }{\mathop{\int }}\,\mathbf{w}\ \cdot \ \text{Td}\Gamma $ (2)

式中,$\Gamma $代表边界,T是包括断层面边界或者人工边界以及自有界面的引力矢量。波动方程的弱化形式(2)与波动方程强化形式(1)是等价的。Komatitsch等(1999)在三维地震波谱元法的研究中提到,采用弱化形式不仅能很好地模拟波动方程在线弹性介质中的传播,而且对非弹性介质、表面地形和面波传播的模拟也更为精确。

在谱元法的应用中,一般先将二维问题的计算域划分成离散四边形单元网格,而三维问题的计算域则被划分为离散六边形网格,单元网格用${{\text{ }\!\!\Omega\!\!\text{ }}_{\text{e}}}$表示。再将所有单元网格投射到局部坐标的参考区域[-1,1]中,选取GLL积分节点,并根据GLL插值积分计算法则进行全局聚合计算,最后按照一定顺序进行叠加,即可得到公式(2)的矩阵形式:

$[\mathbf{M}]\left\{ {\mathbf{\ddot{U}}} \right\}+[\mathbf{K}]\left\{ \mathbf{U} \right\}=\left\{ \mathbf{F} \right\}$ (3)

式中,M为全局质量矩阵,U为全局位移向量,K为全局刚度矩阵,F为震源项。

该方法的位移场采用GLL插值点的高阶单元。在本文所使用的程序中,将地表设置为自由表面,计算区域的其他边界选择廖式透射边界条件。

1.2 断层面位移公式

破裂和波的传播过程可以用边界条件来限制。我们先将一个断层简化成两个分开的破裂面,换言之,该断层上的同一计算节点被分配给这个断层的两个破裂面。每一个破裂节点的位移表示成${{u}^{+}}$${{u}^{-}}$。因此,滑移与滑移速率可以被写作:

$\left\{ \begin{matrix} D={{u}^{+}}({{X}^{+}},t)-{{u}^{-}}({{X}^{-}},t) \\ \dot{D}={{{\dot{u}}}^{+}}({{X}^{+}},t)-{{{\dot{u}}}^{-}}({{X}^{-}},t) \\ \end{matrix} \right.$ (4)

这里的${{X}^{+}}$${{X}^{-}}$表示破裂面两端任意节点的位置。

然后利用公式(5)所示的边界关系计算断层两侧的位移,为了方便描述,用各节点在走向、倾向和法向的应力和位移来代替公式(2)中的xyz向:

$\left\{ \begin{matrix} a_{\gamma ijk}^{+}(\ddot{u})_{\gamma }^{+}=S_{\gamma ijk}^{+}+a_{\gamma ijk}^{\prime +}[{{T}_{\gamma ijk}}+{{T}_{\gamma 0ijk}}] \\ a_{\gamma ijk}^{-}(\ddot{u})_{\gamma }^{-}=S_{\gamma ijk}^{-}+a_{\gamma ijk}^{-}[{{T}_{\gamma ijk}}+{{T}_{\gamma 0ijk}}] \\ \end{matrix} \right.$ (5)

这里的脚标$\gamma $代表断层的三分量(走向、倾向和法向),ijk代表节点号,${{T}_{\gamma ijk}}$${{T}_{\gamma 0ijk}}$分别是总牵引力和节点的预加应力,${{a}_{\gamma ijk}}$为有关位移的体积系数,$a_{\gamma ijk}^{\prime }$为有关节点应力的表面集成系数。${{S}_{\gamma ijk}}$为应力、位移以及震源项的体积集成。为了获得正应力,将(5)式同除${{a}_{\text{n}ijk}}$(下标n代表法向),得到下式:

$\left\{ \begin{matrix} \left( {\ddot{u}} \right)_{\text{n}ijk}^{+}=\frac{S_{\text{n}ijk}^{+}}{a_{nijk}^{+}}+\frac{a_{\text{n}ijk}^{\prime +}}{a_{nijk}^{+}}[{{T}_{\text{n}ijk}}+{{T}_{\text{n}0ijk}}] \\ \left( {\ddot{u}} \right)_{\text{n}ijk}^{-}=\frac{S_{\text{n}ijk}^{-}}{a_{\text{n}ijk}^{-}}+\frac{a_{\text{n}ijk}^{-}}{a_{\text{n}ijk}^{-}}[{{T}_{\text{n}ijk}}+{{T}_{\text{n}0ijk}}] \\ \end{matrix} \right.$ (6)

由于断层面上法向位移连续而得到法向加速度也为连续,因此上式在法向上的差值可得法向应力:

${{T}_{\text{n}}}=\frac{S_{\text{n}}^{-}/a_{\text{n}}^{-}-S_{\text{n}}^{+}/a_{\text{n}}^{+}\ \ }{a_{\text{n}}^{\prime +}/a_{\text{n}}^{+}-a_{\text{n}}^{\prime -}/a_{\text{n}}^{-}\ \ }+{{T}_{\text{n}0}}$ (7)

Tn0表示节点的法向预加应力,式中省略了节点号。如果Tn>0,则为了避免断层沿法向破裂,使Tn=0。为了得到剪切应力Tv,我们运用滑移弱化摩擦公式求得抗剪强度${{\tau }_{\text{c}}}$

${{\tau }_{\text{c}}}=-{{T}_{\text{n}}}{{\mu }_{\text{f}}}$ (8)

其中:

${{\mu }_{\text{f}}}=\left\{ \begin{array}{*{35}{l}} {{\mu }_{\text{s}}}-\left( {{\mu }_{\text{s}}}-{{\mu }_{\text{d}}} \right)\frac{D}{{{d}_{0}}},D\le {{d}_{0}} \\ {{\mu }_{\text{d}}},D>{{d}_{0}} \\ \end{array} \right.$ (9)

上述τc是抗剪强度,${{\mu }_{\text{s}}}$${{\mu }_{\text{d}}}$分别代表静摩擦和动摩擦系数,${{d}_{0}}$表示临界滑移弱化距离。在求解动态破裂信息之前,需要校正剪应力Tv,这里的下标v表示走向(s)或倾向(d)。可以利用Day等(2005)的方法引入${{\tilde{T}}_{\text{v}}}$来校正Tv${{\tilde{T}}_{\text{v}}}$可以利用切向速度的连续性从(5)式中计算得到:

$\dot{u}_{\text{v}}^{+}(n+1)-\dot{u}_{\text{v}}^{-}(n+1)=0$ (10)

根据(6)式的步骤并省去节点号,就得到${{\tilde{T}}_{\text{v}}}$

${{{\tilde{T}}}_{\text{v}}}(n+1)=\frac{\left[ \left( S_{\text{v}}^{-}/a_{\text{v}}^{-}-S_{\text{v}}^{+}/a_{\text{v}}^{+}\ \ \right)-\dot{D}(n)/\Delta t\ \right]}{a_{\text{v}}^{\prime +}/a_{\text{v}}^{+}-a_{\text{v}}^{\prime -}/a_{\text{v}}^{-}\ \ }+{{T}_{\text{v}0}}(n+1)$ (11)

而第(n+1)个步长的Tv由下式获得:

$\begin{align} & {{T}_{\text{v}}}(n+1)= \\ & \left\{ \begin{matrix} {{{\tilde{T}}}_{\text{v}}}(n+1) & \text{v}=\text{s}、\text{d},\sqrt{\tilde{T}_{\text{s}}^{2}(n+1)+\tilde{T}_{\text{d}}^{2}(n+1)}\le {{\tau }_{\text{c}}} \\ {{\tau }_{\text{c}}}\frac{{{{\tilde{T}}}_{\text{v}}}(n+1)}{\tilde{T}_{\text{s}}^{2}(n+1)+\tilde{T}_{\text{d}}^{2}(n+1)} & \text{v}=\text{s}、\text{d},\sqrt{\tilde{T}_{\text{s}}^{2}(n+1)+\tilde{T}_{\text{d}}^{2}(n+1)}>{{\tau }_{\text{c}}} \\ \end{matrix} \right. \\ \end{align}$ (12)

将计算得到的正应力Tn和剪应力Tv代入(5)式中,即可得到每一个节点第(n+1)时间步长的加速度、速度和位移。通过坐标变换,求得相应O-XYZ坐标系中的动态破裂信息。

1.3 加权速度Newmark方法(The Weighted Velocity Newmark Scheme)

为了抑制谱元法在模拟动态破裂问题中出现的伪震荡,本文对于时间域的离散方式采用了Zhou等(2015)提出的加权速度Newmark方法。Newmark方法是在预测方法(The Prediction Scheme)的基础上改进而成。在预测方法中,tn+1时位移和速度的预测值可写作:

$\left\{ \begin{array}{*{35}{l}} {{{\tilde{u}}}_{n+1}}={{u}_{n}}+\frac{1}{2}\Delta t\cdot \ {{{\dot{u}}}_{n}} \\ {{{\dot{\tilde{u}}}}_{n+1}}={{{\dot{u}}}_{n}} \\ \end{array} \right.$ (13)

式中,${{u}_{n}}$${{\dot{u}}_{n}}$分别表示tn时的位移和速度。${{\tilde{u}}_{n+1}}$表示tn+1时位移的预测值。将预测位移代入(3)式,则可以计算得到tn+1时的加速度预测值${{\ddot{\tilde{u}}}_{n+1}}$。最终,tn+1时的位移un+1和速度${{\dot{u}}_{n+1}}$可由下式得到:

$\left\{ \begin{array}{*{35}{l}} {{u}_{n+1}}={{{\tilde{u}}}_{n}}+\frac{1}{2}\Delta t\ \cdot \ {{{\dot{\tilde{u}}}}_{n+1}}+\frac{1}{2}\Delta {{t}^{2}}\ \cdot \ {{{\ddot{\tilde{u}}}}_{n+1}}={{{\tilde{u}}}_{n+1}}+\frac{1}{2}\Delta t\ \cdot \ {{{\dot{u}}}_{n}}+\frac{1}{2}\Delta {{t}^{2}}\ \cdot \ {{{\ddot{\tilde{u}}}}_{n+1}} \\ {{{\dot{u}}}_{n+1}}={{{\dot{\tilde{u}}}}_{n+1}}+\Delta t\ \cdot \ {{{\ddot{\tilde{u}}}}_{n+1}} \\ \end{array} \right.$ (14)

而加权速度Newmark方法以$[\theta {{\dot{u}}_{n}}+(1-\theta ){{\dot{u}}_{n+1}}]$代替公式(8)中的${{\dot{u}}_{n}}$,其中不同的$\theta $值意味着第n个步长和第n+1个步长的速度在求解第n+1步长的位移时对其数值大小有着不同程度的影响。运用该方法能够有效规避破裂模拟中出现的虚假震荡问题,使得模拟数据真实可信。

1.4 介质选取的依据

根据波动方程可以推导出P波和S波的速度:

$\left\{ \begin{matrix} {{v}_{\text{p}}}=\sqrt{\frac{2\mu +\lambda }{\rho }} \\ {{v}_{\text{s}}}=\sqrt{\frac{\mu }{\rho }} \\ \end{matrix} \right.$ (15)

在不同介质中,波速显然是不同的,因此在讨论层状介质时,其介质参数需要通过地层速度结构来确定。本文的速度结构参考了时华星等(2004)对第三系岩石地震波传播速度的分析,取均匀介质中vP=6000m/s,vS=3464m/s;取2km深处层状介质的vP=2720m/s,vS=1600m/s。

1.5 断层模型

本文选用矩形空间的动态破裂断层模型(SCEC/USGS开展的自发破裂代码验证项目中的TPV5模型),该模型是一个垂直右旋走滑断层,破裂到达地表,并且破裂是在一个长30km、深15km的矩形区域M内传播。初始破裂发生在矩形区域中心的一个3km×3km的正方形区域内。破裂中心与断层右端的中点区域(P2)初始剪应力较低,而破裂中心与断裂左端的中点区域(P1)初始剪应力较高。这两个区域均为3km×3km大小,且沿走滑方向距破裂中心7.5km,深7.5km。模型如图 1所示。


图 1 断层三维模型 Fig. 1 The three-dimension fault model
2 数据分析
2.1 层状介质对于断层面上破裂传播的影响
2.1.1 断层面位错峰值

断层破裂面上的位错值决定了滑动时间函数的高度,其直接影响地震动峰值,位错的值越大,地震动峰值越大。因此对断层面位错的考察可以定量地看出层状介质的存在对地震破裂传播的影响。

图 2(a)为均匀介质模型的位错等值线分布,其具有如下特点:① 位错峰值中心位于破裂中心,位错值为0.55cm;② 位错峰值次级中心为破裂中心在地表的垂直投影,位错值为0.52cm;③ 位错等值线分布以两个中心向外辐射状递减,向地下衰减速度较快,但向地表传播时衰减速度较慢;④ 由于断层模型中两个异常初始剪应力区域的影响,在破裂中心的两侧存在明显的等值线异常区域,即左侧初始剪应力高的区域的位错值比同一水平位置的其他值明显增大,而右侧剪应力低的区域位错值则比同一水平位置的其他区域的位错值小,因此造成了等值线的外凸和内凹,即初始剪应力高的地方能量释放剧烈。


(a)与埋深2km的层状介质模型(b)中断层面上位错等值线分布对比图
图 2 均匀介质模型 Fig. 2 The contour map of dislocation on the fault plane in the uniform medium (a) and in the layered medium (b)

图 2(b)为层状介质模型的位错等值线分布,从图中可以看出,各等值线均在深度2km处闭合,呈现出一个明显的分界面。将异常带的密集等值线放大可看到(图 3),2km深的层状介质分异面将断层面分为两个等值线区域:① 2km以下,即存在破裂中心的区域。其位错峰值中心位于深2.1—2.4km的破裂中心垂直投影处,中心值为0.64cm,次级中心在破裂中心,中心值为0.62cm(类似于均匀介质模型的表现);② 2km处至地表区域。该区域的位错普遍比下层区域大1倍,峰值中心位于深1.9km处的两个剪应力异常投影区,中心值为1.2cm,次级中心出现在地表,但地表位错普遍偏大,均值在1.04cm左右。值得注意的是,两区域的峰值中心均靠近2km处的层状介质分异面。


图 3 深度为2km的异常带位错放大图 Fig. 3 The magnification effect of anomaly zone at the depth of 2km

通过对图 2(a)(b)两图的比较,可以发现:

(1)层状介质的存在使得断层面上各点的位错值均有增加,因此层状介质可以看做是能量释放的“加强带”,致使破裂通过层状介质向地表传播时地震动更为剧烈。

(2)断层面上的位错峰值中心有明显变化,初步观测是层状介质与初始剪应力的共同作用。

(3)地表位错变化明显。

(4)介质分异面附近地震动响应强烈。

2.1.2 地震矩及震级

考虑到文中模型为弹性体系,品质因子无穷大,因此在本模型基础上引入粘弹性介质进行对比试验,以此验证本文模型计算得到的地震矩的可靠性。

地层的吸收衰减作用往往发生在地震波在地下介质传播时,通常将这种地下介质称为粘弹性介质。本文利用粘弹性介质基本模型——Kelvin-Voigt模型中弹性系数和粘滞系数的关系(孙成禹等,2010)进行了实验,同时参考中国东部地区品质因子范围(李光品等,2000),当vS为3464m/s和1600m/s时介质品质因子分别取500和200。在考虑品质因子的条件下计算得到均匀介质模型合成地震矩MM约为1.77×1024dyn/cm,反演出的震级约为MW5.43;而层状介质模型中合成地震矩MM为2.37×1024dyn/cm,震级为MW5.52。而在品质因子无穷大的模型中计算得到的均匀介质模型合成地震矩MM约为1.78×1024dyn/cm,震级约为MW5.43;层状介质模型合成地震矩MM为2.34×1024dyn/cm,震级为MW 5.51,两者相差不大。这是因为本文模型中的波场为低频,而杨祥森(2009)指出:粘弹性介质对地震波具有高频吸收和能量衰减作用,而对低频波的吸收较小。因此,我们通过设置品质因子而验证了本文线弹性模型在计算地震矩方面的合理性。

同时,通过对比可以发现,不管是线弹性模型还是粘弹性模型,层状介质的地震矩均比均匀介质地震矩大,前者是后者的约1.3倍,可以看出层状介质对破裂所释放的能量有放大作用。

2.1.3 破裂传播时间

图 4(a)所示的均匀介质模型破裂传播时间等值线图中可以看出,破裂传播时间变化较为均匀,以破裂中心为原点,传播时间向外呈类椭圆状辐射递增。等值线的异常外凸和内凹则是由于破裂中心两侧初始应力高低不同而产生。将上述结果与图 4(b)所示的层状介质模型破裂传播时间等值线图对比,发现大部分研究区域的传播时间没有明显差异,只是由于层状介质分异面的存在,导致2km深度处等值线收敛闭合加剧,因此与均匀介质模型相比,层状介质模型的地表处破裂传播时间较长。在均匀介质模型中破裂传播至地表的时间为3.3—5.7s,而在层状介质模型中破裂传播至地表的时间延长为4—7s。由此可知,层状介质对于破裂传播到地表的时间有延迟作用,导致破裂传播速度较在均匀介质中慢。


(a)与埋深为2km的层状介质模型(b)中破裂传播时间等值线对比图
图 4 均匀介质模型 Fig. 4 The contour map of wave propagation time in the uniform medium (a) and in the layered medium (b)
2.1.4 上升时间

上升时间指断层破裂面上某一质点发生的错动达到一个稳定值所经历的时间,它直接影响地震动的峰值和频率成分。上升时间越短,地震动的峰值越大,高频成分越多。

图 5中两模型进行对比,层状介质模型中的上升时间存在以下特点:


(a)与埋深为2km厚的层状介质模型(b)的上升时间等值线对比图
图 5 均匀介质模型 Fig. 5 The contour map of rising time in the uniform medium (a) and in the layered medium (b)

(1)在断层面上深度大于2km的区域,各点的上升时间略有增加,其中心值比均匀模型延长0.9s,中心位置也比均匀模型更靠近破裂中心,同时破裂中心附近的等值线也更为规则。

(2)等值线在深2km的层状介质分异处出现弯折,地表上升时间为0—4.5s不等,地表中心区域上升时间在1—4s内,而不同于均匀介质下地表上升时间为0—6.3s、地表中心区域多为5s左右的情况。

综合对位错、地震矩、破裂传播时间以及上升时间这几个震源参数的分析,层状介质作为能量释放的“加强带”,除了对断层面造成影响之外,其对地表的地震动响应具有主要且明显的增强效果,加剧了地表破裂,同时地表破裂范围也较均匀介质情况下更为广泛。

2.2 层状介质对地表位移的影响

从前文所述的研究结果可以看出层状介质对地表影响显著,由于篇幅的限制,本文只展示地表x向位移分布对比图。在针对地表研究时我们放大了观察区域,即对地表的观察范围沿xy方向展布48×44km2。这与谱元法并行计算设定的cpu大小有关。在并行计算过程中,每个cpu在x方向计算尺度为12km,在y方向计算尺度为11km,一共沿xy方向上有4×4个cpu,因此观察区域可扩大到48×44km2

图 6所示,图(a)中最大峰值位于破裂中心在地表的投影区域,值为0.255cm;图(b)中最大峰值位于水平位置15km以及35km的破裂面处,值为0.525cm。层状介质模型水平x向位移场的峰值最大区域由水平x轴方向的中心区域扩散至水平x方向15km以及35km处两个区域,且最大峰值为均匀模型的2倍,同时峰值区域也相对突出集中,与均匀模型相比能量释放区域的划分更为明显和尖锐。


(a)均匀介质模型,(b)层状介质模型
图 6 地表水平向位移峰值分布对比 Fig. 6 The contour map of surface peak displacement in the uniform medium (a) and in the layered medium (b)

通过对比地表水平位移场和垂直位移场的分布等值线,可以发现结论都是相似的,即均匀或者不均匀模型中的位移等值线都以地表y方向中心y=22km处的破裂断层面为对称轴,以破裂中心投射到地表的点为原点,大致呈蝴蝶状辐射对称分布,在断裂的两端部(x方向9km及39km处)变化较强烈,但可明显看出层状介质模型整体位移值均比均匀模型大,即层状介质对地表位移具有放大作用。

水平位移的峰值分布的变化与破裂面位错峰值图中的地表形态相吻合,说明之前推测的水平位移场的变化不仅由层状介质决定,还伴随初始剪应力的影响是有依据的。

3 结论

通过对两个模型下的断层面位错、上升时间、破裂传播时间、地震矩以及地表位错的分析,可以初步得出以下结论:

(1)层状介质的不均匀性对于破裂过程的能量释放有放大作用。图 2中层状介质模型的位错值较均匀模型整体偏大(最大值相差0.63cm),同时层状介质模型中的地震矩也是均匀模型的1.3倍。因此,可以说层状介质是断层破裂能量释放的“加强带”,可使破裂位错以及地震矩出现明显增加。

(2)将异常带放大后(图 2)可发现介质分异面附近地震动强烈,同时可以看出在介质分异面上有两个位错峰值中心,地表也有3个较大的位错峰值中心。通过观察其位置,推测这些峰值中心除了跟层状介质有关,还与初始应力异常以及破裂中心有关。

(3)层状介质的存在使得破裂传播至地表的速度减慢,并缩短了地表区域的上升时间,增强了地表地震动响应,而对于2km深的层状介质以下的断层面区域来说,层状介质延长了其上升时间,减弱地震动响应。

(4)层状介质对于地表位移值有着明显的增加作用,同时协同破裂面上的初始应力异常区域,对位移峰值中心的改变有显著影响。

综上所述,可将层状介质看成一个分界线和作用带。对于介质分异面以上区域,该作用带在破裂能量的释放以及地震动响应方面有加强作用;对于介质分异面以下区域,作用带延长了上升时间,使得该范围地震动响应减弱。因此,在具有地下层状介质地区,需要充分考虑层状介质产生的场地效应,否则可能会低估该地区的地震危险性。

致谢: 本文结果的完成得到了国家超级计算天津中心的大力帮助,在此表示衷心感谢。
参考文献
白玉柱, 徐杰, 周本刚, 2009. 2008年汶川8.0级地震地表垂直和水平位移场模拟计算分析[J]. 震灾防御技术, 4(3): 256-265. DOI:10.11899/zzfy20090302
陈立春, 王虎, 冉勇康, 等, 2010. 玉树MS 7.1级地震地表破裂与历史大地震[J]. 科学通报, 55(13): 1200-1205.
杜晨晓, 谢富仁, 史保平, 2008. 隐伏和出露地表断层近断层地表运动特征的研究进展[J]. 震灾防御技术, 3(2): 172-181. DOI:10.11899/zzfy20080208
李光品, 徐果明, 高尔根, 等, 2000. 中国东部地壳、上地幔横波品质因子的三维层析成像[J]. 地震学报, 22(1): 73-81.
李孝波, 薄景山, 齐文浩, 等, 2014. 地震动模拟中的谱元法[J]. 地球物理学进展, 29(5): 2029-2039. DOI:10.6038/pg20140506
刘成利, 郑勇, 葛粲, 等, 2013. 2013年芦山7.0级地震的动态破裂过程[J]. 中国科学:地球科学, 43(6): 1020-1026.
刘启方, 2005. 基于运动学和动力学震源模型的近断层地震动研究. 哈尔滨: 中国地震局工程力学研究所.
时华星, 曹建军, 伍向阳, 等, 2004. 济阳坳陷下第三系岩石地震波传播速度分析[J]. 油气地质与采收率, 11(4): 28-30.
孙成禹, 肖云飞, 印兴耀, 等, 2010. 黏弹介质波动方程有限差分解的稳定性研究[J]. 地震学报, 32(2): 147-156.
王童奎, 李瑞华, 李小凡, 等, 2007. 谱元法数值模拟地震波传播[J]. 防灾减灾工程学报, 27(4): 470-477.
杨祥森, 2009. 粘弹介质中有限差分法地震波场数值模拟. 见: 中国地球物理2009. 北京: 中国地球物理学会.
张冬丽, 陶夏新, 周正华, 2004. 近场地震动格林函数的解析法与数值法对比研究[J]. 地震工程学报, 26(3): 199-205.
张勇, 冯万鹏, 许力生, 等, 2008. 2008年汶川大地震的时空破裂过程[J]. 中国科学D辑:地球科学, 38(10): 1186-1194.
Andrews D. J., Ben-Zion Y., 1997. Wrinkle-like slip pulse on a fault between different materials[J]. Journal of Geophysical Research, 102(B1): 553-571. DOI:10.1029/96JB02856
Burridge R., Levy C., 1974. Self-similar circular shear cracks lacking cohesion[J]. Bulletin of the Seismological Society of America, 64(6): 1789-1808.
Day S. M., Dalguer L. A., Lapusta N., et al, 2005. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture[J]. Journal of Geophysical Research, 110(B12): B12307. DOI:10.1029/2005JB003813
Duan B. C., Oglesby D. D., 2006. Heterogeneous fault stresses from previous earthquakes and the effect on dynamics of parallel strike-slip faults[J]. Journal of Geophysical Research, 111(B5): B05309.
Fletcher J. B., Spudich P., 1998. Rupture characteristics of the three M~4.7 (1992-1994) Parkfield earthquakes[J]. Journal of Geophysical Research, 103(B1): 835-854. DOI:10.1029/97JB01797
Harris R. A., Archuleta R. J., Day S. M., 1991. Fault steps and the dynamic rupture process:2-D numerical simulations of a spontaneously propagating shear fracture[J]. Geophysical Research Letters, 18(5): 893-896. DOI:10.1029/91GL01061
Harris R. A., 2004. Numerical simulations of large earthquakes:dynamic rupture propagation on heterogeneous faults[J]. Pure and Applied Geophysics, 161(11-12): 2171-2181.
Komatitsch D., Tromp J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation[J]. Geophysical Journal International, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x
Komatitsch D., Barnes C., Tromp J., 2000. Simulation of anisotropic wave propagation based upon a spectral element method[J]. Geophysics, 65(4): 1251-1260. DOI:10.1190/1.1444816
Komatitsch D., Tromp J., 2002. Spectral-element simulations of global seismic wave propagation-Ⅱ. Three-dimensional models, oceans, rotation and self-gravitation[J]. Geophysical Journal International, 150(1): 303-318.
Komatitsch D., Liu Q., Tromp J., et al, 2004. Simulations of ground motion in the Los Angeles basin based upon the spectral-element method[J]. Bulletin of the Seismological Society of America, 94(1): 187-206. DOI:10.1785/0120030077
Madariaga R., Olsen K. B., 2000. Criticality of rupture dynamics in 3-D[J]. Pure and Applied Geophysics, 157(11): 1981-2001. DOI:10.1007/PL00001071
Oglesby D. D., Archuleta R. J., Nielsen S. B., 2000a. Correction to "dynamics of dip-slip faulting:explorations in two dimensions" by David D. Oglesby, Ralph J. Archuleta, and Stefan B. Nielsen[J]. Journal of Geophysical Research, 105(B8): 19149-19150.
Oglesby D. D., Archuleta R. J., Nielsen S. B., 2000b. The three-dimensional dynamics of dipping faults[J]. Bulletin of the Seismological Society of America, 90(3): 616-628. DOI:10.1785/0119990113
Oglesby D. D., Day S. M., 2001. The effect of fault geometry on the 1999 Chi-Chi (Taiwan) earthquake[J]. Geophysical Research Letters, 28(9): 1831-1834. DOI:10.1029/2000GL012043
Patera A. T., 1984. A spectral element method for fluid dynamics:laminar flow in a channel expansion[J]. Journal of Computational Physics, 54(3): 468-488. DOI:10.1016/0021-9991(84)90128-1
Richards P. G., 1976. Dynamic motions near an earthquake fault a three-dimensional solution[J]. Bulletin of the Seismological Society of America, 66(1): 1-32.
Ronquist E. M., Patera A. T., 1988. A Legendre spectral element method for the incompressible Navier-Stokes equations. The 7th GAMM-Conference on Numerical Methods in Fluid Mechanics.
Seriani G., Priolo E., Carcione J., et al., 1992. High-order spectral element method for elastic wave modeling. In:62nd Annual International Meeting, SEG Technical Program Expanded Abstracts. Society of Exploration Geophysicists, 1285-1288.
Zhang H. M., Chen X. F., 2006. Dynamic rupture on a planar fault in three-dimensional half space-Ⅰ. Theory[J]. Geophysical Journal International, 164(3): 633-652. DOI:10.1111/gji.2006.164.issue-3
Zhang W. B., Tomotaka I., Kojiro I., et al, 2004. Dynamic rupture process of the 1999 Chi-Chi, Taiwan, earthquake[J]. Geophysical Research Letters, 31(10): L10605.
Zhang Z. G., Zhang W., Chen X. F., 2014. Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics[J]. Geophysical Journal International, 199(2): 860-879. DOI:10.1093/gji/ggu308
Zhou H., Jiang H., 2015. A new time-marching scheme that suppresses spurious oscillations in the dynamic rupture problem of the spectral element method:the weighted velocity Newmark scheme[J]. Geophysical Journal International, 203(2): 927-942. DOI:10.1093/gji/ggv341


The Effect of Layered Medium on Rupture Process
Qi Mengxue*, Zhou Hong*
(Institute of Geophysics, China Earthquake Administration, Beijing 100081, China)
Abstract

The propagation simulation of ground motion caused by dynamic rupture is one of the most frontier topics in the field of earthquake and earthquake engineering research. Through three-dimensional models based on the spectral element method, which combines the flexibility in terms of mesh just like the finite element method with the accuracy like the pseudo spectral method as well as the computation efficiency, seismic rupture process in layered medium is simulated with Weighted Velocity Newmark Scheme and multiple transmitting artificial boundary condition in TPV5 model of the SCEC/USGS Spontaneous Rupture Code Verification Project. The results of layered medium and uniform medium are compared in terms of dislocation of fracture surface, earthquake moment, time of rupture propagation, rising time and surface peak displacement. The analysis shows that the layered medium can exert a great influence on earthquake rupture propagation:①The dislocation of fracture surface in layered medium is larger than that in the uniform counterpart, and the earthquake moment in layered medium is about 1.3 times as that of uniform medium, which indicates that layered medium intensifies the release of energy during rupture processes. ②Layered medium slows down the propagation of rupture to ground surface, shortens the rising time of various points on the ground surface and amplifies surface ground motion responses. ③Layered medium results in the magnification of surface displacements, and changes the location of peak surface displacements together with stress anomaly on faults.④There are strong ground motion near the surface of layered medium. Our results suggested that the site effects of layered medium should be carefully considered in order to avoid underestimating the seismic hazard in areas with layered medium. These conclusions mentioned above can provide implications for the seismic hazard mitigation of the region with layered earth medium.



主办单位:中国地震台网中心
版权所有:《震灾防御技术》编辑部
地址:北京西城区三里河南横街5号,   邮编:100045
邮箱:zzfy2006@126.com   电话:59959251;59959139;59959462
访问人数:1161615
层状介质对破裂过程的影响
齐梦雪*, 周红*
《震灾防御技术》, DOI:10.11899/zzfy20170316