珠江口伶仃洋风暴潮传播特性及影响因素数值模拟研究——以台风“山竹”为例

陈鹏，张卓,，宋志尧，章卫胜，叶荣辉，李玉婷

（1.南京师范大学 虚拟地理环境教育部重点实验室，江苏 南京 210023；2.江苏省地理信息资源开发与利用协同创新中心，江苏 南京 210023；3.水利部水旱灾害防御重点实验室，江苏 南京 210029；4.珠江水利委员会珠江水利科学研究院，广东 广州 510611；5.水利部珠江河口动力学及伴生过程调控重点实验室，广东 广州 510611）

风暴潮是一种重要的天气现象，会导致沿海地区洪水泛滥、屋舍冲毁、农田淹没等灾害，威胁人民的生命安全且造成巨大的财产损失[1]。珠江口面临南海，是台风的多发区，珠江三角洲又是中国重要的经济区域，风暴潮灾害受到普遍的关注和重视。据统计，仅2018 年台风“山竹”就造成了广东、广西、海南、湖南、贵州5 省（自治区）近300 万人受灾，直接经济损失超50 亿元[2]。因此，风暴潮传播规律的研究对于灾害区域的防灾减灾及应急响应措施的制定是非常必要的[3]。

通过理论推导，Proudman 认为风暴潮增水最大值应出现在天文潮潮位较低时刻[4]；Prandle 等在研究泰晤士河风暴潮增水时发现，风暴潮增水在天文潮涨潮阶段较落潮时偏大[5]。国内，很多学者对我国沿海地区的风暴潮开展了大量研究。在观测研究方面，王永信等利用40 年的风暴潮增水资料，分析了风暴潮沿珠江河道上溯运动问题[6]；史键辉等根据赤湾等站点的实测资料，研究了风暴潮和天文潮之间的相互作用[7]；陈杏文等利用浮标观测数据分析了2021 年台风期间海洋动力学和热力学的响应特征[8]；罗志发等利用1970－2018年共49 年的珠江口实测数据，结合数值模拟结果分析得到各测点历史最高潮位以0.02～0.03 m/a 的平均速率增加[9]；刘媛媛等将风速、风向、气压以及前时序的潮位数据与LSTM 模型结合，建立了风暴潮预报模型[10]。在数值模拟方面，于斌等对风暴潮沿珠江河道上溯进行模拟，发现风暴潮位和风暴强度和路径有很大关系[11]；李杰等发现当台风中心距离验潮站最短时该站通常会出现最大增水[12]；Du 等对珠江口不同台风路径及台风中心移动速度对风暴潮增水和海岸淹没的影响进行比较，发现当台风移动速度越慢时引发的海岸淹没越严重[13]；严枫等对双台风作用下的增水、流场的相互影响及影响区域进行了探讨[14]。近十年来，珠江口风暴潮数值也在逐步完善中，这方面有傅赐福等对风场适应性研究[15]；叶荣辉等的珠江三角洲风暴潮模型[16]；李心雨等研究了WRF 与台风经验模型在台风“山竹”汇总的应用对比，认为台风经验模型在台风模拟中仍是可接受且便捷的方法[17]。在未来气候变化影响研究方面，Yin 等研究了未来海平面上升和台风强度增强对珠江口风暴潮的影响[18]。

总体来说，目前从台风登陆位置和移动方向对珠江口风暴增水时空分布影响的研究不少，但多是实测点的风暴潮增水分析，而采样点数和台风案例有限，无法对影响因素实行单因子定量分析。数值模拟多针对整个珠江口，而针对西进型台风在广东近海改变登陆地点、登陆角度以及登陆时间条件下，对伶仃洋内风暴潮增水的时空分布特征变化进行影响因子量化分析的相关研究较少，更缺乏伶仃洋内天文潮与风暴潮的相互作用对风暴潮增水影响的典型案例。因此，研究选取近年来对广东地区影响较大的一次西进型台风——“山竹”作为典型案例进行数值模拟，分析西进型台风在伶仃洋内的风暴潮分布及变化特征，探讨西进型台风改变条件（如登陆时间、中心气压、登陆位置和移动方向）对伶仃洋风暴潮增水的响应规律，为有关部门的决策制定提供依据。

