引用本文
徐志国, 王君成, 王宗辰, 梁姗姗, 史健宇. 2019年11月14日印尼马鲁古海7.1级地震的震源机制及海啸数值模拟[J]. 地震地质, 2020,42(6): 1417-1431.
XU Zhi-guo, WANG Jun-cheng, WANG Zong-chen, LIANG Shan-shan, SHI Jian-yu. FOCAL MECHANISM AND TSUNAMI NUMERICAL SIMULATION OF THE NOVEMBER 14, 2019 MOLUCCA SEA MW7.1 EARTHQUAKE [J]. SEISMOLOGY AND GEOLOGY, 2020,42(6): 1417-1431.  
DOI:10.3969/j.issn.0253-4967.2020.06.010
Permissions
Copyright©2020, 《地震地质》编辑部
《地震地质》编辑部 所有
2019年11月14日印尼马鲁古海7.1级地震的震源机制及海啸数值模拟
徐志国1,2), 王君成1), 王宗辰1), 梁姗姗3)),*, 史健宇1)
1)国家海洋环境预报中心, 北京 100081
2)中国科学院大学, 计算地球动力学重点实验室, 北京 100049
3)中国地震台网中心, 北京 100045
*通讯作者: 梁姗姗, 女, 工程师, E-mail: liangshanshan@seis.ac.cn。

〔作者简介〕 徐志国, 男, 1979年生, 2007年于中国地震局地球物理研究所获地球探测与信息专业硕士学位, 高级工程师, 主要从事实时地震学与海啸预警预报应用研究, 电话: 010-62105791, E-mail: xuzhg04@sina.com。

收稿日期: 2020-02-26
基金项目: 自然资源部 “全球变化与海气相互作用”(GASI-GEOGE-05)和同济大学海洋地质国家重点实验室开放课题(MGK20204)共同资助
摘要

2019年11月14日16时17分(UTC), 印尼马鲁古附近海域发生了 MW7.1地震, 并引发小规模海啸。 为深入理解和认识该地震的震源参数和发震构造特征, 评估其引发海啸的危险性, 文中初步分析了该地震的区域构造背景、 震源机制以及海啸数值模拟等内容。 W震相快速矩张量解反演结果表明马鲁古海地震是一次浅源、 以高倾角右旋斜向逆冲为主的地震事件, 余震呈SSW-NNE向带状分布, 推测此次地震是在马鲁古海东、 西两侧哈马黑拉弧-桑义赫弧碰撞所产生的区域应力场的作用下发生的以挤压破裂为主的地震事件。 海啸的数值模拟结果表明, 在高倾角发震断层近垂直倾滑的逆冲过程中, 海底地形产生了同震垂直位移, 使得震源上方的水体突然抬升, 从而产生了小规模的局地海啸; 震源周边大部分潮位站记录的海啸首波理论波形和观测波形的到达时间和波形幅度都较为一致, 说明结合W震相反演所得的地震断层面的几何参数能够用于海啸早期预警, 可有效地预测海啸地震产生的海啸波高, 并对于认识海啸成灾过程及灾害分布具有重要的现实意义。

关键词: 马鲁古海地震; 海啸; 弧-弧碰撞; 震源机制; 数值模拟
中图分类号:P315.2 文献标志码:A 文章编号:0253-4967(2020)06-1417-15
FOCAL MECHANISM AND TSUNAMI NUMERICAL SIMULATION OF THE NOVEMBER 14, 2019 MOLUCCA SEA MW7.1 EARTHQUAKE
XU Zhi-guo1,2), WANG Jun-cheng1), WANG Zong-chen1), LIANG Shan-shan3), SHI Jian-yu1)
1)National Marine Environmental Forecasting Center, Beijing 100081, China
2)Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China
3)China Earthquake Networks Center, Beijing 100045, China
Abstract

A strong earthquake with magnitude MW=7.1 occurred in the area of Molucca Sea, Indonesia on November 14, 2019(Coordinated Universal Time, UTC), and then generated a small-scale local tsunami. In order to better understand the earthquake source characteristics and seismogenic structure, as well as to assess the hazard of tsunami caused by earthquake, this paper mainly focuses on the regional tectonic background, the focal mechanism, and tsunami numerical simulation for the Molucca Sea MW7.1 earthquake. The broadband seismic waveforms from IRIS Data Management Center are used to estimate the moment tensor solution of this earthquake by W phase method. The result shows that the Molucca Sea earthquake occurred at a shallow depth on a high dip-angle, right-lateral reverse fault, the aftershocks were distributed along the SSW-NNE direction and concentrated near the main shock. These results indicate the Molucca Sea earthquake with characteristic of compressional rupture occurred in the complex plate boundary region of eastern Indonesia, which is dominated mostly by the collision interaction of the Halmahera slab and the Sangihe slab in the east and west sides of Molucca Sea under control of current regional stress field. The coseismic displacements of Molucca Sea MW7.1 earthquake calculated using Okada's model of rectangular dislocation in a uniform elastic half-space show that the Molucca Sea earthquake generated vertical coseismic deformation with a maximum uplift of 0.15m when the rupture occurred along the high dip-angle reverse fault. The synthetic tsunami waveforms are provided by COMCOT tsunami modelling package solving the nonlinear shallow water wave equations based on the determined fault geometry from W phase inversion. These studies indicate the vertical coseismic deformation resulting in the sudden uplift of water volume above the earthquake source, and finally inducing a small-scale local tsunami. The energy of tsunami mainly propagates to both side of the fault, and part of energy propagates to Sula Islands of Indonesia along the fault dislocation direction; and compared with the first cycle of tsunami records observed by tide gauges deployed along the coastal line of earthquake source region, the observed tsunami head wave fits well with the synthetic wave, both are consistent in amplitude and tsunami arrival time, but the follow-up waveforms are quite different. The numerical simulation of tsunami shows that, in combination with the fault geometry parameters obtained by W phase fast inversion, the tsunami numerical model can be used for tsunami early warning, and it provides sufficient accuracy for forecasting tsunami wave height, thus, having great practical significance for understanding the propagation process and disaster distribution of tsunami.

