2016年全球地震海啸监测预警与数值模拟研究

任智源, 原 野, 3, 赵联大, 3, 王培涛, 3, 侯京明, 3, 徐志国, 3, 高 义, 3, 李宏伟, 闪 迪, 王君成, 范婷婷, 王宗辰



2016年全球地震海啸监测预警与数值模拟研究

任智源1, 2, 原 野1, 2, 3, 赵联大1, 2, 3, 王培涛1, 2, 3, 侯京明1, 2, 3, 徐志国1, 2, 3, 高 义1, 2, 3, 李宏伟1, 2, 闪 迪1, 2, 王君成1, 2, 范婷婷1, 2, 王宗辰1, 2

(1. 国家海洋环境预报中心, 北京100081; 2. 国家海洋局海啸预警中心, 北京 100081; 3. 国家海洋局海洋灾害预报技术研究重点实验室, 北京 100081)

回顾了国家海洋环境预报中心(国家海洋局海啸预警中心) 2016年全球地震海啸监测预警的总体状况, 并基于震源生成模型和海啸传播数值模型的计算结果详细介绍了几次主要海啸事件及其影响特性。2016年全年国家海洋环境预报中心总共对全球6.5级(中国近海5.5级)以上海底地震响应了45次,发布海啸信息81期, 没有发生对我国有明显影响的海啸。结合精细化的数值模拟结果和浮标监测数据, 重点介绍了苏门达腊7.8级地震海啸、厄瓜多尔7.8级地震海啸、新西兰7.1级和7.8级地震海啸, 以及所罗门7.8级地震海啸的波动特征和传播规律, 模拟结果与实测海啸波符合较好。针对厄瓜多尔7.8级地震海啸事件, 本文比较分析了均匀断层模型和多源有限断层模型对模拟结果的影响; 针对新西兰7.1级地震海啸, 探讨了色散效应对海啸波在大水深、远距离传播过程的影响规律。

地震海啸; 海啸预警; 数值模拟; 浅水波方程; Boussinesq方程

海啸一般是由海底地震、海底滑坡及火山喷发等地质灾害所引起, 作为一种突发性、且破坏力极大的海洋灾害, 近几年来受到高度重视[1]。据统计, 超过90%的海啸是由地震引发[2-3], 且海底地震是可实时监测的, 因此对地震海啸的监测预警及研究成为了主要方面。2004年12月26日, 在印度尼西亚苏门答腊北部地区发生的震级为M=9.2级的地震并引发了特大海啸, 对印度洋周边国家造成了约30万的人员伤亡[4-6]。许多涉海国家开始重视海啸预警系统的建设与研究工作。而2010年智利的M=8.8级地震海啸和2011年日本东北M=9.0级地震海啸对近场区域所造成的严重灾害和人员伤亡对局地海啸的近岸预警提出了新的挑战[7]。

海啸的全生命周期包括了海啸的生成、传播, 以及爬高和淹没[8]。地震海啸的生成阶段主要是由于地震引发的海底的突然变形进而引起海水的大范围波动。Mansinha和Lmylie[9]基于弹性变形理论, 根据Volterra公式推导出了由板块运动所造成的海底变形量的计算方法, 被COMCOT模型采用。Okada[10-11]采用有限域积分模型由震源参数位错形态, 为海啸形成的初始形态提供海底边界的变化过程。在得到海底的变形之后, 再将海底地形的静态形变直接投影到水面作为初始形态处理[12]。该方法被称为瞬时响应模型。对于强度较大的地震, 板块运动形变历时较短, 该方法是适用的。对于运动时间较长的地震或者滑坡, 可以通过求解Cauchy-Poisson问题将海底变形由空间域的Fourier变换和时间域的Laplace变换得到水面位移和海底位移的关系, 即通过传递函数联系海底地形变化和水面波动,该方法可以称为动态响应模型[1]。Lynette和Liu[13]基于弱非线性Boussinesq方程研究了海底滑坡诱发的海啸。王本龙[14]根据Madsen等[15]提供的建议速度值, 改写了将Euler方程在自由面和底面的投影方程, 可以较好地模拟由非均匀海底引起的大范围非线性色散波浪的传播问题。赵曦[16]将该模型进一步推广为海底运动激发海啸的动态响应模型。在得到海啸的初始波动形态后, 海啸波将以自由表面重力波的形式向外传播。

目前, 海啸的数值传播模型主要有浅水波数值模型和Boussinesq色散模型。基于浅水波方程的数值模型主要有MOST[17]、GeoClaw[18]、COMCOT[19]和TUNAMI[20]等。浅水波方程是从N-S(Navier-Stokes)方程出发, 忽略了速度在垂直方向的变化而采用平均速度得到的。因此, 当流体垂直分量的速度较小, 即动压相对较小而以静压为主的条件下, 浅水波方程是适用的。Boussinesq方程是将流场的垂直分布以多项式级数形式近似表示, 并将连续性方程和动量方程进行沿水深的垂向积分, 将空间三维流动问题简化为平面二维问题。Boussinesq方程考虑了波浪的非线性特性和垂直加速度所导致的非静压效应, 即色散效应, 主要模型有FUNWAVE[21]和CULWAVE[22]。作为二十世纪九十年代海岸工程水动力学的重要进展, 通过提高模型的色散性使得Boussinesq类方程的适用范围被拓展至深水区域, 并通过改善模型的色散性和非线性等特性进一步扩大了其应用的范围[8, 23-25]。此外, 在近些年得到发展的非静压模型也值得关注[26-27]。

