地震活动性显示深源地震多发生在板块俯冲区域，与俯冲带的热结构及物质组成密切相关.对于深源地震的形成机制，目前仍存在诸多争议，国内外学者的观点主要包括脱水脆化(Meade and Jealoz，1991; Dobson et al.，2002; Zhang et al.，2004)、剪切失稳(Griggs and Baker，1969; Ogawa，1987; 叶正仁等，1991)和反裂纹破裂机制(Green and Burnley，1990; Green et al.，1990; 宁杰远和臧绍先，1999)等.其中，反裂纹破裂机制利用高温高压实验将亚稳态相变与深源地震分布相对应，较好地解释了深源地震的特征，被认为是较为合理的一种深源地震机 制.实验表明，当岩石圈俯冲板块内部温度低于600~700 ℃时，橄榄石向尖晶石的相变过程受到抑制，低压相的橄榄石以亚稳态的形式存在于板块内部(Chien-Min and Burns，1976; Rubie and Ross，1994).由相变引发的密度及温度变化，促进了亚稳态橄榄石中反裂纹的生长与断裂的形成，进而诱发地震.

大量的俯冲带地震观测资料与热力学数值模拟研究表明与亚橄榄石相变有关的剪切失稳可能是诱发深源地震的原因(Green and Burnley，1989; Iidaka and Suetsugu，1992; Iidaka and Furukawa，1994; Kirby et al.，1996; Yoshioka et al.，1997; Chen and Brudzinski，2003; Carminati et al.，2005; Kaneshima et al.，2007). 亚稳态橄榄石的存在一定程度上决定了深源地震能够发生的区域，对地震的 活动性有一定的指示作用.Iidaka和Suetsugu(1992)利用深源地震的走时数据研究太平洋板块内亚稳态橄榄石的存在，结果表明由亚稳态模型计算出的走时残差与观测结果更为符合；Carminati等(2005)计算了Calabria俯冲带的温度分布，认为俯冲带热模拟所预测的亚稳态橄榄石存在范围与深源地震分布相符；Kaneshima等(2007)通过对马里亚纳板块深源地震后续P波波形的分析，发现地震波低速异常区域与亚稳态橄榄石楔状带的形态十分接近，以此来解释震源位置的分布能得到较为理想的结果.

由于亚稳态橄榄石的形成受到俯冲带温度的控制，因此通过俯冲带热结构的模拟定量地讨论俯冲带各参数的影响，有助于深入研究亚稳态橄榄石与深源地震活动性的关系，探讨深源地震的形成机制.目前俯冲带热模拟研究表明，俯冲带的年龄、俯冲角度和速度等参数会影响俯冲带内亚稳态橄榄石的分布(臧绍先和宁杰远，1994，2001; Negredo et al.，2004，2007).但关于俯冲带形态对亚稳态橄榄石及其相变边界的影响研究较少.通过俯冲带热结构的形态确定亚稳态橄榄石相变边界是目前最通用的研究方法.

俯冲带的形态和地幔对流模式有一定的关联性.目前地幔中存在的对流模式主要分为全地幔对流模式和双层地幔对流模式.O′Connell(1977)最早提出在整个地幔范围内，地幔对流从核幔边界上升至岩石圈底部，形成全地幔大尺度对流.在这种模式下，俯冲带能够穿透660 km相变界面到达下地幔(Van der Hilst et al.，1997; Monnereau and Quéré，2001). 而双层地幔对流模式主张上地幔和下地幔各自形成独立的对流系统，上、下地幔物质化学成分的差异会阻碍地幔对流，使俯冲带停留在660 km相变界面附近(Jeanloz and Richter，1979; Allègre et al.，1987).

本文利用有限差分方法计算日本海沟俯冲带热结构、浮力及P波速度异常分布，结合绝热压缩、相变潜热和剪切生热等多热源的综合影响，比较两种地幔对流模式下亚稳态橄榄石的存在范围，探讨其与日本海沟俯冲带深源地震的关系. 2 日本海地区地震活动特征

