翁輝輝　黃金水

1　中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院地震和地球內(nèi)部物理實(shí)驗(yàn)室，合肥市金寨路96號(hào)，230026　2　蒙城地球物理野外科學(xué)觀測(cè)研究站，合肥市金寨路96號(hào)，230026

?

俯沖帶地震循環(huán)的數(shù)值模擬
——以日本Tohoku MW9.0地震為例

翁輝輝1,2黃金水1,2

1中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院地震和地球內(nèi)部物理實(shí)驗(yàn)室，合肥市金寨路96號(hào)，2300262蒙城地球物理野外科學(xué)觀測(cè)研究站，合肥市金寨路96號(hào)，230026

摘要：基于滑移弱化摩擦準(zhǔn)則，以日本2011年Tohoku MW9.0地震為例，建立一個(gè)以物理規(guī)律控制地震循環(huán)過程的二維有限元數(shù)值模型。結(jié)果顯示，參考模型在1 000 a間的6次大地震表現(xiàn)出特征地震的規(guī)律，地震重復(fù)周期約161 a,單位破裂長(zhǎng)度地震矩為1.13×1020Nm/km，在兩次大地震中間會(huì)發(fā)生一次5.62×1018Nm/km的小地震。參考模型的數(shù)值結(jié)果與地表同震GPS位移、震間GPS速度分布具有較好的一致性。模型彈性參數(shù)的不確定性對(duì)同震以及震間形變的影響有限，模型粘性參數(shù)主要影響震間形變場(chǎng)。數(shù)值計(jì)算也顯示，假設(shè)震間形變僅由斷層運(yùn)動(dòng)規(guī)律所決定，那么在一個(gè)地震周期內(nèi)，模型空間重力異?；旧想S時(shí)間均勻變化，在大陸一側(cè)距海溝100 km處可達(dá)-370 μGal；速度場(chǎng)的變化主要發(fā)生在震后約5 a的時(shí)間內(nèi)，此后基本保持穩(wěn)定增加。

關(guān)鍵詞：特征地震；有限元數(shù)值模擬；滑移弱化摩擦準(zhǔn)則；Tohoku地震

歷史地震資料[1]和地震波形分析[2]都顯示，無論大地震還是小微地震都會(huì)在一個(gè)地方有規(guī)律地重復(fù)發(fā)生，即所謂的特征地震特性。地質(zhì)和地球物理研究顯示，Cascadia俯沖帶在過去的6 700 a內(nèi)發(fā)生過12次大地震，平均地震周期為570～590 a[1]；加利福尼亞San Andreas斷層Parkfield區(qū)域平均每22 a就會(huì)發(fā)生一次6級(jí)左右的大地震[3]；在Parkfield和日本東北俯沖區(qū)域，許多小微地震群也具有特征地震的規(guī)律[2]。但是，當(dāng)利用這種特征地震規(guī)律來進(jìn)行地震預(yù)測(cè)時(shí)往往會(huì)遇到困難。一個(gè)最典型的案例是在加利福尼亞San Andreas斷層Parkfield區(qū)域所作的地震預(yù)報(bào)實(shí)驗(yàn)[3]，實(shí)際地震時(shí)間比預(yù)測(cè)時(shí)間晚了10余a[3]，原因可能是該區(qū)應(yīng)力場(chǎng)受附近地震影響[4]，改變了特征地震規(guī)律。

地震循環(huán)的數(shù)值模擬是了解地震規(guī)律的重要手段[5-8]。翁輝輝等[6]的數(shù)值模擬研究顯示，基于摩擦弱化準(zhǔn)則，在剪切應(yīng)力的作用下，斷層的粘滑行為顯示出典型的特征地震規(guī)律，且附近應(yīng)力擾動(dòng)確實(shí)會(huì)對(duì)該規(guī)律產(chǎn)生影響。如果應(yīng)力擾動(dòng)導(dǎo)致斷層面上壓應(yīng)力增加，將會(huì)明顯延遲地震的發(fā)生，并相應(yīng)增加地震釋放的能量；反之，如果應(yīng)力擾動(dòng)導(dǎo)致斷層面上壓應(yīng)力減小，則會(huì)導(dǎo)致地震提前發(fā)生或引發(fā)較小規(guī)模的地震。但其重點(diǎn)是關(guān)注應(yīng)力變化對(duì)地震規(guī)律的影響，并沒有針對(duì)地球進(jìn)行實(shí)際參數(shù)選擇，模型參數(shù)具有較大的不確定性。因此，要利用這種數(shù)值模型來了解實(shí)際地震發(fā)生的規(guī)律，需要結(jié)合地震發(fā)生處的具體地質(zhì)和地球物理環(huán)境來建立和檢驗(yàn)數(shù)值模型。

現(xiàn)代空間大地測(cè)量技術(shù)提供的高精度形變觀測(cè)提高了人們對(duì)地震規(guī)律的認(rèn)識(shí)，如GPS觀測(cè)使科學(xué)家們能夠發(fā)現(xiàn)并研究俯沖斷層帶的慢地震和粘滑現(xiàn)象[9]、震后余滑[10]和下地殼粘滯性[11]，GPS觀測(cè)也幫助科學(xué)家確定地震同震形變[12]和震間形變[11]等重要地震活動(dòng)特征的具體數(shù)據(jù)。