针对波浪爬高问题, Carrier和Greenspan[28]最早将非线性浅水波方程经过无量纲化变化为一组线性方程, 给出了简谐波爬高问题的解析解。Synolakis[29]利用合成孤立波谱的方法得当了孤立波爬高的非线性解析解, 并发现线性理论和非线性理论得到的波浪爬高结果一致。Carrier[30]基于新的坐标系统, 将非线性浅水波方程转化为线性双曲方程, 推导出了任意形态的初始波动在斜坡爬高的解析解。Tadepalli和Synolakis[31-32]基于浅水方程研究了N波(由两个孤立波组成, 类似字母“N”)在斜坡上的最大爬高值, 发现波谷在前的N波比波峰在前的N波产生更大的爬高。Madsen和Schffer[33]给出了单峰波和N波更普遍的表达式, 并通过非线性浅水方程研究了它们在斜坡上的爬高。Zhao等[34]研究了单峰波和N波在在斜坡上的爬高特性, 并进一步揭示了在爬高过程中的能量转化规律。

目前, 针对地震海啸的预警技术主要基于地震实时监测、浮标实测数据监测、快速数值模拟技术, 以及海啸预警数据库。Satake[35]最早提出了海啸数据库的构建方法, 并通过格林函数反问题方法计算断层滑移量。Wei等[36]提出了基于最小二乘法的反问题预报模型, 并应用于分析阿拉斯加潜在震源对夏威夷的影响分析。Liu等[37]将南海马尼拉潜在断层划分为39个单元板块, 采用COMCOT 模型建立了海啸预警系统。Ren等[38]将该方法拓展为多个浮标共同参与反演, 提高了预警的准确性, 并提出了一种优化布置浮标的方法。近几年发生的几次大海啸, 均被经过相应海面的卫星高度计捕获到, 发展卫星实时监测技术也为海啸预警提供了有益的参考[39-41]。此外, Wang和Liu[42-43]根据海啸数值模拟和基于海啸波影响的地磁场麦克斯韦方程, 建立了根据局部地磁场异常变化探测海啸特性的方法, 是近年来海啸预警研究的一个亮点。

我国的海啸警报信息是由国家海洋环境预报中心(国家海洋局海啸预警中心)发布, 所采用的新一代海啸预警系统包括了全球地震监测系统、60个海啸浮标与800余个潮位站组成的水位实测系统、基于高性能计算机的数值模拟系统, 以及警报信息发布平台。2016年全年, 国家海洋环境预报中心共针对45次地震发布了81次海啸信息(参见表1), 其中震级较大且监测到明显海啸波的分别是苏门答腊7.8级地震海啸、厄瓜多尔7.8级地震海啸、新西兰7.1级和7.8级地震海啸, 以及所罗门7.8级地震海啸。震源位置如图1所示, 其中全球板块边界采用Bird[44]提出的模型PB2002, 图中显示了18个较大板块的边界。可以发现, 环太平洋地震带依然是最活跃的地震带。在太平洋板块与周边其他板块的相互作用下, 发生了多次海底地震, 尤其是太平洋板块的西边界与欧亚板块、鄂霍次克板块、菲律宾板块、澳大利亚板块的相互作用非常频繁。此外, 菲律宾板块与琉球板块、扬子板块, 以及巽他板块之间的作用也比较明显, 因此也有可能诱发对中国沿海有较大威胁的海啸。本文将基于数值模型与实测数据分别对2016年的典型海啸事件进行深入地分析。

表1 2010~2016年地震海啸预警响应频次(中国海洋灾害公报2010~2016)

AF: 非洲板块; AM: 阿默尔板块; AN: 南极板块; AR: 阿拉伯板块; AU: 澳大利亚板块; CA: 加勒比板块; CO: 科库斯板块; EU: 欧亚板块; IN: 印度板块; NZ: 纳斯卡板块; NA: 北美板块; OK: 鄂霍次克板块; PA: 太平洋板块; PS: 菲律宾板块; SA: 南美板块; SO: 索马里板块; SU: 巽他板块; YA: 扬子板块

AF: plates of Africa; AM: plates of Amur; AN: plates of Antarctica; AR: plates of Arabia; AU: plates of Australia; CA: plates of Caribbean; CO: plates of Cocos; EU: plates of Eurasia; IN: plates of India; NZ: plates of Nazca; NA: plates of North America; OK: plates of Okhotsk; PA: plates of Pacific; PS: plates of Philippines; SA: plates of South America; SO: plates of Somalia; SU: plates of Sunda; YA: plates of Yangtze

1 研究方法

1.1 实测数据

DART (Deep-ocean Assessment and Reporting of Tsunamis)是由美国太平洋海洋环境实验室(PMEL, Pacific Marine Environmental Laboratory)开发的一套用于全球海啸监测浮标系统。原始的DART实测数据包括了15 s、1 min、15 min等不同的采样频率的数据。在一次海啸事件中, 浮标将从15 min一次的采样频率转换为15 s一次的采样频率并持续一段事件, 随后再次转换为1 min一次的频率直至事件结束。此外, 全球潮位站的实测数据以及我国近岸的潮位站可以提供海啸到达近岸附近的实测数据。原始的水位数据中, 包含了潮位的水位波动数据和海啸波的水位波动序列, 但由于二者频率不同, 因此可以通过滤波方法将海啸波分离出来, 得到海啸波的实测波形时间序列。

