香港海啸数值模型发展——以日本3.11海啸为主要案例对模型的验证

刘廸森， 莫庆炎

（香港天文台， 中国 香港）

0 引言

海啸可以由海底地震、 海底巨大山泥倾泻、海底火山爆发和陨石撞击海洋而引发， 其中海底地震是引发海啸最常见的原因， 占所有海啸个案90%以上。 过去10年， 全球发生的大海啸引起了世界关注， 如2004年12月印度尼西亚苏门答腊北部以西的印度洋海底9 级地震引发跨海洋海啸，及2011年3月在日本东北的太平洋发生9 级地震所产生的海啸均造成严重的人命伤亡及破坏。 香港位处珠江口， 面对南海， 来自太平洋的海啸有时会进入南海， 而在马尼拉海沟所产生的海啸也会对华南沿岸地区以至香港构成威胁。

香港天文台于2005年底通过政府间海洋学委员会（IOC）的海啸淹浸模型交换计划（TIME）引进了近场海啸数值模型， 利用南海及香港的水深和地形数据， 模拟海啸抵达香港时香港沿岸的海啸高度①黄永德， 杨继兴， 李健威. 香港海啸警报服务的新发展和数值模式的应用. 第二届粤港澳地震科技研讨会论文集.中国澳门， 2006.。 受当时技术限制， 该数值模型模拟海啸所需的运算长达数十小时， 只能用作建立预设地震情景之下的模拟海啸结果数据库， 在需要时提取最接近地震情景的海啸模拟数据以作参考， 而不能在发生地震后利用实时分析所得的地震参数进行运算和提供模拟海啸的结果。

随着海啸数值模型技术不断改进及计算机运算速度的不断提升， 香港天文台在2012年引进了一套可应用平行运算技术操作的海啸数值模型-Cornell Multigrid COupled Tsunami（COMCOT）model[1]。 该模型由美国康乃尔大学土木与环境工程学系刘立方教授所领导的团队研发， 大大缩短模拟海啸的运算时间， 并可与香港天文台在2010年开始运作的地震数据处理及分析系统结合， 实时把系统探测和分析所得的太平洋或南海地震参数， 输入模型以模拟海啸传播过程， 并计算海啸在海面上和抵岸时的情况， 为应急决策提供及时和贴近实际情景的模拟结果。 本文利用日本311地震所引发的海啸和过去香港曾经录得的海啸记录以验证COMCOT 模拟海啸的能力， 并探讨COMCOT 在香港天文台海啸预警工作上的应用。

1 COMCOT 海啸数值模型概述

COMCOT 是基于浅水波理论方程而发展的数值模型， 现时在全球多个海啸预报中心及学术机构的业务运作或学术研究上被广泛地应用。

COMCOT 采用平行运算技术以加快运算速度，在发生地震后可以利用实时地震参数进行运算，并在短时间内提供运算结果。 此外， COMCOT 可以根据波高与水深的比值以及模拟范围的大小，适当地选择使用球坐标系统或卡式坐标系统、 网格大小、 线性或非线性波方程进行模拟运算， 在基本上不影响计算准确度的情况下尽量缩短运算的时间。 在大范围运算的深水区中， 由于海啸在深水传播时， 网格分辨率对模式的准确度影响不大， 模式采用较大的网格及使用线性方程， 因应模拟范围大和需要考虑地球面弧度而采用球坐标系统。 到近岸及浅水地区， 因需要考虑海洋底部磨擦项的影响及模拟海啸造成的淹浸情况， 便采用非线性浅水波方程及高分辨率的网格以获取更详尽的海啸高度。 由于模拟范围缩小， 地球曲率效应影响相对地轻微， 模式改为采用卡式坐标系统以减少运算时间[2]。

网格方面， COMCOT 采用嵌套方法， 网格层的数目视乎最外层与最内层分辨率的差别而定。要确保模式运行的稳定性， 时间步长的选取、 嵌套网格区的设定和水深的变化必须配合， 因此依据目标地的地理情况融合水深和地形资料是成功运行模式的重要工作。 根据香港的地理情况， 天文台在南海制作了一套五层的嵌套网格（图1）， 最外层的网格覆盖南海（5° ～25°N 及105° ～125°E），而最内层的网格包括整个香港（22°10′ ～22°35′ N及113°50′ ～114°30′ E）， 各层网格的分辨率从内至外分别为1、 3、 9、 27 和81 弧秒（1 弧秒约等于30 m)。 各网格层之间的分辨率比例为1： 3。至于南海以外的太平洋则采用美国 NOAA 的National Geophysical Data Centre（NGDC）所提供的ETOP1 地形及水深数据[3]， 分辨率为2 弧分。

