地球物理学报  2012, Vol. 55 Issue (1): 117-125   PDF    
马尼拉海沟俯冲带热结构的模拟研究
高翔1, 张健1 , 孙玉军1, 吴时国2     
1. 中国科学院研究生院计算地球动力学重点实验室, 北京 100049;
2. 中国科学院海洋研究所海洋地质与环境重点实验室, 青岛 266071
摘要: 俯冲带热结构的数值模拟研究是对地表观测研究的重要补充,也是验证地球动力学模型的重要方法.本文沿马尼拉海沟俯冲带东火山链(EVC)和西火山链(WVC)各取一条剖面,依据地质、地球物理条件,进行了有限元热模拟计算.计算过程中,分析了摩擦和剪切热对俯冲带热结构的影响,模拟了EVC和WVC两条测线下俯冲带的热结构,并结合岩石学实验结果预测了俯冲板块发生脱水和部分熔融的位置.模拟结果表明,在100 km深度处,考虑摩擦和剪切热时,俯冲板块表面的温度约为865 ℃;而不考虑摩擦和剪切时,俯冲板块表面的温度仅为770 ℃,二者温差可达95 ℃.在相同深度处,考虑摩擦和剪切热时,在EVC和WVC测线下俯冲板块表面的温度分别为865 ℃和895 ℃,俯冲洋壳底部温度分别为560 ℃和605 ℃.俯冲板块表面少量矿物开始脱水的深度小于50 km,但大量脱水和部分熔融主要发生在深度100 km左右,这与地表观测的火山活动位置一致.
关键词: 马尼拉海沟      俯冲带      热模拟      摩擦和剪切热      有限元     
A simulation study on the thermal structure of Manila trench subduction zone
GAO Xiang1, ZHANG Jian1, SUN Yu-Jun1, WU Shi-Guo2     
1. Key Laboratory of Computational Geodynamics, Graduate University of Chinese Academy of Sciences, Beijing 100049, China;
2. Key Laboratory of Marine Geology and Environment in the Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China
Abstract: Numerical simulation study of thermal structure is a supplement to field observations and an important method to verify geodynamic models. In this study, based on geological and geophysical conditions, we construct two thermal models across the East Volcano Chain (EVC) and the West Volcano Chain (WVC) along the Manila trench subduction zone and do the thermal simulations using finite-element numerical method. During these simulations, we study the thermal structure of the subduction zones affected by the frictional and shear heating and simulate the thermal structure of the subduction zones beneath the seismic lines at EVC and WVC. We also employ the results of rock experiments to predict the regions of dehydration and partial melting in the subducting slab. The simulation results show that, at 100 km depth, at the slab interface with frictional and shear heating the temperature is 865 ℃, about 95 ℃ warmer than temperature without frictional and shear heating. Considering the frictional and shear heating, the temperatures at the slab interface are 865 ℃ and 895 ℃ and at the base of subducting oceanic crust are 560 ℃ and 605 ℃ at 100 km beneath the seismic lines of EVC and WVC. At slab surface, the depth of dehydration onset is less than 50 km. However, partial melting and dehydration may mainly occur at about 100 km depth, which is well associated to the observed location of the volcanic activity.
Key words: Manila trench      Subduction zone      Thermal simulation      Frictional and shearing heat      Finite-element     
1 引 言

俯冲带的热结构与俯冲板块上覆地幔楔的部分熔融、岛弧火山的形成及地震发生机制有着紧密的联系[1-2].摩擦和剪切热是影响俯冲带热结构的主要热源之一,它能提高俯冲板块与上覆板块界面周边的温度,从而改变俯冲板块内含水矿物的脱水深度和板块界面周边岩石的强度.近年来,断层间的摩擦和剪切热与断层带上地震活动的关系得到了深入的研究[3-5];而海洋俯冲带板块界面间的摩擦和剪切热对俯冲带热结构的影响并没有单独的分析,在计算中基本都采用简单近似的处理方式,有的模型将摩擦和剪切热视为一个不变的热源加到俯冲界面之上[1],还有一些模型在计算中将其忽略[6-7].对俯冲带的热结构,特别是俯冲带界面附近的热结构的研究,在计算上采用不同的处理方式会对计算结果产生很大的差异,因此,选用合理的处理方式尤为重要.