1 模型介绍

1.1 研究区域和计算范围

本文中，为了获取热带气旋引起风暴潮的整个过程，模型计算采用大、小区双层嵌套的方法。大区范围包括南海和部分西北太平洋区域，模型边界采用主要4 个半日潮（M2，S2，K2，N2）和4 个全日潮（K1，O1，Q1，P1）来估算外海边界上的天文潮。小区范围（图1）包括珠江三角洲河网、八大入海口门以及伶仃洋和黄茅海，简称珠江口模型，模型采用无结构三角网格，网格分辨率在外海为5 km×5 km，河口海湾内加密为100 m×100 m。该模型在北江、西江和东江的上游设置流量边界，外海开边界由大范围计算的水位结果提供。研究的主要区域位于伶仃洋内，该区域口门外受到香港离岛及珠海东澳岛等众多岛屿的掩护，风浪对于潮位的影响相对较小[19]，因此本文不考虑风浪对风暴潮的影响。

图1 珠江口模型网格和地形

1.2 台风模型

选取登陆华南地区的最强台风—“山竹”作为典型案例进行模拟（图2），台风“山竹”在西北太平洋洋面生成，并于2018年9月16日17时在广东台山登陆，登陆时仍有14 级风力。“山竹”引发的大风、风暴潮以及暴雨洪水等对广东、香港和澳门等地造成了相当大的影响，据媒体报道，至少有70人死亡，超过252万人被转移安置。

图2 台风“山竹”行进路径及在登陆时刻图

本文采用Holland[20]模型作为风场驱动，与其他台风气压参数模型相比，Holland 模型引入气压剖面特征参数B 来描述相同的强度和半径上的风速分布。其值可由Holland 模型和梯度风公式结合日本气象厅公布的50 节（25.722 m/s）风速圈半径估算得到。台风“山竹”在移动和发展过程中，估算B的变化范围在0.9～1.0之间，本文取B=1.0。

1.3 风暴潮模型

FVCOM 是一个非结构网格、有限体积的海岸和海洋模型，最初由陈长胜团队开发[21]，并由UMASDSD 和WHOI 研 究 团 队 改 进[22-23]。基 于FVCOM 主程序构建的风暴潮模型，采用干湿处理和通用湍流模块，其中驱动力有外海潮边界，上游径流以及海表面风应力。本文构建中国海-珠江口风暴增水计算模型，该模型既能模拟天文潮和风暴潮的相互作用，又可以模拟波流相互作用，能够得到比前人研究更合理的结果[24]。珠江口模型的外海边界可由中国海模型的结果提供，上游径流洪水为西江的高要、北江的石角和东江的博罗等3个水文站当日的实测径流量。

风应力是海气能量交互的关键，在数值模拟中风应力一般按海表面10 m风速的2次方来计算：

式中：τ为海表面切应力；ρa为空气密度；Cd为拖曳力系数；U10为海面10 m风速。其中风拖曳力系数采用Wu[25]基于实测资料得到的经验公式计算得到，公式如下：

2 模型验证与结果分析

2.1 风场验证

对风场的验证采用收集到的台风登陆前后的3个实测站位（W1、W2和W3）的风速风向资料，站位分布见图3。其中W1 站点位于澳门机场东南侧，W2 号站点位于磨刀门口门处，W3 号站点位于黄茅海湾口处。由验证结果可知，“山竹”台风期间风速风向的模拟结果较好（图4），说明本文选取的台风风场计算方案合理。

图3 珠江口模型风速风向及水位测站分布图

图4 台风登陆前后风速风向计算值（线）与实测值（点）验证结果

2.2 水位验证