Keyword: Molucca Sea earthquake; tsunami; arc-arc collision; focal mechanism; numerical simulation
0 引言

2019年11月14日16时17分(UTC), 印尼马鲁古附近海域(1.621°N, 126.416°E)发生了MW7.1地震, 震源深度为33km。 地震发生后, 位于北京的联合国教科文组织政府间海洋学委员会南中国海区域海啸预警中心(UNESCO-IOC/SCSTAC)快速向南中国海区域周边国家发布了海啸预警信息, 预计在震源周边约300km范围内可能发生局地海啸。 震源区附近的潮位站水位观测数据表明, 该地震引发了高约10cm的海啸波。 马鲁古海MW7.1地震的震感强烈, 最大烈度达Ⅵ 度①(① https://earthquake.usgs.gov/earthquakes/eventpage/us60006bjldyfiintensity。), 但由于此次地震震中位于马鲁古海中, 距离沿岸最近处约130km, 故并未给马鲁古海沿岸地区造成人员伤亡和损失。

从震中位置来看(图 1), 马鲁古海MW7.1地震发生于印尼东部地区, 处于太平洋板块、 印度-澳大利亚板块、 菲律宾海板块和巽他板块的会聚交界带。 GPS测量结果显示, 菲律宾海板块相对巽他板块以90~100mm/a的速度向NW向俯冲会聚(Bird, 2003)。 太平洋板块、 印度-澳大利亚板块和欧亚板块的相互作用及其长期演化造就了马鲁古海地区复杂的地质构造背景, 形成了极其复杂的俯冲、 碰撞、 增生以及弧后扩张等构造过程。 通过马鲁古海历史地震的震源深度剖面AA'(图1b)可知, 俯冲于哈马黑拉弧之下的岩石圈板片长达200~300km, 形成贝尼奥夫带; 而在桑义赫弧之下的贝尼奥夫带在深度上至少可识别至约600km, 其历史地震的分布特征与该区域整体的构造背景一致。 由于马鲁古海区域构造复杂, 因此对该区域的研究一直是地震学与地球动力学的热点之一。 大量研究结果表明, 马鲁古海板块向E下插于哈马黑拉板块之下, 向W下插于桑义赫板块之下, 形成反向、 非对称的倒U字形双向俯冲会聚带(McCaffrey, 1982), 也被称为离散型双俯冲带(Divergent Double Subduction)(Zhang et al., 2017), 这使得马鲁古海板块几乎完全俯冲消减在上覆的桑义赫板块和哈马黑拉板块之下, 形成了桑义赫弧和哈马黑拉火山弧。 马鲁古海两侧的哈马黑拉弧和桑义赫弧以10cm/a的速度相向俯冲会聚(Hinschberger et al., 2005), 形成典型的弧-弧碰撞带, 也是现今全球构造上弧-弧碰撞造山的惟一实例(潘桂棠等, 2004)。 马鲁古弧后洋盆向西侧的桑义赫火山弧和东侧的哈马黑拉弧俯冲, 两者的弧前带已相互碰撞, 极性相反的两弧之间均为俯冲碰撞杂岩(Collision complex), 碰撞杂岩向外分别逆冲至哈马黑拉弧和桑义赫弧这2个面对面的弧前区上(Silver et al., 1978), 两者的弧前区均向后逆冲于各自依附的弧体之上, 并形成塔劳海岭(Talaud Ridge)。 沿碰撞杂岩的冲断层系向S追踪可达到苏朗断层(Sorong Fault), 向N则进入棉兰老岛。 在塔劳群岛与棉兰老岛之间, 并无马鲁古海板块在菲律宾海沟与哥达巴托海沟之间发生俯冲的迹象(Hall, 1987; Hall et al., 1990)。

图 1 研究区域的构造背景
a 马鲁古海及周缘地区主要的俯冲带(锯齿状线)和历史地震(彩色圆)分布, 黑色五角星为马鲁古海MW7.1主震震中, 黑色三角形表示火山; b 研究区域AA'剖面的震源深度分布, 红色火山符号表示桑义赫弧和哈马黑拉弧Fig. 1 Regional tectonic setting of the study areas.

受板块俯冲作用的影响, 马鲁古海地区地震、 海啸、 火山喷发等灾害频发, 是全球地质灾害最为活跃的区域之一。 据USGS历史地震目录统计, 自1900年以来, 在震源区150km范围内发生6级以上地震115次, 其中7级以上地震17次, 震级最大的一次为1968年8月10日发生的7.6级地震, 位于此次地震(即2019年7.1级地震)的SW侧, 距离约为23km。 美国全球历史海啸数据库(NGDC)资料表明, 自2000BC至今, 由海底地震或火山喷发等因素引发的确定性历史海啸事件有7次, 其中最大的一次海啸灾害由1871年3月3日位于印尼桑义赫群岛最南端的鲁昂火山喷发所致, 观测到的最大海啸波幅高达25m, 死亡人数超过600人, 该次灾难给沿岸人民的生产生活带来严重影响; 1859年6月28日马鲁古海附近发生的7.0级地震也引发了大规模区域海啸, 最大海啸波幅达10m, 但并未造成人员伤亡; 时间相隔最近的一次为2014年11月15日马鲁古海发生的7.1级逆冲型地震所引发的小规模海啸, 震源区附近印尼Jailo潮位站于此次灾难中观测到9cm高的海啸波。

尽管震源区附近地震频繁多发, 但由于这些地震远离海岸, 单纯的地震灾害并未给沿岸地区造成严重的人员伤亡和经济损失, 但由地震、 火山喷发所诱发的海啸灾害却不容小觑。 马鲁古海MW7.1地震发生后, 美国地质调查局(USGS)和哈佛大学(Harvard)等相关地震机构给出了此次地震的震源机制解, 结果均显示此次地震为一次高倾角的逆断层型事件, 断层面走向为NE或SW向。 此次地震所处地理位置较为特殊——马鲁古海板块完全俯冲并埋藏于碰撞杂岩之下, 无地表迹象可依。 为了深入理解和认识此次马鲁古海地震的发震机制和发震断层特征, 进一步解释此次地震产生海啸的原因, 本文首先利用W震相方法(Kanamori et al., 2008; Hayes et al., 2009; Zacharie et al., 2012)快速得到此次地震的断层面几何参数、 地震矩、 矩心位置和深度以及震源区应力状态等震源特征参数, 并结合震源区的地质构造背景探讨地震的发震构造; 随后, 采用本文中获得的点源震源模型, 应用美国Cornell大学开发的COMCOT(Cornell Multi-grid Coupled Tsunami Model)模型(Liu et al., 1995; Wang et al., 2006)进行海啸数值模拟, 通过模拟结果分析马鲁古地震引发的海啸波的能量分布特征, 并比较近岸潮位站观测海啸波高与模拟海啸波高之间的差别, 从而验证实际断层参数对模拟结果的影响。

1 数据和方法

本研究采用美国地震学研究联合会数据管理中心(IRIS/DMC)提供下载的全球数字地震台网(GSN)记录的宽频带(BH)三分量地震波形数据, 选取了信噪比较高且相对震源具有均匀方位角覆盖的数十个台站的波形资料, 用于研究马鲁古海MW7.1地震的震源特征(图 3)。 在海啸模拟过程中, 采用国际海道测量组织(IHO)和政府间海洋学委员会(IOC)联合编制的GEBCO_2014海底地形和沿岸高程数据集进行数值模拟, 其地形分辨率为30arc-sec(1/120°)。 同时, 为验证海啸数值模拟结果的可靠性, 我们收集了印尼地理信息局(Geospatial Information Agency, BIG)提供的马鲁古海周边地区沿岸的4个验潮站(图7a)记录的海啸波形数据①(①http://tides.big.go.id。), 并对观测值与模拟结果进行比较, 以验证海啸数值模拟的可靠性。

图 3 部分参加反演的台站记录到的观测波形(黑色)与W震相获得的矩张量解理论波形(红色)的对比
图中蓝色五角星为地震震中, 黄色和橙色实心圆为反演计算所用台站的分布, 橙色实心圆为理论地震图(红线)与观测地震图(黑线)对应的台站位置。 波形上的红点表示W震相的时间窗; 波形上侧的字母分别表示台网代码、 台站代码、 通道分量和台站位号Fig. 3 A part of waveform comparison between the observed(black)and its synthetic by W-phase moment tensor inversion(red).

首先, 我们采用基于W震相求取矩张量解的方法反演确定此次地震的震源特征参数。 由于W震相是P波和S波之间斜坡状的长周期震相(100~1 000s), 其能量为P、 PP、 SP和S等多个体波震相的叠加, 该震相的群速度为4.5~9km/s, 比传统的面波速度快, 且振幅不容易出现限幅状态, 适用于快速测定较大震级地震的震源参数, 可为大地震应急救灾和海啸预警提供服务(郭志等, 2018; 梁姗姗等, 2018)。 Kanamori等(2008)采用W震相方法反演了MW≥ 7.5地震的震源参数, 分析结果表明, W震相能够快速反演得到可靠的震源参数, 能够很好地确定大地震引发海啸的潜能。

由W震相反演未知震源的参数向量m可由式(1)表示:

m=fηc(1)

其中, 地震矩张量f=[Mrr, Mpp, Mtt, Mrp, Mrt, Mpt]T, η c表征矩心时空坐标, η c=[θ c, φ c, hc, τ c]T, 参数θ c、 φ c、 hc和τ c分别表示矩心余纬度、 矩心经度、 矩心深度及矩心相对发震时刻的时间偏移。

在求解震源未知参量m时, 假设点源位置随着震源过程时间的变化而变化。 若已知矩心位置、 震源持续时间和矩心时间偏移, 则把其作为最佳矩心时空坐标η c的初始值, 震源参数的求解过程可归纳为最小二乘线性意义上的线性方程求解。 为了确定最佳的矩心时空坐标η c, 可通过空间网格搜索法, 使得定义的误差函数最小:

χ(m)=12(SW(m)-dW)(SW(m)-dW)(2)

误差函数主要用于表征理论和实际波形二者之间的拟合程度, 其中SW和dW为W震相理论地震图和实际观测波形。 为了提高震源机制解的反演效率, 需要预先构建格林函数库。 我们基于简振正型叠加方法(Woodhouse, 1988), 采用PREM全球一维速度模型(Dziewonski et al., 1981)计算格林函数。

在地震矩张量反演之初, 我们对5°~85°震中距范围内的宽频带三分量数字波形记录进行去除仪器响应、 去均值和倾斜趋势等处理, 采用0.001~0.01Hz带通滤波器对波形数据进行滤波, 自动提取理论初至P波后15倍震中距(单位为°)时间窗(单位为s)范围内的W震相。 在反演过程中, 将USGS测定的地震基本参数(包括发震时刻、 震中经纬度、 震源深度和震级)作为初始值, 采用空间网格搜索方法反演得到震源参数的最优解。

随后, 应用美国Cornell大学开发的COMCOT模式对马鲁古海MW7.1地震开展海啸数值计算。 该模式基于球坐标系下线性和非线性浅水长波方程模拟海啸在深海和浅海中的传播。 文中采用球坐标系下非线性浅水方程模拟近岸潮位站的海啸波形, 考虑底摩擦效应, 并对模拟结果和观测结果进行比较。 在进行地震海啸的数值模拟过程时, 通常认为初始海啸波由海底断层的瞬时垂向错动引起, 同时假设海水表面的升降与海底位移变形一致, 忽略断层破裂过程的复杂性对海啸的影响。

基于上述矩张量反演得到的震源机制解参数, 利用Blaser等(2010)提出的大洋俯冲带逆冲型地震震级(MW)与断层长度(L)、 宽度(W)之间的统计关系式计算断层参数:

log10L=a+b×MW, a=-2.81, b=0.62(3)

log10W=a+b×MW, a=-1.78, b=0.45(4)

平均滑移量(D)通过式(5)、 式(6)得到:

M0=μDS(5)

S=W×L(6)

式中, S为断层面积, μ为剪切模量, 取μ=3×1010N/m2。 利用Okada(1985)提出的均匀弹性半空间中的有限矩形位错模型计算震后海底表面的垂直同震形变, 为海啸的生成提供初始条件。 图 2 为有限矩形位错模型示意图, 其中L为断层长度, W为断层宽度, H为震源深度(以断层下底面为参考), θ 为走向角, δ 为倾角, λ 为滑动角。 以断层的走向方向为X轴, 水平地面的垂线方向为Z轴, 在水平面内垂直于X、 Z轴的方向为Y轴。 d为位错矢量, d=[d1, d2, d3]。 其中, d1、 d2、 d3分别表示断层上盘相对下盘任意位错的走滑、 倾滑和张性位错分量, 在整个断层面上各位错分量是一致的。 滑动矢量的绝对值D表示位错量(或平均滑移量)的大小。 n为位错面法线矢量, 其中n·d=0。

图 2 有限矩形位错模型示意图Fig. 2 Sketch of a finite rectangular dislocation model.

2 结果分析
2.1 矩张量反演

经过多次反复迭代反演, 最终得到马鲁古海MW7.1地震的矩张量解和相关震源特征参数, 并通过空间网格搜索法获得了最佳矩心位置和矩心深度。 其中, 断层面节面Ⅰ 的走向、 倾角和滑动角分别为225.5°、 48.5°和109.9°; 节面Ⅱ 的走向、 倾角和滑动角分别为16.8°、 45.2°和69°; 地震矩M0=4.921×1019N·m, 对应矩震级MW=7.1; 所得到的矩心位置为(1.721°N, 126.316°E), 比初始破裂点(1.63°N, 126.4°E)向NW偏移约12km, 矩心深度为35.5km, 与破裂起始深度33km相差不多; 震源区应力主轴的空间取向为:主压力轴P的方位角为305°, 倾角为4°; 主张力轴T的方位角为191°, 倾角为80°。 通过对比参与反演台站记录的W震相的观测波形与理论波形拟合情况(图 3)可以看出, 大部分波形拟合较好, 观测波形和理论波形的拟合均方根残差仅为0.053mm。

为了更好地评价本研究结果的可靠性和稳定性, 我们收集了哈佛大学全球矩心矩张量解(GCMT)和USGS给出的此次地震的相关震源参数(表1), 采用Kagan(1991)提出的最小空间旋转角衡量震源机制解之间的差别。 以本文反演的震源机制解为参考解计算其与其他结果之间的最小旋转角, 结果如表1所示。 由表1可以看出, 所得最小旋转角(Kagan角)均较小, 表明3种结果具有较好的一致性, 均显示马鲁古海MW7.1地震是一次高倾角逆冲型事件。

表1 不同研究机构给出的马鲁古海MW7.1地震震源参数对比 Table1 Comparison of the Molucca Sea MW7.1 earthquake source mechanism from different institutes

马鲁古海MW7.1地震发生后中强余震活动不断(图4b), 截至2019年12月31日共记录到M4以上余震131次, 余震序列沿SSW-NNE向呈块状分布, 与塔劳海岭的走向基本一致, 余震延展长度约为20km。 沿剖面BB'余震震源主要集中分布在马鲁古海板片顶部10~50km深度范围的岩石圈内(图4c), 位于桑义赫弧和哈马黑拉弧碰撞的杂岩带内, 与历史地震的震源分布特征一致。 本次主震的震中位于马鲁古海中, 且远离岸边, 附近台站分布稀疏, 同时由于该区对中小地震监测和定位的能力有限, 故针对本次地震的余震分布无法提供断层面倾向等信息, 尚不能明确此次地震的发震构造。 由主震矩张量反演的震源区应力主轴的空间取向可知, B轴为NNE向, 与塔劳海岭的走向一致, 主压应力P轴近水平, 与塔劳海岭走向垂直, 反映了此区受到了NW-SE向挤压作用力的影响。 马鲁古海地震5次主要中强余震的震源机制与主震相差不多, 都为近NE向逆断型破裂, 反映余震与主震具有相同的震源应力场特征, 其主压应力轴基本为NW-SE向, 且与历史地震的震源应力场特征相似。 相关研究结果表明, 塔劳海岭下的高角度逆断层作用几乎是碰撞带内所有地震产生的原因(McCaffrey, 1982)。 马鲁古海地震及其强余震是在一个稳定的、 主压应力P轴以NW-SE向为主的震源应力场控制下发生的, 与弧-弧逆冲挤压方向相一致。

图 4 印尼马鲁古海地震的震源机制和余震空间分布
a 研究区域Slab2.0(Hayes et al., 2018)板片的俯冲深度; b 马鲁古海地震震源区余震序列的震中(红色实心圆)分布和主震(黄色沙滩球)及中强余震(红色沙滩球)的震源机制解, 锯齿状线分别为桑义赫弧和哈马黑拉弧; 黑色三角形为火山, 黑色虚线为马鲁古海板块的俯冲深度; c 剖面BB'震源深度分布, 紫色实线为马鲁古海板块沿BB'的俯冲深度Fig. 4 Focal mechanisms and spatial distribution of mainshock and its main aftershocks.

欧亚板块和菲律宾海板块碰撞导致马鲁古洋盆向西侧的桑义赫弧和东侧的哈马黑拉弧俯冲, 且两者的弧前带已经碰撞, 桑义赫弧前带逆冲于哈马黑拉弧前带之上, 形成塔劳海岭, 两者的弧前区均向后逆冲于各自依附的弧体上。 两弧之间的浅层地壳由弧前盆地(forearcbasin)和增生杂岩(accretionalcomplex)组成(图 5)。 马鲁古海地震的反射、 折射和重力剖面显示, 增生杂岩广泛逆冲于邻近的桑义赫弧和哈马黑拉弧的弧前盆地之上(Silver et al., 1978; Widiwijayanti et al., 2004; Zhang et al., 2017)。 在板块运动产生的强大的挤压作用下, 增生杂岩体是地震最为活跃的区域, 以强烈的浅源地震活动为特点。 因此, 结合震中区域附近的地质构造特征, 初步推断震源机制节面I为可能的发震断层面。 该地震可能由于桑义赫弧前逆冲至哈马黑拉弧前并发生强烈形变所致, 其发生受欧亚板块和太平洋板块相向俯冲作用所影响, 与欧亚板块和太平洋板块强烈推挤下产生EW向缩短的双弧构造有关。

图 5 马鲁古海地区弧-弧碰撞构造的示意图(修改自Zhang et al., 2017)Fig. 5 Sketch of arc-arc collision zone of cross-section BB' in Molucca Sea(adapted after Zhang et al., 2017).

2.2 海啸数值模拟

采用马鲁古海矩张量反演所得实际断层面几何参数(节面I, 走向为225.5°、 倾角为48.5°、 滑动角为109.9°)构建断层模型, 由Blaser等(2010)推导的统计关系式估算断层长44km、 宽25km, 平均滑动量为1.2m。 采用Okada弹性半无限空间位错模型计算同震垂直位移, 并作为马鲁古海MW7.1地震海啸的初始位移场(图6a)。 利用COMCOT模式非线性浅水方程模拟海啸的3h传播过程, 分析海啸波的能量分布以及近岸潮位站的实际观测和理论合成海啸波形的特征。

图 6 马鲁古海MW7.1地震的海啸数值模拟结果
a 海啸初始位移场; b 最大海啸波幅场Fig. 6 Tsunami simulation of Molucca Sea MW7.1 earthquake.

由图6a可知, 海啸初始位移场以沉降-抬升-沉降相间为特征, 以垂直抬升为主, 最大抬升位移为0.15m。 从最大海啸波幅场分布(图6b)来看, 最大海啸波幅主要集中分布在震源近岸区域, 海啸波以扩散的方式向各个方向传播。 由于实际发震断层的走向以SW为主, 海啸波的能量主要向断层两侧方向传播扩散, 部分能量沿断层错动方向传播至印尼苏拉群岛地区。

在近岸潮位站记录的波形序列中, 除海啸引起的海平面变化外, 还包括潮汐、 异常扰动等波动成分。 为了更好地识别海啸波形的时间序列特征, 我们采用2~120min的带通滤波器对水位波形进行滤波, 并去均值和倾斜趋势影响。 监测数据显示, 马鲁古海MW7.1地震发生后, 在震源附近区域引发了小规模的局地海啸。 其中, 位于震中西侧的Bitung潮位站在地震发生后18min首先记录到海啸波, 随后震中东侧的Ternate、 Jailolo和Tidore潮位站也相继监测到海啸波, 最大海啸波达0.15m。

为定量化评估海啸数值模拟效果及其可靠性, 我们采用方差减小量VR(Variance reduction)值表示潮位站观测波形和理论波形的拟合程度:

VR=1-∑(o-s)∑o22(7)

式中, o和s分别表示潮位站的观测波形和理论波形。 VR值越大则两者残差越小、 波形的拟合程度越高。

当海啸波传播至近岸边并在近岸爬高时, 由于水深变浅, 受近海水深、 海岸局地地形等影响, 其波动的非线性作用将增强, 而这往往会导致模拟波形和观测波形之间存在显著差异, 特别是更难以对后续的海啸波形序列进行准确模拟。 因此, 我们仅利用模拟和记录海啸波形的第一波周期, 对各潮位站记录的波形与理论波形进行比较, 并计算其所对应的VR值, 如图7b所示。 通过对比近岸潮位站记录到的观测波形与理论波形可知, Tidore潮位站的拟合程度最高(VR=0.83), Bitung潮位站的拟合程度较低(VR=-0.57)。 由图7b可以看出, 除Bitung潮位站外, 其他3个潮位站的VR值均在0.6以上, 其初至海啸波的到达时间和波幅都较为接近。 首波记录之后的后续实际波形振幅随时间增加逐渐衰减, 而模拟结果的衰减效应却不明显。 另外, 4个潮位站的后续观测波形均具有明显的高频成分, 与理论波形序列的形态差异较大。 究其原因, 可能受以下几方面影响: 1)模拟区域的海底地形和水深数据精度有限, 并不能完全真实反映地形及海岸线特征, 从而无法真实描述海啸波与近岸地形的相互作用; 2)在计算海底同震垂直位移形变时, 仅采用简化的均匀滑动位错模型, 没有考虑实际断层破裂的复杂性, 包括断层的非平面性、 滑移分布的非均匀性和方向性以及震源的动态破裂过程; 3)此次地震所处的地理位置特殊, 位于马鲁古海浅层地壳的增生杂岩体中, 增生杂岩主要由弧-弧碰撞刮削下来的海底沉积物和洋壳碎片堆积而成, 地震可能诱发海底滑坡, 导致海面出现剧烈变动, 从而产生高频的海啸波成分; 4)除此之外, 地球自转产生的科氏力、 频散效应、 底摩擦项的设置也是影响模拟结果的重要因素。

