马尼拉俯冲带热结构数值模拟与地震意义

陈爱华， 许鹤华

（1. 国土资源部海底矿产资源重点实验室 广州海洋地质调查局， 广州 510760； 2. 中国科学院边缘海地质重点实验室 中国科学院南海海洋研究所， 广州 510301）

0 引言

通过分析俯冲板块年龄和热流值， 马尼拉俯冲带属于热俯冲[5-11]。 陈爱华等[12]尝试利用热俯冲探讨马尼拉俯冲带的地震分布规律。 但马尼拉俯冲带的俯冲速度和俯冲角度比其它热俯冲大得多，会对俯冲洋壳的热结构分布产生很大影响[13]。 马尼拉俯冲带经历了复杂的地质演化过程， 俯冲板块深度较深， 同时存在陆壳、 洋陆过渡带和洋壳，且不同区域的地质条件相差较大[14]。 高翔等[15]选取马尼拉俯冲带南、 北两条测线， 根据两处俯冲板块的角度、年龄和运动特征的差异， 分析了摩擦热和剪切热对俯冲带热结构的影响， 讨论俯冲板块发生脱水和部分熔融的位置。 由于研究的问题不同， 建立的模型和考虑的因素也不同， 必然会导致模拟结果产生差异。 近年来， 通过对马尼拉俯冲带以及吕宋岛弧地球物理测量， 获得了马尼拉俯冲带各构造单元莫霍面的形态、 走向、 梯度和地壳性质， 为本文研究马尼拉俯冲带的热结构提供了最新的地球物理资料[16]。 本文综合马尼拉俯冲带的地震、 古地磁、 地热与层析成像等多种最新地球物理资料， 选取3 条典型剖面（位置见图1）， 构建地质模型和二维数值模型， 根据剖面板块俯冲角度、 速度等， 应用有限元方法求解非线性瞬态热方程， 得到马尼拉俯冲带的热结构特征，为进一步探讨马尼拉俯冲带地震成因机制提供热力学依据。

1 区域地质背景

马尼拉俯冲带是菲律宾海板块向南海板块仰冲所形成的汇聚边界， 南起民都洛岛西南陆架的海底峡谷， 北部延伸至北纬21°15′， 并且逐渐与台湾造山带西部的变形相连接， 空间上呈南北向弧形展布， 断面呈不对称的“V”形， 整体表现为向西“凸”形， 形态为向东倾的贝尼奥夫带（图1）[17-18]。

图1 马尼拉俯冲带地形图Fig.1 Topographies maps of Manila subduction zone

马尼拉俯冲带主要受近东西向的挤压， 大多数地壳内地震的主压应力轴接近E-W 方向； 菲律宾断层以北以NW 向挤压应力为主， 显示北部以挤压逆冲为特征； 菲律宾海槽及海沟地区， 主压应力方向为NWW 到NEE 之间； 南部较为复杂，存在NW、 NE 和近NS 向挤压应力， NW 向的挤压应力主要来自菲律宾板块继续向西北的运动， 以顺时针旋转为特征， 且存在左旋走滑断层[19-20]。

根据南海海底磁异常对比分析， 多数学者认为， 南海中央次海盆具有东西向对称分布的磁条带(11～5 d)[6]。 大陆裂谷及早期海底扩张发生于32～26 Ma， 期间伸展应力场由SE 向转为SSE 向；24 Ma 后， 转为NW-SE 向扩张， 最后终止于15 Ma[8,21]。 在中新世， 菲律宾海板块沿马尼拉海沟以7.0 cm/a 的速度、 NW55°方向斜向仰冲[22]。 Bautista等[23]的研究也表明， 菲律宾板块在吕宋岛北端以约7 cm/a， 向南逐渐增大， 在棉兰老岛的东南部最大， 达到9 cm/a， 欧亚板块则以1 cm/a 的速度沿相近方向做同向运动， 由于菲律宾板块速度比欧亚板块运动速度快， 二者表现出强烈会聚。 随着俯冲作用继续和菲律宾海微板块向西北方向移动， 吕宋岛弧和欧亚板块发生碰撞[22]。 菲律宾海板块持续地向西北运动， 碰撞作用加剧并向南扩展， 马尼拉海沟的北端逐渐进入碰撞造山阶段， 形成了非火山弧（增生楔）-弧前盆地（北吕宋海槽和西吕宋海槽）-火山弧（吕宋火山弧）构造组合[24，25]。

