地球物理学报  2013, Vol. 56 Issue (7): 2303-2312   PDF    
日本俯冲带地震发生过程的数值模拟研究
姜辉1 , 邓志辉2 , 高祥林1 , 陈梅花3 , 杨竹转1     
1. 中国地震局地质研究所, 北京 100029;
2. 中国地震局地震监测与减灾技术重点实验室, 广州 510070;
3. 浙江师范大学, 浙江 金华 321004
摘要: 了解地震发生的动力学机制是研究地震发震原因的关键, 而数值模拟的方法是高速、有效的手段.2011年3月11日日本东北部宫城县发生9.0级大地震, 文中以该次大地震所在的日本俯冲带为研究对象, 通过使用黏弹性有限元数值模拟, 并引用接触对, 建立了研究区二维数值模型, 模拟俯冲带与上覆板片之间的滑动、黏滞到再滑动的过程, 亦即断层失稳发生地震的过程.模拟结果显示, 随着太平洋板块不断俯冲, 在俯冲带上自发出现了断层闭锁、解锁到再闭锁的黏滑过程, 且这种过程呈现一定的准周期性, 大事件主要集中分布在20~30 km的深度范围内.根据俯冲带可能在俯冲过程中角度的变化, 建立了不同的模型, 进行模拟对比研究, 结果表明, 俯冲带的几何形态, 以及俯冲角度变化所在的不同深度, 对模拟的结果有不同的影响.
关键词: 日本俯冲带      数值模拟      黏滑      地震     
Numerical modeling of earthquake generating processes in the Japan subduction slab
JIANG Hui1, DENG Zhi-Hui2, GAO Xiang-Lin1, CHEN Mei-Hua3, YANG Zhu-Zhuan1     
1. Institute of Geology, China Earthquake Administration, Beijing 100029, China;
2. Key Laboratory of Earthquake Monitoring and Hazard Reduction Technology, China Earthquake Administration, Guangzhou 510070, China;
3. Zhejiang Normal University, Jinhua Zhejiang 321004, China
Abstract: Understanding of earthquake dynamic mechanism is the key of studying the causes of earthquakes, especially the large events. Numerical simulation is a fast and efficient approach to solve this issue. The Mw9.0 earthquake of Miyagi Prefecture, northeastern Japan on 11th March, 2011 occurred at the Japan subduction slab. In this article, we take the Japan subduction slab as the study object, and simulate the stick-slip process between the slab and overlying plate as well as the instability processes of the fault by a two-dimensional viscoelastic finite element model. The results show that the underthrusting of the Pacific plate leads to the spontaneous stick-slip processes, characterized by fault locking, fault unlocking and back to fault locking again, which present certain quasi-periodicity in the slab, and big events mainly concentrated in the depth range of 20~30 km. Using different models with varied dip angles of the subduction slab, we make a comparative study of simulations. The results reveal that the geometry of the subduction slab and different depths where the subduction angle changes have different effects on the simulation results..
Key words: Japan subduction slab      Numerical simulation      Stick-slip      Earthquakes     
1 引言

根据中国国家地震台网测定,北京时间2011年3月11日13时46分,日本东北地区宫城县北部发生里氏9.0级大地震.USGS的观测数据表明,这次大震震中位于北纬38.1°,东经142.6°,震源深度约24.4km.此次大震引发了巨大的海啸,海啸波于震后15min抵达日本沿岸[1],最高达40.5m.截至当地时间4月12日19时,此次地震及其引发的海啸已造成14063人死亡,13691人失踪.

日本大地震发生之后,引起了世界各地地震学家的关注.Ozawa等[2]分析了该次大地震发生后产生的同震及震后滑移.Kerr[3]分析了大震发生的可能原因,认为震前的慢滑移对大地震的发生有一定的促进作用.Kato[4]通过分析大量9.0级大地震的前震数据,发现在大地震发生前,有两个明显的前震序列以2~10km/d的速率向大震震中移动.经过论证认为,第一个序列的慢地震(慢滑移)使得震中处的板块交界处变弱,而第二序列慢地震造成9.0级大地震震中位置应力大量加载,这些慢地震产生的效应可能是导致破裂产生的最初动力.徐彦[5]分析了9.0级大地震的震源特征.随着有限元数值模拟方法的成熟运用,大量专家学者对此次大震进行了模拟分析[6-7].