1.2 数值模型

针对地震产生的海底地形形变, 采用弹性半空间位错理论Okada模型[10-11]。在各项同性的条件下, 根据弹性半空间位错理论, 某一矩形几何面发生滑移引发地表某点的动力响应所产生的位移, 与滑移面的滑动量成正比。其中, 比例系数由滑移面的深度、倾角、几何尺寸, 以及该点与滑移面的相对位置确定。

在得到海底地形初始形变之后, 再用水波模型模拟海啸波的传播特征。平面二维形式的非线性浅水波方程[18]如下:

(2)

(3)

其中,表示时间,()为总水深,()用来描述真实的地形,()和()分别表示和方向的垂向平均速度,表示重力加速度。和表示底摩擦项:

表示曼宁系数, 采用有限体积法对浅水方程进行数值离散, 加入限制器(limiter)使得数值精度达到二阶精度。

本文建立了基于Wei和Kirby[46]提出的完全非线性Boussinesq方程的海啸传播色散模型, 考虑了波浪的色散性、摩阻系数、边界吸收、移动海岸等附加特性。数值方法采取TVD(Total Variation Diminishing)有限体积方法, 对时间项采用四阶强稳定R-K(Runge-Kutta)方法[21]。该模式将方程进一步扩展, 可以根据需要改写成为Nwogu[47]扩展的Boussinesq方程、Peregrine[48]扩展的Boussinesq方程、弱非线性色散模型、非线性浅水方程等。为了模拟浅水处拍岸浪的动力过程, 带有涡度粘性的动力混合项加入方程中, 模拟波浪破碎引起的能量耗散。在球极坐标下的弱非线性Boussinesq方程为

(5)

(6)

(8)

这里,表示总水深,r表示地球半径, (u,v)表示经度和纬度方程的垂向速度,表示自由水面高度,表示底摩擦系数,和分别表示经度方向和纬度方向,表示科氏力。这里取–0.39, 对应。

该模型在模拟波浪浅化、破碎、波浪爬高以及近岸环流等方面均有较好表现。本文将基于浅水波方程和Boussinesq方程分别建立数值模型并针对海啸传播过程进行研究。

2 2016年典型地震海啸模拟研究

2.1 苏门达腊7.8级地震海啸

2016年3月2日12时49分48秒(UTC时间), 苏门达腊附近海域发生了M=7.8级地震海啸。震中位于94.2°E, 4.9°S, 距离苏门达腊西海岸约800 km, 震源深度达到24 km。根据美国地质勘探局(USGS, United States Geological Survey)的评估, 此次地震产生的地震矩达到6.0×1020N·m。震中距离印度-澳大利亚板块和巽他板块交界处约600 km。在该处, 印度-澳大利亚板块每年以55 mm的速度向巽他板块移动。根据USGS提供的震源参数, 此次地震属于走滑型地震, 即断层的上下盘沿着断层面走向相对滑动, 且以水平运动为主。因此, 这种形式的板块运动不易激发海啸。2012年4月11日在苏门达腊海域曾先后发生的8.6级和8.2级地震, 以及2000年6月18日的7.9级地震都属于走滑型地震。2012年苏门达腊8.6级地震, 在浮标DART-56001、DART-23401, 和DART-23227监测到的海啸波的最大波幅分别为1、7和8 cm。而近岸潮位站分别监测到了15 cm (Hanimadhoo)和40 cm(Sabang)的海啸波动。

根据DART深水浮标监测到的数据, 此次地震引发的海啸波幅在1 cm左右。本文结合了USGS和NOAA提供震源参数, 震源参数为: 走向角96°, 倾角84°, 滑动角–170°, 长度100 km, 宽度50 km, 滑移量3.0 m。根据Okada模型, 计算得到海底变形位移量, 进一步利用浅水波模型模拟海啸的传播过程。图2为数值模拟结果, 图2a为最大波幅分布和走时图(传播时间)的等值线。由于震级较小以及板块运动特性, 此次地震海啸相对较弱, 但影响范围依然较大, 海啸波影响到了印度洋的大部分地区, 北至孟加拉湾、东至澳大利亚、西到马达加斯加岛, 对阿拉伯海影响较小。海啸在生成之后约3 h抵达印度东南沿海, 在4 h左右到达澳大利亚。图2b图是数值模拟与实测的海啸波对比结果, 黑线表示实测结果, 红线表示模拟结果。本文选取了4个浮标分别为DART-56001 (109.944°E, 14.016°S)、DART-56003 (117.942°E, 15.023°S)、DART-23401 (88.537°E, 8.905°N),和DART-23227 (88.792°E, 6.255°N)。整体来看, 数值模拟结果对海啸波的首波比较准确。在地震初期, 4个深水浮标处均发生由于地震引起的DART浮标海底模块的震动而产生类似海啸波动的“假性海啸”, 此处应与真正的海啸波相区分开。

2.2 厄瓜多尔7.8级地震海啸

2016年4月16日23时58分48秒(UTC时间), 厄瓜多尔沿岸发生了M=7.8级地震。根据USGS测定的地震基本信息, 震中位于79.94°W, 0.37°N, 距离厄瓜多尔首都基多约170 km, 震源深度为20 km。此次地震产生的地震矩达到7.1×1020N·m。此次地震属于浅源逆冲型, 在该处纳斯卡板块以每年61 mm的速度向南美洲板块之下俯冲。这两个板块的相互作用也曾造成了1960年智利9.5级地震、2010年智利8.8级地震、2014年和2015年智利8.3级地震。

