0 引言

台风风暴潮是威胁沿海低地城市安全的重要气象和海洋灾害^[1-3]。台风可诱发风灾、强降雨、洪水和泥石流等次生灾害，其中极端风暴潮引起的增水效应与风浪效应对沿海低地最具危险性。通常海岸防护工程设计、建设与安全评估等都需要考虑增水量和风浪等级这两个致灾因子^[4]。而且，有效的台风风暴潮增水和风浪模拟对沿岸沉积物输运和污染物传播十分重要^[5]。因此，建立高精度的数值模拟系统实现特定风暴情景下海洋增水和风浪模拟具有重要意义。

台风风暴潮模拟可追溯至19世纪30年代^[6]。Emanuel K.A.(1987)基于实测数据建立了台风风暴潮增水与台风等级之间定量统计关系，后续又有大量研究跟进和深化^[7-8]，但统计模型的最大问题是缺乏物理基础，无法实现区域动态和过程模拟。随着流体力学和计算数学的迅速发展，基于流体静压假设的一维、二维和三维数值模拟技术得到长足发展，例如20世纪80—90年代由美国国家气象中心(ANMB)研发的SLOSH模型、荷兰三角洲研究院(DELTAS)研发的DELFT系列模型、丹麦丹华水力公司(DHI)开发的MIKE模拟软件、荷兰代尔夫特理工大学开发的SWAN模型、法国Electricite de France (EDF)的Research and Development (R&D)实验室开发的TELEMAC-MASCARET数值模拟软件等，被广泛应用于水动力和风暴潮数值模拟。近期研究表明台风风暴潮有效模拟除需要考虑台风路径、中心气压、海表风场和海-气耦合外^[9-10]，还需考虑波-流耦合过程及天文潮、增水和波浪的多重耦合^{[4, 11-13]}。未来，因全球气候变暖导致的海平面上升和风暴增强等问题将使风暴潮灾害进一步恶化^[3]，这也是国际关注的热点问题，目前这些问题已在全球各大海域做过相关研究。Shen等^[14]采用数值模拟方法研究长江口及其附近海域海平面上升对风暴潮增水的影响。Karim和Mimura^[3]研究印度洋海域气候变化与海平面上升对风暴潮洪水的作用。Ferrarin等^[15]采用数值模拟方法研究地中海风暴潮增水对意大利沿岸的影响。Ebersole等^[16]基于ADCIRC模型研究墨西哥湾“卡特里娜”台风(2005年)期间风暴增水、风浪效应及漫堤淹没灾害。

北部湾位于中国南海北部海域，是由雷州半岛、中南半岛和广西海岸围成的半封闭湾区(图 1)。北部湾是台风风暴潮影响重灾区，平均每年受到2~3次台风袭击^[7]。近期，采用数值模拟方法研究北部湾台风风暴潮的文献逐渐增加^{[9, 17-18]}，但采用大气-波浪-潮流耦合方法研究风暴潮增水和风浪模拟的研究较少。本文基于TELEMAC和TOMAWAC数值模拟软件开发大气-波浪-潮流耦合风暴潮数值模拟系统，研究北部湾附近风暴增水和风浪效应。模型采用2012年台风“山神”和台风“启德”实测水位数据做率定和验证，最后基于模拟所得增水和波浪强度，探讨台风风暴潮对北部湾天文潮与风暴潮增水交互作用的非线性叠加效应。

1 材料与方法 1.1 模型耦合

台风风暴潮增水和风浪模拟采用大气-潮汐-波浪模型耦合而成：台风模型采用Holland模式^[19]，潮汐模型采用TELEMAC-3D模式^[20]，波浪模型采用TOMAWAC模式^[21]。模拟采用三级嵌套方式实现耦合。台风模型分辨率为12.5 km，覆盖南海和菲律宾部分地区(图 1)；潮汐模型采用可变分辨率三角网，覆盖南海大部分海域(图 2a)；增水模型采用TELEMAC-TOMAWAC耦合实现，覆盖整个北部湾(图 2c)。