自1973年以来,日本海沟的俯冲带已经发生过9次7级以上的地震事件.大多数近海地震都发生在日本同一个俯冲带上.其中1611年、1896年和1933年的近海地震都在日本东北的太平洋沿岸引起毁灭性的海啸.日本俯冲带地震发生十分频繁,且有很多震级很大的事件,可为研究俯冲断层地震孕育和发生机理提供数据支持.同时,理解俯冲带上逆断层地震发生过程也对以后的地震预测研究提供理论依据.已有的研究结果也表明,日本俯冲带的地震活动对我国的地震活动有重要的影响.数值模拟方法近年来发展迅速,并已大量用于地学领域.因此,本文结合地震地质资料和地球物理观测数据,使用有限单元数值模拟方法模拟地震黏滑失稳过程,探讨俯冲带地震发生的黏滑机制和构造变形特征具有重要的理论意义和应用价值.

2 日本俯冲带地震动力学数值模拟 2.1 日本俯冲带构造背景

在新生代,西太平洋岛弧和边缘海地区已成为强烈的构造活动区.至晚第三纪以来,在西太平洋岛弧带上还发生了强烈的火山作用.在西太平洋岛弧的外侧发育着地球表面上最深的沟谷带,即深海沟带.西太平洋板块向西俯冲,在日本海形成俯冲带,称为日本海俯冲带(图 1).日本俯冲带位于日本岛弧东侧,北起千岛,南至菲律宾海,是西太平洋边缘俯冲带的重要组成部分.该俯冲带处于环太平洋-火山地震带上,地质构造复杂,地震活动频繁,地震发生数量占全球20%.沿千岛群岛和日本本州岛的东侧,太平洋板块正以不同倾角向欧亚板块俯冲.在本州东部的日本海沟,震源分布显示的消减板片平均厚度大约80~100km[8],倾角29°,下插的最大深度接近600km.太平洋板块的北西向运动和俯冲,对从北海道到本州的日本东北部产生SEE-NWW向挤压.地震断层面解显示,日本海的东边缘也同样是SEE-NWW向挤压变形[9],可能是阿穆尔板块的东向移动产生的.由于日本俯冲带处于这样的挤压环境下,导致该带发生的大部分地震都具有逆断层性质.西太平洋地区震源分布图显示,沿着堪察加半岛-千岛群岛-日本本州岛北部,分布有大量浅源逆断层型地震,而中深源地震多分布在岛弧下面的俯冲带中.这种地震震源机制也印证了太平洋板块和欧亚板块的相聚,意味着该区域处于压缩应力状态.2011年3月11日日本本州东海岸附近海域发生9.0级特大地震,是日本有历史记载以来最大的一次地震.震后数十分钟内日本东海岸又遭受最高浪达10m的海啸袭击.鉴于该区域地震活动剧烈,破坏力大,并且对中国大陆,特别是华北、东北地区的现今构造变形有较大影响,本研究参考前人的研究成果[10-11],选取垂直俯冲带走向,切过9.0级大地震震源区的剖面(图 1)作为研究对象,建立日本俯冲带二维有限元数值模型(图 2),研究俯冲带的断层黏滑过程,模拟日本俯冲带的地震活动.