a: 最大波幅分布和走时图等值线分布; b: 海啸波实测与数值模拟结果对比

a: Distribution of the maximum tsunami wave amplitude and traveling time; b: Comparison of the tsunami wave series at DART locations between the measurement (black line) and numerical results (red line)

实时监测数据显示, 距离震中较近的La Libertad站监测到海啸波的首波约为0.12 m, 最大波幅为0.24 m。附近的浮标也捕获到了轻微的海水波动信息。本文根据USGS提供的震源参数(长度160 km, 宽度60 km, 走向角29°, 倾角15°, 滑移角123°), 基于Okada模型分别对均一板块和多个单元板块组成的有限断层模型(Finite fault model, 240个单元板块)反演得到的多板块震源模拟得到地震激发的海底地形变化, 作为海啸波传播的初始条件。进一步采用非线性浅水波模型模拟研究海啸波的传播演化特征。全场的地形和网格分辨率均为4 min, 近岸的地形和网格分辨率分别为0.5 min和0.25 min。

图3为数值模拟结果, 其中图3a为均一化震源模型计算得到的最大波幅分布和走时图的等值线。此次地震海啸主要对近场存在一定影响, 在厄瓜多尔附近海域, 包括巴拿马、哥伦比亚、哥斯达黎加、尼加拉瓜等均产生一定影响, 但并没有在近岸产生海啸灾害。尽管海啸波较小, 但对整个太平洋地区都有一定影响。太平洋水深比较均匀, 因此传播时间均匀分布, 抵达中国沿海至少需要22 h。图3b为数值模拟和监测得到的海啸波时间序列的比较图, 图中黑色线为实测值, 红色线和绿色线分别为均一化震源模型和有限断层多个单元板块模拟的结果。两个深水浮标32411(90.850°W, 4.995°N)和32413(93.477°W, 7.393°S)监测到海啸的最大波幅均不足1 cm。模拟的结果和实测海啸波在首波基本一致, 但后续的几个波存在一定偏差。在距离震中最近的DART-32067 (81.262°W, 0.641°N), 多板块震源可以模拟得到初始时刻达到5 cm的海啸波。在近岸的La Libertad站(89.319°W, 13.4851°N), 受限于地形分辨率和震源机制解的影响, 两种震源均可以模拟得到海啸波首波, 但无法精细描述海啸波在近岸的持续震荡过程。在夏威夷附近的DART51407没有监测到海啸波, 而此处的模拟结果显示仅有1 mm的海水波动。

a: 最大波幅分布和走时图等值线; b: 海啸波实测与数值模拟结果对比

a: Distribution of the maximum tsunami wave amplitude and traveling time; b: Comparison of the tsunami wave series at DART locations between the measurement (black line), numerical results of the uniform earthquake source (red line), and the finite fault model (green line)

2.3 新西兰7.1级地震海啸

2016年9月1日16时37分57秒(UTC时间), 新西兰附近海域发生了M=7.1级地震。震中位于179.055°E, 37.401°S, 位于新西兰的吉斯伯恩港东北部, 距离新西兰约50 km, 震源深度为20 km。此次地震产生的地震矩为4.4×1019N·m。此次地震属于浅源型地震, 主要是由于太平洋板块和澳大利亚板块相互作用而产生的。太平洋板块以47 mm/a的速度向澳大利亚板块之下俯冲, 使得这一区域发生6级以上地震超过28次, 包括了1995年2月一次7.1级地震。

新西兰的两个潮位站East Cape(178.159°E, 37.55°S)和Tauranga(176.18°E, 37.64°S)监测到了海啸波动, 最大波幅分别为23 cm和5 cm。采取USGS提供的震源参数: 长度为60 km, 宽度为40 km, 震源深度19 km, 走向角351°, 倾角26°, 滑动角–145°, 滑移量0.6 m。在震源附近选取了一个包含潮位站的小区域, 网格为0.5分, 采用非线性浅水方程进行模拟。图4a给出了在潮位站位置, 实测海啸波与数值模拟结果的比较。数值模拟结果与实测结果比较吻合。尽管在East Cape监测到的海啸波周期约为30 min。但由于近岸的反射所造成的边缘波及其俘获效应, 海啸在近岸的波动持续了4个小时以上。

图4b为浅水波方程模拟计算得到整个太平洋范围的能量分布和传播时间图。此次地震震级相对较小, 影响的范围也相对较小。但由于海底山脊的导向作用, 部分海啸波也传播到了南美洲和南极地区。震中位于南太平洋中部地区, 海啸传播至太平洋西海岸和美洲大陆均需要13 h左右。本文针对此次海啸事件, 建立了太平洋范围内的基于弱色散boussinesq方程模拟海底变形激发海水大范围扰动的传播模型。图4c为分别在考虑色散效应和不考虑色散效应时海啸波引起的最大波幅的差异。尽管海啸波相对较小, 但依然存在一定的色散效应, 与震源存在一定距离, 并沿海啸的最大能量方向分布于东太平洋。

2.4 新西兰7.8级地震海啸