南海水深和地形资料以美国国家航空航天局（NASA）和美国影像与制图局（NIMA）合作通过航天飞机雷达地形测量任务制作的SRTM30 PLUS 数位高程资料（DEM）为基础[4]， 并在近岸处根据较精确的区域资料作出修定。 SRTM30 PLUS 的最高分辨率为30 弧秒， 即约1 km。 而最内层网格则采用香港地政总署的香港地形资料及香港海事处的香港境内水深资料， 最高分辨率达30 m。

图1 迭加五层嵌套网格Fig 1 Overlap of 5 layer nested grids

本文利用COMCOT 在模拟海啸对香港近岸的影响时， 网格分辨率达30 m 及90 m 的最内两层网格层将采用非线性方程运算， 而其余的网格层则采用线性浅方程。 一般在水深少于100 m 的情况下， 浅水非线性效应比较明显， 利用非线性浅方程运算可以得到较佳的计算精度。 而天文台在香港境内制作的最内两层网格亦刚好在水深少于100 m 的范围内而分辨率也是最高。 但这两层在运算的过程中所需要的时间会相对较长。 在水深大于100 m 的海洋上， 利用线性方程可节省运算时间。

以一个南海马尼拉海沟发生的地震所引发的海啸为例， 使用全部五层嵌套网格作运算的时间为45 min。 但如果减少最内一层而只使用四层嵌套网格， 运算时间则可缩短至7 min。

2 模拟海啸及验证

COMCOT 在一台内有两颗八核中央处理器的服务器运行， 运算速度为每秒 294 兆浮点（GFLOPS）。 模拟海啸之前， 需要把地震数据处理及分析系统所计算的各项参数和模拟时间的长度、嵌套网格等数据输入COMCOT 内作好准备。

图2 断层模型Fig 2 Fault model

初始条件方面， 模型假设水面初始垂直位移幅度与海底因断层而造成的变形幅度相等， 并采用了Manshinha and Smylie 断层模型（1971）[5]及配合日本气象厅提供的经验方程计算断层在海底的变形幅度， 例如长度（L）、 宽度（W）及滑动距离（D）（图2）。 而震中位置、 矩震级（Mw）、 震源深度、 走向角、 倾角及滑动角可透过香港天文台的地震数据处理及分析系统在地震发生后短时间内自动计算出来。

我们利用日本311 地震在日本及香港录得的海啸资料和其它香港曾经录得的海啸个案以验证COMCOT 模拟海啸的能力。

2.1 2011年3月11日在日本东北发生九级地震所引发的海啸

根据美国地质勘察局（USGS）的地震参数（表1）及利用日本气象厅提供的经验方程求得的断层模型变形幅度参数（表2）， COMCOT 模拟海啸影响日本近岸的情况， 模拟时间长度为20 h（图3）。

表1 2011年3月11日在日本东北发生9 级地震的地震参数Table 1 Earthquake parameters for M9 Tohoku earthquake on 11 Mar, 2011

表2 2011年3月11日在日本东北发生九级地震的断层模型变形幅度参数Table 2 Fault model parameters for M9 Tohoku earthquake on 11 Mar，2011

模拟结果显示影响相马及石卷的第一个海啸波高度分别为9.5 m 和7.0 m（图4）， 与当地两个验潮站分别录得的9.3 m 和7.6 m 接近（图5）。 在模拟最大海啸波高度图显示（图6）， 其他日本东北沿岸的海啸波高度跟验潮站录得的海啸波高度亦大致相约。

图3 利用COMCOT 模拟2011年3月11日在日本东北发生的海啸Fig.3 Simulation of Tohoku tsunami on 11 March，2011 by using COMCOT

图4 COMCOT 模拟日本沿岸相马及石卷的海啸波的高度时间序列Fig.4 Time series of tsunami wave height of Sōma， Fukushima and Ishinomaki， Miyagi of the coastal area of Japan simulated by COMCOT

海啸波在地震发生后约8 h 到达香港（图7）。香港三个验潮站（图8）所录得的海啸到达时间及第一个海啸波的高度和COMCOT 模拟结果的比较列于表3。

结果显示各潮站录得的第一个海啸波的高度与COMCOT 的模拟颇为接近， 相差不多于0.1 m，而录得的海啸到达时间与模拟所得也大致相约，差距少于15 min。 此外， 在验潮站的时间序列图中显示（图9、 10 及11)， 海啸到达后5 h 内的海啸波幅变化与模拟的波幅变化接近同步， 显示模式模拟海啸拍岸时的变化表现相当理想。

2.2 香港录得的其它海啸记录

除日本311 地震所引发的越洋海啸， 香港自1950年开始运作自动验潮站以来另外曾录得6 次海啸， 全部属于轻微水位变化（少于0.5 m）。 我们利用USGS 的地震参数模拟这6 次海啸事件， 模拟结果与香港鲗鱼涌验潮站的观测数据作比较（表4）。 结果显示其中5 次的模拟结果与实测的数值颇为接近， 只有1988年6月24日吕宋海峡所发生的5.7 级地震， COMCOT 没有显示任何海啸到达香港， 似乎模式的海啸波在抵达香港前已完全减弱。