风暴潮模拟计算从2018 年9 月14 日0 时至2018 年9 月18 日0 时，历时4 天，涵盖了台风“山竹”从临近、登陆和消亡的全过程。采用6 个潮位测站（图3）的实测资料对珠江口风暴潮模型进行验证，通过对糙率和台风参数的率定，除西侧万顷沙站的风暴增水退潮阶段相对误差较大外，其他站位模拟结果与实测值吻合较好（图5）。万顷沙站退潮阶段误差较大的原因可能是没有考虑珠江三角区内台风暴雨引发的汇流洪水，也可能是该站位上游的河道地形有误所致。

图5 台风期间水位过程计算值（线）与实测值（点）验证结果

2.3 台风期间的风暴增水分析

伶仃洋南部湾口宽约35 km，北部湾顶宽仅有5 km，呈口外宽里湾窄的喇叭状。在喇叭状海湾中，风暴潮从口外深海向浅水湾顶传播时，受到河口断面收缩和底摩擦的影响，潮波能力向湾顶聚集，风暴增水被放大，从而在湾顶形成极端高水[26]。珠江口模型计算的风暴潮总水位减去天文潮可以得到风暴增水几个关键时刻的空间分布（图6）。台风登陆前8小时，东北风扫过海湾，吹向西南，导致伶仃洋湾顶和虎门的风暴增水为负。随着台风中心向珠江口靠近，风向逐渐沿顺时针旋转，16 日13－16 时，在风力驱动下水体向伶仃洋西海岸堆积，伶仃洋出现西高东低的增水分布。登陆时风力达到最大，西北风向几乎与海湾走向平行，海岸及河口附近的多个观测站的增水均达到了最大值。登陆后台风强度明显减弱，风暴增水也随之减弱。

图6 台风“山竹”登陆期间的增水过程

2.4 风暴潮水位过程线分析

为研究台风登陆前后风暴增水过程的变化，在伶仃洋的湾口、中部、湾顶及虎门河口选取12个特征点（图7（a）），记录特征点逐时的风暴潮位变化。从特征点的风暴潮水位过程线图（7（b））可以看出，12 个特征点的风暴潮过程曲线相似，可分为初振段、主振段和余振段。特征点所在位置的不同，各曲线显示出不同的特征。

图7 伶仃洋和虎门口12个特征点在各阶段增水过程（台风在第65 h登陆）

初振段之前，一系列较弱增水波传入伶仃洋，持续时间约30 h（第17～48 h），振幅约为0.2 m（图7（c）、（f））。这种先兆增水是由大洋或外海远程风引发的边缘波形成的，在靠近海湾时会略有增强。在第48～60 h 的增水上升段，湾口附近的1～5点（图7（c））和湾中至湾顶6～12点（图7（f））的水位线特征不同，主要区别在于后者上升段叠加了更强的短期振荡，表明在湾内天文潮对风暴增水的影响更大。在增水的峰值阶段（图7（d）、（g）），由于伶仃洋为喇叭状，风暴引起的风暴增水可以无阻挡地向湾内堆积，从湾口至湾顶的峰值增水逐渐加大。同时，湾内不同点之间的增水峰值出现的时间也有明显差异。因伶仃洋是开敞式湾口，增水水位沿横向坡降较小，如1 和2、6和8、9和10点，而4点和5点由于受到岛屿遮挡，此两点的峰值增水小于3 点。主增水阶段过后，湾口附近（1～5）和湾内（6～12）之间的风暴增水下降段有所不同（图7（e）、（h））。与上升段相似，特征点在湾内的短期波动比海湾口附近的更强烈。在峰值增水出现约10 h后，第75～77 h，湾内各点出现明显的增水次高峰，而湾口附近各点则不明显。

从对图6、图7 的分析中，可以得到珠江口风暴潮的时空分布受台风风场和海岸几何形态及水深地形的共同影响。由于海湾开口宽阔，湾内外水量交换通畅，故12 个特征点的水位过程线形态趋于一致，而伶仃洋湾内比湾口附近点的水位过程线的波动现象更强烈，这是地形浅化效应和天文潮与风暴潮相互作用的结果。

3 珠江口伶仃洋风暴潮影响因子分析

3.1 台风登陆时间对风暴增水的影响