2 计算模型

2.1 剖面特征与初始条件

马尼拉俯冲带是南海东部的俯冲边界， 北部受台湾岛碰撞带的影响， 南部由于巴拉望地块和菲律宾板块的作用， 使得地质构造复杂。

AA′剖面位于台湾岛南部附近， 沿20°N 的纵波速度结构剖面显示， 东沙隆起（116°～118°E）处地壳厚度约32 km， 岩石圈厚度约95 km； 往东至120°10′E 处， 地壳减薄至20 km， 此段西部为陆坡区， 东部南海东北部海盆是洋壳[26]。 由此可以推断， 剖面AA′俯冲板块为陆壳， 地壳厚度取为25 km， 岩石圈厚度取95 km。

BB′剖面处于洋陆过渡带， 往北逐渐为陆壳，往南则为南海中央海盆。 沿18°N 的P 波速度结构剖面表明， 海南岛之下的岩石圈厚度在80～85 km左右， 从这里开始至南海中央海盆东部边缘（118°E）， 岩石圈厚90～95 km。 在南海西北海盆至中央海盆之下（114°～119°E）， 地壳厚度减薄到5～6 km， 莫霍面深度10 km， 在马尼拉海沟（119°30′E）之下， 莫霍面深度从10 km 向40 km 下降， 但岩石圈厚度未见变化[26]。

尾矿资源化和减量化是指将尾矿直接作为资源进行回收利用或者对尾矿进行再生利用，其方法包括尾矿再选和尾矿制备材料等。尾矿的资源化和减量化利用不仅可提高我国矿业开发综合利用水平，同时还可减少尾矿和废水等矿业固体废弃物的排放，减少矿业开发过程对环境的影响。

CC′剖面总体属于中央海盆区域， 也处于古俯冲带的洋中脊附近， 俯冲板块为洋壳。 根据上地幔三维S 波速度结构结果：

（1）沿北纬16°的速度剖面上， 在东经114°附近， 软流层顶界面深度最浅， 在60 km 左右， 向东逐步变深， 在东经120°最深， 达到70 km 左右，然后又变浅， 达到60 km 左右。

（2）东经115°的速度剖面上， 在北纬16°附近， 出现一个隆起， 此处软流层顶界面深度最浅，在60 km 左右， 北纬16°以北部分， 向北逐渐变深， 最深达到100 km 左右， 北纬16°以南部分，向南先是变深， 在65 km 左右， 然后又逐步变浅，最浅达到60 km 左右[27]。

南海海盆的岩石圈厚度为65 km 左右， 最浅可达到60 km， 向东到菲律宾海盆平均地壳厚度25.7 km， 平均岩石圈厚度80 km[28]。 BB′和CC′两条剖面处的磁条带异常对比大概为12（约35 Ma）和6（约20 Ma）， 根据半空间冷却模型进行计算，两处剖面的岩石圈厚度约为77 km 和58 km [7,29]。 海沟两侧的台湾-吕宋岛地区的地壳厚度约为34 km， 岩石圈厚度为70～80 km[15，30]。 根据层析成像结果可以推测， 三条剖面处的板块俯冲深度均超过300 km[31-32]。

根据南海及邻区莫霍面深度反演计算， 南海板块在BB′剖面处地壳中的俯冲角度很平缓， 大约3°， 中央海盆莫霍面深度较厚， 约为19～22 km，马尼拉俯冲带处地壳厚度为15～20 km； 由于CC′剖面处靠近古俯冲洋中脊， 中央海盆莫霍面深度减薄至10 km 左右， 马尼拉俯冲带地壳厚度约为10～15 km， 地壳中的俯冲角度不大， 约4°左右[16]。 B autista 等[23]根据地震分布和震源机制解， 推断南海俯冲板片表面形态， 马尼拉俯冲带南海俯冲板片的形态复杂， 底部的俯冲角度变化非常大， AA′剖面俯冲板块在浅部的俯冲角度约5°， 深部约60°，至18°逐渐变缓， 向南俯冲角度又逐渐变陡， 至最南端的民都洛岛处， 俯冲角度近乎垂直。

