2011年MW9.0東日本大地震動(dòng)力學(xué)破裂過程的數(shù)值模擬

劉敦宇, 胡才博, 蔡永恩

1 北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系， 北京 100871 2 中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049 3 中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院， 北京 100049

?

2011年MW9.0東日本大地震動(dòng)力學(xué)破裂過程的數(shù)值模擬

劉敦宇1, 胡才博2,3, 蔡永恩1

1 北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系， 北京 100871 2 中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049 3 中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院， 北京 100049

力學(xué)上，地震可以看作在應(yīng)力場(chǎng)作用下由于斷層帶介質(zhì)的突然損傷或軟化導(dǎo)致的斷層帶失穩(wěn)事件.本文基于這個(gè)地震動(dòng)力學(xué)模型，利用一種可以模擬斷層大位錯(cuò)的有限元方法，研究了2011年MW9.0東日本大地震(Tohoku-Oki)的動(dòng)力學(xué)破裂過程.比較了無障礙體和具有不同剛度障礙體的斷層帶模型產(chǎn)生的斷層位移、位錯(cuò)和應(yīng)力降.主要結(jié)果表明，障礙體的存在并不明顯地改變障礙體區(qū)域的初始構(gòu)造應(yīng)力場(chǎng).對(duì)有障礙體情形，準(zhǔn)靜態(tài)結(jié)果顯示斷層上盤最大逆沖位移和最大剪切位錯(cuò)分別為51 m和58 m，均發(fā)生在海底表面海溝處，與無障礙體的結(jié)果(最大剪切位錯(cuò)約55 m)相比差別不大；下盤最大傾向位移(-10 m)并不與上盤最大值出現(xiàn)在同一位置，而是在障礙體處.障礙體處剪應(yīng)力降(約11 MPa)大于周圍非障礙體區(qū)域.障礙體處正應(yīng)力降的最大值約為3 MPa.模擬結(jié)果似乎不支持海山是導(dǎo)致本次地震異乎尋常大位錯(cuò)的原因，而傾向于斷層帶剪切剛度在地震過程中極度損傷或軟化.

東日本大地震； 斷層破裂過程； 斷層位移； 應(yīng)力降； 有限元方法

1 引言

發(fā)生于2011年3月11日的MW9.0東日本大地震(Tohoku-Oki)產(chǎn)生的巨大位錯(cuò)和海嘯，大大出乎地震研究者的預(yù)料.海內(nèi)外學(xué)者利用各種觀測(cè)資料對(duì)其同震位錯(cuò)以及破裂過程進(jìn)行了反演或聯(lián)合反演(Lay et al.，2011；Yamazaki et al.，2011；Yue and Lay，2011；Simons et al.，2011；Ide et al.，2011；Yoshida et al.，2011；Zhang et al.，2012)，給出了本次地震準(zhǔn)靜態(tài)位錯(cuò)分布以及震源破裂過程等重要信息.然而如Lay等(2011)指出：利用地震波、GPS或者海嘯等資料反演得到的同震位錯(cuò)差異較大，特別是分辨海溝附近斷層面上位錯(cuò)量的大小和分布，這是由于數(shù)據(jù)對(duì)于近海溝的位錯(cuò)分辨能力有限.因此靠近海溝的觀測(cè)數(shù)據(jù)就顯得尤為重要.重復(fù)多激光束跨海溝測(cè)深數(shù)據(jù)表明本次地震在海溝附近產(chǎn)生了巨大的位移量：上盤介質(zhì)朝向海溝水平位移量高達(dá)約50 m，垂直隆升約7到10 m(Fujiwara et al.，2011).準(zhǔn)靜態(tài)同震位錯(cuò)反演和震源破裂過程的運(yùn)動(dòng)學(xué)反演得到的最大位錯(cuò)從30 m到60 m不等.Lay等(2011)利用遠(yuǎn)震P波最小二乘反演得到斷層?xùn)|北部海溝附近最大位錯(cuò)達(dá)62 m.Yamazaki等人(2011)利用近場(chǎng)海嘯數(shù)據(jù)約束反演的結(jié)果顯示斷層面上最大位錯(cuò)約60 m，集中在震源附近以及海溝附近.Simons等(2011)利用GPS和深海壓力測(cè)量數(shù)據(jù)聯(lián)合反演了本次地震的同震位錯(cuò)，其結(jié)果顯示震源附近位錯(cuò)達(dá)到60 m.震源破裂過程的運(yùn)動(dòng)學(xué)反演中，Yue和Lay(2011)利用高頻GPS資料反演破裂過程，結(jié)果顯示位錯(cuò)較大的區(qū)域靠近海溝和震源，最大位錯(cuò)約60 m.Ide等(2011)利用寬頻帶地震波資料反演了本次地震的破裂過程，得到的最大累積位錯(cuò)值約30 m，靠近海溝.Yoshida等(2011)利用長(zhǎng)周期(20～200 s)強(qiáng)地面震動(dòng)數(shù)據(jù)反演破裂過程，結(jié)果顯示本次地震分3個(gè)階段，第一階段破裂在斷層北側(cè)深部發(fā)展；第二階段，淺部發(fā)生持續(xù)時(shí)間較長(zhǎng)且位錯(cuò)量較大的破裂；第三階段破裂在斷層南段發(fā)展.最大位錯(cuò)位于斷層北側(cè)靠近海溝處，約47 m.Zhang等(2012)利用遠(yuǎn)震P波垂直分量和GPS數(shù)據(jù)聯(lián)合反演了本次地震的破裂過程，其結(jié)果顯示最大位錯(cuò)位于海溝處，為30 m.最大應(yīng)力降達(dá)111 MPa，平均應(yīng)力降46 MPa.

利用數(shù)值方法對(duì)地震斷層破裂過程進(jìn)行動(dòng)力學(xué)正演能夠提供震源行為的原因.Duan(2012)利用TSN (Traction-at-Split-Nodes)方法結(jié)合摩擦破裂準(zhǔn)則正演了本次地震的自發(fā)破裂和地震波傳播過程，結(jié)果表明位于震源上方的大小約70 km×23 km的海山對(duì)整個(gè)地震破裂過程起著重要的作用.海山處位錯(cuò)值最大，約47 m，且產(chǎn)生了超過50 MPa的剪切應(yīng)力降.