日本本土部署的GPS觀測(cè)站記錄了2011年MW9.0 Tohoku大地震以及震前的大量形變資料[12-13]。受太平洋板塊俯沖和日本海擴(kuò)張影響，日本主要表現(xiàn)為東西向擠壓。這一特征與Tohoku大地震前(1997～2000年間)該區(qū)陸地相對(duì)于穩(wěn)定的歐亞板塊的GPS水平速度[13]一致。震前GPS速度顯示，該區(qū)在震前表現(xiàn)為西向運(yùn)動(dòng)，運(yùn)動(dòng)速度從東向西逐漸減小。Tohoku大地震的同震GPS 位移顯示，地震導(dǎo)致該區(qū)東向運(yùn)動(dòng)，位移值從東向西逐漸減小，最大水平位移約為5.3 m[12]。

除了同震地表位移，日本本土部署的地震觀測(cè)網(wǎng)絡(luò)也獲得高精度的斷層同震滑移[12]。這無疑是一個(gè)重要的約束和檢驗(yàn)數(shù)據(jù)。同時(shí)，日本地球物理學(xué)界對(duì)該區(qū)大量的研究使我們對(duì)其俯沖斷層界面的幾何形態(tài)、殼幔物質(zhì)粘彈性特征等有了充分了解[14]。因而，本文試圖利用該區(qū)的具體地質(zhì)和地球物理特征，以該區(qū)高精度的GPS位移和速度以及同震滑移等為約束，建立一個(gè)地震循環(huán)數(shù)值模型，以期進(jìn)一步了解該區(qū)地震發(fā)生的規(guī)律。

1　數(shù)值模型的建立

本文采用二維模型模擬地震的循環(huán)過程。選擇一個(gè)過地震同震滑移最大值位置，且垂直于海溝的剖面建立模型，見圖1。模型長(zhǎng)1 600 km，深340 km。坐標(biāo)原點(diǎn)設(shè)在表面海溝處，距右邊界600 km，X軸向右為正，Y軸向上為正。俯沖海洋巖石圈上表面按USGS給出的等深線數(shù)據(jù)確定，然后依據(jù)100 km厚度通過樣條插值生成海洋巖石圈下邊界。大陸地殼和地幔巖石圈的厚度分別設(shè)為40 km和60 km。忽略地表界面起伏和自重力作用的影響。

圖1　模型設(shè)置示意圖Fig.1　The 2D model

模型假定大陸地殼是完全彈性體，而其他部分都為麥克斯韋粘彈性體。模型各部分的彈性參數(shù)根據(jù)Miyamachi[14]的地震速度模型給出：大陸地殼、地幔巖石圈以及軟流圈的彈性參數(shù)(楊氏模量)分別取為90 GPa、160 GPa和200 GPa，海洋巖石圈與軟流圈的彈性參數(shù)分別取130 GPa和190 GPa，模型的泊松比都取0.25。大陸地殼、地幔巖石圈以及軟流圈的密度分別取為2 670 kg/m3、3 200 kg/m3和3 300 kg/m3，海洋巖石圈與軟流圈的密度為3 500 kg/m3和3 300 kg/m3。海洋軟流圈的粘性設(shè)為5×1020Pa·s，海洋巖石圈粘性設(shè)為1023Pa·s，大陸地幔巖石圈和軟流圈的粘性值都設(shè)置為1020Pa·s?？紤]到模型參數(shù)選取的不確定性，同時(shí)選取不同的彈性和粘性參數(shù)對(duì)模型進(jìn)行參數(shù)敏感性測(cè)試。測(cè)試中彈性參數(shù)的變化范圍約20%，粘性參數(shù)的變化范圍為1～2個(gè)量級(jí)。大陸地幔巖石圈和軟流圈粘性取值相同，是考慮到大陸的地幔在俯沖帶上方由于脫水和區(qū)域?qū)α鞯挠绊?，粘性可能較小[15]。為進(jìn)一步考察粘性變化可能產(chǎn)生的影響，考慮并計(jì)算了較小的軟流圈粘性(如2.5×1019Pa·s[10])。

為模擬俯沖帶發(fā)震過程，模型將斷層分為具有不同特性的兩段：深度40 km以下的斷層持續(xù)滑動(dòng)，不產(chǎn)生地震，稱為運(yùn)動(dòng)斷層；深度40 km以上的斷層是發(fā)震區(qū)域，發(fā)震和滑移規(guī)律由滑移弱化摩擦準(zhǔn)則控制，稱為動(dòng)力斷層(圖1)。動(dòng)力斷層根據(jù)應(yīng)力變化可處于閉鎖狀態(tài)或產(chǎn)生滑動(dòng)。參考Zhao等[16]的模型，將海洋巖石圈底部邊界面也設(shè)置為運(yùn)動(dòng)斷層(圖1)。動(dòng)力斷層深度設(shè)為40 km，一方面考慮到俯沖帶構(gòu)造地震一般發(fā)生在一定深度范圍上方(如30～50 km[17])，另一方面Tohoku地震的斷層同震滑移[12]也主要發(fā)生在深度40 km以上。其他邊界條件設(shè)置如下：上表面都為自由邊界條件；除海洋巖石圈右側(cè)和下側(cè)外，其他邊界都為自由滑移邊界條件。為避免邊界條件對(duì)俯沖巖石圈的速度分布產(chǎn)生影響，海洋巖石圈的右側(cè)和下側(cè)設(shè)為自由邊界條件(圖1)。