模型中选取的尺度较大， 洋壳沉积层厚度较薄， 且本文的模拟结果主要针对中深源地震， 为方便建模和计算， 忽略洋壳沉积层。 在BB′剖面菲律宾以下的南海板块俯冲角度达到45°； CC′剖面南部剖面深部约65°， 综合层析成像和地震剖面资料， 本文中CC′剖面取为60°。

综合上述地质资料及分析， 地质模型参数选取如表1 所示。

2.2 模型参数

根据马尼拉俯冲带的地质资料， 可以建立理论数值模型（图2）。 为了研究俯冲板块的热结构及减少边界影响， 本文中设定计算区域马尼拉海沟横向宽度为1 000 km， 模拟底边界深度为600 km，板片俯冲深度为300 km。 陆壳模型的初始温度由稳态热传导方程给出， 洋壳俯冲模型的初始温度根据菲律宾板块的年龄及厚度， 通过半空间冷却模型的计算方法给出。 上边界温度取为0℃， 底边界处温度取为1 600 ℃， 海洋板块岩石圈边界和上覆板块岩石圈边界均采用半空间冷却模型计算出的温度分布作为边界条件， 软流圈地幔左右边界均为自由边界。 模型各参数如下表2。

表1 地质模型参数Table 1 Parameters of geology model

图2 数值模型示意图Fig.2 Schematic diagrams of numerical model

表2 数值模型参数Table 2 Parameters of numerical model

根据臧绍先等[33]的研究结果， 俯冲洋壳的剪切生热以及地幔中的放射性物质的产热量A 相对其他热源项要小很多， 忽略它们作为一种近似是合理的。 由于AA′剖面俯冲板块为陆壳， 考虑到在同一圈层差异很小， 也为方便计算， 放射性生热率取为1*10-6W/m3。

3 模拟结果

对建立的两个数值模型赋予相同的边界条件，为验证本模拟结果， 采用热流作为检验标准。 AA′剖面属陆壳， 地壳厚度较厚， 海沟西侧热流值为66 mW/m2， 最低值约为45 mW/m2， 海沟东侧约为72 mW/m2， 陆壳本身热流值较低， 加上角落流、岩浆活动和火山活动等因素的影响， 上覆板块热流值会相对上升， 使得海沟东侧热流值高于西侧；BB′剖面中， 在俯冲较浅位置， 热流值为68 mW/m2， 海沟最低处热流值为45 mW/m2， 在菲律宾板块， 热流值约为65.8 mW/m2； CC′剖面由于岩石圈厚度较薄， 浅部达到70 mW/m2， 而海沟最低热流值约为20 mW/m2， 与BB′剖面处的俯冲角度和俯冲速度有关， 菲律宾板块热流值与BB′剖面基本一致。 模型中的热流值和热流趋势与南海热流特征一致[9,10,34]， 因此， 此模拟结果可以较好地反应马尼拉俯冲带的热结构。

图3 AA′、BB′和CC′剖面计算热结构图Fig.3 Calculated thermal structure for AA′、BB′ and CC′cross-section