JFAST(Japan Trench Fast Drilling Project)鉆孔的結(jié)果能夠提供關(guān)于斷層帶的材料和結(jié)構(gòu)更多的認(rèn)識(shí)(Chester et al.，2013).Chester等(2013)對(duì)本次地震的鉆孔巖樣進(jìn)行了分析，結(jié)果表明：本次地震斷層帶厚度約5 m，材料區(qū)別于圍巖介質(zhì)，層內(nèi)材料主要是粘土礦物,大位錯(cuò)主要發(fā)生在斷層帶薄層內(nèi).Wang和Kinoshita(2013)總結(jié)了JFAST得到的結(jié)果，認(rèn)為該區(qū)域斷層帶厚度薄且強(qiáng)度弱是這次地震能夠產(chǎn)生如此巨大位錯(cuò)的原因.

關(guān)于薄弱斷層帶對(duì)破裂過程影響的研究還比較少.為討論弱斷層帶的影響，本文將用彈性三維節(jié)理單元(殷有泉和張宏，1982；殷有泉，2014)來模擬斷層帶薄層，并在胡才博等(2009)發(fā)展的模擬地震引起的穩(wěn)態(tài)變形場(chǎng)的有限元方法基礎(chǔ)上，繼續(xù)研發(fā)了模擬斷層破裂動(dòng)力學(xué)過程的方法和程序，并用其模擬2011年東日本大地震的破裂過程和探討它異乎尋常大位錯(cuò)的動(dòng)力學(xué)機(jī)制.

2 有限元方法

野外地質(zhì)考察和地震導(dǎo)波反演表明斷層通常有一定的厚度(Li et al.，2002)，可以分為斷層核、破損區(qū)和原巖三個(gè)主要區(qū)域(Caine et al.，1996；Gudmundsson，2004).從大地震前后的地震波速異?？梢酝茢喟l(fā)震斷層的力學(xué)性質(zhì)在地震過程中發(fā)生了改變(Scholz，1990).因此地震可以看作在應(yīng)力場(chǎng)作用下，由斷層帶介質(zhì)的損傷或軟化導(dǎo)致的斷層帶突然失穩(wěn)事件.在這種地震動(dòng)力學(xué)模型中，損傷或軟化可以通過斷層帶介質(zhì)的彈性參數(shù)降低(損傷力學(xué)本構(gòu))或其內(nèi)摩擦系數(shù)和內(nèi)聚力的減小(塑性力學(xué)本構(gòu))實(shí)現(xiàn).地震斷層的主要滑動(dòng)集中在斷層核(相對(duì)斷層帶的寬度而言非常窄)，因此如果不關(guān)心斷層附近的變形場(chǎng)或應(yīng)力場(chǎng)時(shí)，可以把斷層視為沒有厚度的位錯(cuò)面(地震運(yùn)動(dòng)學(xué)模型)或由摩擦強(qiáng)度控制的接觸面(地震動(dòng)力學(xué)模型或接觸力學(xué)模型).無論是斷層帶模型還是位錯(cuò)或接觸模型，雖然角度不同，它們都可以模擬地震斷層滑動(dòng)現(xiàn)象.本文將使用斷層帶介質(zhì)彈性參數(shù)損傷的思想，研究2011年東日本大地震的斷層破裂動(dòng)力學(xué)過程，并探討其大位錯(cuò)的動(dòng)力學(xué)機(jī)制.

如果不考慮斷層的動(dòng)力學(xué)破裂過程，俯沖帶斷層地震引起的穩(wěn)態(tài)變形場(chǎng)可以用有限元方法模擬(胡才博等，2009)

(1)

(2)

當(dāng)斷層帶非常薄時(shí)，利用連續(xù)變形單元模擬斷層帶將會(huì)導(dǎo)致斷層帶內(nèi)外單元?jiǎng)偠认嗖顟沂?，使總體剛度矩陣條件數(shù)變差，較難模擬斷層大位移滑動(dòng).為了有效地模擬斷層大位移滑動(dòng)，可采用能夠模擬斷層位移間斷(位錯(cuò))的節(jié)理元(殷有泉和張宏，1982；殷有泉，2014)，其相當(dāng)于把斷層兩盤用彈簧連接起來.圖1為節(jié)理元示意圖，其沿傾向x′長(zhǎng)L1，沿走向y′長(zhǎng)L2，厚度為b.當(dāng)L1?b，L2?b時(shí)，節(jié)理元內(nèi)的應(yīng)變可近似表達(dá)為

εx′≈0,εy′≈0,γx′y′≈0,εz′≈Δw/b,

γz′x′≈Δu/b,γz′y′≈Δv/b,

(3)

式中Δu,Δv,Δw分別為節(jié)理元沿?cái)鄬觾A向、走向和法向的位移間斷量或位錯(cuò).

利用斷層坐標(biāo)系x′，y′，z′定義的方向，斷層上下盤的位移間斷量可以表示為

Δu=[ΔuΔvΔw]T

=[u上-u下v上-v下w上-w下]T，

(4)

u上,v上,w上u下,v下,w下分別為節(jié)理元上下盤x′，y′，z′方向的位移，上角標(biāo)“T”表示轉(zhuǎn)置.本研究選用8節(jié)點(diǎn)六面體節(jié)理單元，每個(gè)節(jié)點(diǎn)有3個(gè)自由度，節(jié)點(diǎn)編號(hào)如圖1所示，1,3,5,7號(hào)節(jié)點(diǎn)屬于上盤，對(duì)應(yīng)的下盤節(jié)點(diǎn)編號(hào)為2,4,6,8.

圖1 節(jié)理元示意圖(x′，y′，z′)為斷層滑動(dòng)坐標(biāo)系. 節(jié)理元沿x′方向長(zhǎng)L1，沿y′方向長(zhǎng)L2，厚度b，其中L1?b，L2?b.1,3,5,7為上盤面節(jié)點(diǎn)編號(hào)，對(duì)應(yīng)的下盤面節(jié)點(diǎn)編號(hào)為2,4,6,8.Fig.1 Schematic diagram of joint element in local fault coordinate system (x′, y′, z′)L1, L2 and b are its length, width and thickness, respectively. Nodal numbers of the element on the hanging wall fault surface are 1,3,5,7 and on footwall are 2,4,6,8.

Δu=BUe,

(5)

h5I-h5Ih7I-h7I],

節(jié)理元的應(yīng)力可以由位移間斷量表示為