日本海沟位于太平洋西北部，是太平洋板块向欧亚板块俯冲的汇聚带.太平洋板块沿着日本海沟向西的俯冲始于26 Ma前，以平均9 cm/a的速度向日本岛弧下运动(Creager and Jordan，1986; 臧绍先和宁杰远，1996; 何建坤和刘福田，1998; Peacock and Wang，1999; 孙文斌和和跃时，2004; Huang and Zhao，2006)，由此形成的日本海沟俯冲带是西太平洋地震活动频度最高、强度最大的地区.由1960—2008年国际地震中心(ISC)提供的地震分布(图 1)可知，日本海地区浅源地震活动主要集中在日本岛弧附近，地震震源深度自日本海沟向西逐渐加深.日本海沟俯冲带深源地震活动高峰在350~500 km深度范围内，500~600 km内地震活动性明显降低，发生深源地震的最大深度达到605 km(如图 2所示).

图 1

Fig. 1

图 1 日本海及中国东北地区Mb≥4.0地震分布图. 区间范围34°N—45°N，107°E—147°E; AA′为模型剖面所在位置; 箭头指向为太平洋板块运动方向Fig. 1 Earthquakes distribution in Japan Sea and northeast China with magnitudes greater than 4.0 in the region of 34°N—45°N and 107°E—147°E. Line AA′ shows the position of model profile in this study. The arrow represents the movement direction of the Pacific plate

图 2

Fig. 2

图 2 日本海地区深源地震统计Fig. 2 Statistic of deep-focus earthquakes in Japan Sea

利用反裂纹相变破裂机制来解释深源地震活动性的一个重要条件是需要有亚稳态相存在的证据.那么，日本海沟俯冲带上的深源地震是如何形成的？能否利用反裂纹相变破裂机制来解释日本海沟俯冲带上的深源地震活动性？针对这一问题，本文将通过俯冲带热结构的数值模拟加以验证. 3 模型 3.1 控制方程

俯冲带热结构影响着橄榄石向高压相尖晶石的相变过程、俯冲带内物质密度和浮力分布，进而影响到俯冲板块的动力学过程(臧绍先等，1994; 臧绍先和宁杰远，2001; Bina and Kawakatsu，2010).俯冲带的热模拟一直是研究俯冲带动力学演化，验证俯冲带上地震机制与地震观测结果的重要手段(臧绍先和宁杰远，1996).以往的研究通常采用Minear和Toksöz提出的准动力学计算方案计算俯冲带的热结构(臧绍先和宁杰远，1994; Kirby et al.，1996; Stein and Stein，1996; Bina，1997).Negredo等(2004)在此基础上加入了相变过程的自洽模拟，得到更加精确的结果.

本文基于Negredo等(2004)提出的TEMSPOL模型，利用有限差分方法计算全地幔对流模式和双层地幔对流模式下日本海沟俯冲带的热结构，所采用的控制方程为含平流项的二维对流热传导方程：

计算过程中，考虑橄榄石向尖晶石相变产生的密度变化及由温度引起密度变化的浮力效应，则密度表示为

同时，热源项Q包含因俯冲带在俯冲过程中受挤压做功而产生的热量，即绝热压缩生热(Toksöz et al.，1971)

俯冲带与周围介质摩擦而产生的剪切热与剪应力密切相关.Minear和Toksöz(1970)给出的剪切生热为

橄榄石向尖晶石相变过程中产生的相变热也是影响俯冲带热结构的重要热源之一.其计算公式为

我们选取日本海沟以西39°N，110°E—145°E之间的剖面AA′(如图 1所示)作为计算区域，计算130 Ma，95 km厚的太平洋板块，以29°倾角(臧绍先和宁杰远，1996; Huang and Zhao，2006)、9 cm/a(Peacock and Wang，1999; 孙文斌和和跃时，2004)的速度俯冲至日本岛弧下方产生的热结构.考虑全地幔对流和双层地幔对流模式下日本海沟俯冲带在660 km相变界面上不同的形态特征.图 3为 模型示意图.初始温度根据Stein C A和Stein S(1992)给出的GDH1温度模型计算得出.边界条件在地表温度取0 ℃，模型底部为给定热流边界，两侧边界采用绝热温度约束.基于现有的研究成果(Ringwood and Major，1970; Honda，1985; 臧绍先等，1993，1994; 臧绍先和宁杰远，2001; Bina and Kawakatsu，2010)，计算模型中所使用的参数如表 1所示.

图 3

Fig. 3