图 1 研究区构造背景.图中实线表示板块边界;虚线表示推测的阿穆尔板块与鄂霍茨克海板块边界;实心箭头表示太平洋板块相对于欧亚板块的运动方向;锯齿状曲线代表海沟[12].五角星表示9级大地震震源位置;剖面AA′为研究剖面 Fig. 1 Tectonic setting of the study area.Bold lines denote the plate boundaries; the dashed line shows the estimated boundary between the Amurian and Okhotsk plates; the bold arrow indicates the motion direction of the Pacific plate relative to the Eurasian plate; zigzag curves are trenches[12].Star represents the location of Mw9.0 earthquake; profile AA′is study profile
图 2 日本俯冲带二维几何模型 (a)俯冲带在岩石圈内部弯折;(b)俯冲带在岩石圈底部弯折;(c)俯冲带无弯折,倾角为30°. Fig. 2 Two-dimensional geometric model of Japan subduction slab Two-dimensional geometric model of Japan subduction slab (a) Subduction slab bending in lithosphere; (b) Subduction slab bending at the bottom of lithosphere; (c) Subduction slab without bending, dip angle is 30°.
2.2 黏弹介质断层接触有限元模型

根据地质构造背景及地球物理探测得到的结果,本研究建立了横跨日本俯冲带的二维有限元模型.以往的研究中多使用弹性介质,主要是因为弹性介质变形本构关系简单,参数相对较少,比较容易获得数据,同时弹性介质模型的计算相对容易实现,对计算设备不需要太高的要求.但是使用弹性物性参数模拟出的结果仅能近似地壳浅部的快速变形,由于实际地质体结构的复杂性,使得使用简单的弹性本构关系无法得到精确的结果.采用黏弹性介质,则赋予介质黏性及弹性的属性,在外力的作用下,不但产生弹性变形,而且变形随时间变化,更加贴近真实地质体,反映实际情况.另外,模型中引入了有限元分析中的接触对,模拟断层摩擦滑动.在俯冲带与上覆板块之间设置摩擦接触对,通过摩擦单元的运用,可以计算获得模型黏滑失稳产生的黏着、滑动的现象,模拟分析地震孕育、发生的过程.并且,在计算过程中,考虑了重力的作用,以及初始应力场的影响等.

模型东西宽度1200km,纵向深度为660km.俯冲带两侧地块分为上地壳、下地壳、上地幔上部(地幔岩石圈)和上地幔下部四层.俯冲带板片包括两层:俯冲地壳和俯冲岩石圈地幔.参考前人的研究结果,并综合考虑模型模拟计算的复杂程度,将模型中地壳厚度设为30km,其中上地壳为7km,中地壳为10km,下地壳为13km;岩石圈厚度为100km,上地幔延伸到560km;俯冲板片整体厚度为80km,其中俯冲地壳8km,俯冲上地幔岩石圈72km.俯冲板片向西俯冲,在560~660km处转为水平.俯冲断层活动不仅与动力边界和岩层介质性质有关,还受到俯冲带几何形状很大的影响.因此本文在模拟日本俯冲带黏滑失稳过程的基础上,重点研究了几何形状的变化对黏滑事件发生的影响.

本文建立3种几何结构模型,分别进行断层黏滑运动过程模拟,然后进行对比分析:

1)模型A:俯冲板片在岩石圈内部深度为80km处拐折,浅部倾角设为23°,较深部分设为30°(图 2a);

2)模型B:俯冲板片在岩石圈底部深度为100km处拐折,浅部倾角设为23°,较深部倾角设为30°(图 2b);

3)模型C:俯冲板片平直向下俯冲,俯冲倾角设为30°(图 2c).

2.3 模型参数及边界条件

各个几何模型中相同地质体使用的介质参数相同,即三个模型使用同一套介质参数(表 1).模型介质参数所对应的地质体分别为,上地壳(A1)、下地壳(A2)、上地幔岩石圈(A3)和上地幔下部.在俯冲板片的左右(东西)两侧,上地幔下部的介质参数存在差异,可分为两类:西侧上地幔(WM)和东侧上地幔(EM).俯冲板片虽然都为海洋岩石圈,但深度不同,温度压力条件不同,其介质参数也将存在很大的差异.在本研究中,随深度的增大,俯冲板片地壳部分共分4段:slab1-slab4,分别对应海水、陆地地壳、上地幔岩石圈和上地幔下部的深度段;地幔岩石圈部分分为两段,belt1和belt2,对应岩石圈和上地幔下部的深度段(图 2a).介质参数主要包括密度、杨氏模量、泊松比和黏滞系数.其中上地壳认为是弹性介质,其余各层均为Maxwell体.参考前人使用的参数[13-17],经过反复试算,得到认为较合理的参数(表 1).