動(dòng)力斷層的滑動(dòng)由滑移弱化摩擦準(zhǔn)則控制，模型中滑移弱化距離d0設(shè)為2 m[6]，模型的靜摩擦系數(shù)和動(dòng)摩擦系數(shù)分別設(shè)為0.68和0.5。模型在動(dòng)力斷層面設(shè)置了-50 MPa的初始正應(yīng)力(負(fù)值代表應(yīng)力為壓應(yīng)力)[18]，以表示周圍巖石圈的環(huán)境應(yīng)力狀態(tài)[6]。計(jì)算中對(duì)斷層滑動(dòng)的處理方法以及庫侖應(yīng)力的定義與翁輝輝等[6]一致。

地震矩定義為：

(1)

式中μ為剪切模量，S為斷層破裂面積，A為斷層平均滑移距離。對(duì)于本文二維模型，斷層破裂面積采用斷層面的破裂長(zhǎng)度來表示，利用二維模型的斷層破裂長(zhǎng)度計(jì)算的地震矩相當(dāng)于實(shí)際地震單位走向破裂長(zhǎng)度的地震矩。在后面的表示中，單位走向破裂長(zhǎng)度都取1 km，因而式(1)的S值是1 km與二維模型斷層破裂長(zhǎng)度的乘積。

模型計(jì)算采用有限元方法。為確保程序能夠正確判斷地震的發(fā)生，動(dòng)力斷層附近的有限元網(wǎng)格必須足夠小。有限元網(wǎng)格采用從邊界向斷層區(qū)域逐漸加密的三角形網(wǎng)格,動(dòng)力斷層處分辨率最大，最小網(wǎng)格邊長(zhǎng)為2 km。計(jì)算時(shí)間步長(zhǎng)設(shè)置為1 a。有限元模型采用PYLITH程序解算。

2　數(shù)值模型結(jié)果

根據(jù)上文介紹的模型及參數(shù)，模擬1 000 a內(nèi)的模型形變以及斷層滑動(dòng)過程。根據(jù)特征地震的假設(shè)，斷層在地震過程中會(huì)釋放掉在震間積累的所有應(yīng)變。在動(dòng)力斷層滑動(dòng)前，由于摩擦力的作用，動(dòng)力斷層的主要部分仍處于閉鎖狀態(tài),閉鎖的動(dòng)力斷層使得其附近區(qū)域的應(yīng)力場(chǎng)發(fā)生變化，并在地震發(fā)生前達(dá)到最大值。從斷層面的應(yīng)力分布來看，運(yùn)動(dòng)斷層的加載使得動(dòng)力斷層深部(沿?cái)鄬用婢啾砻婕s150 ～218 km)的壓應(yīng)力略微減小、中部(距表面約100 ～150 km)的壓應(yīng)力略微增加，而動(dòng)力斷層淺部(距表面約100 km)的正應(yīng)力受影響很小(圖2(a)，橫坐標(biāo)表示從地表起算的斷層面長(zhǎng)度，海溝位置為0 km)。除動(dòng)力斷層與運(yùn)動(dòng)斷層連接處的幾個(gè)節(jié)點(diǎn)正應(yīng)力變化略大外，震前斷層正應(yīng)力的變化幅度總體小于3 MPa(圖2(a))。

圖2　模型動(dòng)力斷層面Fig.2　The dynamic fault of the model

動(dòng)力斷層上的剪切應(yīng)力積累主要集中在深部，剪切應(yīng)力的變化控制了動(dòng)力斷層的滑動(dòng)規(guī)律(圖2(b))。當(dāng)動(dòng)力斷層上的剪切應(yīng)力無法克服摩擦力時(shí)，幾乎所有弱化動(dòng)力斷層都保持閉鎖狀態(tài)。隨著運(yùn)動(dòng)斷層滑移量的增加，弱化斷層部分點(diǎn)的剪切應(yīng)力達(dá)到最大靜摩擦力(圖2(b))，或者說庫侖應(yīng)力達(dá)到臨界值(圖2(c))，這些點(diǎn)將產(chǎn)生滑動(dòng)(圖2(d))。但這些小規(guī)模的滑動(dòng)僅為無震滑移。當(dāng)動(dòng)力斷層應(yīng)力積累時(shí)間達(dá)到80 a，弱化動(dòng)力斷層上的剪切應(yīng)力達(dá)到最大靜摩擦力(或庫侖應(yīng)力達(dá)到臨界值)的區(qū)域增大到一定程度時(shí)，弱化斷層上的部分?jǐn)鄬影l(fā)生地震。地震的發(fā)生使得以前在弱化斷層累積的剪切應(yīng)力得到釋放(圖2(b))，庫侖應(yīng)力減小(圖2(c))。但這次地震僅使得斷層深部一小塊區(qū)域產(chǎn)生破裂，釋放的能量不大，為一個(gè)相對(duì)較小的地震，地震矩約為5.62×1018Nm/km。

