李子正, 丘学林, 贺恩远, 张浩宇, 王强
1 中国科学院南海海洋研究所, 边缘海与大洋地质重点实验室, 广州 510301 2 南方海洋科学与工程广东省实验室(广州), 广州 511458 3 中国科学院南海生态环境工程创新研究院, 广州 510301 4 中国科学院大学, 北京 100049 5 中石化石油物探技术研究院有限公司, 南京 211103 6 海洋地质调查研究所, 福建省厦门地质工程勘察院, 福建厦门 361008 7 自然资源部丘陵山地地质灾害防治重点实验室, 福建省地质工程勘察院, 福州 350002
海底高原的俯冲会对俯冲系统产生一系列影响,例如:影响弧后扩张进而改变俯冲带平面展布,加剧对弧前的构造侵蚀,也可能造成弧前地块的抬升变形等(Miura et al., 2004; Bangs et al., 2003; Wallace et al., 2008).除此之外,数值模拟的结果表明,俯冲角度减小、岩浆活动减弱等现象同样与海底高原俯冲有关(Gerya et al., 2009).
在马里亚纳海沟最南段,东端是太平洋地壳俯冲,俯冲角度较大;西端则靠近水深较浅的卡罗琳海底高原,俯冲角度较小;其上覆板片的岩浆活动和构造演化与输入板片存在密切联系,是研究海底高原俯冲动力学和岩浆过程的极好场所(Martinez et al., 2018; Zhu et al., 2019).但是,马里亚纳海沟最南段已有的深地震测线结果较少,靠近卡罗琳海底高原区域的研究相对空白;因此有必要进一步开展卡罗琳海底高原附近跨越海沟地壳结构的研究,得到对马里亚纳俯冲系统新的认识.
广角地震探测和重力模拟是研究海沟区域壳、幔特征的有效方式,前者可以获取海沟地壳和上地幔的精细速度结构,而重力模拟常常作为其补充,可以分辨大尺度的横向密度变化,尤其是在射线覆盖较少的下地壳和上地幔区域(Shulgin et al., 2009; Planert et al., 2010; Scherwath et al., 2010; Contreras-Reyes et al., 2011; Doo et al., 2015).本文基于2016年采集的TS01测线进行走时试算和模拟,结合对卫星重力异常的计算模拟,获取其地壳结构特征.
马里亚纳海沟是太平洋板块向菲律宾海板块俯冲形成的,构造活动剧烈,是研究俯冲带深部构造和俯冲系统演化历史的理想场所(图1a).自~50 Ma初始俯冲以来,马里亚纳海沟中段和北段具有强烈的岛弧岩浆作用和弧后扩张活动,形成了典型的沟-弧-盆体系(Stern, 2004):第一轮弧后扩张持续到15 Ma,这次弧后扩张形成了四国海盆(Shikoku Basin)和帕里西维拉海盆(Parece Vela Basin),残余的九州—帕劳海脊(Kyushu-Palau Ridge)停留在海盆西侧;自~8 Ma起,第二次的岛弧裂解和弧后扩张形成了马里亚纳海槽(Mariana Trough),海槽两侧是西马里亚纳海脊(West Mariana Ridge)和马里亚纳火山岛弧(Mariana Ridge).由于弧后扩张速率的差异,以及南北两端海底高原的俯冲,海沟和岛弧南北展布呈现弓形(Miller et al., 2006).在马里亚纳俯冲带中部和北部区域,日本和美国科学家进行了大量的主动源反射、折射地震以及重力观测,得到了弧前、岛弧和弧后盆地的二维、三维壳幔速度结构,以及俯冲的太平洋板片的地壳特征(Suyehiro et al., 1996; Kitada et al., 2006; Oakley et al., 2008; Kim et al., 2009; Takahashi et al., 2009;Cai et al., 2018).
马里亚纳海沟南段通常指关岛以南区域(图1b),由于卡罗琳海底高原(Caroline Plateau)碰撞对弧后扩张速率的影响,这里的海沟转变为近东西走向,且俯冲板片、弧前区域的特征与海沟中部、北部有显著不同(Fujiwara et al., 2000; Kobayashi, 2004; Martinez et al., 2018).在关岛西南侧的弧前区域,平行于海沟的弧前扩张分裂了始新世—中新世的弧前地块,这里新地壳增生且伴随广泛的构造变形,残余的岛弧地块夹杂其中(Martinez et al., 2018).马里亚纳海沟南段西端,残余的弧前地块相对完整,弧前地块后方的西南马里亚纳裂谷(South West Mariana Rift,SWMR)仍然处于构造张裂阶段(Martinez et al., 2018; Sleeper et al., 2021).在靠近雅浦(Yap)海沟的地方,难以识别到残余的始新世—中新世弧前地块,也几乎不存在沟-弧-盆体系(图1b).马里亚纳海沟南段最深点位于“挑战者深渊”,也是世界的最深点,可达10956 m(Lim et al., 2013);受到卡罗琳海底高原俯冲的影响,海沟和俯冲板片水深在平行海沟方向逐渐变浅(图1c).除此之外,在靠近海沟区域,俯冲板片弯曲产生了大量平行海沟的正断层,平均断距达320 m,远大于马里亚纳海沟中部和北部,造成了更为强烈的地幔蛇纹石化(Zhou and Lin, 2018; Wan et al., 2019; Zhu et al., 2021; He et al., 2022).尽管马里亚纳海沟南段在水深、构造等特征上有诸多特殊性,但是由于这里远离陆地,以往对这里的研究较少.为此,中国科学家在马里亚纳海沟南段进行了一系列主、被动源探测作业,旨在获取马里亚纳海沟南段的精细壳幔结构和俯冲动力学过程等关键信息(Wan et al., 2019; Zhu et al., 2019, 2021; He et al., 2022).本文中的TS01测线位于马里亚纳海沟南段最西端,靠近卡罗琳海底高原和雅浦海沟,马里亚纳海沟在此处变浅并逐渐尖灭.
TS01测线由“探索1号”科考船于2016年7—8月采集,投放并全部回收了9台海底地震仪(OBS),台站布设间距15 km,其中6台是6000 m级别的短周期OBS,3台是9000 m级别的宽频带OBS,均由中国科学院地质与地球物理研究所研发和提供.仪器投放后,科考船进行了主动源气枪作业,采用了4只大容量Bolt气枪组合,工作压力2000 psi,气枪总容量为6000 in3(1500×4).科考船航速5 km,放炮间隔60 s,共计激发963炮,测线总长160 km.本文将9台OBS记录到的原始数据解编为SAC文件后,根据气枪信号激发时间和位置(ukooa导航文件),采用课题组自编的Sac2y程序将SAC文件裁截并分道排列(赵明辉等,2004),得到SEGY格式数据.
位置校正工作也是OBS数据预处理中不可缺少的一环.二维测线的位置校正采用直达水波走时反演的方法,具体做法是采用蒙特卡洛法在海底生成一系列随机点,并计算每个点的直达水波理论走时,再采用最小二乘法的方式,计算理论走时与拾取直达水波震相的均方根误差(Root Mean Square, RMS),最终选取RMS最小的点为台站的实际坐底位置(敖威等,2010;张莉等,2013;李子正等,2021).值得一提的是,拾取直达水波震相前,本文对地震剖面进行了5~20 Hz的零相位滤波,这一参数不会影响直达水波到时(陈瀚等,2019;李子正等,2021).二维位置校正可以确定台站沿测线方向上的位置,在进行走时正演之前,本文将台站位置投影到测线正下方;经过二维位置校正后,直达水波同相轴关于0偏移距对称.在位置校正过程中,反演程序还会根据水声速度、多波束水深和近偏移距直达水波到时,求出地震数据的时间校正量Tadjust,该校正量将应用于所有近偏移距和远偏移距的震相拾取中(李子正等,2021).
数据预处理及前期校正后,利用SU软件对数据进行增益、滤波(滤波方式为零相位滤波,OBS02、09号台站滤波参数为3~6.5 Hz,其他台站为3~8 Hz)等处理,获得信噪比最高的单台地震剖面图,震相类型的识别除了根据震相的时距曲线特征外,还需要通过建立初始模型、走时试算和射线追踪来进行验证.
根据台站所处区域,本文选取了5个台站进行震相的展示.OBS01、02台站位于上覆板片,它们除了记录到上覆板片内的Pg震相,还记录到了广角反射震相:OBS01台站记录到了俯冲板片中地壳的反射震相(PcP)(图2);OBS02台站记录到的反射震相种类丰富,包括上下板片交界面的反射震相(PtopP)、一小段上覆板片内的反射震相和俯冲板片莫霍面的反射震相(PmP),其中PmP震相形态对应海沟轴线处“V”状地形(图3).OBS05和OBS06台站位于海沟轴线附近,这两个台站的Pg震相在正偏移距延伸较远,PmP震相出现于正偏移距Pg震相下方,OBS06台站还记录到了远偏移距的Pn震相;但是这两个台站的Pg震相在负偏移距延伸很短,几乎没有跨越海沟,表明了地震波能量在海沟附近的快速衰减,暗示这里的地壳高度破碎水化(图4、图5).OBS09台站位于测线最南侧,在负偏移距,Pg震相可以延伸到-45 km,-60~-80 km处有模糊的Pn震相,其形态同样受到海沟轴线处地形的影响;Pg震相下方出现了非常清晰且振幅强烈的中地壳反射(PcP)(图6a).
图5 OBS06台站的地震剖面、震相识别和射线追踪结果(a) OBS地震剖面和震相识别(折合速度为6 km·s-1); (b) 走时拾取(Pg震相使用绿色点拾取; PmP震相使用橙色点拾取;Pn震相使用粉色点拾取)和走时计算结果(黑色实线); (c) 射线路径.Fig.5 Seismic records, seismic phases, and raytracing of OBS06(a) Seismic profile of OBS (reduced to 6 km·s-1) with interpreted seismic phases; (b) Picked arrival times (Pg is picked by green dots; PmP is picked by orange dots; Pn is picked by pink dots) and calculated arrival times (black lines); (c) Calculated ray paths.
经过走时试算和射线追踪进一步确认震相的种类后,本文在xzplot软件中拾取震相的初始起跳点:采用绿色点拾取弧前地壳和俯冲板片内的Pg震相,拾取误差为50 ms;用蓝色点拾取来自俯冲板片中地壳界面的反射(PcP),青色点拾取板片交界面反射(PtopP),这两种震相拾取误差均为60 ms;俯冲板片的PmP震相用橙色点拾取,拾取误差同为60 ms;来自上地幔的Pn震相用粉色点拾取,设置的拾取误差为80 ms(图2b、图3b、图4b、图5b、图6b).震相拾取情况见表1.
表1 正演模型的震相参数Table 1 Seismic phases of the forward modeling
全部站位震相确认和走时拾取后,使用RayInvr软件进行走时正演试算和模拟,通过不断调节速度模型,使理论震相和拾取震相的RMS值和卡方值不断降低(Zelt and Smith, 1992; 丘学林等,2011).以本文中5个台站为例,最终计算的理论走时(黑色实线)和拾取震相(彩色点)拟合良好(图2b、图3b、图4b、图5b、图6b);9个台站所有震相走时拟合RMS值为67 ms,卡方值为1.476,详细震相拟合情况见表1.
最终速度模型如图7a所示:俯冲板片一侧,地壳厚度约13 km,在靠近海沟区域,地壳顶部速度为4.5 km·s-1,底部速度约7.0 km·s-1;远离海沟区域上地壳显著变厚,地壳顶部速度更低,约4.0 km·s-1,但是下地壳速度略高,约7.2 km·s-1.弧前地壳呈现楔状,速度上有显著的横向变化:在靠近海沟部分和模型边缘区域速度较低,但是在OBS02台站下方的地形隆起处,地壳顶部速度可达5.6 km·s-1.
图7 最终的正演结果图(a) TS01测线最终速度模型; (b) 模型射线覆盖图.黑色粗实线为广角反射震相控制的界面, 模型的横向和纵向网格大小分别为0.5 km和0.25 km.Fig.7 Final velocity model obtained by forward modeling(a) Final velocity model of profile TS01; (b) Ray coverage of the model. Black thick solid line represents the interfaces controlled by wide-angle reflection phases, the horizontal and vertical grid of velocity model were 0.5 km and 0.25 km, respectively.
在TS01测线9个台站的地震剖面上,本文均没有识别到Ps震相,证实这里沉积层稀薄.除了壳幔速度值和莫霍面之外,9个台站记录的广角反射震相还约束了上下板片交界面、俯冲板片中地壳界面的形态(粗实线标出).
射线覆盖图显示(图7b),模型的上地壳部分射线覆盖良好,大部分区域射线覆盖大于60次;由于PmP和Pn震相较少,模型的下地壳大部分区域射线覆盖在10~30次之间,在上地幔顶部,射线只覆盖了局部区域.
布格重力异常可以反映莫霍面的深度变化,是研究地球深部构造、划分构造单元的有效方式(曾华霖,2005;郝天珧等,2008).根据是否进行地形校正,布格校正可以分为简单布格校正和完全布格校正(吕川川等,2009;郝天珧等,2014);在海沟区域,水深快速变化,仅进行中间层校正会带来很大误差,所以地形校正是必不可少的.
卫星重力数据具有测量范围广、覆盖均匀等特点,精度足以分辨10~20 km波长的异常特征,常被用于海沟区域构造的研究(Kim et al., 2009;Martinez et al., 2018).本文在1′网格卫星重力数据基础上(Sandwell et al., 2014),首先进行中间层校正,再根据多波束水深数据进行地形校正(海水和地壳密度分别取1.03 g·cm-3和2.67 g·cm-3),由此得到完全布格重力异常.经过这两步校正,相当于将海水物质替换为地壳物质,消除了海底地形的影响,在重力均衡区域(远离海沟)可以大致反映莫霍面的起伏;靠近海沟区域因不符合重力均衡假设,所以这里较高的布格重力异常值不能说明此处地壳厚度较薄.
布格重力异常图显示(图8a),俯冲板片一侧,洋壳部分呈现较高的布格重力异常(550~600 mGal),而靠近卡罗琳海底高原区域布格重力异常较低(<500 mGal).平行海沟方向,俯冲板片的布格重力异常快速变化(图8b,A-B测线):随着向西靠近卡罗琳海底高原,布格重力异常值减少了约200 mGal.岛弧一侧,马里亚纳海槽的布格重力异常值在500 mGal左右,较低的重力异常值(<400 mGal)对应着始新世的弧前地块,以及海底高原北部的一处弧前山脊.
图8 马里亚纳海沟南段布格重力异常图(a) 布格重力异常图; (b) 俯冲板片沿平行海沟方向(对应图(a)中实线A-B)的布格重力异常变化.Fig.8 Bouguer gravity anomaly of Southmost Mariana Trench(a) Map of Bouguer gravity anomaly; (b) Trench-parallel variations of Bouguer gravity anomaly at the subducting crust (Corresponding to black solid line A-B in (a)).
重力数据优势在于可以反映壳幔物质大尺度的横向密度变化,也可以粗略估计地壳厚度,除此之外,重力模拟还可以约束地震射线覆盖不足的区域(郝天珧等,2008; Scherwath et al., 2010).因此,本文基于多波束水深数据和卫星重力数据(Lim et al., 2013; Sandwell et al., 2014),使用Mask软件(姚长利等,2003)进行正演模拟,分析TS01测线密度特征.值得说明的是,卫星重力数据虽不及船测重力数据精确,但是以往研究通过对二者的比较认为,卫星重力数据可以分辨10~20 km波长的异常特征,其精度满足海沟区域构造研究的需求(Kim et al., 2009).
图9a展示了计算数据与卫星重力数据的差异(均方差~11.4 mGal).重力模型如图9b,密度参数的设置参考了2.2节的走时正演结果(Brocher, 2005),具体参数设置如表2.
表2 重力模型参数Table 2 Parameters of the gravity model
图9 TS01测线重力模拟结果(a) 自由空气重力异常(蓝色虚线)和计算重力异常(含蛇纹石化上地幔模型的计算异常用黑色实线表示;不含蛇纹石化上地幔模型的计算 异常用黑色虚线表示); (b) TS01测线的重力模型,虚线上方为地震射线约束区域 (每个区域的密度值用数字标出,单位为g·cm-3).Fig.9 Result of gravity modeling of profile TS01(a) Free air gravity anomaly (blue dashed line) and calculated gravity anomaly (The gravity anomaly of a model with serpentine upper mantle is shown by black solid line; The gravity anomaly of a model without serpentine upper mantle is shown by black dashed line); (b) Gravity model of profile TS01, and above the dotted line is the region constrained by seismic rays (Density of each part of the model is numbered. Units: g·cm-3).
在远离海沟区域,俯冲板片的地壳厚度约12~13 km,与走时正演的结果一致(图7).靠近海沟区域(40~95 km),自由空气重力异常大幅度下降了近240 mGal,这里重力模型的调整值得详细说明:2.3节的走时正演结果表明,靠近海沟区域地壳速度仅在地壳顶底部有小幅变化(图7a),故重力模型这一部分的地壳密度也只有小幅降低,此时计算得到重力异常值远高于实际值(图9a,黑色虚线);进一步调低靠近海沟下方和已俯冲板片的上地幔密度后(图9b),计算重力异常和实际吻合良好(图9a,黑色实线).因此本文认为,重力密度模型反映了上地幔区域密度的大幅度降低.
卡罗琳海底高原是形成于卡罗琳火山热点之上的大型火成岩省(Dong et al., 2018; Zhang J and Zhang G L, 2020; Zhang et al., 2021),海底高原形成过程中通常会发生地壳的增厚(Farnetani et al., 1996);TS01测线靠近卡罗琳海底高原,其12~13 km的地壳厚度也接近一些海底高原边缘区域的地壳厚度(Korenaga and Sager, 2012; Zhang et al., 2016).TS01测线距离Wan等(2019)的测线约155 km,俯冲板片一侧较Wan等(2019)测线更浅的水深(~2000 m)和更低的布格重力异常(~100 mGal)也暗示着更厚的地壳(图1c、图8b).除此之外,OBS01、03、08、09台站记录的PcP震相约束了俯冲板片的中地壳,这种震相在洋壳中较为罕见,但常出现在海底高原区域(Shulgin et al., 2009; Scherwath et al., 2010).
速度模型中(图7a),上覆板片的地壳呈现楔状,其顶部速度在模型20~40 km之间较高,地形上对应弧前区域的一个小山脊;而在40~70 km间上覆板片速度值随着靠近海沟逐渐降低:这部分速度值在4.5~5.2 km·s-1之间,高于沉积物组成的增生楔(Bangs et al., 2003; Eakin et al., 2014; Liu et al., 2018).增生楔常出现在陆缘沉积物丰富的洋-陆型俯冲带,如智利海沟、马尼拉海沟等;而马里亚纳海沟属于洋-洋型俯冲带,远离大陆,因此增生楔不发育.构造侵蚀是马里亚纳海沟更为常见的地质过程(Bloomer, 1983; Miura et al., 2004),该过程通常会导致弧前地壳的破碎水化,进而使靠近海沟区域呈现相对较低的地震波速度(Contreras-Reyes et al., 2011).除此之外,在距离TS01测线很近的雅浦海沟,卡罗琳海底高原和正断层的俯冲增强了构造侵蚀,导致了陡峭的海沟内壁和火山岛弧靠近海沟等地质现象(Lee, 2004; Zhang et al., 2019; Zhang J and Zhang G L, 2020),故本文认为TS01测线具备构造侵蚀的有利条件.综上,本文将上覆板片40~70 km之间区域解释为构造侵蚀后高度破碎水化的海沟内壁和弧前地块.以往研究结果表明,破碎水化的弧前地壳对地震波能量有强烈的吸收衰减作用(Christeson et al., 2000; Zhu et al., 2010),本文的地震剖面也显示,Pg震相几乎无法跨越海沟传播.
在缺少射线覆盖的上地幔区域,本文借助重力模拟揭示了海沟下方厚达25 km的上地幔低密度区域(~2.93 g·cm-3):该区域厚度与朱高华等(Zhu et al., 2021)的研究接近;且根据速度-密度的经验公式(Brocher, 2005),2.93 g·cm-3的岩石密度对应的横波速度约4 km·s-1,同样与朱高华等(Zhu et al., 2021)采用瑞雷波层析成像获得的SV波速度范围对应良好(3.6~4.1 km·s-1).板片俯冲后会受到高温高压影响发生脱水,故已俯冲板片的上地幔密度可能略高于海沟正下方区域.以往研究表明,上地幔的蛇纹石化还会导致壳幔分界的模糊,难以观测到清晰的PmP震相(Brocher et al., 2003),而TS01测线的9个台站没有记录到来自海沟正下方区域(模型65~100 km处)的PmP震相,OBS02、05、06台站的PmP震相也不够清晰(图3a、图4a、图5a、图7b);因此,PmP震相的模糊和缺失可能也是上地幔发生蛇纹石化的证据之一.
本文采用RayInvr软件进行走时试算和模拟,获取了TS01测线附近精细地壳结构;重力模拟与速度模型相互验证、共同约束了TS01测线的壳幔结构,揭示了地幔顶部的横向密度变化.由此得出以下几点认识和结论:
(1) OBS走时正演和卫星重力数据模拟的结果均表明,TS01测线附近俯冲板片的地壳厚度约12~13 km,介于洋壳厚度和海底高原地壳厚度之间.
(2) 弧前一侧靠近海沟区域的低速反映了构造侵蚀现象,海底高原的临近和俯冲可能会增强构造侵蚀.
(3) 重力模拟结果和卫星重力数据吻合良好,靠近海沟区域较低的上地幔密度反映了这里的蛇纹石化现象.
致谢感谢中国科学院深渊科考队在航次中的辛苦工作,以及“探索一号”船长与全体船员的努力和配合.数据采集处理和成文过程中,得到中国科学院南海海洋研究所赵明辉研究员、黄海波副研究员、张佳政副研究员,中国海洋大学朱俊江教授等的帮助指导,中国科学院地质与地球物理研究所王元高级工程师全程负责OBS的投放和回收,匿名审稿专家的建设性意见提升了文章质量,在此表示感谢.