珠江河口波浪—风暴潮耦合数值模拟

罗志发，黄本胜，谭 超，邱 静，黄广灵

(1.广东省水利水电科学研究院，广东 广州 510635；2.广东省水动力学应用研究重点实验室，广东 广州 510635；3.广东省流域水环境治理与水生态修复重点实验室，广东 广州 510635)

1 概述

珠江河口位于中国大陆南端，濒临南海，海岸线狭长，沿岸地势较低。特殊的地理位置及气候条件使得珠江河口成为国内台风登陆最频繁、风暴潮灾害最严重的区域之一。近年来，在全球气候变化背景下，风暴潮灾害次数和强度呈增加的趋势，登陆珠江河口的超强台风引发的风暴潮灾害，屡次刷新沿岸潮位站历史最高潮位和最大增水记录。珠江河口沿岸人口稠密、经济发达，每年由于风暴潮造成的经济损失相当巨大，已经成为影响人民生活质量、制约国民经济高质量发展的重要因素。

为了防御风暴潮灾害，珠江河口沿海陆续兴建了江海堤围。尽管如此，强台风往往伴随着强台风浪，在强浪冲击作用下，江海堤围面临着堤身结构损毁的风险，可能导致海水漫过海堤入侵内陆，引发城市内涝。因此，准确模拟风暴潮及台风浪过程，了解正面登陆台风影响下珠江河口近岸海域风暴潮增水及台风浪变化过程，对准确评估风暴潮的灾害风险以及做好风暴潮灾害防御与应对工作、减少风暴潮灾害损失有着重要现实的意义。

珠江河口是多种动力因子协同作用的复杂系统，呈“三江汇流，八口入海”的形势。河网区水网密布，横向支汊发育，其间多种动力因子相互耦合，动力复杂多变。珠江河口风暴潮、波浪数值模拟已开展较多的研究，如二维风暴潮模型的建立[1]，风浪数值模拟[2-3]，路径、风速对风暴潮增水的影响[4]，地形对局部增水的影响[5]等。波浪在近岸传播过程中受到浅水变形的影响发生破碎，产生的辐射应力驱动水面抬升形成增水，在台风期间波浪增水相对较大[6]。波浪与风暴潮耦合计算已开展了相关的研究，在不同的海域，由于岸线地形及水动力条件的不同，波浪对风暴潮增水的影响有所差异[7-9]。针对珠江河口复杂的动力系统，本文基于无结构三角网格的数值模式，优化拟合复杂岸形，以珠江河网—河口作为整体，综合考虑风、径流、潮流、波浪等动力因子的耦合作用，构建珠江河口波浪—风暴潮耦合数值模型，可有效提高风暴潮数值模拟的准确度，为进一步探讨风暴潮—波浪耦合作用提高有效的科学手段。以1822号“山竹”台风为例，对台风过程的风暴潮增水及波浪进行数值模拟，定量分析考虑波浪—风暴潮耦合作用对风暴潮数值模拟结果的影响。

2 模型的建立与验证

2.1 模型简介

本文选用海洋环流模式SELFE建立南海—珠江河口双重嵌套模型(模型网格见图1)，该模式基于非结构化网格，可精细化拟合复杂岸线和地形，采用半隐式的欧拉—拉格朗日有限元算法求解N-S方程组，可有效提高模型计算效率。南海模型计算范围为98°E～126°E，0°N～30°N，涵盖整个南海及西北太平洋海域，南边界至卡里马塔海峡、北边界至浙江省沿岸海域、东边界至48小时警戒线。网格分辨率从近岸的 1 km 逐渐过渡到外海20 km，水深数据采用ETOP01全球1′×1′分辨率的地形资料。珠江河口模型上边界为西江高要，北江石角，东江博罗，流溪河及潭江上游，外海下边界取约100 m等深线处。模型采用无结构网格，拟合复杂河岸边界，对局部进行加密，提高网格分辨率。模型网格共有101 752个节点，173 045个网格单元，网格大小为从河网区10 m逐渐过渡到外海的20 km；在垂向上采用Sigma坐标，均匀分为10层；珠江三角洲网河区采用2005—2008年的大范围实测地形，河口区及近岸海区采用2000—2008年海图地形。南海模型由风场、气压场及8个主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1)[10]进行驱动，计算得到余水位，将余水位以及风场、气压场、8个分潮作为珠江口模型的驱动条件，对珠江河口风暴潮增水进行数值模拟，模型详细介绍见文献[7]。

