海口湾海洋水动力自然灾害评估分析❋

高 峰， 朱川林， 李 昊， 唐友刚， 黄杰斯， 张慈珩， 刘 针， 姜云鹏

(1. 交通运输部天津水运工程科学研究所，港口水工建筑技术国家工程实验室，工程泥沙交通行业重点实验室，天津 300456；2. 天津大学建筑工程学院，天津 300072； 3. 海南金海湾投资开发有限公司，海南 海口 570000)

中国是世界上海洋灾害最严重、最频发的少数国家之一，而四面环海的海南省又是最易受到海洋灾害影响的地区，对于地处海南省经济发展前沿的最大城市——海口市而言，其濒临海口湾水域所面临的海洋灾害形势十分严峻，超强台风、风暴潮、极端大浪、海啸、海平面持续上升、岸滩侵蚀加剧以及海岸生态环境变化加速等，均对沿海经济造成威胁相当严重，已成为沿海开发规划实施和经济建设的制约因素之一。项目用海灾害风险分析主要是对项目灾害性风险的种类、项目用海对海洋功能的风险、项目用海对周边利益相关者的风险进行论证和分析[1]，进而提出避免灾害性风险发生的方案，或是灾害性风险发生后的应急措施。因此，开展沿海工程海洋灾害风险评价分析，是有效提高灾害防御能力、控制风险和减少灾害损失的关键技术方法之一，对该地区近海涉水工程的建设具有重要意义。

目前，一般风险分析的程序是：风险识别、风险估计、风险评价、风险处理和风险决策[2]。来自海洋水动力的风险主要包括两个方面：(1)一方面是用海项目自身引起的突发或缓发事件导致对海域动力环境造成的危害，发生于施工期居多，如波浪增加、流速加快、悬浮物扩散、船舶溢油等；(2)另一方面是由于自然灾害对海域使用项目造成的危害，发生于营运期居多，如遭受热带气旋、风暴潮、海啸等海洋灾害对工程造成的损坏。海南省四面环海、濒临南海，地理位置特殊，海口湾位于海南岛北部的海岸线上，分布于琼北中部岬间海湾内，隔琼州海峡与大陆的雷州半岛相望，背依海口市及澄迈县的马村(老城)工业开发区。其所在海域的海洋水动力灾害主要包括：风暴潮、海浪、海啸等。按照一般风险分析的程序，就近岸海洋动力风险分析研究来说，重点在于上述灾害的风险识别与估计，本文即通过实测资料分析结合数值模拟计算，确定上述各主要风险因子发生的强度大小和破坏情况。

1 评价内容与结果分析

1.1 风暴潮增水

登陆我国的台风和强温带天气过程往往造成风暴潮灾害，其成灾频率高、致灾强度大，造成的人员和经济损失惨重。风暴潮灾害主要是由异常的风暴增水使潮位大幅升高而形成，致灾因子不仅包括风暴潮，还包括天文大潮、近岸浪及其三者之间的耦合作用。影响海南岛北部的热带气旋较多[3]，是受风暴潮影响较为严重的区域之一，海水堆积引发的潮灾是海口湾水域主要海洋灾害之一。根据相关资料统计，海口年均发生风暴潮约3次，而严重或特大风暴潮约每2年1次[4]。海口湾每年6—12月均有可能发生风暴潮，且集中发生在7—10月，尤以10月份最多。在对历史风暴潮资料进行搜集整理的基础上，本文研究基于现有的实测台风资料通过将Monte Carlo法[5](又称为随机统计模型)与水动力模型相结合的方法实现对可能出现的台风风暴潮进行数值模拟研究，从而推求出重现期为50、100以及500 a的极端高水位预测值。

1.1.1 计算方法 采用假定水体上下层均匀一致的垂向平均二维风暴潮数值模型作为近海台风风暴潮的水动力模型，且二维模型突出了风暴潮预报与计算的最重要变量—风暴潮位，这也是预测台风风暴潮引起的极端高水位的核心变量。模型采用大小模型嵌套方式，其中大模型的北侧和西侧以大陆岸线为陆地边界，南面以北纬15.0°为外海边界，东面以东经117.9°为外海边界，模型范围覆盖了北部湾、琼州海峡、琼州海峡以东及南海部分海域，这种规模的大尺度计算网格可以很好地对一次台风影响过程进行数值模拟，其网格尺寸约为1 km，共有420×808个网格点，339 360个单元。同时，小模型南北两侧以大陆岸线为陆地边界，南北向跨度最大处约180.0 km，东西向从108.7°E～110.9°E，跨度约为229.2 km。计算范围及网格划分如图1(a)～(b)所示。