前人曾通过构建理想的数值模型,来分析影响俯冲带热结构的主要因素对俯冲带热结构的影响[7-12],这些因素包括俯冲板块的年龄、板块的相对运动速度、俯冲角度和上覆岩石的结构等.这些研究方式对于阐明单一因素对俯冲带热结构的作用是很有意义的,例如:Peacock 和Wang(1999)[12]通过对日本东北部年龄较老的快速俯冲和西南部年龄较新的慢速俯冲板块的内部温度进行了比较,发现前者温度比后者低了约300℃,并指出俯冲板块的年龄越老,俯冲速度越快,板内温度越低.但是这些影响因素不是单一存在,而是相互作用的,马尼拉海沟俯冲板块的年龄、运动速度和俯冲角度在EVC(东火山链)和WVC(西火山链)下都存在着差异[13-16],因此,结合俯冲带板块界面间的摩擦和剪切热,综合分析这些因素对俯冲带热结构的影响,有利于更好地解释观测结果.

俯冲板块内不同含水矿物的脱水温度和含水量差异很大[17],而俯冲带的热结构直接影响俯冲板块含水矿物的脱水位置和脱水量.一方面,在上覆岩石圈地幔内,脱去的水会使地幔橄榄石发生蛇纹石化,从而改变俯冲界面岩石的强度及俯冲界面的解耦深度[18];另一方面,在岛弧火山下方,脱掉的水会降低俯冲板块上覆高温地幔岩石熔点,引发岩石的部分熔融,进而产生岩浆和火山活动[19].

本文中,我们将建立二维的数值模型,采用有限元数值模拟方法求解非线性瞬态的热方程.首先,通过选取不同的摩擦系数研究摩擦和剪切热对俯冲带热结构的影响;其次,根据马尼拉海沟EVC 和WVC 下俯冲板块的角度、年龄和运动特征的差异,分析两处俯冲带的热结构;最后,结合岩石学实验结果研究俯冲带热结构对俯冲板块内含水矿物脱水的影响.

2 地质条件

马尼拉海沟是菲律宾海板块向南海板块仰冲所形成的一条具有特殊构造意义的汇聚边界,其西侧为南海海盆,东侧为台湾-吕宋火山岛弧,南起民都洛岛西南陆架的海底峡谷,向北水体变浅并一直延伸到台湾碰撞构造带中,空间上呈南北弧形展布.地质学和火山学研究表明,台湾-吕宋岛火山弧由北向南在吕宋岛分叉为WVC 和EVC,WVC 内的岩浆活动停息于4~2 Ma,而EVC 内的岩浆活动几乎完全在第四纪,由于南海古洋中脊的俯冲,在WVC和EVC 之间约17.8°N 出现了宽约50km 的地震和火山活动的空白区域[15](图 1).