风应力和气压场是台风风暴潮模拟的主要大气驱动力，采用Holland台风模型计算^[19]。Holland台风模型是假定双曲线辐射压力场模型，采用解析方式计算气压场和台风风场，具有控制参数少、精度高的特点。解算10 m高风场和气压场所需参数(包括台风路径、中心点坐标、中心气压值、最大风速和最大风速半径)可从中国台风网www.typhoon.org.cn获取。模拟区下垫面地形数据来源于美国地质调查局全球高程数据。大气模型模拟时长为一个月，前5 d设置为模型冷启动时间。

潮汐模拟采用TELEMAC-3D实现。TELEMAC-3D是由法国EDF的R&D实验室开发的非结构有限元三维数值模拟软件，目前由英国HR Wallingford和法国EDF等机构共同维护更新的开源软件系统。TELEMAC-3D水平向网格采用不规则三角网划分法，垂直向采用sigma分层法。可变分辨率三角网划分方案一方面利于贴合复杂多变的岸线，另一方面可对重点区域加密网格而对非重点区域采用粗分辨率表达。模型基于Boussinesq流体静压假设，解算三维不可压缩雷诺方程组。液体自由水面计算过程考虑垂向sigma分层传导效应。分别采用2次龙格库塔法和2次隐式梯形法解算水平传导项和垂直传导项。通过线性梯度重构技术实现2次空间模拟精度。其他关于TELEMAC-3D控制方程组及详细设置方案参考文献[20]。风场和气压场边界由Holland模型提供，波浪边界由TOMAWAC模型提供，潮汐边界主要考虑TPXO的8个分潮调和而成。

风浪模拟采用TOMAWAC实现。TOMAWA- C同样采用不规则三角网格，是EDF开发的第三代风浪模型。TOMAWAC模型基于波能守恒方程，内嵌波能耗散项，解算浅水方程组波浪传播过程，主要考虑风场边界和水位边界^{[11, 22]}。与SWAN模型类似，TOMAWAC模拟波浪传播和波浪破碎项都考虑波-流耦合过程。其中，方向频率增量因子设置为1.2，初始频率设置为0.051 2 Hz，频率段数设置为25，方向频段设置为36。详细控制方程及模型设置方案参考文献[21]。

1.2 模型设置

TELEMAC和TOMAWAC模型耦合提供波-流耦合计算。考虑风暴潮作用下海-气耦合需要考虑大气场与海洋场相互作用过程，本文建立了覆盖整个北部湾和南海部分区域的三角网格，实现大气-海洋-波浪耦合模拟。整个模拟区域覆盖东经105.6°—114.4°E和北纬15.6°—22.8°N海域。北部湾附近分辨率较高，达800 m，而南海区域分辨率较低，为5 000 m(图 2a)。水平向网格数目为124 000个、节点数为278 000个，垂直向分5层，总共620 000个网格。地形水深数据来源于GEBCO全球水深数据库和北部湾实测水深。底摩擦系数设置为空间均一分布值(Nikuradse=0.001 m)。

冷启动初始水位、流速和风场都设置为0。冷启动时段过后加入海-气耦合和波-流耦合实现风暴潮数值模拟。Holland台风场模拟时长为1个月，TELEMAC潮汐模拟时长为20 d，TELEMAC-TOMAWAC波-流耦合模拟时长为10 d。岸线边界设置为吸收潮汐和波浪的封闭边界，开放边界水位变化设置为

$\eta = \frac{{\mathit{\Delta }P}}{{\rho g}} + \sum {{f_i}} {H_\mathit{i}}\mathit{cos}\left( {2\pi \frac{t}{{{T_i}}} + {\vartheta _i} + {u_i} - {\theta _i}} \right), $

(1)

式中，ΔP为台风中心与边缘大气压强差，ρ为海水密度，g为重力加速度，t为时间，i为8个分潮K1、O1、P1、Q1、M2、S2、N2和M4，T_i为各分潮周期，H_i为分潮振幅，θ_i为分潮相位角，f_i为分潮结点因子，ϑ_i为分潮初始相位，u_i为结点矫正角。

1.3 模型验证