图 7 震源周边最大海啸波高(a)、 潮位站(红色三角形)观测波形(黑色实线)和理论波形(红色实线)对比(b)Fig. 7 Maximum tsunami height around the focal area and waveform comparison of the observed (black lines) recorded by tide gauges (red triangles) and its synthetic (red lines).

值得注意的是, 与其他3个潮位站相比, Bitung潮位站所处地理位置较为特殊(图7a)。 该观测站点位于印尼伦贝岛西侧, 紧邻伦贝海峡, 两岸呈喇叭口状向外敞开, 内窄外宽, 与马鲁古海相连。 当海啸波进入到喇叭口形海湾中时, 受特殊复杂地形的影响, 发生了持续的发射或干涉, 对海啸有明显放大作用, 从而导致Bitung潮位站记录的实际观测波幅和周期都明显大于模拟计算结果。

3 讨论与结论

本文使用IRIS提供的全球地震台网记录的数字波形资料, 采用W震相矩张量解反演方法快速得到MW7.1马鲁古海地震的相关震源参数。 反演结果表明, 此次地震是一次以逆断型为主的事件, 矩心深度约为35.5km, 断层面节面I的走向、 倾角和滑动角分别为225.5°、 48.5°和109.9°; 节面Ⅱ 的走向、 倾角和滑动角分别为16.8°、 45.2°和69°, 矩震级为MW7.1。