根据图3 可知， 在俯冲的浅部， 由于俯冲角度平缓， 俯冲板块的温度主要受热传导作用控制，基本与俯冲板块形状保持一致， 随俯冲深度加深而增加。 当深度加深， 俯冲角度逐渐增大， 纵向的分速度也逐渐增大， 较冷的俯冲板块运动至较深的深度。 至软流圈边界时， 开始受软流圈对流影响， 俯冲板块温度变化明显， AA′剖面俯冲板块顶部约为570 ℃， BB′剖面洋壳顶部温度比AA′高，约为680 ℃， 可能是由于BB′剖面厚度较AA′剖面薄， 相对更 “热”。 CC′剖面的俯冲角度和俯冲速度均较BB′剖面大， 温度也较低， 洋壳顶部仅有约450 ℃。 当洋壳至软流圈后， 洋壳温度快速升高，在90 km 处， AA′剖面板块顶部温度达到730 ℃，底部增加较缓， 约为650 ℃； BB′剖面表面温度约为805 ℃， 底部温度约为673 ℃； CC′剖面表面温度约为885 ℃， 底部温度比BB′剖面低的多， 约为210 ℃。 由于CC′剖面的俯冲速度和角度均比AA′和BB′剖面大， 使得软流圈对洋壳内部“加热”作用比AA′和BB′剖面弱的多， 保持较冷的状态到更深的深度。 经过软流圈的加热后， 洋壳表面和底部温度均有较大升高， 在160 km 时， AA′剖面温度约875 ℃， 底部温度约为671 ℃； BB′剖面洋壳顶部温度升至950 ℃， 洋壳底部温度也较高， 达到840℃； CC′剖面洋壳顶部温度更高， 达到1 050 ℃， 洋壳底部温度也升高到500 ℃。 三条剖面的热结构对比也表明， 俯冲带热结构主要受俯冲角度、 俯冲速度和俯冲板块本身地质条件如软流圈深度、厚度、年龄等因素影响。

4 地震意义

Peacock 等[35]和Kirby 等[36]将热流值＞75 mW/m2，且年龄＜15～25 Ma 的洋壳俯冲称为热俯冲， 其它热流较低、年龄较老的俯冲为冷俯冲。 BB′和CC′剖面的地温梯度值比日本西南部热俯冲的地温梯度还高， 因此， 本模拟结果也进一步证明BB′和CC′剖面属于热俯冲[35]。 根据热俯冲的研究， 地震成因机制主要有以下两种观点：

（1）俯冲深度比较浅[37]。

（2）含水相脱水和榴辉岩的形成， 在俯冲板块与软流圈接触之前停止[36]。

一般认为， 浅源地震是由脆性断裂导致的，大洋俯冲带深部地震基本由岩石相变引起[38-39]。 当俯冲洋壳开始俯冲时， 上覆岩石压力和温度增大，大多数的孔隙被挤压消失， 排出大量孔隙水。 在俯冲的早期阶段， 随着变质程度从沸石相的递增，相应的主要含水矿物为沸石、 绿纤石等， 岩石中的含水量可达到8%～9%[40]。 随着进一步俯冲， 洋壳在约20 km 的深度进入绿帘石蓝片岩相区域，在50 km 的深度进入角闪石榴辉岩相[41]。 根据本文模拟结果， 当俯冲至软流圈边界时， 开始受软流圈的影响， 俯冲板块温度迅速升高， 自由水降低200～300 ℃反应温度， 使俯冲洋壳大量脱水， 孔隙压力继续增加， 上覆岩石压力减小， 岩石脱水脆化发生破裂， 生成的水促进周围无水辉长岩和玄武岩相变， 形成大量榴辉岩， 伴随着密集的地震发生[36,41-43]。 榴辉岩的形成也使洋壳密度增大， 板块应力改变， 负浮力增加。 由于BB′剖面的俯冲角度和俯冲速度比CC′剖面小， 使得BB′剖面发生地震的深度更浅。 俯冲洋壳顶部温度较高， 在南海洋壳到达软流圈之前， 含水相脱水和榴辉岩的形成基本停止， 使软流圈之上的洋壳发生密集的地震活动。 洋壳底部温度更低， 需要在更深的深度才能形成榴辉岩， 地震活动也持续到更大的深度，当洋壳全部进入干榴辉岩相时， 地震活动也基本停止， 这种推测与马尼拉俯冲带地震震源深度分布结果吻合。 但进入榴辉岩相的具体深度， 则需要结合实验岩石学结果， 详细分析俯冲洋壳的相变。

5 结论

通过本次研究， 得出以下几点认识和结论：

（1）模拟计算结果表明， 俯冲带热结构主要受俯冲角度、 俯冲速度和俯冲板块本身地质条件等因素影响。

（2）根据模拟结果的地温梯度数据， 进一步证明BB′剖面和CC′剖面属于热俯冲。

（3）当洋壳俯冲至软流圈边界时， 受软流圈的影响， 俯冲板块温度迅速升高， 容易形成大量榴辉岩， 产生地震活动。