隨著運(yùn)動(dòng)斷層的進(jìn)一步加載，動(dòng)力斷層深部剪切應(yīng)力再次增強(qiáng)(圖2(b))，該區(qū)庫侖應(yīng)力也逐漸增大(圖2(c))。當(dāng)弱化動(dòng)力斷層再次達(dá)到小地震破裂條件時(shí)，弱化斷層的破裂就不僅限于前述的部分?jǐn)鄬?，而是?dǎo)致整個(gè)斷層發(fā)生破裂(圖2(d))，即發(fā)生大地震。整個(gè)弱化斷層累積的剪切應(yīng)力得到釋放(圖2(b))、庫侖應(yīng)力減小(圖2(c))，滑動(dòng)釋放的能量大約相當(dāng)于1.13×1020Nm/km，破裂導(dǎo)致地震矩比前述的小地震大約2個(gè)數(shù)量級(jí)。

在1 000 a內(nèi)發(fā)生了6次這樣的完整過程(圖2(d))，大地震的重復(fù)周期約為161 a，小地震的重復(fù)周期約為80 a。從某種意義上講，深部動(dòng)力斷層產(chǎn)生破裂的周期約為80 a，但由于淺部斷層的應(yīng)力加載速度非常慢，約80 a的加載周期不足以導(dǎo)致整個(gè)動(dòng)力斷層產(chǎn)生破裂。第一次小地震將深部斷層的應(yīng)力傳導(dǎo)到上部，使得第二次地震能夠破裂整個(gè)動(dòng)力斷層，從而產(chǎn)生大地震。因此，該模型的一個(gè)完整周期大約是161 a，包含一次小地震和一次大地震。

圖3給出了模型結(jié)果與空間大地測(cè)量結(jié)果等的對(duì)比?？梢钥闯?，水平位移GPS觀測(cè)與模型符合得很好(圖3(a))；垂直位移盡管存在一定差異，但總體仍基本一致(圖3(b))，不過GPS垂直位移精度遠(yuǎn)不如水平位移高[12]。圖3(c)給出了模型斷層滑移量和利用GPS數(shù)據(jù)反演得到的同震斷層滑移分布[12]。模型大地震的平均同震滑移量約為14 m，最大滑移量可以達(dá)到16 m，盡管從曲線形態(tài)上看兩者存在一定差異，但總的滑移量大致相當(dāng)。

圖3　模型結(jié)果與空間大地測(cè)量結(jié)果的對(duì)比Fig.3　Comparison with the model result and GPS result

圖3(d)是模型大地震周期內(nèi)水平速度場(chǎng)以及1997～2000年間GPS觀測(cè)的平均水平速度分量?？梢钥闯?，模型的震間速度場(chǎng)變化可以達(dá)到13 mm/a，反映了模型粘性效應(yīng)的影響。圖3(d)亦顯示，模型大地震發(fā)生前15 a(或者說一個(gè)地震周期的最后15 a)GPS速度場(chǎng)變化很小，GPS觀測(cè)值是震前11 ～14 a的觀測(cè)值，兩者數(shù)據(jù)具有一定程度的一致性。

計(jì)算結(jié)果顯示，彈性參數(shù)對(duì)地表形變的影響包含兩個(gè)方面：地表同震位移和震間水平速度變化。彈性參數(shù)對(duì)同震位移的影響主要集中在近海溝處200 km范圍內(nèi)，且大陸地殼和海洋巖石圈參數(shù)變化的影響相對(duì)較大。大陸地殼或海洋巖石圈的楊氏模量20%的變化會(huì)造成海溝附近最大約15%的水平位移和垂直位移變化，其他塊體20%的彈性參數(shù)變化造成的地表同震位移變化最大小于1%。在遠(yuǎn)離海溝處，楊氏模量20%的變化造成的同震位移變化都小于2%。彈性參數(shù)對(duì)震間水平速度的影響范圍較廣，對(duì)震前15 a的水平速度場(chǎng)影響最小的區(qū)域位于孕震(動(dòng)力)斷層上方，即距海溝約200 km的陸地上；而影響最大的區(qū)域位于該區(qū)兩側(cè)，即距海溝約100 km和約300 km的陸地區(qū)域。和彈性參數(shù)變化對(duì)同震位移的影響一樣，對(duì)震前15 a的水平速度場(chǎng)，大陸地殼或海洋巖石圈參數(shù)變化導(dǎo)致的震間速度場(chǎng)的變化也最大。大陸地殼或海洋巖石圈楊氏模量20%的變化造成的最大震間水平速度變化約為5 mm/a，而其他塊體20%的彈性參數(shù)變化造成的震間水平速度變化最大值約1～2 mm/a。對(duì)震后20 a的水平速度場(chǎng)，彈性參數(shù)變化影響最大的區(qū)域位于近海溝處和模型左側(cè)邊界附近，而且影響最大的是大陸軟流圈和海洋巖石圈的彈性參數(shù)。20%的彈性參數(shù)變化造成的震間水平速度變化最大值也大約5 mm/a。