((a)大模型; (b)小模型; (c)1409号台风模拟; (d)验证结果。(a)Large scale model; (b) Local nested model; (c)Typhoon 1409 simulation; (d)Verification results. )

图1 模型网格嵌套分布与1409号台风模拟验证情况
Fig.1 Distribution of model nested grid and verification result of No.1409 typhoon’s simulation

气压分布模型采用Holland模型，环流风速模型采用梯度风速公式，移行风速模型采用宫崎正卫模型。由气压分布模型和风速模型组合成台风风场的模拟方案，简称Holl模型[6]。典型的实测台风依据美国FEMA使用的台风历史记录数据处理方法—模拟圆法，提取在一定时期内途径以要模拟的城市为模拟点，以设定的距离为半径的圆的所有台风，以这些台风的记录作为研究对象[7-8]。在1949—2015年期间曾有151个台风经过以海口湾为中心，200 km为半径的圆形研究区域，筛选登陆点中心最低气压在0～980 hPa范围内共计60个。在所有的60个实测台风中，2014年7月18日经过的“威马逊”台风强度最大，故利用铁山港海区英罗湾附近气象观测站实测资料对该场台风进行了风场验证，“威马逊”登陆后整个研究区域于2014年7月18日18:00的气压场与风场模拟情况可分别见下图1(c)所示，模拟与实测的最大风速拟合良好图1(d)。进而采用Monte Carlo模型生成100年内可能出现的模拟台风数据。由于台风历史资料采用的是1949—2015年间的60个典型台风数据，则500年内应当产生的台风数为448个。对这些个可能出现台风的参数进行模拟，可得所需的模拟台风数据。采用随机模型(Monte Carlo)随机模型和水动力模型(Delft3D)相结合的方法得到水位过程线，并得出极值水位。首先采用Monte Carlo随机模型，将筛选出的67年里60条台风的气压和路径分别作为独立的元素随机组合得到448个重组台风。将重组台风要素带入已验证的风暴潮水动力模型进行计算，其中台风过程与天文潮位随机组合，最后得到448个水位过程线，统计每条水位过程线中的最大水位值，用于后面的极端水位分析。

1.1.2 计算结果 通过以上过程，统计每条水位过程线中的最大水位值，部分较大水位过程线如图2(a)所示。将统计出的448个最大水位值进行降序排列，计算得秀英站处最大水位序列的最高值达3.75 m。极端水位的回归周期TR可按水位的排序参照式：TR=500/M估计，其中，M=8是448个最高水位的从高到低的排序，TR为重现期。例如，当M=1(对应的是448个极端水位的最大值)时，这个水位的最大值对的重现期是TR1=500/1=500 a。相应的极端水位所对应的重现期可点绘成图，且可由该图确定出任一重现期所对应的极端水位值。第五高水位3.02 m所对应的是100年一遇的极端高水位值，重现期为500、100和50年所对应的秀英站极端高水位预测值分别为3.75、3.02和2.68 m。

另外，通过采用重现期100年一遇气压P0=883 hPa，最大风速半径R为33 km，台风移速为35 km/h，采用假想路径，对可能最大风暴潮增水(PMSS-ProbableMaximum Storm Surge)[9]进行了计算，计算结果表明最大增水达到4.42 m，其对应增水过程见图2(b)。

((a)水位过程线; (b)最大可能增水过程。(a)Stage hydrograph; (b)Maximum possible water increase process.)

图2 各重现期秀英站台风暴潮水位过程线和最大可能风暴增水过程图
Fig.2 Storm surge water level process of different return period for Xiuying and the greatest possible storm surge elevation

1.2 极端大浪影响

海口湾属台风作用活跃区域，受台风和大陆季风作用明显，但该水域及周边缺乏长期观测站，现仅有玉苞角、白沙门两个临时站。影响该海域的台风主要来自南海自东或东南经广东、海南登陆或从海南岛南部海域向偏西方向移动的台风。海口湾水域全年涌浪较少、以风浪为主，风浪频率为76%～85%，涌浪频率占14%～23%，常浪向为ENE，频率为30.1%，次常浪向为NE，频率为22.9%。波浪出现最少的方向为S～WSW。受季风风向影响，工程海域常浪向随季节风而变化，东北季风期(11 月—次年3 月)，常浪向为ENE，次常浪向为NE，季风转换期(4—5 月，9—10 月)常浪向、次常浪向与东北季风期一样，但常浪向出现频率略低于冬季，次常浪向出现频率略高于冬季；西南季风期(6—8 月)常浪向为NE。目前已知最大波浪由1409 号威马逊台风超强台风产生，该超强台风在海口湾的东北侧通过，极有利于西北向波浪的成长和传入海口湾，在海口湾出现最大波浪，形成有效波高6.8 m，周期约7.8 s，该波高为偏西波向的异常特大值，在N-E-S方位的最大有效波高3.9 m，周期约6.3 s。

