俯冲带地震循环的数值模拟
——以日本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速度分布具有较好的一致性。模型弹性参数的不确定性对同震以及震间形变的影响有限，模型粘性参数主要影响震间形变场。数值计算也显示，假设震间形变仅由断层运动规律所决定，那么在一个地震周期内，模型空间重力异常基本上随时间均匀变化，在大陆一侧距海沟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)是地面点重力异常值。空间重力异常由地面重力场作空间改正得到，因而对于地面点，高度越高，重力异常越小；高度越低，重力异常越大。从图中可以看出，在-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)尽管表面各点位移和重力异常变化基本均匀，难以观测到震前突变，但在一个地震周期内各点的位移和重力变化可以达到现有仪器的观测水平。空间重力异常和地面点重力异常在一个地震周期内的最大变化可达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)。