粘性參數(shù)的變化對(duì)同震位移的影響區(qū)域主要集中在海溝附近，粘性2.5倍或1個(gè)量級(jí)的變化導(dǎo)致同震位移變化約3%，且主要變化在距離海溝100 km的范圍內(nèi)。粘性參數(shù)的變化對(duì)震間形變的影響比較顯著，除了海洋巖石圈粘性變化影響較小外，其他塊體1個(gè)量級(jí)的粘性變化最大可以導(dǎo)致超過10 mm/a的速度場(chǎng)變化。變化最大的區(qū)域主要位于大陸一側(cè)，這可能和模型將大陸地殼設(shè)置為完全彈性有關(guān)。圖4給出了其他模型參數(shù)不變，大陸軟流圈粘性取2.5×1019Pa·s，大陸地幔巖石圈粘性分別取2.5×1019Pa ·s、2.5×1020Pa·s和2.5×1021Pa·s時(shí)3個(gè)模型的GPS水平速度場(chǎng)。當(dāng)大陸巖石圈和軟流圈粘性都只有2.5×1019Pa·s，其震間水平速度場(chǎng)的變化可達(dá)38 mm/a，遠(yuǎn)遠(yuǎn)高于大陸巖石圈和軟流圈粘性都為1020Pa·s的模型。隨著大陸巖石圈粘性的增強(qiáng)，震間水平速度場(chǎng)的變化減小。在大陸地幔巖石圈粘性值為2.5×1020Pa·s時(shí)，其震間水平速度場(chǎng)在震前1～15 a的模型值也與觀測(cè)數(shù)據(jù)吻合(圖4(b))。震間速度場(chǎng)變化與圖3(d)所示模型存在差別，這表明震間地表GPS速度取決于巖石圈和軟流圈粘性的組合。

圖4　不同巖石圈粘性模型的震間地表水平速度Fig.4　The interseismic surface horizontal velocity of the models with different values of viscosity of the continental mantle lithosphere

此外，彈性參數(shù)和粘性參數(shù)對(duì)模型結(jié)果影響的測(cè)試還表明，在上述模型參數(shù)變化范圍內(nèi)，參數(shù)變化對(duì)地震周期的影響小于5%，對(duì)地震矩的影響小于2%。應(yīng)該指出的是，模型的地震復(fù)發(fā)周期和地震矩主要受控于模型初始正應(yīng)力、摩擦系數(shù)和動(dòng)力斷層的深度等，增加模型的初始正應(yīng)力、動(dòng)靜摩擦系數(shù)之差(μs-μd)或動(dòng)力斷層的長(zhǎng)度都會(huì)相應(yīng)地增加模型的地震周期以及地震矩。動(dòng)力斷層的深度不同，會(huì)得到不同的同震滑移分布(圖3(c))。此外，模型的同震位移和震間形變也是確定模型參數(shù)的重要約束。應(yīng)該說，模型中的許多參數(shù)是參考實(shí)際地球的結(jié)果，盡管有些參數(shù)的變化可能與參數(shù)組合有關(guān)，如大陸巖石圈和軟流圈的粘性，但計(jì)算結(jié)果顯示，按照實(shí)際地球選取參數(shù)的模型與現(xiàn)有的同震位移、震間速度場(chǎng)等觀測(cè)基本一致。

3　討　論

本文利用一個(gè)數(shù)值模型模擬了日本Tohoku附近區(qū)域的地震循環(huán)過程，該模型與GPS位移和速度場(chǎng)數(shù)據(jù)符合較好。本文模型是二維的，僅能給出單位走向破裂長(zhǎng)度的地震矩。如果假定斷層面上的滑移分布可用于整個(gè)Tohoku 地震400 km破裂面[12]，那么根據(jù)本文模型得到的等效地震矩約為4.52×1022Nm，相當(dāng)于一個(gè)MW=9.1的地震，與觀測(cè)基本相符[12]。如果本文模型正確，就意味著Tohoku大地震的頻率約為每161 a一次。

本文模型表現(xiàn)出的另一個(gè)現(xiàn)象是，在兩次大地震之間會(huì)發(fā)生一次小地震。假設(shè)小地震的走向破裂長(zhǎng)度可達(dá)90 km (以7級(jí)左右地震的走向破裂長(zhǎng)度來估算)，由此估算小地震的等效地震矩約為5.06×1020Nm，相當(dāng)于一個(gè)MW=7.8的地震。值得注意的是，在Tohoku地區(qū)，在2011年大地震之前確實(shí)發(fā)生過幾次7級(jí)左右的大地震[19]。盡管模型只有一次略小的地震，但觀測(cè)卻不止一個(gè)。