1.2.1 计算方法 该海域主要成风原因有两种，一为季风、二为气旋。因此，研究主要通过风场模拟计算波浪场，风场资料需综合考虑季风和台风两种影响。具体研究思路如下：

(1)采用中国气象数据网(www.data.cma.cn)公布的台风路径，结合欧洲中长期在南海海域的后报风场资料，选取每年影响工程海区的台风及季风过程，对风场进行分析。

(2)根据多年风场数据，采用SWAN模式建立大区域波浪计算模型，得到每年工程区附近外海深水区的年波浪极值，进而采用Pearson-Ⅲ型曲线适线法推算不同重现期波要素。

(3)根据外海不同重现期波浪条件，在自然地形条件下即工程建设前，利用SWAN模式建立小区域波浪计算模型，推算工程区设计波浪条件。

1.2.2 主要结果与分析 根据工程海区最近30年(1986—2015年)台风和季风资料，通过台风场模型和季风场数据分析计算工程海区大风过程，进而采用大范围风浪模型计算工程海区风浪场，从中选取主要计算浪向ENE、NE、NNE、N、NNW与NW向的波浪年极值，分析工程外海深水区不同重现期波浪。提取海口湾北侧外海计算点(水深约-70 m)处的主要计算方向波浪重现期结果见下表1，其中前文所提到的60个典型台风中，分别选取1522号“启德”台风与1117号“飞燕”台风过境琼州海峡时波高分布图(见图3)，其中重现期100年有效波高为7.14 m，重现期50年有效波高为6.80 m。

表1 海口湾不同重现期有效波高结果Table 1 Results of different return period of significant wave height in Haikou-Bay /m

另外，依据上述波高结果，通过HUDSON公式对海口湾近岸工程可能采用斜坡式护岸或者防波堤的护面块体重量进行了核算，如目前海口湾内各人工岛如普遍采用30 t护面块体时，其对应破坏波高可达8.78 m，当于重现期为1 400年波浪(年概率为0.071%)，若以设计年限100年考虑，这一等级的块体稳定性保障率可达93%。同时，根据爬高和越浪的计算分析，极端高水位条件下有效波高超过7.9 m就会出现危害堤顶铺面结构安全的越浪量，该波浪条件对应的重现期约为400年，既年出现频率为0.25%。

图3 典型台风工程海域作用过程Fig.3 Functional process of typical typhoon in engineering sea area

随着波浪由外海向近岸传播，充分考虑了湾内自然水深地形条件下波浪分布。采用SWAN模型对工程区的波浪分布情况进行计算分析。计算方向为ENE、NE、NNE、N、NNW与NW向。包含100年一遇高水位、设计高水位与设计低水位。其中，ENE最大有效波高为6.7 m(见图4(a))，NE向最大有效波高5.8 m(见图4(b))，NNE向为5.5 m，N向为3.8 m，NNW向为4.4 m，NW向为4.9 m。对海口湾近岸影响较大的2个方向为ENE与NE向，重现期100年一遇高水位时H13%波高分布见图4所示。

图4 100年一遇高水位重现期100年有效波高分布结果Fig.4 Distribution of significant wave height for 100 years return period by a high water of 100 years return period

2.3 海啸波

依据美国地质调查局(USGS)提供的一份潜在海啸地震断裂带风险图[10]，我国南海海域附近有三大俯冲带被认为具有地震激发海啸的潜在危险，分别为马尼拉俯冲带、琉球群岛俯冲带和北苏拉威西岛俯冲带。其中，马尼拉俯冲带(又称为马尼拉海沟)对南海周围的国家和地区会带来很大的海啸威胁，其位于亚欧大陆板块和菲律宾海板块交界处，南起菲律宾巴拉望岛北端，沿菲律宾吕宋西部边缘向北发展，北至台湾岛，总长度约1 000 km[11]。USGS根据马尼拉海沟的方位特征和断层几何学将其划分为6个断裂带[12]，沿马尼拉海沟分布如图5(b)所示。