σ=DIIΔu,

(6)

σ=[τz′x′τz′y′σz′]T,

(7)

其中kz′x′,kz′y′,kz′z′分別為節(jié)理元傾向和走向(x′,y′)剪切面剛度以及法向(z′)面剛度，Gz′x′,Gz′y′,Ez′z′分別為節(jié)理元相應(yīng)方向的剪切模量和楊氏模量.

由有限元方法(殷有泉，2014)，可以得到面積為s，厚度為b的節(jié)理元的單元?jiǎng)偠染仃嚍?/p>

Ke=b?sBTDBds.

(8)

(9)

M,C分別為有限元模型的總體質(zhì)量矩陣及總體阻尼矩陣，Kt為時(shí)刻t總體剛度矩陣，由于斷層破裂區(qū)域隨時(shí)間變化，所以Kt也是隨時(shí)間變化.對(duì)于阻尼矩陣取為Rayleigh阻尼形式，即C=αM+βK，其中α主要控制低頻成分的衰減，β主要控制高頻成分的衰減.在實(shí)際巖石破裂過程中，大量的能量消耗在摩擦和破裂過程中，因此在利用彈性損傷模型時(shí)，需要在斷層帶介質(zhì)上施加較大阻尼來模擬該能量損耗過程.

在獲得構(gòu)造運(yùn)動(dòng)加載引起的斷層震前變形場(chǎng)U0或與其對(duì)應(yīng)的背景應(yīng)力場(chǎng)σ0后，利用公式(9)可以計(jì)算出在此應(yīng)力場(chǎng)控制下，地震(斷層損傷或軟化)過程中的動(dòng)態(tài)位移場(chǎng)Ut，利用(5)和(6)式可以求出斷層上的位錯(cuò)和應(yīng)力演化，由公式(10)可以得到地震應(yīng)力降隨時(shí)間的變化

Δσ(t)=σ(t)-σ0.

(10)

3 2011年MW 9.0東日本大地震的數(shù)值模擬

3.1 有限元模型及初始構(gòu)造應(yīng)力場(chǎng)

本文利用上述有限元方法模擬了2011年3月11日的MW9.0東日本大地震的破裂過程.構(gòu)建的有限元模型如圖2a所示.研究區(qū)域500 km×650 km×90 km，x軸方向?yàn)闁|偏北67°，y軸方向?yàn)闁|偏南23°，z軸方向?yàn)榇怪钡乇硐蛏?俯沖帶(圖2a中淺藍(lán)色介質(zhì)層)界面的幾何模型參考日本海溝深度和歷史地震剖面(Rhea et al.，2010)，設(shè)海溝深度7 km，在海底0～40 km深度，俯沖帶傾角為11°；海底深度40 km以下俯沖帶傾角為30°.斷層帶(圖2a紅框區(qū)域)設(shè)為夾在上盤介質(zhì)和俯沖板片之間的5 m厚度的薄層(Chester et al.，2013).本次地震的斷層破裂區(qū)如圖2a黑框區(qū)域所示，其沿走向長(zhǎng)250 km，沿傾向?qū)?30 km.只有在破裂區(qū)內(nèi)的斷層帶材料在地震過程中產(chǎn)生損傷或軟化.

采用六面體網(wǎng)格剖分研究區(qū)域，并在靠近斷層區(qū)域進(jìn)行網(wǎng)格加密.斷層帶單元在斷層面內(nèi)的尺寸為1 km×1 km.模型上表面(地表和洋底)為自由表面；模擬初始構(gòu)造應(yīng)力場(chǎng)時(shí)，東南側(cè)邊界施加推擠位移邊界條件，其他邊界面施加滾筒約束.圖2a中，紅星處為震源，深度26 km.測(cè)線AB位于模型上表面且垂直于斷層跡線(海溝)；測(cè)線CD位于斷層面(地表和洋底)且垂直于斷層跡線.

有學(xué)者把本次大地震異常大的滑動(dòng)歸因于海山障礙體的破裂(Yoshida et al.，2011；Duan，2012).從斷層帶介質(zhì)的彈性參數(shù)的角度看，障礙體(可能為海山)可視為剛度較大巖石材料，為討論其對(duì)斷層位移、位錯(cuò)以及應(yīng)力降的影響，我們?cè)O(shè)計(jì)了不存在障礙體(模型A)和兩個(gè)存在障礙體(模型B和C)的三種斷層模型.對(duì)于模型B和C，在斷層帶破裂區(qū)設(shè)置了一個(gè)尺度為70 km×31 km的障礙體(如圖2c斷層破裂區(qū)中間的淺藍(lán)色區(qū)域)，走向方向從290～360 km，傾向從-26.2～-57.6 km.

這三個(gè)模型的斷層帶外部的材料參數(shù)均相同(PREM模型)，如表1所示.鑒于JFAST鉆孔得到的淺部(海底下812 m)斷層帶主要由粘土類礦物組成(Chester et al.，2013)，我們參考Mondol等(2007)對(duì)于粘土聚合物壓縮試驗(yàn)的研究結(jié)果設(shè)置斷層帶的材料參數(shù).為區(qū)別斷層帶的破裂和未破裂區(qū)，本文假設(shè)未破裂斷層帶的材料參數(shù)高于破裂區(qū).根據(jù)Mondol等人的研究結(jié)果，如果考慮孔隙流體影響，斷層帶內(nèi)的彈性常數(shù)大小與有效應(yīng)力有關(guān)：對(duì)于模型A，選取破裂區(qū)的材料為垂直有效應(yīng)力10 MPa時(shí)的楊氏模量2.0 GPa和剪切模量0.7 GPa；未破裂區(qū)的楊氏模量和剪切模量分別為20 GPa和7 GPa.參考Duan(2012)利用比非障礙體斷層區(qū)域較高的靜摩擦系數(shù)、較低的孔隙壓力以及較高的初始應(yīng)力來模擬障礙體.本研究的模型B中，假設(shè)障礙體仍為粘土聚合物，但其材料參數(shù)選擇垂直有效應(yīng)力較高時(shí)(50 MPa)的剪切模量和楊氏模量，分別為1.75 GPa和5 GPa，其剪切模量為破裂區(qū)非障礙體的2.5倍；在模型C中，假設(shè)障礙體為海山，物質(zhì)來源為海洋地幔，其楊氏模量和剪切模量分別取為151 GPa和59.3 GPa(參見表1)，其中剪切模量為破裂區(qū)非障礙體的84倍.除了障礙體材料參數(shù)的差別，模型B和C的區(qū)別還在于模型B的障礙體在地震中剪切面剛度損傷更嚴(yán)重.