表 1 日本俯冲带有限元模型介质参数 Table 1 Material parameters of blocks in the Japan subduction slab

本研究采用黏弹性有限元方法模拟接触问题.初始应力场对断层面的接触状态有很大影响,主要与接触问题的高度非线性有关.同时黏弹性的介质属性也决定了模型计算的时间需要一定的尺度.鉴于此,初始应力场在模拟中是必须要考虑的.目前初始应力场仍无法完全确定,本研究采用如下方法还原初始应力场.

在模型建立后只施加重力作用,使得在垂向上达到重力均衡[18],之后再在水平方向上施加自东向西的水平位移.太平洋板块相对于北美板块以70~90mm/a的速率向西运动,并向西俯冲到日本海沟下方.模型中按照80mm/a的速率计算,时间为200万年,相当于从第四纪早期开始,此时的模型状态作为初始状态,应力场为初始应力场.在此状态下再施加以下边界条件(图 3):模型底边界垂直方向固定,西边界岩石圈部分水平方向固定,东边界施加自东向西80mm/a的水平运动速率,作用时间为10万年,累计水平位移为8000 m.上边界陆地顶面认为是自由边界,上边界海洋部分由于受海水重力作用,在slab1板片及西侧相邻地块界面上施加垂向及侧向海水压力.同时考虑到上地幔深部岩石的围压作用,在西边界岩石圈之下到俯冲带接触面施加自西向东随深度递增的法向静岩压力(西边界递增箭头).

图 3 日本俯冲带模型边界条件示意图 Fig. 3 Schematic diagram showing boundary constrains of Japan subduction slab model
3 结果

为了使计算结果准确可靠,本研究在ANSYS并行计算平台上,模拟计算了西太平洋俯冲带在外载边界条件作用下10万年的变化情况.各个模型均计算2万步,步长为5年,得到了各节点随时间过程的应力、应变、位移、位错、应变能、接触力、摩擦力、接触状态等海量数据结果(每次计算结果都为180G以上),并在大型图形工作站上进行多视角的后处理,挖掘蕴含在数据中断层活动信息.为了更直观地对比分析各个模型发生黏滑事件的特点与异同,分别绘制了不同模型黏滑事件滑移量与深度随时间分布图(图 4-6),每个模型得到两张图(a、b).因接触面上的各节点在同一计算步中有不同的滑动量,所以这里挑选所发生的黏滑事件中滑移量最大的值作为该时刻目标事件的断层位错量,同时提取滑移量最大值所对应的节点深度作为事件发生深度.观察所得目标事件滑移量随时间分布规律,即第一张图(a).根据黏滑事件数据提取的坐标数据绘制事件发生深度随时间分布图,即第二张图(b),以便分析计算时长内黏滑事件发生深度的特征,从而可以判断日本俯冲带地震发生深度分布规律.

图 4 俯冲带在岩石圈内部弯折模型黏滑事件滑移量(a)和发生深度(b)随时间变化分布图 Fig. 4 Slip (a) and depth distribution (b) with time of stick-slip events in subduction slab with bending in lithosphere
图 5 俯冲带在岩石圈底部弯折模型黏滑事件滑移量(a)与发生深度(b)随时间变化分布图 Fig. 5 Slip (a) and depth distribution (b) with time of stick-slip events in subduction slab with bending at the bottom of lithosphere
图 6 俯冲带无弯折,倾角30°模型事件滑移量(a)与发生深度(b)随时间变化分布图 Fig. 6 Slip (a) and depth distribution (b) with time of stick-slip events in subduction slab without bending and with dip angle 30°