结合震源区的区域构造背景, 我们推断节面I为此次地震的发震断层面, 即NNE向的高倾角断层面, 并认为此次地震的发生与弧-弧碰撞的内部变形有关, 推测该地震可能受欧亚板块和菲律宾海板块相互俯冲作用影响, 由桑义赫弧和哈马黑拉弧逆冲碰撞产生的挤压变形所致。 在高倾角断层近垂直倾滑的逆冲过程中, 海底地形发生同震垂直位移, 使得震源上方的水体突然抬升, 从而产生了小规模的局地海啸。 潮位站波形序列记录到明显的高频海啸波成分, 不排除震源区引发海底滑坡而产生海啸, 本文中暂不做讨论。

文中基于反演获得的地震断层面几何参数构建断层模型, 采用COMCOT模型进行海啸数值模拟。 研究结果表明, 海啸波能量主要向断层两侧传播扩散, 部分能量沿断层错动方向传播至印尼苏拉群岛地区, 海啸的最大波幅主要集中分布在震源近岸区域, 震源周边大部分潮位站记录的海啸首波理论波形和观测波形基本相似, 但后续波形相差较大。 海啸数值模拟结果表明, 结合W震相快速反演的地震断层面几何参数和数值模式能够有效地预测海啸地震产生的海啸波高, 可为海啸预警提供定量化的数值参考。