致谢：感谢中国科学院南海海洋研究所张佳政博士给予的帮助。

[1] 张健， 汪集旸. 南海北部大陆边缘深部地热特征[J]. 科学通报， 2000， 45 (10)： 1095-1100.

[2] 宋茂双， 谢洪森， 徐有生. 俯冲带的热结构研究进展[J].地质地球化学， 1996（2）： 18～22.

[3] 臧绍先， 宁杰远， 景志成. 俯冲带流变性质的研究[J].中国科学：D 辑， 2001， 31（09）： 705-711.

[4] 周蒂， 王万银， 庞雄， 等. 地球物理资料所揭示的南海东北部中生代俯冲增生带[J]. 中国科学： D 辑, 2006,36 (03)： 209-218.

[5] Taylor B, Hayes D E. The tectonic evolution of the South China Basin[J]. The tectonic and geologic evolution of Southeast Asian seas and islands， 1980， 1： 89-104.

[6] Taylor B, Hayes D. Origin and history of the South China Sea basin [J]. The tectonic and geologic evolution of Southeast Asian seas and islands， Part， 1983， 2： 23-56.

[7] Briais A， Patriat P， Tapponnier P. Updated interpretation of magnetic anomalies and seafloor spreading stages in the South China Sea： implications for the Tertiary tectonics of Southeast Asia [J]. Journal of Geophysical Research，1993， 98 (B4)： 6299-6328.

[8] 孙珍， 孙龙涛， 周蒂， 等. 南海岩石圈破裂方式与扩张过程的三维物理模拟[J]. 地球科学(中国地质大学学报)， 2009 (03)： 435-447.

[9] He L， Wang K， Xiong L， et al. Heat flow and thermal history of the South China Sea[J]. Physics of the Earth and Planetary Interiors， 2001， 126 (3-4)： 211-220.

[10] Shi X， Qiu X， Xia K， et al. Characteristics of surface heat flow in the South China Sea[J]. Journal of Asian Earth Sciences, 2003， 22（3）： 265-277.

[11] 刘绍文， 施小斌， 王良书， 等. 南海成因机制及北部岩石圈热-流变结构研究进展[J]. 海洋地质与第四纪地质， 2006， (04)： 117-124.

[12] 陈爱华， 许鹤华， 马辉， 等. 马尼拉俯冲带缺失中深源地震成因初探[J]. 华南地震， 2011， 31（04）： 98-107.

[13] Yang T F， Lee T， Chen C H， et al. A double island arc between Taiwan and Luzon: consequence of ridge subduction[J]. Tectonophysics， 1996， 258（1-4）： 85-101.

[14] Chemenda A, Yang R, Hsieh C, et al. Evolutionary model for the Taiwan collision based on physical modelling[J].Tectonophysics， 1997, 274 (1-3)： 253-274.

[15] 高翔， 张健， 孙玉军， 等. 马尼拉海沟俯冲带热结构的模拟研究[J]. 地球物理学报， 2012， 55（1）： 117-125.

[16] 秦静欣， 郝天珧， 徐亚， 等. 南海及邻区莫霍面深度分布特征及其与各构造单元的关系[J]. 地球物理学报，2011， 54 (12)： 3171-3183.

[17] Hsu S K， Yeh Y， Doo W B， et al. New bathymetry and magnetic lineations identifications in the northernmost South China Sea and their tectonic implications[J]. Marine Geophysical Research, 2004, 25 (1)： 29-44.

[18] Pautot G， Rangin C. Subduction of the south China sea axial ridge below Luzon ( Philippines) [J]. Earth and Planetary Science Letters， 1989， 92 (1)： 57-69.

[19] 臧绍先， 陈奇志， 黄金水. 台湾南部-菲律宾地区的地震分布、 应力状态及板块的相互作用[J]. 地震地质，1994， 16（01）： 29-37.

[20] 朱俊江， 丘学林， 詹文欢， 等. 南海东部海沟的震源机制解及其构造意义[J]. 地震学报， 2005 (03)： 260-268.