地震发生是由于断层解锁发生错动所致,断层无滑动时处于闭锁阶段,积累能量,能量达到一定程度,断层发生错动,产生滑移,形成一次地震.所以,根据每一子步时刻和该时间点最大滑移量对应节点坐标,可以绘制俯冲带上震源深度随时间变化分布图.随着计算时间的不断加长,模型在外载边界条件的作用下,俯冲带上的接触状态也发生复杂的变化.当俯冲板片与上覆板块之间无滑动时,断层处于闭锁阶段,积累应变能.能量超过断层所能承受的静摩擦力时,发生破裂,突然错动,断层解锁,地震发生.因此,认为当有滑移事件发生时可以近似为一次地震事件.模型中提取的滑移量是断面上的最大错动量,一般情况下大于地表断层位错.

图 4为模型A黏滑事件滑移量与发生深度随时间变化分布图.图(a)为滑移量随时间分布图.从图中可以看出,滑移量20 m左右的较大黏滑事件有3次,分别为24.6m,22.8m和18.6m.从发生时间上看,三次事件分别发生在29226年,54770年和94066年,时间间隔大致相当.滑移量5m及以上事件分布也较均匀,较规律地出现在三次大事件的间隔中.这种现象与地震发生存在平静期、活跃期比较相似.此外还有大量不足1 m的小事件,成簇地出现.而且从图中清楚地发现,这些小事件大多情况下紧随在一次较大事件发生之后成簇出现.大部分情况下,大震发生与大量的小事件之间时间衔接比较近,有可能是大事件的前震或余震,但个别事件,例如滑移量为22.8m的事件在发生前后均无小事件发生,表现为孤立型地震特征.也有部分小事件在两次较大事件之间,属震间背景地震活动.实际地震类型存在主震-余震型和孤立型,模拟结果显示的现象与地震发生类型有很好的相似性.

图(b)显示的是断层错动时刻最大滑移发生的深度.由于深部地震发震机理不明确,故本文只讨论深度小于300km的区域.图中圆圈的大小代表事件量级的大小(图 5b6b同).图中显示三处非常大的圆,发生深度分别为114.5km,20.3km和115.4km.对照图a可知是滑移量20m左右的三次大事件.随着滑移量逐渐变小,圆的直径也逐渐变小.从图中看出,除了三次大事件之外,较大一些的圆多分布在30km深度范围.此外,图中明显看出,横向上主要有两个地震密集条带,分别位于20~30km和80km左右.80km处存在最大滑动密集条带,有可能与模型在80km处存在角度变化有关,在几何拐点处容易产生应力集中,也可能造成较多滑动.在20~30km处,事件颇为密集,且在整个计算过程中分布都比较连贯,只在大约68000~88000年出现间断.从该阶段对应的黏滑事件滑移量分布图(图 4a)可以看出,该深度上此时段出现在两次巨大的滑移事件发生间隔的中段,表现出清晰的震间平静的特征.

图 5为模型B黏滑事件滑移量与发生深度随时间变化分布图.图 5a为滑移量随时间分布图,图中显示,滑移量20m及以上的事件有两个活跃期,且活跃期大事件的发生深度都在30km附近.图中显示滑移量5 m及以上的事件数量比俯冲带在岩石圈内部弯折的模型A(图 4a)多很多.且5m以下的事件更多,几乎分布于整个时间过程.从图中还可以明显看出,在一次大的事件发生前,有大量较小事件发生,且小事件的滑移量大小随时间呈增长趋势.这种现象在20000~23000年,70000~80000年和80000~85000年三处最突出.还有的在一次大事件发生后,紧接着出现大量小事件,且小事件的滑移量具有随着时间增加呈衰减趋势,如29000~37000年、43000~50000年和89000~92000年处.还有的事件在发生前后都有大量小事件发生,这些现象都与前震-主震-余震发生模式有很好的对应.

图 5b是黏滑事件发生深度随时间变化分布图.图中整体显示出,在深度100km以上区域的事件数量要远远多于其以下区域.图中显示了清晰的密集条带,分别位于20~30km,100km和125km附近三个区域.20~30km事件密集带与浅源地震对应,100km和125km事件密集带是俯冲带在岩石圈底部弯折和力学性质变换所致.角度变化的区域容易导致应力集中,故而该区域事件数量亦较多.从滑移量大小看,较大滑移量的事件多发生在30km附近区域.而100km和125km深度的密集带事件中5m以下滑移量事件占绝大多数.