图5日本沿岸验潮站在2011年3月11日录得的最大海啸波高度图（数据来源：日本气象厅）Fig.5 Maximum tsunami wave height observed by tide gauging system along the coast of Japan on 11 March，2011

图6 COMCOT 模拟日本沿岸在2011年3月11日的最大海啸波高度图Fig.6 Maximum tsunami wave height simulated by COMCOT along the coast of Japan on 11 March，2011

图7 COMCOT 模拟日本311 海啸到达香港的情况(时间序列：左至右，上至下)Fig.7 Arrival of tsunami at Hong Kong for 311 Japan tsunami simulated by COMCOT（Time sequence from left to right,up to down）

图8 香港其中三个验潮站横澜岛、鲗鱼涌和石壁的位置Fig.8 Positions of three tide gauge stations of Hong Kong (Wagland Island，Quarry Bay, Shek Pik)

表3 香港验潮站在日本311 地震录得的海啸到达时间及海啸高度和COMCOT 模拟结果的比较Table 3 Comparison of arrival time and wave height recorded by Hong Kong tide gauge stations and simulated by COMCOT for 311 Japan earthquake

图9 石壁验潮站的实测及模拟海啸波的时间序列Fig.9 Time series of measured and simulated tsunami wave at Shek Pik tide gauge station

图11 横澜岛验潮站的实测及模拟海啸波的时间序列Fig.11 Time series of measured and simulated tsunami wave at Waglan Island tide gauge station

表4 鲗鱼涌验潮站曾录得的六次海啸与COMCOT 的模拟比较Table 4 Comparison of the six recorded tsunamis at Quarry Bay and COMCOT simulation

3 COMCOT 的应用

应用平行运算技术， COMCOT 在一台内有两颗八核中央处理器， 运算速度为每秒294 兆浮点（GFLOPS）的服务器运行， 模拟日本311 地震后从太平洋到香港长约八小时的海啸过程， 运算时间只需约75 min。 而一个在南海马尼拉海沟发生的地震所引发的海啸， 到达香港的时间约需3 h，COMCOT 的运算时间也只不过是45 min 左右， 可见利用COMCOT 配合香港天文台的实时地震数据处理及分析系统， 足以支持香港天文台海啸预警系统的业务应用。

香港天文台的地震数据处理及分析系统一般可以在6 级或以上的地震发生后大约6～8 min 计算出地震发生的时间、 位置和震级。 而地震的震源机制亦可在发震后大约20 min 透过香港天文台的地震数据处理及分析系统计算出来， 或经由USGS 提供。 如果马尼拉海沟发生地震， COMCOT 在地震发生后大约一个小时便可以提供数值模式的计算结果， 为海啸预警的应急工作带来莫大裨益。

4 总结

本文利用日本311 海啸个案及其他六次香港曾录得的海啸记录验证了COMCOT 模拟海啸的能力。 在311 海啸个案的模拟结果显示， 影响日本东北沿岸的海啸波高度跟验潮站录得的海啸波高度大致相约， 而香港各潮站录得的第一个海啸波的高度与模拟颇为接近， 相差不多于0.1 m， 而录得的海啸到达时间与模拟所得也大致相约， 差距少于15 min。 此外， 由于COMCOT 应用平行运算技术， 缩短模拟所需的运算时间， 使其可以有效地应用于海啸应急的业务运作， 大大提升香港天文台海啸预警方面的能力。 再者， 一个在南海马尼拉海沟发生的地震所引发的海啸， 由于到达香港的时间约需3 h， 如果只使用四层嵌套网格模拟， 运算的时间则大大缩短为只需7 min， 对初步评估海啸的影响起了很大的参考价值。

[1] Xiaoming Wang. Technical document for COMCOT v.1.7[M]. New Zealand： Institute of Geological & Nuclear Science， 2009.

[2] Philip L.F. Liu. Seung -Buhm Woo， Yong -Sik Cho.Computer Programs for Tsunami Propagation and Inundation [M]. Ithaca， NY： School of Civil and Environmental Engineering Cornell University， 1998.

[3] Amante, C. B. W. Eakins, ETOPO1 1 Arc-Minute Global Relief Model： Procedures, Data Sources and Analysis [C].NOAA Technical Memorandum NESDIS NGDC-24, 19 pp,March 2009.

[4] Farr， T.G.， M. Kobrick. Shuttle Radar Topography Mission produces a wealth of data [J]. Amer. Geophys. Union Eos， 2000， 81： 583-585.

[5] Manshinha， L.， and D. Smylie. The displacement field of inclined faults [J]. Bulletin of the Seismological Society of America， 1971， 61（5）： 1433-1440.