新西兰再次于2016年11月13日11时02分59秒(UTC时间), 发生了M=7.8级地震。震中位于173.065°E, 42.725°S, 位于新西兰南岛的陆地上, 距离新西兰最近的城市安伯利约56 km, 震源深度为22 km(USGS)。此次地震产生的地震矩为7.4×1020N·m, 属于浅源型地震, 同样是是由于太平洋板块和澳大利亚板块相互作用而产生的。

新西兰的两个潮位站Kaikoura(173.7028°E, 42.4129°S)和Wellington(174.78°E, 41.28°S)监测到了最大波幅达到了233 cm和46 cm的海啸波。采取的震源参数为: 长度为100 km, 宽度为50 km, 震源深度15 km, 走向角219°, 倾角38°, 滑动角128°, 滑移量3.0 m。在震源附近选取了包含潮位站的小区域, 采用高精度网格和地形数据, 以及非线性浅水方程进行模拟。图5a给出了在潮位站位置, 实测结果与数值模拟结果的比较。数值模拟结果与实测结果比较吻合。海啸波在近岸周期约为30 min, 但由于近岸的边缘波及其俘获效应, 使得海啸在近岸的波动达到了4个小时。图5b是海啸波在太平洋传播过程中的最大波幅分布和走时图。可以发现, 此次海啸事件主要是对新西兰南部200 km范围以内有明显的影响, 但同时也对西南太平洋存在一定影响。根据海啸传播的走时图, 海啸约在14 h到达美洲, 18 h到达中国沿海。

2.5 所罗门群岛7.8级地震海啸

2016年12月8日17时38分46秒(UTC时间), 发生了M=7.8级地震。震中位于161.330°E, 10.676°S, 位于所罗门群岛南部海底中, 距离所罗门的基拉基拉(Kirakira)约69.1 km, 震源深度为41.0 km (USGS)。此次地震产生的地震矩为6.4×1020N·m, 属于浅源型地震。在震中位置, 澳大利亚板块以每年96 mm的速度向太平洋板块进行俯冲。根据USGS的震源机制解表明破裂发生在倾斜的逆断层上。

本文选取了接近震源的两个海啸浮标55012 (158.453°E, 15.664°S)和55023(153.537°E, 14.715°S)。板块长度为100 km, 宽度为50 km。采取USGS提供的震源参数: 震源深度41 km, 走向角290°, 倾角44°, 滑动角67°, 滑移量3.0 m。采用非线性浅水方程进行模拟。图6a给出了在浮标位置, 实测结果与数值模拟结果的比较, 相对比较吻合。两个浮标监测到的海啸波均不超过5 cm, 相对较远的浮标没有监测到明显的海啸。主要是由于震源位于所罗门群岛南部, 受周边地形的影响, 激发的海啸难以对较远的区域产生较大的影响, 但影响范围依然是整个太平洋区域。

3 讨论

2016年全年针对地震事件的响应频次与往年持平, 没有发生影响较大的地震海啸, 但依然有五次可以监测到明显的水面波动, 并且对整个太平洋或者印度洋都产生了一定的影响。近十多年来, 先后发生了2004年印度洋9.2级海啸、2010年智利8.8级地震海啸、2011年日本9.0级地震海啸、2014年8.2级和2015年8.3级智利地震海啸, 尤其是极端地震海啸均造成了不可忽视的人员伤亡与财产损失。因此, 对于极端特大地震海啸的监测预警必须保持高度关注。

对于一个典型的海啸过程, 震源的机制对海啸的传播和灾害影响有重要影响。对于今年苏门答腊7.8级地震, 震源机制属于走滑型地震, 断层运动以平移运动为主, 因而产生的海啸波较小。尽管地震所引起的板块运动及海底变形较为复杂, 利用反演得到多源有限断层震源机制解可以对震源的复杂性进行较为准确的描述。准确的复杂震源有助于计算得到更为精确的海啸波动[49]。但对于此次厄瓜多尔7.8级地震, 由于部分震源位于陆地上, 并且地形精度有限, 因此复杂高精度震源的优势并没有很好的体现出来。而大部分海啸事件所采用的均一化震源模型也可以计算得到比较准确的结果。

a: 近岸测点的数值模拟结果与实测结果比较; b: 最大波幅分布和走时图等值线; c: 海啸传播过程中的色散效应

a: Comparison of the tsunami wave series at coastal locations between the measurement and numerical results; b: Distribution of the maximum tsunami wave amplitude and traveling time; c: Difference of maximum amplitude between the shallow water equation and the Boussinesq equation

a: 近岸测点的数值模拟结果与实测结果比较; b: 最大波幅分布和走时图等值线

a: Comparison of the tsunami wave series at coastal locations between the measurement (black line) and numerical results (red line); b: Distribution of the maximum tsunami wave amplitude and traveling time

海啸波在近岸和大洋的传播过程中, 主要表现非线性和色散性[50-51]:

(10)

在公式10中色散性主要表现为频散特性的积分效应, 它与波长、传播距离和传播水深相关。传播距离越远、水深越大、或者波长越短, 色散性表现的越显著。假设在常水深条件下(4 000 m), 海啸传播过程中波长不变(100 km), 在传播距离为1 000 km和5 000 km时得到色散指数分别为0.096和0.48。而根据Glimsdal等[51]的结论, 而当大于0.1时, 色散性变得比较显著。也就是说, 对于太平洋或印度洋平均水深接近4 000 m的深水范围, 若初始波长为100 km时, 当海啸的传播距离达到1 000 km以上时, 色散效应逐渐表现出来, 与图4c基本符合。这也与2011年日本地震海啸所表现出的色散特性基本一致[21]。