图 3 模型示意图.全地幔对流模式: 俯冲角度不变; 双层地幔对流模式：俯冲角度发生水平偏转Fig. 3 Schematic diagram of the model. Whole mantle convection: no subducted angle change; layered mantle convection: horizontal deflection of subducted angle

表 1(Table 1)

表 1 计算模型中所使用的参数 Table 1 Parameters used for the model 参数 数值 参考密度ρ0 3400 kg·m-3 相变密度变化Δρ 181 kg·m-3 比热c 1300 J·K-1·kg-1 热传导系数k 3.0 W·m-1·K-1→变化 热膨胀系数α 3.6×10-5 K-1→变化

表 1 计算模型中所使用的参数 Table 1 Parameters used for the model

综合考虑绝热压缩生热、剪切生热和相变潜热的影响，图 4给出了不同地幔对流模式下日本海沟俯冲带热结构.根据Schmeling等(1999)给出的动力学相变图(图 5)，以600 ℃作为判定亚稳态橄榄石存在范围的极限温度值，计算了高压相尖晶石占1%和99%的深度范围，以确定橄榄石的相变边界.

图 4

Fig. 4

图 4 全地幔对流模式(a)和双层地幔对流模式(b)下日本海沟俯冲带热结构(单位℃).上下两条虚线分别表示尖晶石相占1%和99%的相变界面，其间为相变区Fig. 4 Thermal structures of Japan Trench subduction zone under the whole mantle convection(a) and layered mantle convection(b). Position of 1% and 99% olivine- spinel transformation contours are shown by dashed lines

图 5

Fig. 5

图 5 橄榄石向尖晶石相变的动力学相变图 实线和虚线分别表示1%和99%的尖晶石分布(修改自Schmeling等(1999)).Fig. 5 Kinetic phase diagram of olivine-spinel transition The solid and dashed lines represent 1% and 99% spinel fraction respectively(After reference(Schmeling et al.，1999)).

由图 4可以看出，当太平洋板块俯冲进入地幔时，板块内部的温度明显低于周围地幔温度，橄榄石相变深度随板块俯冲而加深.对于全地幔对流模式(图 4a)和双层地幔对流模式(图 4b)，橄榄石相变深度约为720 km、590 km，达到了亚稳态橄榄石的存在条件.当俯冲带穿越相变区域时，由于410 km处的克拉伯龙斜率为正，低温的俯冲带通过相变界面时，相变深度会相应抬升，且当温度小于600 ℃时，会在板块内部形成亚稳态橄榄石，相变边界表现为两侧升高中间降低的形态.亚稳态橄榄石的存在范围呈现楔形带状分布，随着深度的加深逐渐变薄. 图 4b中在双层地幔对流模式下，俯冲板块到达660 km 相变界面后发生水平偏转，在相变界面上形成一个约1200 km的水平低温区域，与Fukao等(1992)所获得的成果相近.在板块俯冲26 Ma后，橄榄石的相变边界也随之向水平方向延伸了约400 km. 4.2 日本海沟俯冲带浮力分布

俯冲带所受浮力与密度密切相关.俯冲带的岩石密度与地幔围岩的密度差越大，所受的浮力作用也越强.我们定义负浮力指向下，利用以下公式计算单位体积俯冲带的浮力(臧绍先和宁杰远，1994，2001; Bina and Kawakatsu，2010)

图 6

Fig. 6

图 6 全地幔对流模式(a)和双层地幔对流模式(b)下日本海沟俯冲带浮力(单位103 N·m-3)Fig. 6 Buoyancies of Japan Trench subduction zone under the whole mantle convection(a) and layered mantle convection(b)

Birch(1952)最早提出了P波速度与介质密度的线性关系.根据已得的俯冲带密度变化，即可计算出日本海沟俯冲带P波速度异常分布(Birch，1952; 宁杰远和臧绍先，2001)

图 7

Fig. 7

图 7 日本海沟俯冲带P波速度异常(单位%). 全地幔对流模式(a)和双层地幔对流模式(b)Fig. 7 P-wave velocity anomalies of Japan Trench subduction zone under the whole mantle convection(a) and layered mantle convection(b)