图 6为模型C黏滑事件滑移量与发生深度随时间变化分布图.图 6a为滑移量随时间分布图,滑移量20m以上事件共8次,是三个模型中滑移量在20~30m的黏滑事件数量最多的模型.但滑移量在5~20m范围内的黏滑事件数量却是三个模型中最少的.10m及以上,20m以下滑移量黏滑事件数目与俯冲带在岩石圈弯折的模型A相当.5 m及以上,10m以下滑移量的黏滑事件数量远小于其它两个模型.滑移量在5m以下的小事件,量级都较小,基本都是1m以下的较小事件,与模型A相似,只是在20000~30000年之间出现滑移量稍微大些且较成簇出现的小事件.对于小事件,在整个计算时长中随机分布.

观察黏滑事件发生深度随时间变化分布图 6b看出,该模型黏滑事件基本上都分布在100km以上的深度内.尤以20~30km和100km的地方较为突出,呈现明显的条带状.100km以下至300km的区域内也有少量分布.

4 结论与讨论 4.1 初步认识

根据上述模拟结果及分析,得到以下初步认识:

(1)基于黏弹介质和断层接触模型,利用有限元数值模拟方法可以模拟俯冲带地震发生的过程(断层闭锁、解锁再到闭锁的黏滑失稳过程),创建人工地震目录.

(2)模拟得到的断层错动滑移量与实际地震断层错动量相当,并且滑移量越大的事件,其数量越少.黏滑事件的大小、时间间隔、发生深度是随机性和有序性的统一,宏观上有序(准周期、特征震级和多震层密集),微观上随机(时间、深度、滑移量等都不完全相同).

(3)从黏滑事件发生深度看(图 7),三个模型模拟得到的黏滑事件,大多发生在100km以内的深度区域,岩石力学性质变化和断层几何形态变化的深度是黏滑事件多发的区域,主要表现为存在两个明显的条带状事件密集区.这两个密集区分别为20~30km范围和100km附近深度.其中深度为20~30km范围内地震数量最多,且事件量级也较大;100km附近深度虽然也有大量事件发生,但总体数量较前者少,且以相对小事件为主,较大量级事件比较少.

图 7 各个模型黏滑事件发生深度数量统计图 Fig. 7 Statistics of events in different depth intervals of 4 models

(4)俯冲带角度变化及发生变化的深度不同对模型模拟结果有影响.模型B模拟得到的黏滑事件数量最多.模型C模拟得到的黏滑事件数量最少.但值得注意的是,前者多为滑移量5m以下的相对小事件.而模型C得到的滑移量20 m及以上黏滑事件比例最高,可能是因为模型C相比存在角度变化的模型A和模型B来说,模型C断层无弯折,断层长度变长,上覆板片与俯冲带发生错动时,能够同时错动的断层面积更大.模型A和模型B滑移量20m及以上黏滑事件数量所占各自事件总量比例相当,对于滑移量5 m及以下相对小事件来说,模型A发生的此类黏滑事件占其总量比例最高,而模型C则是比例最小的.也就是说模型A更倾向于发生量级较小的事件.

4.2 存在的问题

计算过程中,摩擦系数的选取十分重要.本文模型中考虑到接触摩擦问题,使用的是速度状态摩擦本构关系[19].根据已有对俯冲带摩擦系数的研究结果[20-26],模型中静摩擦系数取值0.2,动摩擦系数取值0.15.摩擦系数选取太大,那么在整个计算过程中,有可能俯冲带一直处在闭锁状态,不能产生“地震”现象;但如果摩擦系数取得太小,能量积累较少,则摩擦滑动的幅度减小,滑移之间的时间间隔也相应较短,即产生不了“大地震”现象.只有选取合适的摩擦关系,才能产生类似于天然地震的“地震”现象.