表1 有限元模型材料參數(shù)Table 1 Material parameters in the finite element model

本文將主要討論模型B的結(jié)果.為方便起見，如無特殊說明，后文中的斷層或斷層帶專指斷層帶破裂區(qū).

為了模擬地震斷層的破裂過程，需要模擬震前的構(gòu)造應(yīng)力場(chǎng).假設(shè)地震遵循彈性回跳學(xué)說，并且本次大地震把構(gòu)造運(yùn)動(dòng)積累起的彈性能全部釋放(Ide et al.，2011)，那么斷層帶上初始構(gòu)造剪應(yīng)力應(yīng)該與剪應(yīng)力降大小相近.根據(jù)Simons等(2011)估算的本次地震應(yīng)力降約為2到10 MPa，當(dāng)在本模型東南側(cè)邊界施加累積推擠位移240 m時(shí)(如圖2a所示)，斷層帶初始構(gòu)造剪應(yīng)力在10 MPa左右.假設(shè)太平洋板塊一直以80～85 mm·a-1(Simons et al.，2011)的速率向日本島俯沖，要在斷層帶積累起這么大的初始剪應(yīng)力需要約2900年.因?yàn)楹Ｋ静豢蓧嚎s，在海底表面，地震引起的海水載荷變化與地震斷層應(yīng)力降相比可以忽略，因此，地震動(dòng)力載荷在此處可視為零(表面自由).

圖2c和2d分別為由邊界推擠240 m得到的斷層上盤面初始剪應(yīng)力和正應(yīng)力(正應(yīng)力和剪應(yīng)力的正方向如圖2b所示).由圖2c可見，在構(gòu)造運(yùn)動(dòng)加載下，斷層破裂區(qū)上積累起的剪應(yīng)力隨深度增加而增大.其中海底海溝處(水深度7 km)剪應(yīng)力最小，約為9 MPa；最大值位于斷層深部，約為15 MPa.障礙體上的初始剪應(yīng)力略大于周圍斷層帶，特別是在障礙體邊界上，剪應(yīng)力(13.1 MPa)高于周圍非障礙體斷層帶0.5 MPa，這是由于該處材料不連續(xù)所致.從圖2d上可看到，斷層面法向呈擠壓狀態(tài)，在障礙體處正應(yīng)力較周圍大，最大值約3.3 MPa.傾向位于-20～-60 km的非障礙體區(qū)域的正應(yīng)力較大.此外，從圖中可以看到材料分界引起的應(yīng)力集中.

兩種障礙體模型(B和C)的初始剪應(yīng)力僅在障礙體邊界處與模型A區(qū)別較大.在障礙體邊界處，模型C的初始剪應(yīng)力比模型A約大1 MPa，障礙體內(nèi)部初始應(yīng)力狀態(tài)與模型A相差不多，這是由于在本研究所采用的構(gòu)造應(yīng)力場(chǎng)的加載方式下，障礙體雖然剛度較大，但其應(yīng)變較小.

3.2 模擬結(jié)果及分析

地震過程通過逐步下降斷層帶的剪切剛度模擬.剛度下降多少可以利用觀測(cè)資料或反演得到的斷層位錯(cuò)量來約束.Yoshida等(2011)指出破裂傳播速度在此次地震過程中約3～4 km·s-1.本文假設(shè)斷層破裂以圓盤形式傳播且破裂傳播速率為3 km·s-1.模擬所選時(shí)間步長(zhǎng)為0.5 s，以保證在每一個(gè)時(shí)間步長(zhǎng)內(nèi)均有斷層帶單元破裂.在初始構(gòu)造應(yīng)力場(chǎng)作用下，當(dāng)破裂傳至某一斷層帶單元時(shí)，即降低其傾向剪切面剛度，降低后不再恢復(fù).通過給斷層帶上下盤質(zhì)點(diǎn)施加大阻尼來模擬地震過程中的能量耗散.

數(shù)值試驗(yàn)表明：對(duì)于模型B，若產(chǎn)生約58 m的最大剪切位錯(cuò)(Lay et al.，2011；Yue and Lay，2011)，地震斷層帶非障礙體單元的傾向剪切面剛度需損傷至0.1 MPa/m(剪切模量損傷至0.5 MPa)，障礙體傾向剪切面剛度需損傷至0.02 MPa/m(剪切模量損傷至0.1 MPa).在模型A和C中，所有地震斷層帶單元傾向剪切面剛度均損傷至0.1 MPa/m.非斷層帶介質(zhì)阻尼系數(shù)α=0.1，β=0.03.由于斷層帶單元很薄，所以需施加較大阻尼模擬破裂和摩擦過程中的耗散，數(shù)值試驗(yàn)給出斷層帶單元阻尼系數(shù)α=1000.0，β=0.03.

圖3a給出了破裂開始后六個(gè)不同時(shí)刻的傾向剪切位錯(cuò)快照(正值表示上盤介質(zhì)相對(duì)下盤沿傾向向上運(yùn)動(dòng)).在33 s時(shí)，破裂剛剛到達(dá)海溝；在150 s左右，破裂趨于穩(wěn)定，這個(gè)結(jié)果與運(yùn)動(dòng)學(xué)反演結(jié)果一致(Zhang et al.，2012).對(duì)比34～70 s的快照?qǐng)D，可以看到海溝處的位錯(cuò)在36 s內(nèi)從不到10 m增大到40多米；在隨后70～150 s內(nèi)，海溝處位錯(cuò)逐漸增加到58 m.從70～150 s的快照?qǐng)D中可以看到，最大值剪切位錯(cuò)為58 m，位于海溝處；障礙體處及位于障礙體前方的斷層帶淺部的位錯(cuò)值均比較大.如此大的位錯(cuò)量可能是因?yàn)檎鸷髷鄬訋Ш苋?傾向剪切面剛度只有0.1 MPa/m)和斷層傾角很低，當(dāng)破裂破穿海溝后，下盤介質(zhì)無法有效地約束上盤介質(zhì).