地震层析成像可以在一定程度上揭示出俯冲带的几何形态特征.从日本海地区的层析成像结果来看(Huang and Zhao，2006)，俯冲到上地幔的西太平洋板块表现为P波高速异常带，板块俯冲至地幔转换带内发生滞留并平卧于660 km相变界面之上(如图 8所示).我们的模拟结果表明双层地幔模式下的日本海沟俯冲带P波速度异常结果与层析成像结果有较好的一致性.然而，与层析成像结果不同的是，由于俯冲带内部存在亚稳态橄榄石，因而在模拟结果中会在俯冲带内部形成一个楔形的低速区，这与Bina和Kawakatsu(2010)的研究结果相同.Iidaka和Suetsugu(1992)及Kaneshima等(2007)利用区域台网观测资料研究伊豆—小笠原俯冲带和马里亚纳俯冲带的波速异常时也验证了这一结果.地震层析成像的结果并没有反映出亚稳态橄榄石的低速异常可能与其分辨率有关(宁杰远和臧绍先，2001).

图 8

Fig. 8

图 8 沿39°N剖面的层析成像结果(修改自Huang和Zhao(2006))Fig. 8 Cross section of seismic tomography along latitude 39°N(After reference(Huang and Zhao，2006))

亚稳态橄榄石的存在范围决定了深源地震能够 发生的最大深度(Green and Burnley，1989; Negredo et al.，2004). 但不同模型给出的亚稳态橄榄石存在范围相差较大.Jiang和Zhao(2011)利用正演模拟法研究日本海地区太平洋俯冲带内的橄榄石亚稳态楔结构，认为亚稳态橄榄石楔的顶端和底端深度分别为400 km和570 km，大部分深源地震发生在亚稳态楔内部.而我们的模拟结果显示亚稳态橄榄石的存在范围在300~600 km之间，与日本海地区深源地震深度分布更为相符，因此我们推断日本海沟俯冲带上深源地震的发生可能是由橄榄石的亚稳态相变引起的，与亚稳态橄榄石的存在有关，受反裂纹机制的控制. 5.2 热传导系数k和热膨胀系数α的影响

热传导系数k和热膨胀系数α是制约俯冲带热结构的两个重要参数.由于地球内部的不均匀性，热参数也应是随深度变化的.但在目前的数值模拟研究中，通常将热传导系数和热膨胀系数取为常数作 近似简化处理.然而对于这两个参数的具体取值，目前仍存 在许多争议.本文基于Xu等(2004)和Bouhifd等(1996)的实验结果(如图 9所示)，选取5种不同的热参数模型(如表 2所示)，计算热传导系数和热膨胀系数对亚稳态橄榄石分布的影响.热传导系数k最小值 取1.8 W·m^-1·K^-1，最大值取4.2 W·m^-1·K^-1； 热膨胀系数α最小值取2.7×10^-5 K^-1，最大值取4.5×10^-5 K^-1.

图 9

Fig. 9

图 9 随温度变化的热传导系数(实线)和热膨胀系数(虚线) Xu等(2004)研究表明，橄榄石的热传导系数随温度升高呈指数递减；Bouhifd等(1996)对镁橄榄石热膨胀系数的实验结果表明 热膨胀系数随温度升高线性增加.Fig. 9 The variation of the thermal conductivity(solid line) and thermal expansion(dashed line)with temperature The result of Xu et al.(2004)shows the thermal conductivity coefficient of olivine decreases exponentially with the temperature，while the experimental data of Bouhifd et al.(1996)indicate a linear increase of the thermal expansion coefficient of forsterite with temperature increment.

表2(Table 2)

表2 不同热参数模型所使用的参数 Table 2 Parameters used for different thermal models Model k/ (W·m-1·K-1) α/ (10-5 K-1) c/ (J·K-1·kg-1) ρ/ (kg·m-3) 1 3.0 3.6 1300 3400 2 1.8 3.6 1300 3400 3 4.2 3.6 1300 3400 4 3.0 2.7 1300 3400 5 3.0 4.5 1300 3400

表2 不同热参数模型所使用的参数 Table 2 Parameters used for different thermal models