计算中,时间步长的选择也非常重要.步长太大,有可能漏掉中间大量的细节,滑动位移的突然变化被平滑或放大,则可能无法产生突变位移(即“大地震”现象)或产生超大事件;时间步长选择太小,则需要大量的计算时间,一般的机器很难满足,没有超级计算设施是不可能实现的.同时,计算中产生的海量数据对存储媒介及显示设备都有很高的要求.对于模拟地震发生的具体过程来说,时间步长以“秒”或“十秒”的量级为佳,然而地震孕育过程需要几百年甚至更长的时间,本文模拟的时间尺度更是达到了十万年之久,以5年为步长计算,数据量已经相当庞大,以秒为量级计算目前无法实现.因此,新的计算方法及应用新的技术是今后工作中需要解决的重要问题.

初始应力场的确定目前仍没有统一的数据或方法.文中虽然在理论上恢复了初始应力场,再在此基础上施加外载进行计算,但真实的初始应力场对结果的准确性和真实性仍是有影响的.如何能更准确地逼近初始应力场也是今后研究工作中需要考虑的问题.

致谢

在本研究工作中,张培震、王海涛、马胜利、陆远忠、冉勇康、刘杰、周永胜、甘卫军研究员提出了宝贵的建议,在此表示感谢.

参考文献
[1] 温瑞智, 任叶飞, 李小军. 日本Mw9.0级地震海啸数值模拟与启示. 国际地震动态 , 2011(4): 22–27. Wen R Z, Ren Y F, Li X J. The tsunami simulation for off the Pacific coast of Tohoku earthquake and disaster mitigation in China. Recent Developments in World Seismology (in Chinese) , 2011(4): 22-27.
[2] Ozawa S, Nishimura T, Suito H, et al. Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake. Nature , 2011, 475(7356): 373-376. DOI:10.1038/nature10227
[3] Kerr R A. A tantalizing view of what set off Japan's killer quake. Science , 2012, 335(6066): 272. DOI:10.1126/science.335.6066.272
[4] Kato A, Obara K, Lgarashi T, et al. Propagation of slow slip leading up to the 2011 Mw9.0 Tohoku-Oki earthquake. Science , 2012, 335(6069): 705-708. DOI:10.1126/science.1215141
[5] 徐彦. 2011年日本9.0级及7.5级地震震源破裂反投影初步结果. 国际地震动态 , 2011(4): 34–37, 5. Xu Y. Back projection results of the 2011 Japan Mw9.0 and M7.5 earthquakes. Recent Developments in World Seismology (in Chinese) , 2011(4): 34-37, 5.
[6] Shinzaburo O, Takuya N, Hisashi S, et al. Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake. Nature , 2011, 475: 373-376. DOI:10.1038/nature10227
[7] Kelin W L, Yan H, He J H. Deformation cycles of subduction earthquakes in a viscoelastic Earth. Nature , 2012, 484(7394): 327-332. DOI:10.1038/nature11032
[8] Fukao Y, Obayashi M, Nakakuki T, et al. Stagnant slab: A review. Annu Rev Earth Planet Sci , 2009, 37(1): 19-46. DOI:10.1146/annurev.earth.36.031207.124224
[9] DeMets C, Gorden R G, Argus D F, et al. Current plate motions. Geophys. J. Int. , 1990, 101(2): 425-478. DOI:10.1111/gji.1990.101.issue-2
[10] Yoshi T. A detailed cross-section of the deep seismic zone beneath northeastern Honshu, Japan. Tectonophysics , 1979, 55(3-4): 349-360. DOI:10.1016/0040-1951(79)90183-5
[11] Yoshioka S, Loo H Y, Mikumo T, et al. A model of post-seismic recovery induced by a deep-focus earthquake. Physics of the Earth and Planetary Interiors , 1992, 72(1-2): 83-98. DOI:10.1016/0031-9201(92)90051-V
[12] Huang Z C, Zhao D P, Wang L S. Shear wave anisotropy in the crust, mantle wedge, and subducting Pacific slab under northeast Japan. Geochem., Geophy. Geosyst. , 2011, 12(1): Q01002. DOI:10.1029/2010GC003343
[13] Loo H Y, Gao X L, Sun J X, et al. Three-dimensional numerical modeling of earthquake migration along a northwestern Pacific subduction slab. Geophys. Res. Lett. , 1992, 19(3): 313-316. DOI:10.1029/92GL00021
[14] 张建, 汪集旸. 南海北部陆缘带构造扩张的深部地球动力学特征. 中国科学D辑 , 2001, 44(5): 437–445. Zhang J, Wang J Y. Geodynamic characteristics of tectonic extension in the northern margin of South China Sea. Science in China (Series D) (in Chinese) , 2001, 44(5): 437-445. DOI:10.1007/BF02909782
[15] 臧绍先, 宁杰远, 景志成. 俯冲带流变性质的研究. 中国科学D辑 , 2001, 44(12): 1119–1127. Zang S X, Ning J Y, Jing Z C. Study on the rheology of subducting slabs. Science in China (Series D) (in Chinese) , 2001, 44(12): 1119-1127. DOI:10.1007/BF02906868
[16] Zhao D P. Global tomographic images of mantle plumes and subducting slabs: insight into deep Earth dynamics. Phys. Earth. Planet. Inter. , 2004, 146(1-2): 3-34. DOI:10.1016/j.pepi.2003.07.032
[17] Zhu G Z, Shi Y L, Tackley P. Subduction of the Western Pacific Plate underneath Northeast China: Implications of numerical studies. Physics of the Earth and Planetary Interiors , 2010, 178(1-2): 92-99. DOI:10.1016/j.pepi.2009.10.008
[18] 朱守彪, 张培震. 2008年汶川Ms8.0地震发生过程的动力学机制研究. 地球物理学报 , 2009, 52(2): 418–427. Zhu S B, Zhang P Z. A study on the dynamical mechanisms of the Wenchuan Ms8.0 earthquake, 2008. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 418-427.
[19] 朱守彪, 邢会林, 谢富仁, 等. 地震发生过程的有限单元法模拟--以苏门答腊俯冲带上的大地震为例. 地球物理学报 , 2008, 51(2): 460–468. Zhu S B, Xing H L, Xie F R, et al. Simulation of earthquake processes by finite element method: The case of megathrust earthquakes on the Sumatra subduction zone. Chinese J. Geophys. (in Chinese) , 2008, 51(2): 460-468.
[20] 石耀霖, 王其允. 斜俯冲板块边界变形分配的力学分析. 地球物理学报 , 1994, 37(5): 606–612. Shi Y L, Wang Q Y. Mechanics of deformation partitioning at plate boundaries of oblique subduction. Chinese J. Geophys. (in Chinese) , 1994, 37(5): 606-612.
[21] 李大鹏, 陈岳龙, 靳野. 板块俯冲带研究中的数值实验. 地球科学进展 , 2010, 25(6): 582–596. Li D P, Chen Y L, Jin Y. Numerical simulation in subduction zone study. Advances in Earth Science (in Chinese) , 2010, 25(6): 582-596.
[22] Gorczy W, Willner A P, Gerya T V, et al. Physical controls of magmatic productivity at Pacific-type convergent margins: Numerical modelling. Physics of the Earth and Planetary Interiors , 2007, 163(1): 209-232.
[23] Hall C E, Gurnis M, Sdrolias M, et al. Catastrophic initiation of subduction following forced convergence across fracture zones. Earth and Planetary Science Letters , 2003, 212(1-2): 15-30. DOI:10.1016/S0012-821X(03)00242-5
[24] Hassani R, Jongmans D, Chery J. Study of plate deformation and stress in subduction processes using two-dimensional numerical models. Journal of Geophysical Research , 1997, 102(B8): 17951-17965. DOI:10.1029/97JB01354
[25] Sobolev S V, Babey A Y. What drives orogeny in the Andes?. Geology , 2005, 33(8): 617-620. DOI:10.1130/G21557.1
[26] Tagawa M, Nakakuki T, Kameyama M, et al. The role of history-dependent rheology in plate boundary lubrication for generating one-sided subduction. Pure and Applied Geophysics , 2007, 164(5): 879-907. DOI:10.1007/s00024-007-0197-4