Orginal Article

Numerical Simulation and Analysis of Tsunami Impacts on Hainan Island and Beibu Gulf Zone

  • LI Zhiguang , 1, 2, 3 ,
  • XIE Shunping , 1, 2, 3, * ,
  • DU Jinkang 1, 2, 3 ,
  • HUANG Yang 4 ,
  • ZUO Tianhui 4 ,
  • ZHENG Wenlong 1, 2, 3
Expand
  • 1. Department of Geographic and Oceanographic Sciences, Nanjing University, Nanjing 210046, China
  • 2. Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, Nanjing University, Nanjing 210046, China
  • 3. Collaborative Innovation Center of South China Sea Studies, Nanjing 210093, China
  • 4. Guangxi Earthquake Administration, Nanning 530022, China
*Corresponding author: XIE Shunping, E-mail:

Received date: 2015-01-04

  Request revised date: 2015-03-02

  Online published: 2015-08-05

Copyright

《地球信息科学学报》编辑部 所有

Abstract

The impact of tsunami hazard to the coastal countries and regions is a hot topic within the field of tsunami research. One of the main reasons that cause a tsunami is submarine earthquake. Since the Manila Trench and a part of China´s territory cover across the South China Sea, and the trench has frequent seismic activity, it becomes a potential tsunami source which may influence the south China. A hypothetical tsunami is simulated by adopting the COMCOT numerical tsunami model, to research its effects on the area of Hainan Island and Beibu Gulf. The tsunami is triggered by an earthquake in the Manila Trench, and the entire source is discretized into 33 rectangular elements, which would produce an earthquake with a magnitude of Mw=9.0. The tsunami wave characteristics and its impacts on this region are analyzed with the help of records from the virtual tidal stations installed in the study area. The majority of tsunami energy is directed northward to the coast of China. The results indicate that the tsunami waves propagating to the Beibu Gulf Zone have distinct nonlinear characteristics, and spread more slowly than in the deep ocean. After 2 hours and 10 minutes when the earthquake occurs, the head wave propagates to the east coast of Hainan Island with the amplitude of 2.6 meters, which could bring great impacts or hazards to this area. However, because of sea bottom friction and the blocking of Hainan Island, it takes another six hours for the head wave to reach Beihai offshore, and the wave amplitude then is reduced to be less than 50 centimeters. The energy of tsunami undergoes a great loss, and it may have little impacts on the southern coast of Guangxi Province and the Beibu Gulf area. And yet, the additive effects of the tsunami and the tidal currents in the gulf area require further research.

Cite this article

LI Zhiguang , XIE Shunping , DU Jinkang , HUANG Yang , ZUO Tianhui , ZHENG Wenlong . Numerical Simulation and Analysis of Tsunami Impacts on Hainan Island and Beibu Gulf Zone[J]. Journal of Geo-information Science, 2015 , 17(8) : 937 -944 . DOI: 10.3724/SP.J.1047.2015.00937

1 引言