圖3b給出了破裂開始后六個(gè)不同時(shí)刻的法向位錯(cuò)快照(正值表示拉張).動(dòng)態(tài)最大幅值約0.07 m，與斷層傾向位錯(cuò)相比非常小.法向位錯(cuò)穩(wěn)定后的幅值不超過0.02 m.斷層上半部分震后法向位錯(cuò)為拉張，而下半部分為壓縮.障礙體處由于法向剛度較大，所以位錯(cuò)較小，但仍為拉張狀態(tài).

根據(jù)位錯(cuò)分布以及斷層圍巖剪切模量可以計(jì)算地震矩，其計(jì)算公式(Aki，1966)為

M0=μAD,

其中μ為斷層圍巖剪切模量，D為斷層平均位錯(cuò)，A為斷層面積.對(duì)于位錯(cuò)分布不均勻且斷層圍巖材料不均勻的情況，可計(jì)算地震矩密度m0=μD，并對(duì)整個(gè)斷層面積分得到總地震矩.在本研究中，在一些斷層區(qū)域上下盤圍巖的剪切模量有可能不同，此時(shí)取其平均值.圖4給出了地震矩密度分布.對(duì)斷層面積分得到本次地震的地震矩M0=3.98×1022N·m.由矩震級(jí)計(jì)算公式(Kanamori, 1977)

得到矩震級(jí)MW=9.0.該結(jié)果與USGS給出的矩震級(jí)一致.

地震在150s時(shí)已經(jīng)基本停止，圖5a—5d展示了震后200s時(shí)斷層帶的位移(位移正方向如圖2b所示)和位錯(cuò)分布.圖5a為斷層上盤面逆沖位移分布圖，最大值51m，位于海溝處.圖5b為斷層下盤面傾向位移分布圖，最大值-10m，位于障礙體處；斷層邊界處位移值為正，說明該處介質(zhì)被上盤拖曳向海溝方向運(yùn)動(dòng).圖5c和5d分別為斷層上盤面和下盤面法向位移分布，最大值2.6m，位于斷層底邊界；斷層淺部沿?cái)鄬酉卤P面外法線方向平均位移了約2m，斷層深部沿?cái)鄬由媳P面外法線方向平均位移了約2m.

圖2 (a)有限元模型; (b)斷層上下盤面應(yīng)力和位移的正方向示意圖; (c)斷層上盤面初始傾向剪應(yīng)力; (d)斷層上盤面初始法向正應(yīng)力Fig.2 (a) Finite element model; (b) Definition of the positive directions of displacements and stresses on the fault surfaces; (c) Initial shear stress on the hanging wall fault surface; (d) Initial normal stress on the hanging wall fault surface

圖3 斷層位錯(cuò)快照(a)傾向剪切位錯(cuò); (b)法向位錯(cuò).Fig.3 Snapshots of dislocation(a) Shear dislocation; (b) Normal dislocation.

圖5 震后200 s斷層面上的位移和位錯(cuò)分布(a)和(b)分別是斷層上盤和下盤剪切位移；(c)和(d)分別是斷層上盤和下盤法向位移；(e)、(f)和(g)分別是斷層上盤、下盤AB測(cè)線上的6個(gè)測(cè)點(diǎn)的剪切位移和剪切位錯(cuò)隨時(shí)間變化.Fig.5 Displacement and dislocation distributions on the fault surfaces at 200 s(a) and (b) are shear displacements of the hanging wall and footwall, respectively; (c) and (d), normal displacements of the hanging wall and footwall, respectively; (e), (f) and (g), shear displacements of the hanging wall, footwall and shear dislocation at six observation sites on Profile AB, respectively.

圖6 CD測(cè)線在不同時(shí)刻的位移(a) 水平位移; (b) 垂向位移. 橫坐標(biāo)‘0’和負(fù)值分別表示斷層在海溝的出露點(diǎn)和斷層上盤.Fig.6 Displacements on the profile CD at different times(a) Horizontal displacements; (b) Vertical displacements. ‘0’ and negative coordinate value represent the trench and the hanging wall, respectively.

圖4 地震矩密度分布Fig.4 Seismic moment density distribution

圖5e—5g展示了圖2a中位于斷層上的AB測(cè)線上6個(gè)測(cè)點(diǎn)的位移和位錯(cuò)隨時(shí)間變化(位移正方向如圖2b所示).0號(hào)測(cè)點(diǎn)為震源，深度26km；1號(hào)測(cè)點(diǎn)深度22km；3、2號(hào)測(cè)點(diǎn)分別為障礙體的上下邊界所在深度，分別為12km和18km；4號(hào)測(cè)點(diǎn)深度8.9km；5號(hào)測(cè)點(diǎn)位于海溝處，深度7km.圖5e為斷層上盤面逆沖位移，海溝處位移值最大，為51m.圖5f為斷層下盤面傾向剪切位移，6個(gè)測(cè)點(diǎn)中3號(hào)(深度為-12km)測(cè)點(diǎn)處位移值最大，略小于-10m.圖5g為斷層傾向剪切位錯(cuò)，海溝處測(cè)點(diǎn)位錯(cuò)最大，為58m.比較六條曲線的斜率可以看出，斷層初始滑動(dòng)速率(斷層上一點(diǎn)在單位時(shí)間內(nèi)位錯(cuò)的大小)隨著破裂向海溝傳播逐漸增大：破裂初始階段為0.7m·s-1(震源，0號(hào)測(cè)點(diǎn))，逐漸增大到2m·s-1(3號(hào)測(cè)點(diǎn))，當(dāng)破裂達(dá)到海溝時(shí)，高達(dá)6m·s-1(5號(hào)測(cè)點(diǎn))，這是由于破裂前端動(dòng)力積累的應(yīng)力集中導(dǎo)致的雪崩式的結(jié)果.

圖6為位于地表和洋底的CD測(cè)線(圖2a)在不同時(shí)刻水平方向和垂向位移.從圖中可以看到非常明顯的彈性回跳的動(dòng)態(tài)過程：隨著破裂向海溝擴(kuò)展，上盤地表最大水平位移和垂直位移的位置向海溝移動(dòng)；在33s時(shí)，破裂即將傳至海溝，該處介質(zhì)變形達(dá)到最大；在34s時(shí)海溝破裂，上下盤介質(zhì)即刻產(chǎn)生顯著相對(duì)位移并隨時(shí)間迅速增大.從圖6a可以看到，地震引起的主要地表變形來自上盤：在海溝處，上盤在150s時(shí)垂直于海溝向下盤水平位移了近50m，而下盤的位移量?jī)H約8m.從圖6b可以看到：在海溝處，上盤在150s時(shí)垂直隆升達(dá)12m，足以產(chǎn)生巨大的海嘯.本文的結(jié)果與重復(fù)多激光束跨海溝測(cè)深(Fujiwaraetal.，2011)得到的結(jié)果比較一致.