图 10显示了5种不同热参数模型对亚稳态橄榄石的影响.模型1中热传导系数k和热膨胀系数α分别取均值3.0 W·m^-1·K^-1和3.6×10^-5 K^-1时，亚稳态橄榄石楔的范围约为38841 km²，其最大深度达到595 km.保持α值不变，将k减小至1.8 W·m^-1·K^-1(模型2)，亚稳态橄榄石楔的范围扩大至99115 km²，是均值状态下的1.55倍；当k增大至4.2 W·m^-1·K^-1(模型3)时，其范围减小了64%，约为14106 km².同时，保持k值不变，α分别取2.7×10^-5 K^-1(模型4)和4.5×10^-5 K^-1(模型5)时，亚稳态橄榄石的范围分别略大于或小于模型1的结果，其变化范围均在15%左右.橄榄石相变深度与亚稳态橄榄石的存在范围有较为一致的变化趋势，其中k值的变化对相变深度影响较大，由最浅520 km到最深620 km，其余深度均在600 km左右.

图 10

Fig. 10

图 10 亚稳态橄榄石存在范围及相变深度Fig. 10 Area and depth of metastable olivine wedge

本文在俯冲带热模拟中，考虑俯冲带绝热压缩生热、橄榄石向尖晶石相变潜热及剪切生热这三种热源对俯冲带热结构的影响.臧绍先等(1993)曾经基于准动力学的计算方案，利用有限元法对日本海沟俯冲带的热结构及热源的影响进行了一些研究，认为橄榄石向尖晶石的相变生热是控制400 km深度以下温度分布的主要热源因素.本文的研究结果 也表明，相变潜热和剪切生热对俯冲带的温度及亚稳态橄榄石的分布有较大影响.若忽略热源，亚稳态橄榄石的存在范围约为48277 km²，与考虑相变潜热和剪切生热的情况相比，橄榄石相变边界向前延伸了约240 km，面积增大了24%(如图 11所示).由此可见，忽略这两种热源会造成模型所预测的亚稳态橄榄石楔范围偏大，导致俯冲带上深源地震活动性的错误估计.

图 11

Fig. 11

图 11 热源对亚稳态橄榄石楔的影响Fig. 11 The influence of heat sources on metastable olivine wedge

俯冲带所受负浮力是板块运动的主要动力来源，与俯冲带的形态密切相关(臧绍先和宁杰远，2001; Mahatsente and Ranalli，2004).臧绍先和宁杰远(2001)与Bina和Kawakatsu(2010)的研究认为双层地幔对流模式下的板块形变特征与亚稳态橄榄石相变作用有关.本文对日本海沟俯冲带负浮力的计算结果也表明，当较冷的板块进入高温高压的地幔时，温度和压力的变化使板块内部物质密度增大，俯冲带上产生的负浮力在400 km深度附近达到最大值，负浮力促进了板块向下运动，使俯冲带的形变在410 km相变界面处有向下弯曲的趋势.而在400~600 km深度之间，橄榄石向尖晶石相变过程中形成的亚稳态橄榄石会使板块内部物质密度逐渐减小，甚至小于周围地幔密度，板块内部分区域所受浮力变为正值，俯冲带向下进入下地幔的过程受到阻碍，俯冲带的形变在660 km相变界面处有向上弯曲的趋势.这在一定程度上可以解释日本海沟俯冲带的水平偏转并停留于660 km相变界面之上的现象. 6 结论

本文通过对日本海沟俯冲带热结构、浮力及P波速度异常的计算及结果讨论，可以得到如下结论：

(1)不同的地幔对流模式决定了亚稳态橄榄石的分布特征.双层地幔对流模式下日本海沟俯冲带上P波速度异常与层析成像结果较为相符，由此得到的亚稳态橄榄石存在范围与日本海沟俯冲带深源地震的活动有较好的相关性.

(2)热传导系数对亚稳态橄榄石的影响较大，且热传导系数和热膨胀系数越小，亚稳态橄榄石的存在范围越大，橄榄石相变深度越深.

(3)相变潜热和剪切生热是影响日本海沟俯冲带热结构的重要热源因素，忽略这两种热源会造成亚稳态橄榄石分布及深源地震活动性的错误估计.

(4)亚稳态橄榄石的存在与俯冲带形态有关，俯冲带所受负浮力在400 km深处最大，随深度增加 而逐渐减小，甚至产生正浮力，阻碍俯冲带穿透660 km 相变界面.