海啸是由于突然的海底形变或水体扰动所产生的一列长波和周期极长的海洋重力波。进入21世纪以来,全球已发生了多次灾害性海啸,如2004年苏门答腊海啸、2011年日本海啸,造成了巨大的人员伤亡和财产损失,这给位于板块交界边缘并拥有漫长海岸线的中国敲响了警钟[1-3]。为此,许多国家的科学家和组织,对地震海啸的形成机理及其影响作了深入研究,构建了国际海啸预警系统,共同预防此类自然灾害。
中国海岸线外的大陆架平缓宽阔,外围又有包括印尼诸岛、菲律宾群岛、台湾岛、琉球群岛、日本群岛、千岛群岛等一系列的岛屿或半岛,形成一道天然屏障[4],阻碍了越洋地震海啸波的传播。海底地震是造成海啸的主要原因,有研究统计全球86%的海啸由海底地震引发[5]。地处中国南海中央海盆东侧的马尼拉海沟,北与台湾碰撞带相连,东边邻接吕宋岛弧,是菲律宾海板块和东亚大陆板块的分界线[6]。观测记录显示1970年1月1日到2014年11月23日间,马尼拉海沟发生Mw≥5.0震级的地震达1519次,表明该地区板块运动活跃,是值得关注的一个海啸源。一些学者通过数值模拟方法研究了该区域海底地震引发的潜在海啸,并做了影响分析。Ruangrassame等、Wu等分别分析了马尼拉海沟引发海啸对泰国湾和台湾的影响[7-8],Ha等对南中国海的海啸传播场景做了模拟和分析[9]
数值模拟是研究地震及波浪传播等自然过程和现象的科学方法[10-11],模拟地震海啸过程需综合海底地震和波浪传播2个过程。COMCOT是美国康奈尔大学在前人工作基础上研发的海啸数值模拟模式[12],后经相关学者多次改进趋于成熟,该模式已成功用于1960年智利海啸、2004年苏门答腊海啸、2006年台湾南部海啸、2010年智利海啸、2011年日本海啸等的数值模拟,为该领域较为经典和实用的数值模拟模式之一。本文采用COMCOT模式,对可能发生在马尼拉海沟一场震级为Mw=9.0级的地震活动所引起的海啸过程进行数值模拟,分析海啸传播到海南岛及北部湾海域的过程和特点,评估海啸灾害对该地区可能造成的影响,并进行相关的机理分析,研究结果可为该地区的海啸预警与防灾提供一定的参考。

2 COMCOT模式理论

海啸过程可抽象为激发、传播及上岸3个阶段[13],激发海啸的原因有很多,如海底地震、火山喷发、海底滑坡、陨石撞击等。在COMCOT模式中可通过波浪生成器、滑坡模型和断层模型3种方式模拟激发海啸,COMCOT模式中采用目前国际上较为通用的Mansinha & Smylie模型[14]和Okada模型[15],这2套模型都是以弹性错移理论发展的震源断层模型。本文通过断层模型的参数设置,构造可能诱发海啸的马尼拉海沟海底地震。
根据海啸波在不同水域的传播特点,COMCOT模式分别针对大洋深水传播和近岸浅水传播提供了线性和非线性控制方程,非线性方程添加了近岸海底摩阻等非线性项;另外,大洋传播时需考虑纬度变化和科氏力影响,模式中有球坐标系下的浅水波方程,也有适用于小范围内的直角坐标系方程。式(1)、(2)分别是球坐标系下的线性和非线性方程组。
η t + 1 Rcosφ P φ + φ cosφQ = - h t P t + gh Rcosφ η ψ - fQ = 0 Q t + gh R η φ + fP = 0 (1)
η t + 1 Rcosφ P φ + φ cosϕQ = - h t P t + 1 Rcosφ ψ P 2 H + 1 R φ PQ H + gH Rcosφ η ψ - fQ + F x = 0 Q t + 1 Rcosφ ψ PQ H + 1 R φ P 2 H + gh R η φ + fP + F y = 0 (2)
式中,η为相对于平均海平面的自由表面位移(m);φ为纬度(°),Ψ为经度(°);R为地球半径(km);h为静水深(m); H = η + h 为总水深(m);f为科氏力系数;g为重力加速度(m/s2);(P,Q)为沿纬度和经度单位宽度的通量(m3/s)。FxFy分别为xy方向的底部摩阻,其计算如式(3)所示,其中n为曼宁系数。
F x = g n 2 H 7 3 P P 2 + Q 2 F y = g n 2 H 7 3 Q P 2 + Q 2 (3)
COMCOT模式采用的水边界是开放边界,线性浅水波方程的海陆边界为垂直反射边界,非线性浅水波方程的则为移动边界。该模式用交错式显性蛙跳有限差分法计算浅水波方程,由此可借助差分方程的数值频散来近似代替波浪在浅水中传播所带来的物理频散。针对波浪传播在大洋海域和近岸区域的复杂性差异,并顾及对计算效率与模拟精度的不同需求,可采用不同分辨率的海底地形资料。在COMCOT模式中采用了多层格网嵌套方式,外层使用大格网,覆盖震源到研究区的广大区域;内层嵌套更高分辨率的子格网,聚焦研究区的不同子区域。随着计算区域的水深变化,不同的格网解析度和模拟时间步长可确保数值模拟能适当地表示海啸物理过程中的频散效应[16]