為了探討可能的震前地表觀測(cè)信息，圖5給出了模型在不同時(shí)間的重力異常和垂直位移的空間分布，以及幾個(gè)不同位置點(diǎn)的重力異常和位移的時(shí)間變化。從圖5(a)和5(c)可以看到，空間重力異常和垂直位移正相關(guān)。進(jìn)一步的分析表明，空間重力異常主要來源于模型地表形變?cè)斐傻牡匦纹鸱Ｔ谝粋€(gè)地震周期內(nèi)，海洋巖石圈俯沖導(dǎo)致海溝處附近(-100 km處)地形降低超過3 m、陸地一側(cè)距海溝約200 km處(-200 km)抬升1 m多；與此對(duì)應(yīng)，空間重力異常在-100 km處減小300多μGal、在 -200 km處增加200多μGal(圖5(a)和5(c))。當(dāng)?shù)卣鸢l(fā)生時(shí)，空間重力異常變化與此相反，在 -100 km和 -200 km處空間重力異常分別增加和減小約300μGal(圖5(e))。

圖5　模型在不同時(shí)間的重力異常和垂直位移的空間分布以及不同位置點(diǎn)重力異常和位移的時(shí)間變化Fig.5　The free-air gravity anomaly distribution of the model at different times and the free-air gravity anomaly versus interseismic time on different positions

為探討重力異常的時(shí)間變化規(guī)律，計(jì)算在空間上重力場(chǎng)變化具有代表性的幾個(gè)特征點(diǎn)的重力異常時(shí)間變化曲線(圖5(b)和5(f))。這些點(diǎn)分別位于大陸一側(cè)距海溝900 km、500 km、200 km和100 km以及海洋一側(cè)距海溝400 km處。圖5(d)、5(g)和5(h)分別給出了這些特征點(diǎn)的位移和速度場(chǎng)的時(shí)間變化，圖5(b)是空間重力異常值，圖5(f)是地面點(diǎn)重力異常值。空間重力異常由地面重力場(chǎng)作空間改正得到，因而對(duì)于地面點(diǎn)，高度越高，重力異常越小；高度越低，重力異常越大。從圖中可以看出，在-100 km處，空間重力異常隨時(shí)間減小最多(圖5(b))，但地面點(diǎn)重力異常隨時(shí)間增加最多(圖5(f))，這主要是地面點(diǎn)高程隨時(shí)間減小的緣故(圖5(d))。在-200 km的情況和 -100 km處正好相反，因?yàn)榇颂幐叱屉S時(shí)間增加，空間重力異常和地面點(diǎn)重力異常分別隨時(shí)間增加和減小，-200 km處高程變化小于-100 km處，-200 km處重力異常變化小于-100 km處。其他點(diǎn)變化類似，只是變化值更小。值得注意的是，在-500 km處，空間重力異常隨時(shí)間略微增加，但由于受地形略微降低的影響，地面點(diǎn)重力異常也略微增加。

從圖3(b)可以看到，模型大地震的同震垂直位移可達(dá)2～3 m，空間重力異常的變化大約為300μGal(圖5(e))，小地震的重力變化也可達(dá)約100μGal。但由于圖5(b)、5(d)和5(f)的采樣點(diǎn)不在變化最大的地方，圖中并沒有明顯看到小地震造成的變化。對(duì)同震重力變化，除了靠近海溝處造成的重力增加和減小外，在遠(yuǎn)離海溝的大陸和海洋一側(cè)都有一個(gè)微小的負(fù)空間重力異常(圖5(e))，變化約為11μGal。本文的重力變化與Matsuo等[20]利用彈性形變模型得到的結(jié)果略微不同，其僅在日本Tohoku大地震區(qū)域大陸一側(cè)發(fā)現(xiàn)約-7μGal的空間重力變化。這種不同可能是由于本文的模型斷層滑移基本均勻，且能到達(dá)海溝，而Matsuo等[20]采用的地震同震滑移模型中心極大造成的。

圖5(g)給出幾個(gè)特征點(diǎn)的水平位移在一個(gè)地震周期內(nèi)的時(shí)間變化?？梢钥闯?，水平位移變化大于垂直位移，同時(shí)由于斷層基本處于閉鎖狀態(tài)，海洋巖石圈表面(400 km處)和近海溝處(-100 km處)的位移變化基本一致，越往內(nèi)陸水平位移越小(圖5(g))。圖5(h)顯示了幾個(gè)特征點(diǎn)的水平速度場(chǎng)在一個(gè)地震周期內(nèi)的時(shí)間變化?？梢钥闯?，海洋巖石圈和大陸深處(如400 km和-900 km處)水平速度變化最小，-200 km處速度變化最大，一個(gè)周期內(nèi)最大速度變化可達(dá)約10 mm/a(圖5(h))。除了中間小地震造成速度躍變外，模型的初始階段在-200 km和-500 km處的速度變化也較大，這主要是模型粘彈性效應(yīng)的影響[11]，在一個(gè)地震周期內(nèi)水平速度的變化基本穩(wěn)定(圖5(h))。

從圖5可以看到，在一個(gè)地震周期內(nèi)，重力和位移等表面觀測(cè)都能被儀器檢測(cè)到，但震前都難以觀測(cè)到各量的突然變化。然而，空間重力異常和地面點(diǎn)重力異常在一個(gè)地震周期內(nèi)的最大變化可以分別達(dá)到378μGal和655μGal，垂直位移和水平位移最大變化可以分別達(dá)到3.3 m和16 m，水平速度最大變化可達(dá)13 mm/a。這也許意味著震前觀測(cè)信號(hào)的檢測(cè)應(yīng)該從關(guān)注震前突變轉(zhuǎn)到觀測(cè)量的總體變化上來。