a: 近岸测点的数值模拟结果与实测结果比较; b: 最大波幅分布和走时图等值线分布

a: Comparison of the tsunami wave series at coastal locations between the measurement (black line) and numerical results (red line); b: Distribution of the maximum tsunami wave amplitude and traveling time

此外, 对于考虑地球曲率的控制方程需要考虑科氏力项的影响。而海啸在深水中的传播过程中的能量主要通过波高表现为势能, 而通过流场所表现的动能相对极小, 因此科氏力对结果的影响也可以忽略[18]。海啸在传播过程中, 底摩擦也存在一定影响, 其量级为, 因此对在深水中传播影响很小。但当水深小于100 m, 以及海啸波的爬高和淹没过程中, 海底糙率对结果有较大的影响[52]。

4 结论

本文介绍了2016年全年针对地震海啸的预警的总体情况, 建立了基于震源生成模型和基于浅水波方程及Boussinesq方程的海啸传播模型的数值模拟手段, 并结合实测数据详细介绍了苏门答腊7.8级地震海啸、厄瓜多尔7.8级地震海啸、新西兰7.1级和7.8级地震海啸, 和所罗门群岛7.8级地震海啸的近岸波动规律与在大洋中的传播特征。尽管2016年未发生类似311日本9.0级地震海啸这样极端性的海啸事件, 也未发生8.0级以上的地震海啸, 但7.8级的地震海啸达到了4次, 总体数量偏多, 表明全球板块的活跃性有一定加剧。从数值模拟结果和灾后调查来看, 这些海啸普遍对近岸造成了一定的影响, 使民众产生了恐慌情绪, 并对整个大洋产生了较大范围的影响。总体而言, 地震海啸作为一种传播速度快、释放能量高、影响范围广、物理机制复杂的极端海洋灾害, 需要保持持续地关注与深入地研究。

致谢: 感谢上海交通大学何友声教授和加州理工学院吴耀祖教授(Theodore Yao-Tsu Wu)的鼓励与支持。对上海交通大学刘桦教授和丹麦科技大学Per A. Madsen教授的指导和帮助表示衷心感谢。本文得到了国家重点研发计划(2016YFC1401500), 国家海洋公益性行业科研专项项目(201405026), 海洋工程国家重点实验室开放课题(1604)和留学人员科技活动项目择优资助的支持。

[1] 任智源. 南海海啸数值模拟研究[D].上海: 上海交通大学, 2015. Ren Zhiyuan. Numerical simulation of tsunami in South China Sea[D]. Shanghai: Shanghai Jiao Tong University, 2015.

[2] 王培涛, 于福江, 赵联大, 等. 越洋海啸的数值模拟及其对我国的影响分析[J]. 海洋学报, 2012, 34(2): 39-47. Wang Peitao, Yu Fujiang, Zhao Lianda, et al. Numerical simulation of trans-oceanic tsunami and its impact analysis on Chinese coasts[J]. Acta Oceanologica Sinica, 2012, 34(2): 39-47.

[3] 叶琳, 于福江, 吴玮. 我国海啸灾害及预警现状与建议[J]. 海洋预报, 2005, 22(S1): 147-157. Ye Lin, Yu Fujiang, Wu Wei. The disaster and warning of tsunami in China and the suggestion in future[J]. Marine Forecasts, 2005, 22(S1): 147-157.

[4] Titov V, Rabinovich A B, Mofjeld H O, et al. The Global Reach of the 26 December 2004 Sumatra Tsunami[J]. Science, 2005, 309(5743): 2045-2048.