3 模拟数据和方法

本文以海南岛和北部湾为研究区,引发海啸的震源位于马尼拉海沟,将研究区设置成2个级别的格网,其中,一级格网覆盖从震源到海南岛及北部湾的整个研究区,一个二级格网覆盖海南岛及北部湾的大部分区域。研究使用美国地质调查局(USGS)地理数据库中分辨率为0.5′的USGSDEM地形数据,通过数据预处理和裁剪获得覆盖北部湾大部分区域精度为0.5′的地形测深数据,数据规模为600×540网格单元,用于二级格网运算;通过邻域重采样方法和裁剪处理获得覆盖马尼拉海沟震源到北部湾区域精度为2′的地形测深数据,数据规模为751×601网格单元,用于一级格网运算,如图1所示。表1列出了2层格网的基本参数及所选择的计算模式。
Tab. 1 Parameters configured for two-level grids in the simulated domain

表1 模拟区域格网设置参数

格网 分辨率 范围 行列数 控制方程 坐标系 曼宁系数
一级 2′ 5°~30°N, 105°~125°E 751×601 线性 球坐标系 缺省
二级 0.5′ 17°~22°N, 106.5°~111°E 600×540 非线性 球坐标系 0.013-0.052
Fig. 1 Configuration of two-level grids and virtual tidal stations in the simulated domain

图1 计算区域格网和虚拟验潮站设置

波高及周期是刻画海啸波最基本的2个特征量,尤其当海啸波趋向近岸时,波形发生变化。为了分析北部湾地区海啸波的波高、到达时间、速率等特征,从北部湾湾口向广西北海沿岸依次设置了10个虚拟潮位站(TS1-TS10),所处位置的水深逐渐变浅;海南岛周围设置了4个(TS11-TS14),如图1(b)所示。
目前,抽象马尼拉海沟断层构造的方法主要有2类:(1) 将整个海沟离散化为6个断层块,2006年USGS发布了该区域的6个断层块的构造参数,Wu等通过对6个断层块构造参数的调整,模拟和研究了海啸灾害对台湾的影响[8];(2) 断层块单元离散化的方法,Liu等将马尼拉海沟断层构造分解为39个长70 km、宽35 km的断层平面单元,用于南中国海海啸风险和预警系统研究[17]。Megawati等综合利用马尼拉海沟的多种数据,考虑潜在断裂区域的几何结构,对不同位置俯冲块的沉潜角度变化判断,对12.5°~23.5°N之间的10个地震横截面[18]进行插值,得到破裂模型[19]。在此基础上利用GPS数据,在最大相对运动区域设定40 m的滑移值,并假定滑移方向近似正交于海沟轴线方向,确定地震的滑移模型,如图2(a)所示。为确定垂直海底位移,根据断层带不同位置的构造特征将震源离散为33个断层块,所对应的震级是Mw=9.0,如图2(b)所示。表2给出了断层参数。由于该模型能较为全面和恰当地表达震源区域的构造情况,本文研究采用这种断层单元离散化模型。为模拟带状震源深度的空间变化,根据参数理论区间对各离散单元的震源深度设置如下:1、2、18~20、33为10 km;3、4、17、21、32为20 km;5~8、16,22~24、30、31为30 km;9、10、15、25、29为40 km;11~14、26~28为50 km。模型中的其他参数不变。
Tab. 2 The parameters of the 33 discretization fault units in Manila trench

表2 马尼拉海沟33个离散化断层单元参数