4　結(jié)　語

1)日本Tohoku區(qū)域大地震的重復(fù)周期約為161 a。單位破裂長(zhǎng)度地震矩約為1.13×1020Nm/km，400 km走向破裂長(zhǎng)度的等效地震矩約為4.52×1022Nm，相當(dāng)于矩震級(jí)為9的地震。每?jī)纱未蟮卣痖g，由于斷層區(qū)域應(yīng)力積累不足，斷層會(huì)產(chǎn)生一個(gè)范圍有限的破裂，形成一個(gè)大小約為5.62×1018Nm/km的小地震。

2)模型結(jié)果與地表同震GPS位移、震間GPS速度分布具有較好的一致性，說明本文2D有限元模型能很好地反映日本2011年Tohoku MW9.0大地震以及一個(gè)地震周期內(nèi)的形變特征。模型參數(shù)測(cè)試表明，彈性參數(shù)對(duì)同震以及震間形變的影響有限；粘性參數(shù)對(duì)同震形變的影響較小，但對(duì)震間形變影響較大。模型的震間形變?nèi)Q于各塊體的粘性組合，如大陸巖石圈粘性減小的影響可由軟流圈粘性增加的效應(yīng)補(bǔ)償。

3)盡管表面各點(diǎn)位移和重力異常變化基本均勻，難以觀測(cè)到震前突變，但在一個(gè)地震周期內(nèi)各點(diǎn)的位移和重力變化可以達(dá)到現(xiàn)有儀器的觀測(cè)水平。空間重力異常和地面點(diǎn)重力異常在一個(gè)地震周期內(nèi)的最大變化可達(dá)378μGal和650μGal，垂直位移和水平位移最大變化可達(dá)3.3 m和16 m，水平速度的最大變化可達(dá)13 mm/a。

致謝：感謝美國地質(zhì)調(diào)查所(USGS)提供俯沖板塊接觸面數(shù)據(jù)。感謝CIG研究小組以及Aagaard等在軟件方面提供支持。工作中與楊宏峰博士進(jìn)行過多次有益討論，在此表示感謝。

參考文獻(xiàn)

[1]Witter R C, Kelsey H M, Hemphill-Haley E. Great Cascadia Earthquakes and Tsunamis of the Past 6 700 Years, Coquille River Estuary, Southern Coastal Oregon [J]. Geological Society of America Bulletin, 2003, 115(10):1 289-1 306

[2]Igarashi T, Matsuzawa T, Hasegawa A. Repeating Earthquakes and Interplate Aseismic Slip in the Northeastern Japan Subduction Zone [J]. Journal of Geophysical Research, 2003, 108(B5)

[3]Bakun W H, Lindh A G. The Parkfield, California, Earthquake Prediction Experiment[J]. Science, 1985, 229(4 714): 619-624

[4]Ben-Zion Y, Rice J R, Dmowska R. Interaction of the San Andreas Fault Creeping Segment with Adjacent Great Rupture Zones and Earthquake Recurrence at Parkfield [J]. Journal of Geophysical Research, 1993, 98:2 135-2 144

[5]蔡永恩，何濤，王仁. 1976 年唐山地震震源動(dòng)力過程的數(shù)值模擬[J]. 地震學(xué)報(bào)，1999,21(5): 469-477(Cai Yongen, He Tao, Wang Ren. Modelling the Source Dynamic Process of the 1976 Tangshan Earthquake [J]. Acta Seismologica Sinica, 1999, 21(5): 469-477)

[6]翁輝輝,黃金水. 應(yīng)力擾動(dòng)對(duì)地震周期和地震矩影響的數(shù)值模擬[J]. 地震學(xué)報(bào),2015,37(1): 66-80(Weng Huihui, Huang Jinshui. Numerical Simulations about the Influence of Stress Disturbance on Earthquake Cycle [J]. Acta Seismologica Sinica,2015,37(1): 66-80)

[7]朱守彪,邢會(huì)林,謝富仁, 等. 地震發(fā)生過程的有限單元法模擬———以蘇門答臘俯沖帶上的大地震為例[J]. 地球物理學(xué)報(bào), 2008,51(2): 460-468(Zhu Shoubiao, Xing Huilin, Xie Furen, et al. Simulation of Earthquake Processes by Finite Element Method: the Case of Megathrust Earthquakes on the Sumatra Subduction Zone [J]. Chinese Journal of Geophysics, 2008, 51(2): 460-468)

[8]朱守彪,張培震. 2008 年汶川MS8.0 地震發(fā)生過程的動(dòng)力學(xué)機(jī)制研究[J]. 地球物理學(xué)報(bào), 2009,52(2): 418-427(Zhu Shoubiao, Zhang Peizhen. A Study on the Dynamical Mechanisms of the Wenchuan MS8.0 Earthquake [J]. Chinese Journal of Geophysics, 2009, 52(2): 418-427)

[9]Miller M M, Melbourne T, Johnson D J, et al. Periodic Slow Earthquakes from the Cascadia Subduction Zone [J]. Science, 2002, 295(5 564): 2 423-2 423

