劉敦宇, 胡才博, 蔡永恩

1 北京大學地球與空間科學學院地球物理系， 北京 100871 2 中國科學院計算地球動力學重點實驗室， 北京 100049 3 中國科學院大學地球科學學院， 北京 100049

?

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

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

1 北京大學地球與空間科學學院地球物理系， 北京 100871 2 中國科學院計算地球動力學重點實驗室， 北京 100049 3 中國科學院大學地球科學學院， 北京 100049

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

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

1 引言

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

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

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

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

2 有限元方法

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

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

(1)

(2)

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

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

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

(3)

式中Δu,Δv,Δw分別為節(jié)理元沿斷層傾向、走向和法向的位移間斷量或位錯.

利用斷層坐標系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′方向的位移，上角標“T”表示轉置.本研究選用8節(jié)點六面體節(jié)理單元，每個節(jié)點有3個自由度，節(jié)點編號如圖1所示，1,3,5,7號節(jié)點屬于上盤，對應的下盤節(jié)點編號為2,4,6,8.

圖1 節(jié)理元示意圖(x′，y′，z′)為斷層滑動坐標系. 節(jié)理元沿x′方向長L1，沿y′方向長L2，厚度b，其中L1?b，L2?b.1,3,5,7為上盤面節(jié)點編號，對應的下盤面節(jié)點編號為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é)理元的應力可以由位移間斷量表示為

σ=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é)理元相應方向的剪切模量和楊氏模量.

由有限元方法(殷有泉，2014)，可以得到面積為s，厚度為b的節(jié)理元的單元剛度矩陣為

Ke=b?sBTDBds.

(8)

(9)

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

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

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

(10)

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

3.1 有限元模型及初始構造應力場

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

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

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

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

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

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

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

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

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

3.2 模擬結果及分析

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

數(shù)值試驗表明：對于模型B，若產生約58 m的最大剪切位錯(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.非斷層帶介質阻尼系數(shù)α=0.1，β=0.03.由于斷層帶單元很薄，所以需施加較大阻尼模擬破裂和摩擦過程中的耗散，數(shù)值試驗給出斷層帶單元阻尼系數(shù)α=1000.0，β=0.03.

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

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

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

M0=μAD,

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

得到矩震級MW=9.0.該結果與USGS給出的矩震級一致.

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

圖2 (a)有限元模型; (b)斷層上下盤面應力和位移的正方向示意圖; (c)斷層上盤面初始傾向剪應力; (d)斷層上盤面初始法向正應力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 斷層位錯快照(a)傾向剪切位錯; (b)法向位錯.Fig.3 Snapshots of dislocation(a) Shear dislocation; (b) Normal dislocation.

圖5 震后200 s斷層面上的位移和位錯分布(a)和(b)分別是斷層上盤和下盤剪切位移；(c)和(d)分別是斷層上盤和下盤法向位移；(e)、(f)和(g)分別是斷層上盤、下盤AB測線上的6個測點的剪切位移和剪切位錯隨時間變化.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測線在不同時刻的位移(a) 水平位移; (b) 垂向位移. 橫坐標‘0’和負值分別表示斷層在海溝的出露點和斷層上盤.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測線上6個測點的位移和位錯隨時間變化(位移正方向如圖2b所示).0號測點為震源，深度26km；1號測點深度22km；3、2號測點分別為障礙體的上下邊界所在深度，分別為12km和18km；4號測點深度8.9km；5號測點位于海溝處，深度7km.圖5e為斷層上盤面逆沖位移，海溝處位移值最大，為51m.圖5f為斷層下盤面傾向剪切位移，6個測點中3號(深度為-12km)測點處位移值最大，略小于-10m.圖5g為斷層傾向剪切位錯，海溝處測點位錯最大，為58m.比較六條曲線的斜率可以看出，斷層初始滑動速率(斷層上一點在單位時間內位錯的大小)隨著破裂向海溝傳播逐漸增大：破裂初始階段為0.7m·s-1(震源，0號測點)，逐漸增大到2m·s-1(3號測點)，當破裂達到海溝時，高達6m·s-1(5號測點)，這是由于破裂前端動力積累的應力集中導致的雪崩式的結果.

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

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

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

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

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

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

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

4 結論與討論

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

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

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

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

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

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

圖7 震后200 s斷層上盤面應力變化(a) 剪切應力降；(b) 正應力變化; (c)和(d)分別是AB測線上6個測點剪切應力降和正應力變化隨時間變化曲線.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 有、無障礙體模型的斷層位錯和應力降比較(a)—(d)分別是A模型 (無障礙體)的剪切位錯、法向位錯、斷層面剪應力降和斷層面法向正應力降；(e)—(h)分別是模型B與模型A的剪切位錯之差、法向位錯之差、斷層面剪應力降之差和斷層面法向應力降之差；(i)—(l)分別是模型C與模型A的剪切位錯之差、法向位錯之差、斷層面剪應力降之差和斷層面法向應力降之差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.

附中文參考文獻

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

殷有泉, 張宏. 1982. 模擬地震的應變軟化的數(shù)學模型. 地球物理學報, 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東日本大地震動力學破裂過程的數(shù)值模擬.地球物理學報，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收修定稿

國家自然科學基金(41474080)資助.

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