圖7a和7b展示了震后200s時(shí)斷層上盤面剪應(yīng)力降和正應(yīng)力變化分布(應(yīng)力變化的正方向如圖2b所示).應(yīng)力變化按照公式(10)計(jì)算，正值表示應(yīng)力降，負(fù)值表示應(yīng)力增加.剪應(yīng)力降隨著深度的增大而增大，這主要是因?yàn)閿鄬用嫔铣跏技魬?yīng)力隨深度增加而增大且本文假設(shè)斷層帶材料損傷量相同.障礙體處應(yīng)力降(約11MPa)比周圍非障礙體處(約9MPa)大.由于斷層帶法向剛度沒有損傷，法向應(yīng)力變化趨勢(shì)同法向位錯(cuò)：以傾向深度-80km線為分界，傾向位于0～-80km的斷層帶區(qū)域法向正應(yīng)力釋放；-80～-130km的區(qū)域法向正應(yīng)力增加.障礙體區(qū)域法向正應(yīng)力釋放較大，約2～3MPa.

圖7c和7d分別為AB測(cè)線上6個(gè)測(cè)點(diǎn)處的剪應(yīng)力和正應(yīng)力變化隨時(shí)間變化圖.從圖7c可以看到，動(dòng)態(tài)應(yīng)力降大于準(zhǔn)靜態(tài)應(yīng)力降.當(dāng)破裂傳播到新的破裂部位之前，由于斷層破裂的端部效應(yīng)，該處會(huì)產(chǎn)生較大的應(yīng)力增加(負(fù)值).斷層破裂過程中，正是由于這種迅速增加的端部動(dòng)態(tài)應(yīng)力，使斷層破裂得以持續(xù)，直到碰到高強(qiáng)度的障礙體或者初始應(yīng)力很低的斷層帶介質(zhì)使破裂停止.圖7d表明，動(dòng)態(tài)正應(yīng)力變化的幅值振動(dòng)幅度達(dá)10MPa；地震結(jié)束后，傾向位于0～-80km區(qū)域的測(cè)點(diǎn)法向應(yīng)力增加約1MPa，而位于-80～-130km區(qū)域的測(cè)點(diǎn)法向應(yīng)力降約2MPa.

圖8a—8d分別給出了無障礙體情況下(模型A)的剪切位錯(cuò)、法向位錯(cuò)、斷層上盤面剪切應(yīng)力降以及法向應(yīng)力降(負(fù)值表示應(yīng)力增加).

圖8e—8h分別為有障礙體情況(模型B)與模型A的剪切位錯(cuò)、法向位錯(cuò)、斷層上盤面剪切應(yīng)力降以及法向應(yīng)力降(負(fù)值表示應(yīng)力增加)的差值分布(用模型B減去模型A).結(jié)果顯示：障礙體處(模型B)的剪切位錯(cuò)要比A模型大3～5m，并且障礙體前方(向海溝方向)位錯(cuò)也相應(yīng)增大.由于障礙體處初始剪應(yīng)力沒有明顯高于無障礙體情況，上述3～5m的差位錯(cuò)是由障礙體更多的材料損傷，破裂后上下盤更容易滑動(dòng)導(dǎo)致.相應(yīng)的障礙體處的剪應(yīng)力降要比沒有障礙體的情形約大3～5MPa；相應(yīng)的法向應(yīng)力降在障礙體上半部分及其前方(向海溝方向)比模型A約大1MPa.

圖8i—8l分別為有障礙體情況(模型C)與模型A的剪切位錯(cuò)、法向位錯(cuò)、斷層上盤面剪切應(yīng)力降以及法向應(yīng)力降的差值分布(用模型C減去模型A).結(jié)果顯示：當(dāng)障礙體和非障礙體斷層帶傾向剪切面剛度均損傷至0.1MPa/m時(shí)，障礙體處(模型C)的剪切位錯(cuò)并沒有明顯大于模型A，最大位錯(cuò)差只有0.1m.模型C的剪應(yīng)力降僅在障礙體邊界處比模型A約大2MPa.

綜上所述，當(dāng)斷層帶材料均損傷至同一值時(shí)，即便障礙體的材料參數(shù)為非障礙體材料參數(shù)84倍(模型C)，其震后斷層位錯(cuò)分布并沒有明顯區(qū)別于無障礙體模型(模型A).如果障礙體處斷層帶材料損傷后，剩余剪切剛度更低的話(模型B)，障礙體處及障礙體前方的位錯(cuò)僅比無障礙體模型大3～5m.因此障礙體可能不是導(dǎo)致東日本大地震異乎尋常大滑動(dòng)的主要原因，即使沒有障礙體，當(dāng)斷層帶剪切面剛度損傷至0.1MPa/m時(shí)，仍會(huì)產(chǎn)生高達(dá)55m的最大剪切位錯(cuò).

4 結(jié)論與討論

本文利用能夠模擬斷層大位移間斷的三維彈性節(jié)理單元探討了2011年東日本大地震的斷層破裂動(dòng)力學(xué)過程.對(duì)于有、無障礙體模型，主要認(rèn)識(shí)如下：

這次大地震的斷層初始剪應(yīng)力隨深度增加而增大.對(duì)于障礙體的剪切模量和非障礙體斷層帶的差別不大時(shí)(模型B)，障礙體邊緣處初始剪應(yīng)力約為13.1MPa，略高于同一深度無障礙體斷層帶區(qū)域；障礙體法向初始擠壓應(yīng)力最大值約3.3MPa.對(duì)于障礙體的剪切模量大大高于非障礙體斷層帶時(shí)(模型C)，障礙體邊緣初始剪應(yīng)力只比模型B約大0.5MPa.即不同剛度的障礙體不會(huì)顯著改變斷層面上的初始應(yīng)力狀態(tài).