序号 经度(°E) 纬度(°N) 深度(km) 走向角(°) 倾角(°) 滑移角(°) 长度(km) 宽度(km) 相对位移(m)
1 120.5 12.6 10 324.46 21.78 90 40.41 19.30 5
2 120.1 13.2 10 325.20 11.26 90 54.19 38.68 5
3 119.8 13.5 20 318.66 6.55 90 54.09 66.93 12
4 119.5 14.0 20 332.40 5.79 90 54.00 75.68 12
5 119.3 14.6 30 0.26 6.47 90 53.89 67.52 25
6 119.2 15.1 30 7.39 11.50 90 80.59 37.51 28
7 119.3 15.8 30 5.85 10.01 90 53.56 43.13 28
8 119.4 16.3 30 355.99 8.46 90 53.43 51.06 28
9 119.3 16.8 40 358.34 7.18 90 53.28 60.09 30
10 119.3 17.4 40 2.50 6.16 90 53.14 69.93 30
11 119.3 17.9 50 16.26 6.52 90 52.99 65.93 35
12 119.7 18.5 50 40.34 5.93 90 105.57 72.15 40
13 120.2 19.0 50 35.93 5.36 90 52.60 79.67 40
14 120.5 19.6 50 21.46 5.70 90 78.62 74.61 35
15 120.9 20.3 40 352.23 3.28 90 78.24 129.36 30
16 120.3 21.0 30 332.43 6.25 90 103.81 67.36 25
17 119.8 21.8 20 339.52 7.62 90 51.68 54.86 12
18 119.7 22.2 10 341.26 9.89 90 46.82 41.97 5
19 120.6 13.0 10 326.63 28.87 90 37.50 37.27 5
20 120.5 13.4 10 351.15 25.45 90 54.14 43.11 5
21 120.0 14.1 20 333.5 30.96 90 134.88 34.22 12
22 119.8 15.0 30 357.94 20.98 90 53.76 53.11 28
23 119.8 16.0 30 11.30 24.22 90 133.84 45.09 28
24 119.8 16.8 30 9.30 18.22 90 53.29 61.36 28
25 119.8 17.3 40 10.90 15.85 90 53.15 70.96 30
26 120.2 17.8 50 47.78 12.76 90 105.90 88.63 40
27 121.0 18.8 50 30.91 14.73 90 131.68 75.92 40
28 121.3 19.7 50 37.06 16.31 90 52.41 67.88 35
29 121.7 20.3 40 24.75 23.79 90 104.39 44.85 30
30 121.4 20.8 30 302.02 31.89 90 51.99 31.66 25
31 120.9 21.2 30 340.64 19.63 90 103.67 55.05 25
32 120.2 21.9 20 320.54 16.07 90 77.33 67.83 12
33 120.0 22.4 10 326.04 24.96 90 21.26 41.86 5
Fig. 2 The slip model and the discretized model for computation of seafloor displacement[19]

图2 海底震源的滑移模型和计算海底位移的离散化模型[19]

4 海啸数值模拟与结果分析

4.1 海啸数值模拟

本研究采用COMCOT模式,以2个级别嵌套格网和震源断层模型,对近岸区域曼宁系数分别采用0.013、0.026和0.052,进行地震引发海啸模拟。海啸波的传播速度V与海水深度h的关系为 V = gh [13],即海啸波速仅与海水深度和重力加速度相关,不同的曼宁系数影响的是海啸推进波幅。图3为采用0.026曼宁系数模拟获得的海啸传播时间序列图,从图3可看出,地震发生时震源附近产生了波幅约20 m的海啸,地震发生1 h海啸前导波已向南中国海方向传播约700 km,地震发生2 h海啸前导波逼近海南岛东南沿岸,此时海啸传播了约1070 km,到达海南岛东岸的首轮波高约为3 m。地震发生4 h后海啸首波抵达北部湾中部,8 h后海啸首波到达广西西南沿岸,到达时的首轮波高不足1 m。海啸在外海传播的平均速度约为535 km/h,而在北部湾内的平均速度则降至100 km/h。
Fig. 3 Time series of tsunami propagation

图3 海啸传播时间序列图

4.2 海啸模拟结果分析