图1 模型嵌套计算网格示意

在风暴潮模型的基础上耦合波浪计算模块，建立波浪—风暴潮耦合模型。波浪模块运用Wind Wave Model-III(WWM-III)模式[11]，该模式主要用于模拟大洋—近岸波浪传播过程，模拟的物理过程包括：波浪折射作用，浅水变形作用，风成浪作用，白浪的耗散作用，水深引起的破碎作用，海底摩擦作用等。控制方程为能量谱平衡方程：

(1)

式中N为波作用密度谱;Cgx和Cgy分别为x和y方向的波群传播速度;u和v分别为x和y方向的水流速度;σ为相对波频;θ为波向;Cσ和Cθ分别为σ和θ方向的波浪传播速度;Stot为谱密度表示的源项，是各种物理现象的源函数的叠加形式。

风暴潮模型为波浪模型提供水位、流速、流向数据计算变化水动力条件下的波浪要素。波浪模型输出有效波高、平均波向、有效波周期、波长和底部波轨流速，输出结果作为边界波浪条件输入风暴潮模型，计算波流共同作用下的底切应力及垂向混合系数。

2.2 模型验证

文献[12]以“山竹”台风为算例，对模型计算的风速、风向、天文潮位、风暴潮位等进行了率定验证，验证结果良好，本文在此基础上对波浪计算结果进行进一步验证。珠江口波浪模型外海波浪开边界条件通过南海模型提供，海表面台风场及气压场通过台风模型及气压模型给出[12]，底摩擦系数设置为0.067[13]，波浪破波系数为0.78[11]。波浪计算时间步长设置为1 800 s，每个时间步长计算结束后与风暴潮模型进行数据反馈。模拟时间为2018年9月10日—9月20日，前5 d用于计算的稳定，后5 d对风暴潮及波浪过程进行模拟。

图2为外海浮标站(113.9947°E，21.4948°N)有效波高验证示意，由图2可知，计算结果有效波高变化过程与实测数据有效波高变化趋势一致，相关系数为0.862，平均绝对误差为0.53 m，二者吻合较好。计算有效波高的最大值与实测结果分别为9.4 m、9.9 m，相对误差仅5.3%，表明模型计算结果能较好地反映有效波高的极值情况。由图3可知，考虑波浪作用之后，潮位数值模拟结果更接近观测值。横门站考虑波浪作用后，最高潮位的绝对误差由0.24 m减小为0.06 m；南沙站的模拟精度同样有所提高。由此可见，考虑波浪—风暴潮耦合作用后可提高风暴潮数值预报的精度。综上所述，模型计算结果能够较好地反映“山竹”台风期间波浪及风暴潮增水的变化过程，可为进一步分析波浪与风暴潮耦合提供基础。

图2 外海浮标站有效波高验证示意

图3 横门(a)、南沙(b)站考虑和不考虑波浪作用潮位模拟结果示意

3 结果与讨论

本节利用上文已建立的波浪—风暴潮耦合数值模型，模拟“山竹”台风风暴潮增水及波浪过程，讨论波浪—风暴潮耦合对风暴潮增水数值模拟的影响。

3.1 波浪及风暴潮增水计算结果分析

“山竹”台风过程对珠江口近岸波浪要素分布影响显著，2018年9月16日13：00，台风中心位于珠江口南侧(见图4a)，台风中心风力较小，有效波高较低，约为6～7 m；台风中心右侧风力较大，有效波高超过9 m。此时，伶仃洋海域为东北偏东风，珠江口有效波高相对较小，约为2～4 m。16日17:00，台风于广东台山海宴镇登陆，珠江口盛行东南偏南风，有效波高有显著增大趋势，主要高值区更偏向近岸，有效波高最大值介于4～5 m(见图4b)。选取灯笼山、南沙、万顷沙、横门站位(站位位置见图4a)分析台风期间有效波高变化过程。由图5可知，台风过境期间，各站位有效波高呈现先增大后减小的过程。磨刀门往外延伸濒临南海，因此灯笼山站受波浪影响较大，最大有效波高达3.2 m。南沙、万顷沙、横门等站位均位于伶仃洋湾内，外海波浪传入伶仃洋经底摩擦、折射、绕射等作用，波能消耗，有效波高均小于2 m。