對(duì)于模型B，地震引起的傾向剪切位錯(cuò)和斷層上盤面逆沖位移最大值在海溝處，分別為58m和51m.地震引起的斷層下盤面傾向位移最大值不在海溝，而位于障礙體區(qū)域，為-10m.地震引起的洋底變形主要由上盤產(chǎn)生.在海溝處，上盤介質(zhì)向海溝方向水平位移了近50m，垂向隆升近12m，模擬結(jié)果與重復(fù)多激光束跨海溝測(cè)深數(shù)據(jù)結(jié)果一致.斷層帶破裂后，其殘余剛度很小，下盤介質(zhì)無法有效地約束上盤介質(zhì)，而且本次地震斷層傾角很低，這可能是本次地震能夠產(chǎn)生如此巨大位錯(cuò)的主要原因.障礙體的影響有限.

對(duì)于模型B，障礙體處剪應(yīng)力降大于其周圍非障礙體斷層帶區(qū)域.如果考慮了斷層厚度，震后法向正應(yīng)力在斷層深部增加，而淺部減小.

如果假設(shè)太平洋板塊平均以80～85mm·a-1速率向日本島俯沖以及本次地震的斷層破裂區(qū)構(gòu)造剪應(yīng)力全部釋放，構(gòu)建初始構(gòu)造應(yīng)力場(chǎng)的結(jié)果表明這個(gè)地區(qū)MW9.0大地震的復(fù)發(fā)期約2900年.

利用簡(jiǎn)單的斷層帶彈性材料損傷或軟化模型模擬地震過程還存在著一定的局限性，例如在破裂前端，存在很大的應(yīng)力集中.真實(shí)巖石材料在應(yīng)力達(dá)到其屈服強(qiáng)度后，將會(huì)產(chǎn)生塑性變形，因此利用彈塑性介質(zhì)模擬斷層帶可能更加符合實(shí)際，并能通過破裂準(zhǔn)則判斷和預(yù)測(cè)破裂的發(fā)展，這也將是下一步工作需要考慮的方向.

圖7 震后200 s斷層上盤面應(yīng)力變化(a) 剪切應(yīng)力降；(b) 正應(yīng)力變化; (c)和(d)分別是AB測(cè)線上6個(gè)測(cè)點(diǎn)剪切應(yīng)力降和正應(yīng)力變化隨時(shí)間變化曲線.Fig.7 Stress changes on the hanging wall fault surface at 200 s(a) and (b) are distributions of shear stress drop and normal stress change, respectively; (c) and (d), shear stress drop and normal stress change at six observation sites on Profile AB, respectively.

圖8 有、無障礙體模型的斷層位錯(cuò)和應(yīng)力降比較(a)—(d)分別是A模型 (無障礙體)的剪切位錯(cuò)、法向位錯(cuò)、斷層面剪應(yīng)力降和斷層面法向正應(yīng)力降；(e)—(h)分別是模型B與模型A的剪切位錯(cuò)之差、法向位錯(cuò)之差、斷層面剪應(yīng)力降之差和斷層面法向應(yīng)力降之差；(i)—(l)分別是模型C與模型A的剪切位錯(cuò)之差、法向位錯(cuò)之差、斷層面剪應(yīng)力降之差和斷層面法向應(yīng)力降之差Fig.8 Comparisons of the fault dislocation and stress drop among the models with and without the asperityModel A (without asperity): (a) shear dislocation; (b) normal dislocation; (c) shear stress drop; (d) normal stress drop. Comparisons of differences of the shear dislocation (e), the normal dislocation (f), the shear stress drop (g), and the normal stress drop (h) between Model B and Model A; Comparisons of differences of the shear dislocation (i), the normal dislocation (j), the shear stress drop (k), and the normal stress drop (l) between Model C and Model A.

Aki K. 1966. Generation and propagation of G waves from the Niigata earthquake of June 14, 1964. Part 2. Estimation of earthquake moment, released energy, and stress-strain drop from G-waves spectrum.BulletinoftheEarthquakeResearchInstitute, 44: 73-88.

Caine J S, Evans J P, Forster C B. 1996. Fault zone architecture and permeability structure.Geology, 24(11): 1025-1028.

Chester F M, Rowe C, Ujiie K, et al. 2013. Structure and composition of the plate-boundary slip zone for the 2011 Tohoku-Oki earthquake.Science, 342(6163): 1208-1211, doi: 10.1126/science.1243719.

Duan B C. 2012. Dynamic rupture of the 2011MW9.0 Tohoku-Oki earthquake: Roles of a possible subducting seamount.J.Geophys.Res., 117: B05311, doi: 10.1029/2011JB009124.

Fujiwara T, Kodaira S, No T, et al. 2011. The 2011 Tohoku-Oki earthquake: displacement reaching the trench axis.Science, 334(6060): 1240, doi: 10.1126/science.1211554.

Gudmundsson A. 2004. Effects of Young′s modulus on fault displacement.ComptesRendusGeoscience, 336(1): 85-92.

Hu C B, Zhou Y J, Cai Y E. 2009. A new finite element model in studying earthquake triggering and continuous evolution of stress field.Sci.ChinaSer.D-EarthSci., 52(7): 994-1004, doi: 10.1007/S11430-009-0082-3.

Ide S, Baltay A, Beroza G C. 2011. Shallow dynamic overshoot and energetic deep rupture in the 2011MW9.0 Tohoku-Oki earthquake.Science, 332(6036): 1426-1429, doi: 10.1126/science.1207020.

Kanamori H. 1977. The energy release in great earthquakes.J.Geophys.Res., 82(20): 2981-2876.

Lay T, Ammon C J, Kanamori H, et al. 2011. Possible large near-trench slip during the 2011Mw9.0 off the Pacific coast of Tohoku earthquake.Earth,PlanetsandSpace, 63(7): 687-692.Li Y G, Vidale J E, Day S M, et al. 2002. Study of the 1999M7.1 Hector Mine, California, earthquake fault plane by trapped waves.Bull.Seismol.Soc.Am., 92(4): 1318-1332.

Mondol N H, Bj?rlykke K, Jahren J, et al. 2007. Experimental mechanical compaction of clay mineral aggregates-Changes in physical properties of mudstones during burial.MarineandPetroleumGeology, 24(5): 289-311.

Rhea S, Tarr A C, Hayes G P, et al. 2010. Seismicity of the Earth 1900-2007, Japan and vicinity. U. S. Geological Survey, Open-File Report 2010-1083-D.

Scholz C H. 1990. The Mechanics of Earthquakes and Faulting. New York: Cambridge University Press.