[21] 孙珍， 钟志洪， 周蒂， 等. 南海的发育机制研究： 相似模拟证据[J]. 中国科学： D 辑， 2006（09）： 797-810.

[22] 丁巍伟， 程晓敢， 陈汉林， 等. 台湾增生楔的构造单元划分及其变形特征[J]. 热带海洋学报， 2005， 24（05）： 53-59.

[23] Bautista B C， Bautista M L P， Oike K， et al. A new insight on the geometry of subducting slabs in northern Luzon, Philippines[J]. Tectonophysics， 2001， 339 (3-4)：279-310.

[24] 陈志豪， 李家彪， 吴自银， 等. 马尼拉海沟集合形态特征的构造演化意义[J]. 海洋地质与第四纪地质，2009， 29（2）： 59-65.

[25] Clift P, Vannucchi P. Controls on tectonic accretion versus erosion in subduction zones: Implications for the origin and recycling of the continental crust[J]. Reviews of Geophysics, 2004， 42， RG2001.

[26] 姚伯初， 万玲. 南海岩石圈厚度变化特征及其构造意义[J]. 中国地质， 2010（04）： 888-899.

[27] 曹小林， 朱介寿， 赵连锋， 等. 南海及邻区地壳上地幔三维S 波速度结构的面波波形反演[J]. 地震学报,2001（02）： 113-124+223.

[28] 蔡学林， 朱介寿， 曹家敏， 等. 东亚西太平洋岩石圈三维结构及其地幔动力学[J]. 地学前缘， 2007（03）：21-38.

[29] Dl T，G S，Eds.. Geodynamics[M]. Lambridge： Cambridge University Press， 2002.

[30] 李志伟， 胥颐， 郝天珧， 等. 南海东北部及其邻近地区地壳上地幔P 波速度结构[J]. 地学前缘， 2009， 16（04）： 252-260.

[31] Lallemand S， Font Y， Bijwaard H， et al. New insights on 3-D plates interaction near Taiwan from tomography and tectonic implications[J]. Tectonophysics， 2001， 335(3-4)： 229-253.

[32] 瞿辰， 周蕙兰， 赵大鹏. 使用纵波和横波走时层析成像研究菲律宾海板块西边缘带和南海地区的深部结构[J]. 地球物理学报， 2007， 50（6）： 1757-1768.

[33] 臧绍先， 宋惠珍， 宁杰远.日本海俯冲带的热结构及热源的影响[J]. 地球物理学报， 1993， 36（2）： 164-173.

[34] 施小斌， 丘学林， 夏戡原， 等. 南海热流特征及其构造意义[J]. 热带海洋学报， 2003（2）： 63-73.

[35] Peacock S M， Wang K. Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan[J]. Science， 1999， 286(29)： 937-939.

[36] Kirby S, Engdahl E, Denlinger R. Intermediate-depth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs [J]. Geophysical monograph， 1996（96）： 195-214.

[37] Sugi N, Uyeda S. Subduction of young oceanic plates without deep focus earthquakes [J]. Bull. Soc. Geol.France， 1984， 26： 245-254.

[38] Liu L. Phase transformations, earthquakes and the descending lithosphere [J]. Physics of the Earth and Planetary Interiors， 1983， 32（3）： 226-240.

[39] Dobson D， Meredith P， Boon S. Simulation of subduction zone seismicity by dehydration of serpentine[J]. Science， 2002， 298 (5597)： 1407-1410.

[40] Peacock S M. The importance of blueschist to eclogite dehydration reactions in subducting oceanic crust [J].Geol Soc Am Bull， 1993， 105： 684-694.

[41] 魏春景， 张颖慧. 俯冲大洋岩石圈的相转变与俯冲带岩浆作用[J]. 科学通报， 2008， 53（20）： 2449-2459.

[42] Wang K， Hyndman R D， Yamano M. Thermal regime of the Southwest Japan subduction zone: effects of age history of the subducting plate [J]. Tectonophysics，1995， 248 (1-2）： 53-69.

[43] Hyndman R， Wang K. Thermal constraints on the zone of major thrust earthquake failure: The Cascadia subduction zone [J]. Journal of Geophysical Research，1993， 98（B2）： 2039-2060.