TELEMAC潮位模拟结果精度较高，在之前研究中作者已做过大量验证工作^[23-24]，非台风期间TELEMAC和TOMAWAC耦合模拟的风浪验证见文献[11]，台风期间风暴潮引起的海岸洪水淹没验证见文献[12]。本文仅给出2012年两次袭击广西的台风事件(台风“山神”和台风“启德”)期间实测水位与模拟结果的对比验证。施测单位是广西近海海洋环境科学重点实验室，施测地点在防城港果子山附近(图 2a)。

台风“山神”于2012年10月24日2时在菲律宾东南部西北太平洋洋面生成，10月25日下午进入南海东部海域，10月28日上午减弱为台风，进入北部湾海域，中心附近最大风力为13级，中心最低气压为96 000 Pa，之后台风“山神”以10~15 km/h的速度向西北转偏北方向移动，登陆越南北部地区。台风“启德”于2012年8月13日9时在西北太平洋洋面上生成。8月15日4时在菲律宾吕宋岛登陆，登陆时中心附近最大风力为10级(25 m/s)。8月17日12时在广东省湛江市再次登陆，登陆时中心附近最大风力为13级(38 m/s)，17日21时在中越边境交界处沿海第3次登陆，登陆时中心附近最大风力为12级(33 m/s)。

结果表明采用波-流耦合数值模拟方法能有效再现历史时期台风增水和风浪过程。两次台风期间水位施测时长各为2 d，总共包括96个测量点位。其中误差小于10%的测量点位占据60%以上，误差小于20%的测量点位占据90%以上。图 3为实测水位和相位与模拟结果的相关性比较，红色表示台风“山神”，测量时间段为2012年10月28日21:00—30日22:00；绿色表示台风“启德”，测量时间段为2012年8月17日10:00—19日10:00。1:1线以上部分为高估的模拟点位，1:1线以下部分为低估的模拟点位。总体来说最佳拟合线坡度为0.97和0.95，R²为0.89和0.85，Skill值为0.93和0.91(>0.9)，模拟效果较好。台风“山神”Skill值为0.95，台风“启德”Skill值为0.91，因此台风“山神”模拟效果优于台风“启德”模拟效果。模拟水位略微低估实测峰值水位，可能是因为风浪干扰作用导致测量误差。另外，风场精度将直接影响增水和风浪模拟结果，因此有必要对风场准确性做进一步验证。

2 结果与分析 2.1 水位变化模拟结果

台风“山神”自东南向西北方向横穿南海，当接近北部湾时强风和低气压导致北部湾大规模增水(图 4)；随着台风继续向北移动，水位从外海逐渐向广西沿岸抬升。虽然台风“山神”接近北部湾时西北风最高速度达28 m/s，但最大潮位并不高，仅0.0~0.5 m(图 4a)。10月28日6:00，东偏东南向风速达到最强，为35 m/s，沿广西沿岸潮位达到最高，为1.2 m。整个台风过程中，最高潮位出现在钦州湾和大风江河口附近，分别达到1.0 m和1.2 m。随后台风逐渐向西北方向消减，增水也随着风速减弱而逐渐减小。

图 4模拟结果表明风暴潮增水从广西沿岸向外海逐渐减弱，潮差也向外海逐渐缩小。台风中心附近风速较小，水位也较低，仅1.0 m以内。低水位出现位置与台风中心迁移方向密切相关。最大风速和最大水位出现在台风臂扫过区域附近。影响北部湾最大风速出现在10月28日6:00—12:00，但整体水体较低，仅-0.3 m。最高潮位出现在10月28日18:00左右，整体水位达到0.4 m。因此仅当高潮位与最大风速叠加时易形成超高水位。

2.2 风浪及增水模拟结果

台风作用引起水位变化主要包括风浪和增水两方面。风浪主要受局部风场剪应力驱动形成，增水主要受大范围风场持续驱动下与天文潮交互作用形成。风浪与增水模拟结果如图 5所示。增水高度在0~4 m间变化，台风中心增水较小，最大增水出现在北部湾右侧(图 5g~l)。风浪高度在0~10 m间变化(图 5a~f)，明显高于增水幅度。风浪空间分布与风场空间分布密切相关。台风中心风浪较小，仅1 m以内，而台风臂扫过海域风浪较大，最高可达10 m。

