俯沖帶地震循環的數值模擬
——以日本Tohoku MW9.0地震為例

翁輝輝　黃金水

1　中國科學技術大學地球和空間科學學院地震和地球內部物理實驗室，合肥市金寨路96號，230026　2　蒙城地球物理野外科學觀測研究站，合肥市金寨路96號，230026

?

俯沖帶地震循環的數值模擬
——以日本Tohoku MW9.0地震為例

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

1中國科學技術大學地球和空間科學學院地震和地球內部物理實驗室，合肥市金寨路96號，2300262蒙城地球物理野外科學觀測研究站，合肥市金寨路96號，230026

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

關鍵詞：特征地震；有限元數值模擬；滑移弱化摩擦準則；Tohoku地震

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

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

現代空間大地測量技術提供的高精度形變觀測提高了人們對地震規律的認識，如GPS觀測使科學家們能夠發現并研究俯沖斷層帶的慢地震和粘滑現象[9]、震后余滑[10]和下地殼粘滯性[11]，GPS觀測也幫助科學家確定地震同震形變[12]和震間形變[11]等重要地震活動特征的具體數據。

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

除了同震地表位移，日本本土部署的地震觀測網絡也獲得高精度的斷層同震滑移[12]。這無疑是一個重要的約束和檢驗數據。同時，日本地球物理學界對該區大量的研究使我們對其俯沖斷層界面的幾何形態、殼幔物質粘彈性特征等有了充分了解[14]。因而，本文試圖利用該區的具體地質和地球物理特征，以該區高精度的GPS位移和速度以及同震滑移等為約束，建立一個地震循環數值模型，以期進一步了解該區地震發生的規律。

1　數值模型的建立

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

圖1　模型設置示意圖Fig.1　The 2D model

模型假定大陸地殼是完全彈性體，而其他部分都為麥克斯韋粘彈性體。模型各部分的彈性參數根據Miyamachi[14]的地震速度模型給出：大陸地殼、地幔巖石圈以及軟流圈的彈性參數(楊氏模量)分別取為90 GPa、160 GPa和200 GPa，海洋巖石圈與軟流圈的彈性參數分別取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。海洋軟流圈的粘性設為5×1020Pa·s，海洋巖石圈粘性設為1023Pa·s，大陸地幔巖石圈和軟流圈的粘性值都設置為1020Pa·s?？紤]到模型參數選取的不確定性，同時選取不同的彈性和粘性參數對模型進行參數敏感性測試。測試中彈性參數的變化范圍約20%，粘性參數的變化范圍為1～2個量級。大陸地幔巖石圈和軟流圈粘性取值相同，是考慮到大陸的地幔在俯沖帶上方由于脫水和區域對流的影響，粘性可能較小[15]。為進一步考察粘性變化可能產生的影響，考慮并計算了較小的軟流圈粘性(如2.5×1019Pa·s[10])。

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

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

地震矩定義為：

(1)

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

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

2　數值模型結果

根據上文介紹的模型及參數，模擬1 000 a內的模型形變以及斷層滑動過程。根據特征地震的假設，斷層在地震過程中會釋放掉在震間積累的所有應變。在動力斷層滑動前，由于摩擦力的作用，動力斷層的主要部分仍處于閉鎖狀態,閉鎖的動力斷層使得其附近區域的應力場發生變化，并在地震發生前達到最大值。從斷層面的應力分布來看，運動斷層的加載使得動力斷層深部(沿斷層面距表面約150 ～218 km)的壓應力略微減小、中部(距表面約100 ～150 km)的壓應力略微增加，而動力斷層淺部(距表面約100 km)的正應力受影響很小(圖2(a)，橫坐標表示從地表起算的斷層面長度，海溝位置為0 km)。除動力斷層與運動斷層連接處的幾個節點正應力變化略大外，震前斷層正應力的變化幅度總體小于3 MPa(圖2(a))。

圖2　模型動力斷層面Fig.2　The dynamic fault of the model

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

隨著運動斷層的進一步加載，動力斷層深部剪切應力再次增強(圖2(b))，該區庫侖應力也逐漸增大(圖2(c))。當弱化動力斷層再次達到小地震破裂條件時，弱化斷層的破裂就不僅限于前述的部分斷層，而是導致整個斷層發生破裂(圖2(d))，即發生大地震。整個弱化斷層累積的剪切應力得到釋放(圖2(b))、庫侖應力減小(圖2(c))，滑動釋放的能量大約相當于1.13×1020Nm/km，破裂導致地震矩比前述的小地震大約2個數量級。

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

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

圖3　模型結果與空間大地測量結果的對比Fig.3　Comparison with the model result and GPS result

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

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

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

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

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

3　討　論

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

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

為了探討可能的震前地表觀測信息，圖5給出了模型在不同時間的重力異常和垂直位移的空間分布，以及幾個不同位置點的重力異常和位移的時間變化。從圖5(a)和5(c)可以看到，空間重力異常和垂直位移正相關。進一步的分析表明，空間重力異常主要來源于模型地表形變造成的地形起伏。在一個地震周期內，海洋巖石圈俯沖導致海溝處附近(-100 km處)地形降低超過3 m、陸地一側距海溝約200 km處(-200 km)抬升1 m多；與此對應，空間重力異常在-100 km處減小300多μGal、在 -200 km處增加200多μGal(圖5(a)和5(c))。當地震發生時，空間重力異常變化與此相反，在 -100 km和 -200 km處空間重力異常分別增加和減小約300μGal(圖5(e))。

圖5　模型在不同時間的重力異常和垂直位移的空間分布以及不同位置點重力異常和位移的時間變化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

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

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

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

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

4　結　語

1)日本Tohoku區域大地震的重復周期約為161 a。單位破裂長度地震矩約為1.13×1020Nm/km，400 km走向破裂長度的等效地震矩約為4.52×1022Nm，相當于矩震級為9的地震。每兩次大地震間，由于斷層區域應力積累不足，斷層會產生一個范圍有限的破裂，形成一個大小約為5.62×1018Nm/km的小地震。

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

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

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

參考文獻

[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 年唐山地震震源動力過程的數值模擬[J]. 地震學報，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]翁輝輝,黃金水. 應力擾動對地震周期和地震矩影響的數值模擬[J]. 地震學報,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]朱守彪,邢會林,謝富仁, 等. 地震發生過程的有限單元法模擬———以蘇門答臘俯沖帶上的大地震為例[J]. 地球物理學報, 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 地震發生過程的動力學機制研究[J]. 地球物理學報, 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

第一作者簡介：翁輝輝，博士后，主要從事地震動態破裂過程研究，E-mail：qfkq7850@mail.ustc.edu.cn。 通訊作者：黃金水，教授，博士生導師，主要從事地球動力學數值模擬研究，E-mail：jshhuang@ustc.edu.cn。

DOI：10.14075/j.jgg.2016.08.001

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

中圖分類號：P315

文獻標識碼：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

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