[5] Grilli S T, Ioualalen M, Asavanant J, et al. Source constraints and model simulation of the December 26, 2004, Indian Ocean Tsunami[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2007, 133(6): 414-428.

[6] Wijetunge J J, Wang X, Liu P L-F. Indian Ocean Tsunami on 26 December 2004: numerical modeling of inundation in three cities on the south coast of Sri Lanka[J]. Journal of Earthquake and Tsunami, 2008, 2(2): 133-155.

[7] Melgar D, Allen R M, Riquelme S, et al. Local tsunami warnings: Perspectives from recent large events[J]. Geophysical Research Letters, 2016, 43(3): 1109-1117.

[8] 刘桦, 赵曦, 王本龙, 等. 海啸数值模拟与南海海啸预警方法[J]. 力学季刊, 2015, 3: 351-369. Liu Hua, Zhao Xi, Wang Benlong. Numerical simulation of tsunami and tsunami warning methods for South China Sea region[J]. Chinese Quarterly of Mechanics, 2015, 3: 351-369.

[9] Mansinha L, Smylie D E. The displacement fields of inclined faults[J]. Bulletin of the Seismological Society of America, 1971, 61(5): 1433-1440.

[10] Okada Y. Surface deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154.

[11] Okada Y. Internal deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 1992, 82(2): 1018-1040.

[12] Dao M H, Tkalich P. Tsunami propagation modelling-a sensitivity study[J]. Natural Hazards and Earth System Science, 2007, 7(6): 741-754.

[13] Lynett P, Liu P L F. A numerical study of submarine–landslide–generated waves and run–up[J]. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 2002, 458(2028): 2885-2910.

[14] 王本龙. 基于高阶Boussinesq方程的海岸破波带数学模型研究[D]. 上海: 上海交通大学, 2005. Wang Benlong. Mathematical modelling of surf zone based on the high order Boussinesq equations[D]. Shanghai: Shanghai Jiao Tong University, 2005.

[15] Madsen P A, Bingham H B, Liu H. A new Boussinesq method for fully nonlinear waves from shallow to deep water[J]. Journal of Fluid Mechanics, 2002, 462: 1-30.

[16] 赵曦. 海啸波生成、传播与爬高的数值模拟[D]. 上海: 上海交通大学, 2011. Zhao Xi. Numerical simulation of generation, propagation and runup of tsunamis[D]. Shanghai: Shanghai Jiao Tong University, 2011.

[17] Titov V V, Synolakis C E. Numerical modeling of tidal wave runup[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1998, 124(4): 157-171.

[18] Leveque R J. Tsunami modelling with adaptively refined finite volume methods[J]. Acta Numerica, 2011, 20: 211-289.

[19] 任智源, 赵曦, 刘桦. 海啸预警研究进展与南海海啸预警方法研究[C]//吴有生. 第二十五届全国水动力学研讨会暨第十二届全国水动力学学术会议文集. 北京: 海洋出版社, 2013: 430-438. Ren Zhiyuan, Zhao Xi, Liu Hua. Advances in tsunami warning technology and tsunami warning method in South China Sea[C]// Wu Yousheng. Article Collection of 25th National Symposium on Hydrodynamics. Beijing: China Ocean Press, 2013: 430-438.

[20] Imamura F, Shuto N, Goto C. Numerical simulations of the transoceanic propagation of tsunamis[C]//IAHR, Sixth Congress of the Asian and Pacific Regional Division of the International Association for Hydraulic Research proceedings. Kyoto: Local organizing committee of the congress, 1988: 265-272.

[21] Kirby J T, Shi F, Tehranirad B, et al. Dispersive tsunami waves in the ocean: Model equations and sensitivity to dispersion and Coriolis effects[J]. Ocean Modelling, 2013, 62(62): 39-55.

[22] Lynett P J. Nearshore wave modeling with high-order Boussinesq-type equations[J]. Journal of Waterway, Port,Coastal, and Ocean engineering, 2006, 132(5): 348-357.

[23] Madsen P A, Schäffer H A. Higher–order Boussinesq–type equations for surface gravity waves: derivation and analysis[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1998, 356(1749): 3123-3181.

[24] Gobbi M F, Kirby J T, Wei G. A fully nonlinear Boussinesq model for surface, waves. Part 2. Extension to O(kh)4[J]. Journal of Fluid Mechanics, 2000, 405(4): 181-210.

[25] Madsen P A, Bingham, Liu H. A new Boussinesq method for fully nonlinear waves from shallow to deep water[J]. Journal of Fluid Mechanics, 2002, 462(1): 1-30.

[26] Bai Y, Cheung K F. Dispersion and nonlinearity of multi-layer non-hydrostatic free-surface flow[J]. Journal of Fluid Mechanics, 2013, 726: 226-260.

[27] Fang K, Liu Z, Zou Z. Modelling coastal water waves using a depth-integrated, non-hydrostatic model with shock-capturing ability[J]. Journal of Hydraulic Research, 2015, 53(1): 119-133.

[28] Carrier G F, Greenspan H P. Water waves of finite amplitude on a sloping beach[J]. Journal of Fluid Mechanics, 1958, 4(4): 97-109.

[29] Synolakis C E. The runup of solitary waves[J]. Journal of Fluid Mechanics, 1987, 185: 523-545.

[30] Carrier G F, Wu T T, Yeh H. Tsunami run-up and draw-down on a plane beach[J]. Journal of Fluid Mechanics, 2003, 475: 79-99.

[31] Tadepalli S, Synolakis C E. The run-up of N-waves on sloping beaches[J]. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1994, 445(1923): 99-112.

[32] Tadepalli S, Synolakis C E. Model for the leading waves of tsunamis[J]. Physical Review Letters, 1996, 77(10): 2141-2144.

[33] Madsen P A, Schaeffer H A. Analytical solutions for tsunami runup on a plane beach: single waves, N-waves and transient waves[J]. Journal of Fluid Mechanics, 2010, 645: 27-57.

[34] Zhao X, Wang B, Liu H. Characteristics of tsunami motion and energy budget during runup and rundown processes over a plane beach[J]. Physics of Fluids, 2012, 24(6): 062107.

[35] Satake K. Inversion of tsunami waveforms for the estimation of a fault heterogeneity: Method and numerical experiments[J]. Journal of Physics of the Earth, 1987, 35(3): 241-254.

[36] Wei Y, Cheung K F, Curtis G D, et al. Inverse Algorithm for Tsunami Forecasts[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2003, 129(129): 60-69.

[37] Liu P L-F, Wang X, Salisbury A J. Tsunami hazard and early warning system in South China Sea[J]. Journal of Asian Earth Science, 2009, 36(1): 2-12.

[38] Ren Z Y, Liu H, Wang B L, et al. An investigation on multi-buoy inversion method for Tsunami Warning System in South China Sea[J]. Journal of Earthquake and Tsunami, 2014, 8(3): 1440004.

[39] Okal E A, Piatanesi A, Heinrich P. Tsunami detection by satellite altimetry[J]. Journal of Geophysical Research, 1999, 104(B1): 599-615.

[40] Smith W, Scharroo R, Titov V, et al. Satellite altimeters measure tsunami[J]. Oceanography, 2015, 18(2): 11-13.

[41] Hamlington B D, Leben R R, Godin O A, et al. Could satellite altimetry have improved early detection and warning of the 2011 Tohoku tsunami?[J]. Geophysical Research Letters, 2012, 39(15): L15605.

[42] Wang B, Liu H. Kinematic dynamo by large scale tsunami waves in open ocean[J]. Theoretical and Applied Mechanics Letters, 2013, 3(3): 032003.

[43] Wang B, Liu H. Space-time behaviour of magnetic anomalies induced by tsunami waves in open ocean[J]. Proceedings of the Royal Society A Mathematical Physical & Engineering Sciences, 2013, 469(2157): 20130038.

[44] Bird P. An updated digital model of plate boundaries[J]. Geochemistry, Geophysics, Geosystems, 2003, 4(3): 1027.

[45] Rabinovich A B, Candella R N, Thomson R E. The open ocean energy decay of three recent trans‐Pacific tsunamis[J]. Geophysical Research Letters, 2013, 40(12): 3157-3162.

[46] Wei G, Kirby J T. Time-dependent numerical code for extended Boussinesq equations[J]. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1995, 121(5): 251-261.

[47] Nwogu O. Alternative form of Boussinesq equations for nearshore wave propagation[J]. Journal of Waterway, Port, Coastal, and Ocean engineering, 1993, 119(6): 618-638.

[48] Peregrine D H. Long waves on a beach[J]. Journal of Fluid Mechanics, 1967, 27(4): 815-827.

[49] Ulutas E. Comparison of the seafloor displacement from uniform and non-uniform slip models on tsunami simulation of the 2011 Tohoku–Oki earthquake[J]. Journal of Asian Earth Sciences, 2013, 62(30): 568-585.

[50] 吴耀祖. 水波动力学研究进展[J]. 力学进展, 2001, 31(3): 327-343. Wu T Y. Advances in water wave mechanics[J]. Advances in Mechanics, 2001, 31(3): 327-343.

[51] Glimsdal S, Harbitz C B, Løvholt F. Dispersion of tsunamis: does it really matter?[J]. Natural Hazards and Earth System Sciences, 2013, 13(6): 1507.

[52] Ren Z Y, Wang B L, Fan T T, et al. Numerical analysis of impacts of 2011 Japan Tohoku tsunami on China Coast[J]. Journal of Hydrodynamics, 2013, 25(4): 580- 590.

Monitoring, early warning and numerical study of global tsunamis in 2016

REN Zhi-yuan1, 2, YUAN Ye1, 2, 3, ZHAO Lian-da1, 2, 3, WANG Pei-tao1, 2, 3, HOU Jing-ming1, 2, 3, XU Zhi-guo1, 2, 3, GAO Yi1, 2, 3, LI Hong-wei1, 2, SHAN Di1, 2, WANG Jun-cheng1, 2, FAN Ting-ting1, 2, WANG Zong-chen1, 2

(1. National Marine Environmental Forecasting Center, Beijing 100081, China; 2. State Oceanic AdministrationTsunami Warning Center, Beijing 100081, China; 3. Key Laboratory of Research on Marine Hazards Forecasting, National Marine Environmental Forecasting Center, Beijing 100081, China)

This study reviewed global earthquake tsunami warnings in 2016, monitored by the National Marine Environmental Forecasting Center (National Tsunami Warning Center, SOA), and introduced the five typical tsunami events-basedon tsunami generation and propagation numerical models. We have responded from 45 submarine earthquakes, above a magnitude of 6.5 (M> 5.5 in the coastal region of China), and released 81 times of earthquake information in 2016. No tsunami affected Chinese Coasts. This study has focused on the Sumatran tsunami ofM7.8, the Ecuadorian tsunami ofM7.8, the New Zealand tsunamis ofM7.1 and 7.8, and the Solomon tsunami ofM7.8, and presented maximum tsunami wave distribution and wave fluctuation data based on refined numerical results and measurements. The numerical results could reproduce the tsunami scenarios and match the measured data well. Comparison between the uniform fault model and the finite fault model with the multi-plate for the Ecuadorian tsunami ofM7.8 is presented. The dispersion effect on tsunami propagation is discussed based on the simulation results of the New Zealand tsunami ofM7.1 with Boussinesq equation.

Tsunami; Tsunami early warning; numerical modeling; shallow water equation; Boussinesq equation

(本文编辑: 李晓燕)

Feb., 14, 2017

[The Public Science and Technology Research Funds Projects of Ocean, No. 201405026; the National Key Research and Development Plan, No. 2016YFC1401500; Opening fund of State Key Laboratory of Ocean Engineering, No. 1604; Technology Foundation for Selected Overseas Chinese Scholar]

O352

A

1000-3096(2017)06-0098-13

10.11759/hykx20170214002

2017-02-14;

2017-04-06

国家重点研发计划(2016YFC1401500); 国家海洋公益性行业科研专项项目(201405026); 海洋工程国家重点实验室开放课题(1604); 留学人员科技活动项目择优资助

任智源(1986-), 男, 陕西西安人, 助理研究员, 博士, 研究方向为水动力学和海岸工程, 电话: 010-82102752, E-mail: zhyren@foxmail.com