本次地震的发生反映了马鲁古海地区的主要构造活动过程, 印尼马鲁古海地区位于一个NW-SE向挤压的区域应力场中, 此次地震应是整体应力场作用的产物。 印尼马鲁古海MW7.1主震为逆断性质, 据此认为此次地震是在板块俯冲作用下浅层海洋地壳发生破裂所致。 由于受欧亚板块和太平洋板块相向俯冲作用影响, 该区域弧-弧碰撞的构造运动特征不会发生本质变化, 其东、 西两侧有可能形成新的应力集中区, 未来还具有发生大地震的可能性。 因此, 考虑到海底大地震可能引发的海啸对震源区附近沿岸的影响, 应该进一步加强该区域范围内海啸灾害风险评估和灾害防范能力。

本研究求取了此次马鲁古海地震的震源机制解, 并探讨了其发震构造特征, 进而采用实际震源断层面的几何参数进行了海啸数值模拟, 所得成果可为理解大地震的发震机制、 孕育机理和海啸诱发机理提供科学参考。 海底大地震的震源机制解反演不仅有助于认识发震构造、 震源物理特征, 而且对评估大地震是否能引发海啸灾害潜势具有重要的现实意义。 基于本研究探讨的方法, 在海底大地震发生后, 可根据断层面快速反演结果和已有的地震构造背景, 识别真实发震断层面, 应用海啸模型进行定量化的海啸数值计算, 对海啸可能到达的地区及其波高进行预测, 准确地描述海啸波时空传播过程和海啸灾害分布, 并对可能遭受海啸灾害威胁的区域发布海啸预警, 从而降低海啸对人类造成的损失。