台风“山竹”引发的风暴增水在9月16日17时（台风登陆时刻）到达峰值，此时的天文潮处于落潮阶段。台风的不同登陆时间会对风暴增水造成的影响与天文潮和风暴增水的相互作用有关，也是风暴潮研究的焦点问题。Proudman 对天文潮与风暴潮相互作用对风暴增水的影响进行了理论研究，发现在高潮位附近的风暴增水总是小于低潮附近的风暴增水[4]，为验证珠江口风暴潮是否遵循该规律，设计如下试验：将台风“山竹”的时间序列分别提前1、2、3 和4 h，以+1、+2、+3、+4来表示，和滞后1 和2 h，以-1、-2 来表示，利用修正后的风场进行数值试验。取代表湾口的4点和代表湾内的11 点的总水位线过程线来说明天文潮和风暴潮相互作用的影响。如图8 所示，总水位最大值与潮位成正比，这意味着当风暴潮与潮位同步时可能会出现极大的总水位。而去除掉天文潮后的风暴潮增水却与潮位并非成正比（图8（c）、（d）），而与之相反，在高潮位（+2）登陆的风暴增水小于在低潮（0、-1 和-2）登陆的风暴增水。通过对4 点和11 点的比较，表明湾内的风暴增水符合Proudman 的分析结论[4]，即在高潮位出现风暴增水最小，并且湾顶附近的增水受天文潮的影响更加显著。经过分析认为，这种由于风暴潮和天文潮之间相位差造成的增水差异，主要源于非线性阻力和浅水效应。因此，在对风暴潮增水进行预报时，尤其在湾内浅海区域，必须考虑台风登陆期间天文潮相位的影响。

图8 4点和11点风暴潮总水位（上）和增水（下）的时间变化

3.2 台风中心最低气压的影响

台风强度与中心气压差成正比，为了定量分析不同的中心气压差对增水大小的影响，将“山竹”台风的中心气压原始数据进行了±10 和±20 hPa 的偏移（图9（a）），对四种设计情景的风暴增水过程进行模拟。从4 点（图9（b））和11 点（图9（c））设计试验情景和原始台风之间的总水位对比表明，增大压差（-10，-20）既能增大正、负风暴潮位的峰值，又能延长高水位持续时间。反之，则峰值水位降低，且高水位持续时间缩短。对比4点和11点的水位过程线，中心气压差变化对湾内11点的影响比湾口的4点更为显著，说明风暴潮增水过程对同一台风气压差变化在空间上的分布是不均匀的。一般而言，水深较浅的沿岸和湾内的风暴增水对台风强度（或气压差）的变化更为敏感，这与浅水效应和湾顶潮能放大效应有关。

图9 4点和11点在不同台风中心气压时间序列下总水位对比

3.3 登陆地点的影响

有学者探讨了台风风速强弱、尺度大小以及海岸线的几何形态等因素对伶仃洋风暴潮的影响，但对于西行台风登陆位置的影响探讨较少。下面以台风“山竹”为例，研究登陆位置的变化对伶仃洋风暴潮的影响。实验设计是将“山竹”台风的原路径向左移动4 次、右移动1 次，形成5 条新路径，每条台风路径相差0.5°（图10）。

图10 台风设计路径与原始路径及登陆位置图

利用中国海-珠江风暴潮模型，模拟5 条设计路径的台风登陆时伶仃洋风暴增水的空间分布，并与原始路径进行比较。模拟结果表明，台风路径向西移动0.5°和1°时，风暴潮增水的强度和范围均较大。继续向西移动幅度达1.5°和2°时，风暴潮增水转向逐渐减弱。从图11（b）看出，路径西移0.5°时伶仃洋内产生的增水最大，超过原路径，因为此时台风最大风速半径正好位于伶仃洋口门中心区域，大量海水在风的驱动涌入湾内。相对地，原路径向东移动0.5°时，伶仃洋内的风暴增水强度则急剧减弱，这是因为最大风速半径已经位于伶仃洋湾口的东部，东移后湾口处的风力比原路径更弱。图12 显示了台风路径变化对风暴潮增水过程的影响，可以看出登陆位置的不同对增水过程有很大的影响。具体而言，路径西移使增水上升段提前到来，路径东移则推迟增水上升段的时间，且增水下降段路径东移和西移存在较大差别，分析原因是路径西移使最大风速半径圈提前经过伶仃洋湾口，因此提前引发增水，东移则相反。