风浪高度和增水高度总和为台风引起总水位变化。风浪波周期较短，而增水具有大范围和长时效性。如图 5g~l所示，增水平均水位变化可达1.2 m，北部湾东侧沿岸浅滩破波可导致增水超过2 m，区域性增水最大在3 m以上。风浪叠加增水如与高潮位相遇易形成海岸洪水淹没。即使增水量级比风浪小很多，但其叠加效应，尤其是高能波浪冲击效应，将对堤防造成巨大压力。

3 讨论

本文建立了风暴潮增水数值模拟系统，采用台风模型(Holland)、潮汐模型(TELEMAC)和波浪(TOMAWAC)模型嵌套和耦合技术实现风浪和增水模拟，研究方法简单直接，具有可复制性。研究区域位于大陆架宽广的北部湾地区，具有典型性。虽然世界上局部地区的大陆架狭窄陡峭(例如台湾东海岸、菲律宾东海岸和智利西海岸等)，但坡度大于0.05%的大陆架仅占世界海岸线的2%不到。宽浅的陆架对潮波传播和风暴潮增水产生重要影响，因此本研究结论具有一定的可移植性。

本文模拟水位误差大致为10%~20%，模拟结果具有一定可信度。台风风场和气压场计算采用Holland台风模型，属于解析模式，方法比较简单，难以刻画台风移动过程中形成风圈不对称性变化的现象。因此最终波浪和增水模拟误差可能主要来源于风场驱动力误差。为提高增水和风浪模拟精度，需进一步提高台风大气场模拟精度。因此更多实测数据率定台风风场将有助于提高风浪和极值水位模拟精度。而且，北部湾计算网格分辨率较为粗糙(800~5 000 m)，难以模拟短周期风浪(10~100 m)变化过程，因此可能会对水位造成低估现象。伍志元等^[9]提出采用中尺度大气模式WRF和区域海洋模式相耦合的方法，在一定程度上提高了风场模拟精度。邹怡杰等^[18]进一步提出采用MCT耦合器实现海-气双向耦合，提高模拟精度。但大气模式和海洋模式的耦合将占用大量计算资源，导致模拟效率低下。

最终模拟所得风暴潮总水位高度包含天文潮高度和风暴潮增水相叠加高度。而大潮高潮位与风暴潮增水叠加时易形成超高水位。北部湾最高潮位出现在10月28日18:00左右，滞后于最大风速出现时间10月28日9:00(图 4~5)，因此台风“山神”期间并未出现历史极端高潮位事件。如果台风“山神”提早3 d或推迟5 d登陆，与大潮高潮位叠加，将造成“两碰头”的局面，水位可瞬间提高1 m以上，将严重增加广西沿岸洪水风险。北部湾大陆架宽浅是造成风暴潮增水的另一个重要因素。地形顶托与摩擦作用易形成潮水和增水间的强烈非线性效应。最大非线性效应出现在北部湾东侧靠近琼州海峡一侧，当地平均水深仅25 m，而潮差和增水高度分别在4 m和3 m以上。因此，天文潮与风暴潮增水交互作用非线性叠加过程既受到潮位和波浪间质量和动量传递作用，又受到底摩擦和科式力作用，是一个复杂的迭代解算过程。

4 结论

台风影响下重点区域极值水位模拟、台风登陆后水位变化模拟和地面淹没模拟仍然具有挑战性。本文基于大气-潮汐-波浪耦合数值模拟方法，采用实测北部湾水下地形和全球发布南海地形数据设计了风暴潮增水数值模型系统，实现北部湾及附近海域风浪、增水和水位模拟。结果表明台风期间北部湾水位变化主要受天文潮波动和风暴潮增水变化影响。本研究意在突破北部湾现有模型局限性，为进一步深化北部湾科学研究和工程应用提供基础数据。本文设计模型的优点是简单有效、区域嵌套和模拟精度高，能够满足一般性台风风浪、增水模拟的后报和预报工作需要。