致谢 本研究使用了Duputel等提供的W震相矩张量反演程序和美国Cornell大学开发的COMCOT模型; 所用资料包括美国地震学研究联合会数据管理中心(IRIS/DMC)提供的全球地震台网(GSN)波形数据和印尼地理信息局(Geospatial Information Agency, BIG)提供的印尼沿岸验潮站水位资料; 文中图件使用GMT绘图软件包制作。 在此一并表示感谢!

参考文献
[1] 郭志, 陈立春, 李通, 等. 2018. 2017年8月8日四川九寨沟 M7. 0和9日新疆精河 M6. 6地震震源机制解[J]. 地震地质, 40(6): 119—129. doi: DOI: 103969/j. issn. 0253-4967. 2018. 06. 007.
GUO Zhi, CHEN Li-chun, LI Tong, et al. 2018. Focal mechanism of the August 8, M7. 0 Sichuan Jiuzhaigou and August 9, M6. 6 Xinjiang Jinghe earthquakes of 2017[J]. Seismology and Geology, 40(6): 119—129(in Chinese). [本文引用:1]
[2] 梁姗姗, 雷建设, 徐志国, 等. 2018. 2017年四川九寨沟 MS7. 0强震的余震重定位及主震矩张量解反演[J]. 地球物理学报, 61(5): 2163—2175.
LIANG Shan-shan, LEI Jian-she, XU Zhi-guo, et al. 2018. Relocation of aftershocks of the 2017 Jiuzhaigou, Sichuan, MS7. 0 earthquake and inversion for focal mechanism of the mainshock[J]. Chinese Journal of Geophysics, 61(5): 2163—2175(in Chinese). [本文引用:1]
[3] 潘桂棠, 王立全, 尹福光, 等. 2004. 从多岛弧盆系研究实践看板块构造登陆的魅力[J]. 中国地质, 23(9): 933—939.
PAN Gui-tang, WANG Li-quan, YIN Fu-guang, et al. 2004. Charm of land ing of plate tectonics on the continent as viewed from the study of the archipelagic arc-basin system[J]. Geological Bulletin of China, 23(9): 933—939(in Chinese). [本文引用:1]
[4] Bird P. 2003. An updated digital model of plate boundaries[J]. Geochemistry, Geophysics, Geosystems, 4(3): 1—52. doi: DOI:10.1029/2001GC000252. [本文引用:1]
[5] Blaser L, Krüger F, Ohrnberger M, et al. 2010. Scaling relations of earthquake source parameter estimates with special focus on subduction environment[J]. Bulletin of the Seismological Society of America, 100(6): 2914—2926. [本文引用:2]
[6] Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model[J]. Physics of the Earth and Planetary Interiors, 25(4): 297—356. [本文引用:1]
[7] Hall R. 1987. Plate boundary evolution in the Halmahera region, Indonesia[J]. Tectonophysics, 144(4): 337—352. [本文引用:1]
[8] Hall R, Nichols G J. 1990. Terrane amalgamation in the Philippine Sea margin[J]. Tectonophysics, 181(1-4): 207—222. [本文引用:1]
[9] Hayes G P, Moore G L, Portner D E, et al. 2018. Slab2, a comprehensive subduction zone geometry model[J]. Science, 362(6410): 58—61. https://www.nature.com/articles/s41598-019-47617-3-ref-link-section-d5960e1242. [本文引用:1]
[10] Hayes G P, Rivera L, Kanamori H. 2009. Source inversion of the W-phase: Real-time implementation and extension to low magnitudes[J]. Seismological Research Letters, 80(5): 817—822. [本文引用:1]
[11] Hinschberger F, Malod J A, Réhault J P, et al. 2005. Late Cenozoic geodynamic evolution of eastern Indonesia[J]. Tectonophysics, 404(1-2): 91—118. [本文引用:1]
[12] Kagan Y Y. 1991. 3-D rotation of double-couple earthquake sources[J]. Geophysical Journal International, 106(3): 709—716. [本文引用:1]
[13] Kanamori H, Rivera L. 2008. Source inversion of W phase: Speeding up seismic tsunami warning[J]. Geophysical Journal International, 175(1): 222—238. [本文引用:2]
[14] Liu L F, Cho Y S, Briggs M J, et al. 1995. Runup of solitary waves on a circular Island [J]. Journal of Fluid Mechanics, 32: 259—285. [本文引用:1]
[15] McCaffrey R. 1982. Lithospheric deformation within the Molucca Sea arc-arc collision: Evidence from shallow and intermediate earthquake activity[J]. Journal of Geophysical Research: Solid Earth, 87(B5): 3663—3678. [本文引用:2]
[16] Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 75(4): 1135—1154. [本文引用:1]
[17] Silver E A, Moore J C. 1978. The Molucca Sea collision zone, Indonesia[J]. Journal of Geophysical Research: Solid Earth, 83(B4): 1681—1691. [本文引用:2]
[18] Wang X, Liu P L F. 2006. An analysis of 2004 Sumatra earthquake fault plane mechanisms and Indian Ocean tsunami[J]. Journal of Hydraulic Research, 44(2): 147—154. [本文引用:1]
[19] Widiwijayanti C, Tiberi C, Deplus C, et al. 2004. Geodynamic evolution of the northern Molucca Sea area(Eastern Indonesia)constrained by 3-D gravity field inversion[J]. Tectonophysics, 386(3-4): 203—222. [本文引用:1]
[20] Woodhouse J H. 1988. The calculation of eigenfrequencies and eigenfunctions of the free oscillations of the Earth and the Sun [A]// Doornbos D J(ed). Seismological Algorithms, Computational Methods and Computer Programs. Academic Press, London, UK: 321—370. [本文引用:1]
[21] Zacharie D, Luis R, Hiroo K, et al. 2012. W phase source inversion for moderate to large earthquakes(1990—2010)[J]. Geophysical Journal International, 189(2): 1125—1147. [本文引用:1]
[22] Zhang Q, Guo F, Zhao L, et al. 2017. Geodynamics of divergent double subduction: 3-D numerical modeling of a Cenozoic example in the Molucca Sea region, Indonesia[J]. Journal of Geophysical Research: Solid Earth, 122(5): 3977—3998. doi: DOI:10.1002/2017JB013991. [本文引用:2]