[10]Diao F Q, Xiong X, Wang R J, et al. Overlapping Post-seismic Deformation Processes: Afterslip and Viscoelastic Relaxation Following the 2011 Mw 9.0 Tohoku (Japan) Earthquake[J]. Geophysical Journal International, 2014, 196(1):218-229

[11]Wang K L, Hu Y, He J H. Deformation Cycles of Subduction Earthquakes in a Viscoelastic Earth [J]. Nature, 2012, 484(7 394): 327-332

[12]Ozawa S, Nishimura T, Suito H, et al. Coseismic and Postseismic Slip of the 2011 Magnitude-9 Tohoku-Oki Earthquake [J]. Nature, 2011, 475(7 356): 373-376

[13]Loveless J P, Meade B J. Geodetic Imaging of Plate Motions, Slip Rates, and Partitioning of Deformation in Japan [J]. Journal of Geophysical Research, 2010, 115(B2)

[14]Miyamachi H. Seismic Velocity Structure in the Crust and Upper Mantle Beneath Northern Japan [J]. Journal of Physics of the Earth, 1994, 42(2): 269-301

[15]Hasegawa A, Nakajima J, Umino N, et al. Deep Structure of the Northeastern Japan Arc and Its Implications for Crustal Deformation and Shallow Seismic Activity [J]. Tectonophysics, 2005, 403(1): 59-75

[16]Zhao S R, Takemoto S. Deformation and Stress Change Associated with Plate Interaction at Subduction Zones: a Kinematic Modelling [J]. Geophysical Journal International, 2000, 142(2): 300-318

[17]Hu Y, Wang K, He J, et al. Three-dimensional Viscoelastic Finite Element Model for Postseismic Deformation of the Great 1960 Chile Earthquake [J]. Journal of Geophysical Research, 2004, 109(B12)

[18]Rice J R. Fault Stress States, Pore Pressure Distributions, and the Weakness of the San Andreas Fault [J]. International Geophysics, 1992, 51:475-503

[19]Yamanaka Y, Kikuchi M. Asperity Map along the Subduction Zone in Northeastern Japan Inferred from Regional Seismic Data [J]. Journal of Geophysical Research, 2004, 109(B7)

[20]Matsuo K, Heki K. Coseismic Gravity Changes of the 2011 Tohoku‐Oki Earthquake from Satellite Gravimetry [J]. Geophysical Research Letters, 2011, 38(7)

Foundation support:National Natural Science Foundation of China, No.41474082, 91014005,40774045.

About the first author:WENG Huihui, postdoctor, majors in dynamic rupture simulations,E-mail:qfkq7850@mail.ustc.edu.cn.Corresponding author:HUANG Jinshui, professor, PhD supervisor,majors in geodynamic numerical simulations,E-mail:jshhuang@ustc.edu.cn.

收稿日期：2015-08-03

第一作者簡(jiǎn)介：翁輝輝，博士后，主要從事地震動(dòng)態(tài)破裂過程研究，E-mail：qfkq7850@mail.ustc.edu.cn。 通訊作者：黃金水，教授，博士生導(dǎo)師，主要從事地球動(dòng)力學(xué)數(shù)值模擬研究，E-mail：jshhuang@ustc.edu.cn。

DOI：10.14075/j.jgg.2016.08.001

文章編號(hào)：1671-5942(2016)08-0659-07

中圖分類號(hào)：P315

文獻(xiàn)標(biāo)識(shí)碼：A

Numerical Simulations about Subduction Earthquake Cycles:The Case of Japan Tohoku MW9.0 Earthquake

WENGHuihui1,2HUANGJinshui1,2

1Laboratory for Earthquake and Earth’s Interior, School of Earth and Space Sciences, University of Science and Technology of China, 96 Jinzhai Road, Hefei 230026, China2Mengcheng National Geophysical Observatory, 96 Jinzhai Road, Hefei 230026, China

Abstract:Earthquake cycles based on the site of the 2011 MW9.0 earthquake in Tohoku, Japan was simulated with a 2D finite element model in which motions of the fault is controlled by a slip weakening friction law. The numerical results show that earthquakes occur in a typical pattern of characteristic earthquakes: there are 6 large earthquakes during the 1 000 a simulation, the intervals between two adjacent earthquakes are about 161±4 a, and the seismic moments of per unit length for each earthquake is about 1.13×1020Nm/km. There is a small earthquake with a seismic moment of 5.62×1018Nm/km occurred between each two-adjoin large earthquakes. The coseismic and interseismic surface deformations of the numerical model agree well with the GPS observations. The uncertainty of the elastic parameters has limited effects on the coseismic and interseismic deformations, while variations of the viscosity influence the interseismic deformations. The numerical results also show that if the interseismic deformations were only controlled by the motion of the fault, the gravity anomaly of this model would decrease linearly during the interseismic period and reach about -370 μGal at about 100 km from the trench on the continent. Variations on velocity mainly happened within the first 5 a and velocities change little about 5 a after each earthquake.

Key words:characteristic earthquake; finite element numerical simulation; slip-weakening friction law; the 2011 Tohoku earthquake

項(xiàng)目來源：國家自然科學(xué)基金(41474082，91014005和40774045)。