Simons M, Minson S E, Sladen A, et al. 2011. The 2011 magnitude 9.0 Tohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries.Science, 332(6036): 1421-1425, doi: 10.1126/science.1206731.

Wang K L, Kinoshita M. 2013. Dangers of being thin and weak.Science, 342(6163): 1178-1180, doi: 10.1126/science.1246518. Yamazaki Y, Lay T, Cheung K F, et al. 2011. Modeling near-field tsunami observations to improve finite-fault slip models for the 11 March 2011 Tohoku earthquake.Geophys.Res.Lett., 38(7): L00G15, doi: 10.1029/2011GL049130.

Yin Y Q. 2014. Plastic Mechanics of Rock-Like Materials (in Chinese). Beijing: Peking University Press.

Yin Y Q, Zhang H. 1982. A mathematical model of strain softening in simulating earthquake.ChineseJ.Geophys.(ActaGeophysicaSinica) (in Chinese), 25(5): 414-423. Yoshida K, Miyakoshi K, Irikura K. 2011. Source process of the 2011 off the Pacific coast of Tohoku Earthquake inferred from waveform inversion with long-period strong-motion records.Earth,PlanetsandSpace, 63(7): 577-582.

Yue H, Lay T. 2011. Inversion of high-rate (1 sps) GPS data for rupture process of the 11 March 2011 Tohoku earthquake (MW9.1).Geophys.Res.Lett., 38(7): L00G09, doi: 10.1029/2011GL048700. Zhang Y, Xu L S, Chen Y T. 2012. Rupture process of the 2011 Tohoku earthquake from the joint inversion of teleseismic and GPS data.Earthq.Sci., 25(2): 129-135, doi: 10.1007/s11589-012-0839-1.

附中文參考文獻(xiàn)

胡才博, 周一杰, 蔡永恩. 2009. 如何用有限元新模型研究地震觸發(fā)和應(yīng)力場(chǎng)連續(xù)演化. 中國(guó)科學(xué)D輯: 地球科學(xué), 39(5): 546-555.殷有泉. 2014. 巖石類材料塑性力學(xué). 北京: 北京大學(xué)出版社.

殷有泉, 張宏. 1982. 模擬地震的應(yīng)變軟化的數(shù)學(xué)模型. 地球物理學(xué)報(bào), 25(5): 414-423.

(本文編輯 何燕)

Numerical simulation of the dynamic rupture process of the 2011 Tohoku-Oki, JapanMW9.0 earthquake

LIU Dun-Yu1, HU Cai-Bo2,3, CAI Yong-En1

1InstituteofTheoreticalandAppliedGeophysics,PekingUniversity,Beijing100871,China2KeyLaboratoryofComputationalGeodynamics,ChineseAcademyofSciences,Beijing100049,China3CollegeofEarthScience,UniversityofChineseAcademyofSciences,Beijing100049,China

Earthquakes can be considered as an unstable dynamic process of stress release induced by sudden softening or damage of material on faults. Based on this idea, a finite element method is proposed to simulate the dynamic rupture process of the 2011 Tohoku-Oki, JapanMW9.0 earthquake to find a mechanical explanation for the extremely huge fault slip.According to the results from JFAST (Japan Trench Fast Drilling Project), the ruptured area of this earthquake is modeled as a 5 m-thick layer, which is 250 km along strike and 130 km along dip direction. An asperity of 70 km by 31 km with different material properties is placed at the middle of the ruptured area. A 3-D linear elastic joint element is used to model the fault zone and corresponding equations of the finite element method are derived. The method consists of two parts: (1) The tectonic initial stress field is induced by boundary displacement. (2) Under the initial stress, the rupture process is modeled by gradually reducing the shear modulus of the ruptured area (rupture propagates as expanding circles with a rupture velocity of 3 km·s-1). The amount of reduction of shear modulus in the fault zone is constrained by matching maximum slip from the geodesy observation or inversion results of previous studies.The simulated tectonic stress field shows that shear stress increases with depth. Its minimum value is about 9 MPa at the trench and maximum value is about 15 MPa at depth. The existence of an asperity does not change the initial tectonic stress in the ruptured area remarkably. The displacements, slips and stresses predicted by the models with and without asperities are discussed. If an asperity exists, to match the fault displacements observed at the sea floor, the shear modulus of the asperity and that of non-asperity area have to be reduced to 0.5 MPa and 0.1 MPa, respectively. The maximum quasi-static shear displacement on the hanging wall is 51 m at the trench surface, and that on the footwall is 10 m at the asperity. The maximum fault slips with and without the asperity are 58 m and 55 m, respectively. Both of them appear at the trench surface and their difference is not obvious. At the trench, the hanging wall moves about 50 m horizontally towards the trench and uplifts 12 m, which is consistent with the result of repeated multibeam bathymetric surveys. Shear stress drop at the asperity (about 11 MPa) is larger than that in its surrounding area. Maximum normal stress drop (about 3 MPa) appears at the asperity.Based on the initial tectonic stress simulation, given that the Pacific plate moves towards Japan islands at a rate of 80～85 mm·a-1and that tectonic shear stress is totally released during this earthquake, the recurrence time ofMW9.0 in this area is estimated to be about 2900 years. The predicted results of this paper do not support the inference that the rupture of an asperity is responsible for the unusually large slip of the Tohoku-OkiMW9.0 earthquake, instead the huge damage of fault zone material might be the reason.

Tohoku-Oki earthquake; Fault rupture process; Displacements on fault; Stress drop; Finite element method

劉敦宇, 胡才博, 蔡永恩. 2015. 2011年MW9.0東日本大地震動(dòng)力學(xué)破裂過程的數(shù)值模擬.地球物理學(xué)報(bào)，58(9):3133-3143,

10.6038/cjg20150910.

Liu D Y, Hu C B, Cai Y E. 2015. Numerical simulation of the dynamic rupture process of the 2011 Tohoku-Oki, JapanMW9.0 earthquake.ChineseJ.Geophys. (in Chinese),58(9):3133-3143,doi:10.6038/cjg20150910.

10.6038/cjg20150910

P315

2014-05-30，2015-04-16收修定稿

國(guó)家自然科學(xué)基金(41474080)資助.

劉敦宇，男，1988年生，主要從事震源破裂過程數(shù)值模擬研究.E-mail：dunyuliu@gmail.com