上述区域地震激发的海啸波会使海南等我国南部海域暴露在海啸波威胁之下，虽然历史统计资料显示几次大的跨洋海啸对我国东南沿海地区影响有限，但南海马尼拉海沟地震带形成的海啸仍会对台湾、海南、广东和福建等沿海地区曾造成不同程度的影响，南海曾发生过约15次局地海啸，造成了一定影响[13]。由于近几年全球地质运动较为活跃，特大地震海啸发生频率较高，马尼拉海沟被认为是最有可能引发海啸的潜在地震源。基于上述分析，仍不可忽视海啸波的潜在威胁。

2.3.1 计算方法 本文采用浅水波数值模型GeoClaw，基于有限体积法对非线性浅水方程进行数值离散，考虑了海啸波在近岸传播的非线性作用，同时嵌入了自适应网格技术可追踪海啸波高判断网格加密与否。通过对马尼拉海沟断层特征的分析，Lau等[14]根据历史上曾经在附近区域发生过的地震和经验关系公式，确定了7.0、8.0级地震时E1～E3板块的数值模拟参数(见表2)，并基于Nguyen等[15]在USGS震源参数基础上给出的震级为9.3级的震源参数，此时对应各板块震源参数见表3，从而依此计算相应的传播对海南岛周边海域的影响。

( (a)南海附近俯冲带;(b)马尼拉海沟不同板块分布。(a)Subduction zone near the south China sea。(b)Distribution of different plates in the Manila trench.)

图5 中国南海附近的俯冲带与马尼拉海沟不同板块产生的海底变形及位置
Fig.5 Sea bottom deformation and location of subduction zone and different tectonic plates from the Manila trench nearby China South Sea

利用Okada弹性半空间错位理论模型，根据前文所述六块断裂带参数生成相应的海底变形可分析出海南岛受E1、E2和E3板块的影响较大。据此将对上述各板块(震源深度15 km)分别发生7.0级和8.0级地震形成的海啸进行数值模拟，进而估算对海口湾水域的影响情况，还考虑了整个断裂带(E1～E6)发生破裂形成的特大地震海啸时的极端情况。计算范围从99°E～130°E，1°S～33°N，包含了我国东南沿海至南海海域以及菲律宾、越南等国家。

表2 震级分别在7.8和8.0时马尼拉海沟E1～E3断裂带震源模拟参数Table 2 Source modeling parameters of E1 ～ E3 fault zone in Mw=7.8 and 8.0 respectively in the Manila trench

表3 Mw=9.3时，马尼拉海沟E1～E6断裂带基本参数Table 3 Source modeling parameters of E1 ～ E6 fault zone in Mw=9.3 in the Manila trench

表4 各计算工况下工程区海啸波特征值 Table 4 Characteristic value of tsunami under the different calculation conditions

2.3.2 结果与分析 数值模拟结果如图6所示，分别给出了马尼拉海沟北部3个板块(E1、E2、E3)分别发生里氏7.0、8.0级海底地震形成的海啸波。各工况下工程区的海啸波特征值如表4所示。从走时图中可以看到，海啸波大概经历4～5 h即可到达海口湾附近海域。对于单一板块的一般震级地震，海啸波到达海口附近波幅都小于4 cm，波长在50 km左右，周期1.6～2 h，海啸持续时间(第一个大波的作用时间)约为0.7～1.1 h，海啸波行进速度为8.1 m/s。对于整个断裂带发生9.3级大地震引发的极端海啸，海啸波到达工程区附近海域时最大波幅为2.25 m，波长115 km，周期3.4 h，海啸持续时间1.3 h，海啸波行进速度为9.4 m/s，该情况将对海口湾会产生影响，特别是岬湾形态的海口湾近岸，更易形成能量集中，应予以关注。

图6 各板块分别在7.0和8.0级地震生成的海啸波走势图Fig.6 The tsunami wave distribution of each plate in Mw=7.0 and 8.0 generated by earth quake

3 结论与分析

海洋水动力风险是用海风险的重要组成部分，关键在于风险识别和风险估测。本文根据对海口湾水域可能出现的风暴潮、极端波浪、海啸等问题的研究，对该海域水动力风险因素进行了研究分析，得到风暴潮不同重现期增水水位，并在确定外海重现期100 a有效波高为7.14 m，重现期50 a有效波高为6.8 m的基础上，计算了海口湾近岸不同波向的最大波高，对湾内近岸斜坡式防护工程采用护面块体的重量校核。针对潜在海啸威胁最大的马尼拉海沟，进行了多种工况数值模拟，除9.3级大地震引发极端海啸导致海啸波达2.25 m之外，7、8级海啸对于海口湾所受影响较小。这些计算结果反映海口湾沿岸在不同水位下的漫滩情况，对海口湾沿岸重要的工程项目以及防灾规划的制定，具有参考意义。