通过对模拟获得的波高场时间序列数据定位检索,获得各潮位站的头波波高和最大波高及相应的到达时间数据,图4为海啸到达14个验潮站的时间和首轮波高、最大波高的直方图,其中,TS1-TS7是从湾口向北海地区依次设定的潮位站,TS8、TS9、TS10位于广西近岸地区的3个潮位站,TS11、TS12、TS13、TS14是围绕着海南岛南部的4个站点。从图4(b)可看出,海啸波到达TS12潮位站的首波波高达到2.08 m,最大波为3.10 m,海啸波向海湾内传播过程中,头波波高逐渐减小,到达广西南部近岸3个站点的波高在0.27~0.36 m之间,最大波为0.50 m。设在深海区域或接近深海区域的潮位站的首波和最大波到达时间较为接近,如TS1、TS2、TS3和TS11、TS12、TS13,而相对远离深海区域潮位站的最大波到达时间较首波到达时间明显滞后,如TS4、TS5、TS7、TS8、TS9、TS10、TS14。海啸波在绕过海南岛传播向湾内时,海南岛近岸海啸首波波高也随传播距离逐渐衰减,说明海啸波在向近岸传播的过程中,海啸波呈现明显衰减的趋势。
Fig. 4 The arrival time and wave heights of head-wave and maximum-wave for all tidal stations

图4 各站点海啸首波和最大波到达时间和波高

根据国际通用的判定海啸发生级别的渡边伟夫海啸等级表(表3),本海啸场景下,北部湾湾口的4个验潮站(TS1、TS2、TS11、TS12)为1级,可能会对该区域造成房屋船只损失;其他潮位站为0级,可能造成轻微损失。
Tab. 3 Watanabe Hideo tsunami scales

表3 渡边伟夫海啸等级表

等级 海啸波高(m) 海啸能量(1010 J) 损失程度
-1 <0.5 0.06 能量损失
0 1 0.25 轻微损失
1 2 1 损失房屋船只
2 4~6 4 人员伤亡,房屋倒塌
3 10~20 16 ≤400 km岸段严重受损,人员伤亡大、房屋损毁严重
4 ≥30 64 ≥500 km岸段严重受损,人员伤亡巨大、建筑物尽毁
图5为北部湾从外向内5个验潮站的波高随时间变化曲线图。分析发现,TS1、TS11验潮站在海啸到达初期潮位起伏变化剧烈,随后变化逐渐平缓;位于广西北海沿岸的TS10验潮站的潮位在海啸到达后的1 h内上升至0.5 m,随后逐渐下降。这些验潮站波幅变化没有明显的周期性,但都呈现出逐渐衰减的趋势,表现出明显的非线性特征。另外,当首波到达潮位站时,无论其波高大小,都会在高于平均海水面的位置上维持较长一段时间才出现波谷。
Fig. 5 The time-change curves of wave heights in some virtual tidal stations

图5 不同虚拟站点波高时变特征曲线

从数值模型分析可知,海啸波的非线性特征主要通过非线性波浪控制方程中的水深和曼宁糙率系数表征,即在浅水区海底摩擦是海啸波幅衰减的主要因素之一。由于无法获得北部湾海底平均糙 率,相关研究建议的近岸区域曼宁系数取值参考范围0.025~0.030 m-1/3s[20],本文为非线性模型设置不同曼宁系数进行模拟,并参照线性模型的模拟结果,考察北部湾海域底部不同摩擦糙率对海啸传播的影响(图6)。相比线性模型模拟结果,除TS1外的各站点最大波高均随曼宁系数的增大而减小,TS1站点水深为1249 m,处于深海水域,此时非线性模型退化为线性模型,基本不受海底摩擦影响;其他站点所处水深均在100 m以内,适用于二级格网的非线性模型,最大波高受海底底部摩擦的影响明显。
Fig. 6 Change of maximum wave heights simulated with different Manning coefficient for all virtual tidal stations

图6 不同曼宁系数模拟的各虚拟验潮站最大波高变化

5 结论