图11 伶仃洋湾最大风暴增水分布

图12 4点和11点在不同台风路径下的风暴潮增水过程线

3.4 登陆角度的影响

根据1968－2017 年的统计资料，对珠三角区域影响较大的台风主要是西进型台风，即以相对平直的西北向移动路径靠近海岸并登陆的台风类型。这种台风的特点是路径长且行进速度快，诞生于北太平洋中西部，较长的生命周期使其在海面充分发展，等到临近登陆时往往形成超强台风，2018 年的超强台风“山竹”就是其中一个典型案例。西进型台风的路线往往比较平直，因此可以用一个登陆角度（登陆时台风路径跟海岸线切线方向的夹角）来代表台风的移动方向。为了研究登陆角度变化对于伶仃洋内风暴潮增水的影响，在台风“山竹”原路径的基础上，将整个路径分别沿逆时针（减小登陆角度）和顺时针（增大登陆角度）方向旋转5°，得到两个设计台风路径（图13）并保持其他台风参数不变。从图14 的模拟结果可知，减小登陆角度使伶仃洋风暴增水增加，而增加登陆角度使伶仃洋风暴潮增水明显减小。分析原因，认为广东沿海海岸线基本沿东西向，当台风登陆时与海岸的夹角减小时，台风最大风圈跟伶仃洋之间的距离更小，从而会增大伶仃洋的最大风暴潮影响的面积。从图15 可以看到，当登陆角度变小，使风暴潮增水前期（60 h以前）增水变小，而在峰值及峰值之后增水增加。增大登陆角度，风暴潮前期增水变大，但在峰值之后增水明显变小，这都是由于登陆角度改变伶仃洋风场过程造成的。

图13 原始路径（方案1）与设计路径（方案2、3）分布图

图14 原始路径（a）与设计路径（b）、（c）风暴潮最大增水空间分布

图15 4点和11点不同路径下风暴增水过程线

4 结论

本文采用FVCOM 数值模型，以近十年来登陆广东省的超强台风“山竹”为例，对伶仃洋及其附近海域的风暴潮过程进行模拟和分析。模型计算结果与珠江口6 个台站的实测值吻合较好。伶仃洋风暴潮增水的时间过程及空间分布表明，伶仃洋的风暴增水过程受风向、岸线几何形态和局部水深地形的影响。伶仃洋的12 个特征点的水位过程线表明，风暴潮的时间过程曲线可分为初振段、主振段和余振段，湾内各点的风暴增水表现出比湾口附近和湾外更强的短期振动，表明湾内天文潮与风暴潮相互作用比湾口和湾外更强。

文章研究了台风的登陆时间、中心气压差、登陆位置和移动方向等特征参数对风暴潮增水的影响。结果表明，不同登陆时间变化通过天文潮与风暴潮的相互作用导致了风暴增水的变化，且在湾内变化更为明显，符合低潮位时风暴增水更大的规律。中心最低气压的变化影响增水的峰值大小，当气压差增大时，可以同时增大增水和减水的峰值，但对风暴潮增水的过程线影响较小。台风登陆位置对伶仃洋的风暴增水的时空过程分布有很大影响，当台风在伶仃洋西侧海岸登陆，且风向近似垂直于海岸时，可能会出现极端的风暴增水。对于西北型登陆的台风，当台风登陆方向与海岸夹角减小时，会增大风暴增水的峰值，反之，增大角度则会削弱增水峰值。本文的研究结果有助于对珠江口伶仃洋及其周边海区典型的西北向移动台风风暴潮过程的认识，可为沿海防台、防潮减灾政策制定提供参考。