图 1 马尼拉海沟的区域构造.实心三角代表火山;箭头指向是板块的运动方向;大写字母表示的线段(如:A-A')是图 2剖面所在的位置 Fig. 1 Geotectonic map of Manila Trench. Solidtriangles: volcanoes;arrow : movement direction of the slabs;lines with capital letters (ex.A-A′):the positions of profiles shown in Fig. 2

在菲律宾海板块之下,南海俯冲板片的形态非常复杂,俯冲角度从22°N 到18°N 逐渐变缓,而从18°N 以南开始,板片的俯冲角度陡增并向南逐渐变陡,直到最南端的民都洛岛俯冲角度接近90°.图 2中4条剖面的位置对应于图 1 中的测线A-A′,B-B′,C-C′,D-D′,圆点代表的是震源,实线代表的是Bautista等(2001)[16]根据地震分布和震源机制解推断的板片表面的形态.在增生楔前端,俯冲板块某个界面之上的沉积物被褶皱冲断增生于增生楔中,我们将此界面称为俯冲拆离面,该界面之下的海盆基底和上覆沉积层一起沿该拆离面向海沟俯冲[20].沿着海沟方向,俯冲板块拆离面之上沉积物的厚度一般在1.5~4km 之间,之下的洋壳厚度介于4~6km[13].海沟西侧的台湾-吕宋岛地壳底界埋深约为34km,岩石圈在海沟两侧的底界埋深为70~80km[21].

图 2 地震垂向分布(圆点)和推断的俯冲板块表面形态(实线).水平轴的零点对应的是马尼拉海沟轴线位置,据(Bautista等,2001)[16] Fig. 2 Vertical distribution of seismicity (dots) and inferred slab surface geometry (solid line).The zero distance on the X-axis is the trench-axis of the Manila trench. From (Bautista el al,2001)[16]

李家彪等(2004)[20]通过对海沟增生楔上逆冲断层和褶皱走向的分析,指出沿马尼拉海沟板块的俯冲应力方向为NW55°.Bautista等(2001)[16]依据NUVEL-1模型,指出菲律宾海洋板块的绝对运动速率在吕宋岛最北端约为7cm/a,并向南逐渐增大,在棉兰老岛的东南部达到9cm/a,运动方向为NW55°,与李家彪等(2004)[20]的结论一致;另一方面,欧亚板块在相近的方向以1cm/a的慢速运动,尽管二者在同一方向运动,但由于菲律宾海板块的运动速率比欧亚板块快得多,所以二者表现出强烈的会聚[16].

3 计算模型

计算过程中,我们首先确立研究马尼拉海沟俯冲带热结构的计算方程;其次,根据研究区的地质和地球物理特征构建计算所需的模型,并对模型添加材料参数、初始条件和边界条件,利用有限元的方法对模型进行三角形网格剖分;最后,运用有限元方法求解计算方程.

3.1 计算方程

计算马尼拉海沟俯冲带地区的热结构,我们采用的二维热方程如下:

(1)

其中ρ为岩石密度,c为岩石比热,k为热导率,珔μ 为板块间的相对运动速度,T为开尔文温度,Q为热源.

计算中热源项Q主要包括摩擦剪切生热Qf(T)和岩石放射性产热Qr.摩擦和剪切生热Qf(T)在俯冲界面解耦深度之上可表示为[22]

(2)

当俯冲界面岩石温度低,表现为脆性时,τ(T)大小由摩擦定律控制,可表示为

(3)

当俯冲界面岩石温度高,表现为韧性时,τ(T)大小与温度有关,可表示为[23]

(4)

深度和温度决定着俯冲界面岩石的脆性和韧性状态,也决定着τ(T)的大小,计算中取数值小的一项,即

(5)

(2)~(5)式中:τ(T)为板块间的摩擦和剪切力,σn为正应力,p为孔隙水压,g为重力加速度,H为俯冲板块接触面深度,ρ 为俯冲界面上覆岩石密度,μ为俯冲界面岩石间的摩擦系数,μ′ 为等效摩擦系数,An是与岩石性质相关的常数, 为剪切应变率,E为激活能,R为气体常数,T为Kelvin温度.我们采用将摩擦和剪切热添加到俯冲板块顶部100m厚的薄层的方法,来处理计算得到的板块界面每一深度的摩擦和剪切热.

摩擦和剪切生热Qf(T)在软流圈地幔楔处只存在剪切生热,计算方法为[18]

(6)

式中的τ(T)为软流圈内部岩石间的剪应力,计算方法与式(4)相同,软流圈的应变率 由角落流的速度场得出,速度场的计算公式为[24-25]

(7)

其中u是水平向速度,v是垂向速度,A= 0,B=-Cθ是俯冲角度.

岩石内部的放射性产热Qr 在地下不同深度和位置其值的大小不同,考虑到Qr 在同一圈层差异很小,因此我们在计算放射性产热率时取本圈层的一个平均值.对于板块在俯冲过程中发生弯曲和破裂导致的板内生热,因生热量与前两项相比非常少,在计算中没有考虑.

3.2 初始条件和边界条件

我们设定计算区域包括从马尼拉海沟向台湾-吕宋岛弧方向350km 的范围,深度为200km.顶层为3km 厚的沉积层,俯冲洋壳厚度为5km,上覆板块上地壳厚度为10km,下地壳厚度为14km,岩石圈底边界深度在俯冲板块和上覆板块之下分别为70km 和80km.坐标系水平向零点为海沟轴线位置,垂直向零点为大地水准面,图 3 为模型示意图.计算区域的初始温度采用半空间冷却模型的计算方法[26],根据上覆板块的年龄计算得出.上边界温度取为0 ℃,下边界200km 深处温度取为1350 ℃,海洋板块岩石圈边界和上覆板块岩石圈边界均采用半空间冷却模型计算出的温度分布作为边界条件,软流圈地幔左右边界均为自由边界.

图 3 计算模型示意图 Fig. 3 Schematic diagram of the computational model
3.3 模型参数

等效摩擦系数μ′ 的大小是决定俯冲界面产热量多少的重要因素.Morrow 等(2000)[27]对控制俯冲界面几种岩石的摩擦系数进行分析,得出起主要作用的滑石在潮湿环境下摩擦系数为0.2,水镁石、绿泥石和叶蛇纹石的摩擦系数分别是0.3、0.4 和0.5.Wang 等(1995)[28],Wada等(2008)[18]根据对Cascadia俯冲带和其他俯冲带热模拟结果的分析,并结合地表观测结果,得出俯冲界面的等效摩擦系数为0.03.基于前人对Cascadia俯冲带的研究已经很成熟[29],而马尼拉海沟俯冲带研究区又缺乏实测热流数据[30],并且Cascadia俯冲板块与马尼拉海沟俯冲板块都具有年龄新、沉积物厚和俯冲速度居中等特征[13, 31].因此,我们在计算中采用此数值.

板块在俯冲过程中会带动周边较软的物质一起运动,在俯冲板块上覆的地幔楔处引起角落流,Wada等(2008)[18]通过对全球多处俯冲带的模拟研究发现,虽然不同俯冲带俯冲板块的年龄、速度、角度及沉积物厚度不同,但板块界面间的解耦深度一般都介于70~80km,而LiZhiwei等[21]根据P 波速度成像结果得出马尼拉海沟周边地区岩石圈深度在70~80km 之间.因此,我们将角落流速度场的边界定为岩石圈底边界与俯冲板块的表面边界.

本文计算模型中所使用的物性参数(表 1)和岩石流变相关参数(表 2)主要采用前人成果.

表 1 计算模型中所使用的物性参数[13, 32] Table 1 Physical parameters used for the model[13, 32]
表 2 岩石流变相关参数[33] Table 2 Parameters of dislocation creep[33]
4 模拟结果及分析 4.1 摩擦和剪切生热对模拟结果的影响

首先以测线B-B′之下的俯冲板块为例,分析摩擦和剪切热对马尼拉海沟俯冲带热结构的影响,并把计算得到的地表热流值与前人的结果相比较,来约束我们的计算结果.

等效摩擦系数的大小直接影响摩擦和剪切生热量,因此,我们在计算中分别取等效摩擦系数为0、0.03和0.3三个不同数值,其中0表示计算中不考虑摩擦和剪切热;0.03 是本文计算马尼拉海沟俯冲带热结构所采取的等效摩擦系数;0.3 代表由实验测得的俯冲板块表面岩石矿物摩擦系数的平均值[27].

图 4是在取不同的摩擦系数时俯冲界面的压力(P)-温度(T)图.俯冲界面的温度随着摩擦系数的增大而升高,在深度100km 的温度分别为770 ℃、860 ℃和890 ℃.地质学研究表明马尼拉海沟俯冲带火山活动发生的深度约为100km[16],而俯冲板片在深度100km 处发生大量脱水熔融的温度至少需要825 ℃[34],而不考虑摩擦和剪切热时,温度仅为770℃.因此,摩擦和剪切热是引起马尼拉海沟俯冲带火山活动不可缺少的因素.

图 4 摩擦系数取不同值时,计算得到的俯冲界面的P-T Fig. 4 Lines are P-T paths for subduction interfacecalculated by different frictional coefficients

图 5是取不同摩擦系数时的地表热流分布情况.0km 所对应的位置是海沟轴线处,当摩擦系数为0和0.03时,热流变化趋势是一致的,表现为从海沟向岛弧方向逐渐降低,到达岛弧附近热流升高,这与Wada等(2008)[18]收集的全球多处俯冲带的实测热流的分布情况一致.而当等效摩擦系数为0.3时,从海沟向岛弧方向的热流不但没有降低,反而迅速增加,这与实际的观测结果严重不符,因此,将实验所得的摩擦系数直接应用于模拟计算中是不可取的.

图 5 摩擦系数取不同值时计算得到的地表热流 Fig. 5 The lines are surface heat flows calculated by different frictional coefficients

计算结果必须要以实测资料作为约束,以证明其正确性与可靠性.Shi 等(2003)[35] 与He 等(2001)[30]通过对实测数据的拟合,得出测线B-B′自海沟向岛弧方向150km 范围内的热流介于40~50mW·m-2之间,且两边高中间低,这与等效摩擦系数取为0.03时热流的计算结果相一致,如图 5所示.

4.2 马尼拉海沟俯冲带的热结构

图 6为测线B-B′和D-D′之下俯冲带的热结构,图中曲线是俯冲板块表面及距俯冲板块表面相同深度层面的P-T曲线,相邻层面间的距离为1km.在火山弧之下深度100km 处,俯冲板块界面的温度分别为865 ℃和895 ℃,俯冲板块洋壳底部温度分别为560 ℃和605 ℃,见图 6,洋壳底部的温度差异明显大于板块界面的温度差异,这主要是由于横向上的热传导引起的.而板块界面两侧存着很高的温度差异造成了俯冲板块表面区域有着较高的温度梯度.

图 6 测线B-B'(a)和D-D'(b)之下俯冲板块内的P-T图.曲线分别是板块界面和距离界面相同深度处间隔1 km的P-T Fig. 6 P-T paths for subducting slab beneath lines B-B' (a) andD-D' (b). Paths are located at slab interface and at 1km intervals in the subducting slab down to 6 km below the interface

测线B-B′和D-D′处的俯冲板块的特征存在很大差异,计算中板块现今年龄分别为35 Ma和20Ma[36],马尼拉海沟开始俯冲的时间为15 Ma[13],因此,开始俯冲时板块的年龄分别为20 Ma和5 Ma;板块相对运动速率分别为6cm/a和7cm/a;俯冲角度分别为42°和60°[16].而正是这些因素使得俯冲带热结构存在差异,我们将简要地把这些因素对俯冲带热结构的影响进行逐一分析.

图 7a是俯冲板块年龄对俯冲带热结构的影响,可以看出,年龄越老俯冲带的温度越低.引起这一现象的主要原因是俯冲板块自身的温度随着年龄的增加而降低,由于初始温度低,也必然会引起俯冲带内较低的温度.板块的俯冲速度对俯冲带热结构的影响则是速度越快温度越低,如图 7b所示,这主要是由于俯冲板块相对冷的物质快速下插,周边相对热的物质来不及将其加热,而下插的速度越快,加热的时间就越短,这样温度也就越低.图 7b中俯冲板块界面间温度差异小是速度取值接近和摩擦剪切热等因素引起的.图 7c中俯冲板块角度的变化对俯冲带深部的热结构没有明显的影响,而在相对较浅深度,俯冲角度越大,温度越低,在俯冲速率相同的情况下,俯冲角度越大,垂向速度就越快,因此,俯冲带的温度就会偏低.

图 7 俯冲板块的年龄(a)、俯冲速度(b)和俯冲倾角(c)对俯冲带热结构的影响 Fig. 7 The thermal structure affected by the subducting slab age (a) , convergence velocity (b) and dip angle
4.3 俯冲板块内矿物脱水与岩石部分熔融

岩石的高温高压实验结果为我们分析板块内含水矿物的脱水位置和岩石发生部分熔融的深度提供了很好的依据.Hacker等(2003)[17]通过分析洋中脊玄武岩和贫化的方辉橄榄岩的高温高压实验结果,得出了俯冲板块内不同矿物脱水及岩石发生部分熔融的温压条件,并构建了图 8所示的相变图.图中虚线为含水矿物变质脱水的临界状态,右侧的实线则是俯冲板块岩石在有水的条件下发生部分熔融的临界状态.我们将把计算得到的P-T图与实验相变图相结合,来分析马尼拉海沟俯冲带板块脱水与部分熔融的情况.

图 8 计算的P-T图与实验相图结合.(a)和(b)为玄武岩的相变图;(c)和(d)为贫化的方辉橄榄岩的相变图(据Hacker等,2003)[13].字母代表的含义:A为角闪岩,B为蓝片岩,E为榴辉岩,a为角闪石,e为绿帘石,j为硬玉,l为硬柱石,z为黝帘石.测线B-B'和D-D'之下俯冲板块内的P-T图分别用实线和虚线表示 Fig. 8 Calculated P-T paths for subducting slab superimposed on ultramafic phase diagram. (a) and (b) are the metamorphic facies for basalt,(c) and (d) are the metamorphic facies for depleted harzburgite (Hacker et a1.,2003)[13]. Abbreviations: A,amphibolite; B,blueschist; E,eclogite; a,amphibole; e,epidote;j,jadeite; l,lawsonite; z,zoisite. The solid lines and dot lines represent the P-T paths beneath lines B-B' and D-D' separately

图 8(a,d)中,我们发现俯冲洋壳表面少量矿物开始脱水的深度约为50km(约1.5GPa),俯冲岩石圈地幔开始脱水的深度约为100km(约3.2GPa).俯冲洋壳在浅部脱掉的水分,有助于上覆地幔楔岩石发生蛇纹石化,在深部脱掉的水则有利于板块上方高温岩石发生部分熔融.图 8(a,b)显示俯冲洋壳表面岩石在80km(约2.5GPa)已经开始发生部分熔融,距离俯冲洋壳表面2km 的洋壳岩石在100km也开始发生部分熔融.图 8c中可以看出俯冲洋壳上覆地幔岩石开始发生熔融的深度约为85km(约2.7GPa),而在100km 深度处俯冲界面的温度已高于地幔岩石熔融温度几十摄氏度,有助于俯冲板块上方岩石的大量熔融,进而引起火山活动.由于在相同深度俯冲板块的温度在测线B-B′要低于D-D′,并随着深度的增加温度差异也增大,因此,D-D′之下板块的脱水和部分熔融的位置更浅.

5 结 论

(1) 通过计算分析,我们发现摩擦和剪切热是模拟马尼拉海沟俯冲带热演化过程中重要的热源.前人通过大量模拟研究给出的俯冲板块界面间的等效摩擦系数值(0.03)适用于马尼拉海沟俯冲带.

(2) 依据俯冲板块的年龄、速度和角度等参数,在100km 深度处,计算得到马尼拉海沟俯冲带测线B-B′之下板块界面温度为865 ℃,俯冲洋壳底部温度为560℃;测线D-D′之下板块界面温度为895 ℃,俯冲洋壳底部温度为605 ℃.板块界面的温度差异小于俯冲洋壳底部的温度差异.

(3) 结合岩石学实验结果,我们计算得出马尼拉海沟俯冲带洋壳开始发生熔融的深度约为80km;俯冲板块上覆地幔岩石开始发生熔融的深度为约85km,而更多的地幔岩石熔融可能发生在更深处.

致谢

感谢加拿大地质调查局王克林教授在俯冲带热模拟研究过程中给予的指导,感谢中国科学院研究生院王多君教授和中国地震局地质研究所周永胜研究员在岩石实验知识方面给予的指导,感谢两位匿名审稿人提出的建设性意见.

参考文献
[1] Peacock S M, van Keken P E, Holloway S D, et al. Thermal structure of the Costa Rica-Nicaragua subduction zone. Phys. Earth Planet. Inter. , 2005, 149(1-2): 187-200. DOI:10.1016/j.pepi.2004.08.030
[2] Syracuse E M, van Keken P E, Abers G A. The global range of subduction zone thermal models. Phys. Earth Planet. Inter. , 2010, 183(1-2): 73-90. DOI:10.1016/j.pepi.2010.02.004
[3] Mase C W, Smith L. Effects of frictional heating on the thermal, hydrologic, and mechanical response of a fault. J. Geophys. Res. , 1987, 92(B7): 6249-6272. DOI:10.1029/JB092iB07p06249
[4] Rice J R. Heating and weakening of faults during earthquake slip. J. Geophys. Res. , 2006, 111(B05311): 1-29.
[5] Segall P, Rice J R. Does shear heating of pore fluid contribute to earthquake nucleation?. J. Geophys. Res. , 2006, 11: B09316. DOI:10.1029/2005JB004129.
[6] Arcay D, Tric E, Doin M P. Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics. Phys. Earth Planet. Inter. , 2005, 149(1-2): 133-153. DOI:10.1016/j.pepi.2004.08.020
[7] Rüpke L H, Morgan J P, Hort M, et al. Serpentine and the subduction zone water cycle. Earth Planet. Sci. Lett. , 2004, 223(1-2): 17-34. DOI:10.1016/j.epsl.2004.04.018
[8] Wang K L, Hyndman R D, Yamano M. Thermal regime of the Southwest Japan subduction zone: effects of age history of the subducting plate. Tectonophysics , 1995, 248(1-2): 53-59. DOI:10.1016/0040-1951(95)00028-L
[9] Abers G A, van Keken P E, Kneller E A, et al. The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow. Earth Planet. Sci. Lett. , 2006, 241(3-4): 387-397. DOI:10.1016/j.epsl.2005.11.055
[10] Cagnioncle A M, Parmentier E M, Elkins-Tanton L T. Effect of solid flow above a subducting slab on water distribution and melting at convergent plate boundaries. J. Geophys. Res. , 2007, 112: B09402. DOI:10.1029/2007JB004934.
[11] Conder J A, Wiens D A, Morris J. On the decompression melting structure at volcanic arcs and back-arc spreading centers. Geophys. Res. Lett. , 2002, 29(15): 1727. DOI:10.1029/2002GL015390.
[12] Peacock S M, Wang K L. Seismic consequences of warm versus cool subduction metamorphism: examples from southwest and northeast Japan. Science , 1999, 286(5441): 937-939. DOI:10.1126/science.286.5441.937
[13] Hayes D E, Lewis S D. A geophysical study of the manila trench, Luzon, Philippines 1. Crustal structure, gravity, and regional tectonic evolution. J. Geophys. Res. , 1984, 89(B11): 9171-9195. DOI:10.1029/JB089iB11p09171
[14] Lewis S D, Hayes D E. A geophysical study of the Manila trench, Luzon, Philippines 2. Fore Arc Basin structural and Stratigraphic evolution. J. Geophys. Res. , 1984, 89(B11): 9196-9214. DOI:10.1029/JB089iB11p09196
[15] Yang T F, Lee T, Chen C H, et al. A double island arc between Taiwan and Luzon: consequence of ridge subduction. Tectonophysics , 1996, 258(1-4): 85-101. DOI:10.1016/0040-1951(95)00180-8
[16] Bautista B C, Bautista M L P, Olike K, et al. A new insight on the geometry of subducting slabs in northern Luzon, Philippines. Tectonophysics , 2001, 339(3-4): 279-310. DOI:10.1016/S0040-1951(01)00120-2
[17] Hacker B R, Abers G A, Peacock S M. Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H2O contents. J. Geophys. Res. , 2003, 108(B1): 2029. DOI:10.1029/2001JB001127.
[18] Wada I, Wang K L, He J H, et al. Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization. J. Geophys. Res. , 2008, 113: B04402. DOI:10.1029/2007JB005190.
[19] Wyllie P J. Constraints imposed by experimental petrology on possible and impossible magma sources and products. Phil. Trans. R. Soc. Lond. A , 1984, 310(1514): 43.9-456. DOI:10.1098/rsta.1984.0003
[20] 李家彪, 金翔龙, 阮爱国, 等. 马尼拉海沟增生楔中段的挤入构造. 科学通报 , 2004, 49(12): 1279–1288. Li J B, Jin X L, Ruan A G, et al. Indentation tectonics in the accretionary wedge of middle Manila trench. Chinese Science Bulletin (in Chinese) , 2004, 49(12): 1279-1288.
[21] Li Z W, Xu Y, Hao T Y, et al. P wave velocity structure in the crust and upper mantle beneath Northeastern South China Sea and Surrounding regions. Earth Science Frontiers , 2009, 16(4): 252-260. DOI:10.1016/S1872-5791(08)60085-8
[22] Bird P. Initiation of intracontinental subduction in the Himalaya. J. Geophys. Res. , 1978, 83(B10): 4975-4987. DOI:10.1029/JB083iB10p04975
[23] Schubert G, Yuan D A. Shear heating instability in the earth's upper mantle. Tectonophysics , 1978, 50(2-3): 197-205. DOI:10.1016/0040-1951(78)90135-X
[24] Turcotte D L, Schubert G. Geodynamics (second edition). Cambridge: Cambridge University Press , 2002: 242-244.
[25] 张斯奇, 张怀, 石耀霖. 青藏高原深部热结构模拟的探讨. 中国科学院研究生院学报 , 2009, 26(3): 357–363. Zhang S Q, Zhang H, Shi Y L. Thermal simulation of deep Tibetan-Plateau by FD method. Journal of the Graduate School of the Chinese Academy of Science (in Chinese) , 2009, 26(3): 357-363.
[26] Davis E E, Lister C R B. Fundamentals of ridge crest topography. Earth Planet. Sci. Lett. , 1974, 21(4): 405-413. DOI:10.1016/0012-821X(74)90180-0
[27] Morrow C A, Moore D E, Lockner D A. The effect of mineral bond strength and adsorbed water on fault gouge frictional strength. Geophys. Res. Lett. , 2000, 27(6): 815-818. DOI:10.1029/1999GL008401
[28] Wang K L, Mulder T, Rogers G C, et al. Case for very low coupling stress on the Cascadia subduction fault. J. Geophys. Res. , 1995, 100(B7): 12907-12918. DOI:10.1029/95JB00516
[29] Currie C A, Wang K, Hyndman R D, et al. The thermal effects of steady-state slab-driven mantle flow above a subducting plate: the Cascadia subduction zone and backarc. Earth Planet. Sci. Lett. , 2004, 223(1-2): 35-48. DOI:10.1016/j.epsl.2004.04.020
[30] He L J, Wang K L, Xiong L P, et al. Heat flow and thermal history of the South China Sea. Phys. Earth Planet. Inter. , 2001, 126(3-4): 211-220. DOI:10.1016/S0031-9201(01)00256-4
[31] Hyndman R D, Wang K. Thermal constraints on the zone of major thrust earthquake failure: the Cascadia subduction zone. J. Geophys. Res. , 1993, 98(B2): 2039-2060. DOI:10.1029/92JB02279
[32] Zhou D, Yu H S, Xu H H, et al. Modeling of thermo-rheological structure of lithosphere under the foreland basin and mountain belt of Taiwan. Tectonophysics , 2003, 374(3-4): 115-134. DOI:10.1016/S0040-1951(03)00236-1
[33] Mouthereau F, Petit C. Rheology and strength of the Eurasian continental lithosphere in the foreland of the Taiwan collision belt: Constraints from seismicity, flexure, and structural styles. J. Geophys. Res. , 2003, 108(B11): 2512. DOI:10.1029/2002JB002098.
[34] 石耀霖, 张健. 活动海岭俯冲与岛弧火山活动的热模拟研究. 地球物理学报 , 1998, 41(2): 174–181. Shi Y L, Zhang J. Thermal modeling of active ridge subduction and arc volcanism. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) , 1998, 41(2): 174-181.
[35] Shi X B, Qiu X L, Xia K Y, et al. Characteristics of surface heat flow in the South China Sea. Journal of Asian Earth Sciences , 2003, 22(3): 265-277. DOI:10.1016/S1367-9120(03)00059-2
[36] Hsu S K, Yeh Y C, Doo W B, et al. New bathymetry and magnetic lineations identifications in the northernmost South China Sea and their tectonic implications. Marine Geophysical Researches , 2004, 25(1-2): 29-44. DOI:10.1007/s11001-005-0731-7