图5 各站位台风期间有效波高变化

台风登陆前4 h，伶仃洋海域为东北偏东风，伶仃洋东岸在离岸风的驱动下，产生减水，深圳湾区域最大减水值约为-1.5 m。伶仃洋西岸由于水体的横向堆积，产生0.5～1.0 m的增水(见图4c)。磨刀门、黄茅海区域为偏北风，有利于水体离岸输运，形成-0.5～-1.0 m的减水。台风登陆后(见图4d)，珠江口海域普遍为东南偏南风，珠江河口东南向的开口方向有利于水体向岸堆积并沿河道向上游输运，此时珠江口海域普遍达到增水的最大值。伶仃洋河口湾顶增水值较大，约为3.1 m，是由于伶仃洋河口湾喇叭状的形态有利于水体向湾顶聚集。其余口门区域最大增水值普遍为2.6～3.1 m。

图4 “山竹”台风登陆前后波浪场及风暴潮增水分布(a、c：登陆前，b、d：登陆后)

3.2 波浪对风暴潮影响分析

利用风暴潮与波浪耦合计算的总水位减去未考虑波浪作用的风暴潮水位可得到波浪增水，波浪增水包含了由波浪辐射应力产生的增水和波流相互作用产生的增水。图6为台风登陆前后波浪增水分布示意。

台风登陆前(见图6a)波浪增水影响范围较广，波浪增水由外海向近岸逐渐增大，外海波浪增水约为10 cm，近岸海域波浪增水约为15～20 cm，与相关研究结果较为一致[9]。近岸海域中，伶仃洋、香港岛附近海域波浪增水较大，黄茅海等海域波浪增水相对较小，其空间分布与有效波高空间分布大体一致。台风登录后(见图6b)，波浪增水集中在近岸海域，主要分布在伶仃洋、磨刀门及黄茅海等海域近岸区域，波浪增水最大值约20 cm。

图6 “山竹”台风登陆前后波浪增水分布示意

选取泗盛围、南沙、横门、灯笼山等站位分析风暴潮增水及波浪增水过程，由图7可知，各站位风暴潮增水呈单峰型，增水快速增大至最大值而后迅速下降，最大增水值均超过3 m。波浪增水同样呈现先增大后减小的过程，最大波浪增水提前于最大风暴潮增水。各站位最大波浪增水均超过0.2 m，约占风暴潮增水的6%。

图7 各站位风暴潮增水与波浪增水过程示意

波浪增水占风暴潮增水的比例虽然不大，但对于设计潮位具有一定的影响。表1列举了珠江口潮位站不同频率下的设计潮位值，由表1可知，同一测站其不同频率的设计潮位仅相差20～30 cm，如西炮台、南沙、灯笼山200年一遇潮位与100年一遇设计潮位值分别相差0.17、0.25、0.22 cm，与最大波浪增水值相当。考虑波浪—风暴潮耦合计算的方法，对推算沿岸设计潮位具有一定的参考意义。

表1 各测站潮位设计值

4 结语

1)本文在已建立风暴潮数值模型的基础上，耦合波浪模块，建立了珠江河口波浪—风暴潮耦合数值模型。采用实测波浪资料对“山竹”台风波浪过程进行模拟验证，计算有效波高平均绝对误差为0.52 m，最大有效波高相对误差仅5.3%，模型计算结果较好地反映了台风期间的波浪过程，可用于珠江河口波浪—风暴潮耦合数值模拟研究。

2)基于本文建立的波浪—风暴潮耦合数值模型，以1822号台风“山竹”为例，讨论了波浪对风暴潮增水的影响，结果表明波浪对风暴潮增水的影响仅限于近岸海域，最大值为10～20 cm，占风暴潮增水的6.5%。考虑波浪—风暴潮耦合作用后，计算结果更接近观测值，可提高风暴潮数值预报的精度，同时对沿岸设计潮位的确定具有一定的参考意义。