本文采用COMCOT模式,以及双层嵌套格网结构,对可能发生在马尼拉海沟一场Mw=9.0级的地震活动所引起的海啸过程进行数值模拟,结论如下:
(1)COMCOT模式可有效模拟大洋海底地震及其引发的海啸传播过程,一级格网线性模型和嵌套其中的二级格网非线性模型的结合可分别以不同离散尺度模拟海啸波在深海区域和浅海区域的传播机理。
(2)地震发生2 h海啸前导波到达海南岛东部沿岸,首轮波高约为2.0 m,最大波高2.7 m;地震发生2.1 h海啸前导波到达海南岛东南角,首轮波高约为2.1 m,最大波高3.1 m,根据渡边伟夫海啸等级可划定为1级,可能造成一定的房屋和船只损失。而对海南岛西岸影响轻微,首轮波高减弱至不足1.0 m。
(3)在北部湾内,由于受海南岛的部分阻隔、水深变浅和海底摩擦的综合影响,海啸传播速度急剧衰减,地震发生8 h后海啸首波到达广西西南沿岸,到达时的首轮波高不足0.3 m,最大波高0.5 m,根据渡边伟夫海啸等级可划定为0级以下,基本不会构成灾害。数值模拟结果说明平均水深较浅的北部湾地区总体上受马尼拉海沟海底地震引发海啸的影响较小,但也需注意海啸波与环流、潮汐等的叠加效应造成的破坏。
本研究仅对海南岛及北部湾地震海啸波进行数值模拟和影响分析,未考虑其他并发因素的叠加效应,后续工作将着重研究海啸波与潮汐波的耦合作用及其影响。

The authors have declared that no competing interests exist.

[1]
柯长青. 印度洋地震海啸(2004-12-26)及其对中国的警示[J].中国地质灾害与防治学报,2006,17(4):91-96.

[2]
温瑞智,公茂盛,谢礼立.海啸预警系统及我国海啸减灾任务[J].自然灾害学报,2006,15(3):1-7.

[3]
于福江,原野,赵联大,等.2010年2月27日智利8.8级地震海啸对我国影响分析[J].科学通报,2011,56(3):239-246.

[4]
潘文亮,王盛安,蔡树群.南海潜在海啸灾害的模拟[J].热带海洋学报,2009,28(6):7-14.

[5]
王培涛,赵联大,侯京明,等.中国的海啸灾害危险性及海啸预警系统——解读2011年太平洋海啸演习[J].海洋学报,2012,29(5):9-16.

[6]
陈志豪,李家彪,吴自银,等.马尼拉海沟几何形态特征的构造演化意义[J].海洋地质与第四纪地质,2009,29(2):59-65.

[7]
Ruangrassamee A, Saelem N.Effect of Tsunamis generated in the Manila Trench on the Gulf of Thailand[J]. Journal of Asian Earth Sciences, 2009,36(1):56-66.

[8]
Wu T R, Huang H C.Modeling tsunami hazards from Manila Trench to Taiwan[J]. Journal of Asian Earth Sciences, 2009,36(1):21-28.

[9]
Dao M H, Tkalich P, Chan E S, et al.Tsunami propagation scenarios in the South China Sea[J]. Journal of Asian Earth Sciences, 2009,36(1):67-73.

[10]
黄金,高星,王伟.地震勘探全波形反演的应用与发展分析[J].地球信息科学学报,2014,16(3):396-401.

[11]
姚远,蔡树群,王盛安.海啸波数值模拟的研究现状[J].海洋科学进展,2007,25(4):487-494.

[12]
Liu P L F, Wang X M. COMCOT: A tsunami modeling package[DB/OL]. 2009.

[13]
陈运泰,杨智娴,许力生.海啸、地震海啸与海啸地震[J].物理,2005,34(12):864-872.

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

[15]
Okada M.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.

[16]
潘文亮,王盛安.COMCOT数值模式的介绍和应用[J].海洋预报,2009,26(3):45-52.

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

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

[19]
Megawati K,Shaw F, Sieh K, et al.Tsunami hazard from the subduction megathrust of the South China Sea.Part I.Source characterization and the resulting tsunami[J]. Journal of Asian Earth Sciences, 2009,36(1):13-20.

[20]
刘双庆. 海啸近场特征与海啸触发源之间的定性定量关系分析[D].兰州:中国地震局兰州地震所,2008.

Outlines

/