粤港澳大湾区风暴潮数值模型的建立与应用

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

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

粤港澳大湾区是我国开放程度最高、经济活力最强的区域之一，以不足全国1%的土地面积和不足全国5%的人口，贡献了全国经济总量的17%，在国家发展大局中具有重要战略地位。粤港澳大湾区位于中国大陆南端，濒临南海，地处珠江流域下游，河网水系发达。由于特殊的地理位置和气候，导致大湾区台风暴潮灾害易发频发，由此造成的人员伤亡、经济损失相当巨大，已经成为影响人民生活质量、制约国民经济高质量发展的重要因素。开展粤港澳大湾区风暴潮研究既是重要科学问题也是风暴潮预警预报关键技术研发的基础，对提升大湾区水安全保障能力具有重要的意义。

粤港澳大湾区所在的珠江河口是多种动力因子协同作用的复杂系统，呈“三江汇流，八口入海”的形势。河网区水网密布，横向支汊发育，其间多种动力因子相互耦合，动力复杂多变。珠江河口风暴潮数值模拟已开展较多的研究，如二维风暴潮模型的建立[1]，路径、风速对风暴潮增水的影响[2]，地形对局部增水的影响等[3]。针对珠江河口复杂的动力系统，本文基于无结构三角网格的数值模式，优化拟合复杂岸形，通过二重嵌套的方式合理提供外海余水位边界，构造台风风场和气压场作为模型的海表面边界条件，构建珠江河口三维风暴潮数值模型。对多组风暴潮增水过程进行数值模拟，验证结果较好。选取1822号“山竹”台风的风暴潮增水过程进行数值模拟，研究天文潮与风暴潮非线性作用对风暴潮增水的影响。

1 模型的建立

1.1 台风场的构建

台风是风暴潮模拟的重要因子，台风风场的质量决定了风暴潮水位模拟的准确性。利用台风风场模型可以较好地模拟台风中心区域的气旋风，但是离台风中心较远的区域，由于气旋风场的减小，背景风场的影响逐渐增大，需要利用再分析风场予以弥补。再分析风场能很好的反映大范围的风场结构特征，但由于其空间分辨率较低，难以刻画台风中心的风场结构，其风速值显著小于实际风速。针对台风风场模型和再分析风场资料存在的不足，本文结合两者优势进行台风风场的构造，即克服了经验台风公式外围风场计算结果较小的问题，也弥补了再分析风场台风中心风强不足的缺陷。

台风风场模型根据梯度风原理，由台风气压场计算出风场。本文采用国内外应用较广泛的Fujita[4]、Ueno[5]公式构建气压场和台风风场，其式如下：

(1)

·[(x-x0)sinθ+(y-y0)cosθ]

(2)

·[(x-x0)cosθ-(y-y0)sinθ]

(3)

式中P(r)为距离台风中心r处的气压值，hPa；P∞为台风外围无穷远处的大气压，取1 010 hPa；P0为台风中心气压，hPa；R为台风最大风速半径，km，根据经验公式计算：R=Rk-0.4×(P0-900)+0.01×(P0-900)2，Rk为经验常数，介于30～60，本文取40；r为计算点离台风中心的距离，km；Vdx、Vdy分别为台风移动速度在x、y的分量，m/s；f为科氏力参数；ρα为空气密度，取值1.292 9 g/m3；ΔP=P∞-P0为台风中心气压示度，hPa；x0、y0为台风中心坐标；θ为台风流入角，取20°；C1、C2为订正系数，本文取0.8。文中构建台风风场模型需要的台风相关数据来自于中央气象台台风网发布的信息。

外围背景风场采用欧洲中期数值预报中心的ERA50再分析风场，分辨率为0.125°×0.125°。台风模型风场与ERA50再分析风场合成方法如下：

(4)

1.2 风暴潮模型的建立

风暴潮增水是大尺度的动力过程，模拟好风暴潮需要建立大范围的模型，同时，珠江口区域精细化模拟需要构建高分辨率的计算网格，为了解决模拟范围、网格分辨率、计算效率的问题，本文采用两重嵌套的方法进行模拟计算(见图1所示)。选用海洋环流模式SELFE建立风暴潮模型，该模式基于非结构化网格，可精细化拟合复杂岸线和地形，采用半隐式的欧拉—拉格朗日有限元算法求解N-S方程组。模式最大的特点是减小CFL条件的限制，在保证计算结果准确的前提下，可以适当放大时间步长，以达到计算精度和计算效率的双赢。

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

大范围的南海模型计算范围为98°E～126°E，0°N～30°N，涵盖整个南海及西北太平洋海域，南边界至卡里马塔海峡、北边界至浙江省沿岸海域、东边界至48 h警戒线。网格分辨率从近岸的1 km逐渐过渡到外海20 km，水深数据采用ETOP01全球1′×1′分辨率的地形资料。珠江河口计算范围为111°E～ 116°E， 21°N～ 23.7°N，模型上边界为西江高要，北江石角，东江博罗，流溪河及潭江上游。外海下边界取约100 m等深线处。模型采用无结构网格，拟合复杂河岸边界，对局部进行加密，提高网格分辨率。模型网格共有101 752个节点，173 045个网格单元，网格大小为从河网区10 m逐渐过渡到外海的20 km。模型在垂向上采用Sigma坐标，均匀分为10层。模型中珠江三角洲网河区采用2005—2008年的大范围实测地形，河口区及近岸海区采用2000—2008年海图地形，外海海域采用ETOP01全球海洋地形。

南海模型由风场、气压场及8个主要分潮(M2、S2、N2、K2、K1、O1、P1、Q1)进行驱动，计算得到余水位，将余水位以及风场、气压场、8个分潮作为珠江口模型的驱动条件，对珠江河口风暴潮增水进行数值模拟。

2 模型验证

2.1 台风场验证

采用澳门气象站(站位分布见图2)“山竹”台风实测风速风向资料对台风模型计算的风速风向结果进行验证，结果如图3所示。

图2 水文气象站位分布示意

图3 澳门气象站“山竹”台风过程风速风向验证结果示意

计算的台风风速、风向其变化过程与实测数据相符合，最大风速绝对误差仅1 m/s。说明计算的风场较好地反映了台风过程及最大风速。图4为采用本文台风场构建方法计算的台风场图，由图4可知，在经验模型的基础上融合了ERA50的再分析数据，对台风场的计算避免了经验公式关于台风外围风速计算明显偏小的问题，构建的风场结构较为合理，可用于风暴潮数值模拟。

图4 2018年9月16日6：00“山竹”台风风场示意

2.2 天文潮验证

模拟风暴潮和天文潮耦合的水位过程时，天文潮的模拟计算是基础。在不考虑风应力的情况下，对南海天文潮进行数值模拟，以检验模型对天文潮计算的可靠性。计算时间为2018年7月1—31日，选用沿海大万山、珠海港、上川岛、马鞭洲等潮位站(站位分布见图2)实测资料对模型计算结果进行验证，比较结果见图5所示。

由图5可知，各潮位站计算的天文潮水位过程线，无论是高低潮位还是相位与实测数据基本一致，潮位计算绝对误差介于8.8～18.4 cm，验证结果良好。表明模型开边界条件基本正确，底摩擦参数设置合理，可在此基础上进一步开发风暴潮模式，为珠江口风暴潮模型提供合理的余水位开边界条件。

图5 南海模型天文潮验证示意

珠江口模型已经过2001年、2005年、2008年、2009年等多组水文观测资料的验证，并成功应用于珠江河口三维水动力及盐度输移的数值模拟，模式简介及天文潮、水动力验证结果详见文献[6]。

2.3 风暴潮水位过程验证

以“山竹”台风为例，对珠江河口风暴潮模型进行验证。模型上游高要、石角、博罗、人和给定实测径流数据，石咀因缺乏实测流量数据而给定当月份的多年平均值。外海开边界水位通过天文潮水位+余水位的形式给定，天文潮水位由调和常数计算得到，余水位由南海模型提供。海表面台风场及气压场通过台风模型及气压模型给出。模型时间步长设置为200 s，模拟时间为2018年9月10—20日，前5 d用于天文潮计算的稳定，后5 d对风暴潮过程进行模拟。

选用珠江口门及河网区水文站位(泗盛围、南沙、万顷沙、横门、灯笼山、黄金、西炮台，站位分布见图2)实测水位资料对风暴潮计算结果进行验证。由于篇幅有限，仅给出部分验证结果图(见图6)。由图6可知，计算结果水位变化过程与实测数据水位变化趋势一致，计算结果与实测结果吻合较好。统计了模型计算最高潮位的绝对误差以及相位误差(见表1)，由表1可知，“山竹”台风计算结果的最高潮位绝对误差介于0.15～0.49 m，相对误差控制在25%以内，最高潮位的相位误差均在1 h以内。综上分析，模型计算结果较为合理，所建立的模型能够较好地反映珠江河口径流、潮流、风等多种动力因子相互作用的结果，能够较为准确地模拟珠江河口的风暴潮过程。

图6 “山竹”台风过程水位验证示意

表1 误差统计

3 结果分析

本节利用上文已建立的风暴潮模型，模拟“山竹”台风风暴潮增水过程，讨论天文潮与风暴潮非线性相互作用的特征。

3.1 风暴潮增减水特征分析

2018年9月7日20:00，台风“山竹”在西北太平洋洋面上生成；9月15日，台风“山竹”从菲律宾北部登陆；16日17:00在广东台山海宴镇登陆，登陆时中心附近最大风力14级，中心最低气压为955 hPa。9月16日13:00，台风中心位于珠江口南侧(见图7a)，伶仃洋海域为东北偏东风，伶仃洋东岸在离岸风的驱动下，产生减水，深圳湾区域最大减水值约-1.5 m。伶仃洋西岸由于水体的横向堆积，产生0.5～1.0 m的增水。磨刀门、黄茅海区域为偏北风，有利于水体离岸输运，形成-0.5～-1.0 m的减水。

台风登陆后(见图7b)，珠江口海域普遍为东南偏南风，珠江河口东南向的开口方向有利于水体向岸堆积并沿河道向上游输运，此时珠江口海域普遍达到增水的最大值。伶仃洋河口湾顶增水值较大约3.1 m，是由于伶仃洋河口湾喇叭状的形态有利于水体向湾顶聚集。其余口门区域最大增水值普遍为2.6～3.1 m。

注：a、b分别为2018年9月16日13时以及2018年9月16日17时珠江口风暴潮增减水分布，红色线条为台风路径，红色圆点为台风中心。

3.2 天文潮风暴潮非线性相互作用分析

对于风暴潮增水而言，除了气象因子之外，天文潮与风的非线性作用对风暴潮增水也是重要的影响因子。设置3个数值试验：考虑天文潮风的耦合作用、仅考虑潮的作用、仅考虑风的作用。利用天文潮与风耦合计算的总水位减去天文潮位可得到风暴潮增水，风暴潮增水包含了纯风生增水(仅考虑风的作用)和非线性增水。利用风暴潮增水减去纯风生增水后可得到非线性增水。在珠江口区域设置4个计算点(见图7a)，分别位于虎门(P1)、伶仃洋西岸(P2)、磨刀门(P3)、黄茅海河口湾顶(P4)，分析各计算点不同水位的变化特征。

由风暴潮水位的时间序列(见图8)可以看出，由于风暴潮最大增水发生在小潮的低水位阶段，因此总水位最大值主要是由风暴潮最大增水决定。除了P2点外，其余各计算点均呈现先减水后增水的特征，其中最大增水值位于P1点约3.1 m，与伶仃洋河口湾喇叭状的形态有利于水体聚集有关；最小减水值位于P3点约-2 m，与黄茅海河口湾的开口方向和台风登陆前的风向有关。最大减水发生在台风登陆前1～2 h，最大增水值发生在台风登陆后1～4 h。由于P2点位于伶仃洋西岸，台风登陆前后的风向均有利于水体的向岸堆积，因而在整个台风过程中均处于增水的状态，也是最早出现最大增水的位置。

图8 各代表站不同水位过程

各计算点非线性增水约1.19～1.39 m 对于风暴潮增水的贡献介于37.5%～50.0 %之间。由天文潮水位曲线与非线性增水曲线可以看出，非线性增水极小值一般出现在高潮位阶段不利于增水，极大值一般出现在低潮位阶段有利于增水。计算结果与前人研究结果基本一致[7]，从动力学上的解释认为：耦合作用中风与水位效应的非线性项作用影响的强弱与潮汐、纯风暴潮和水深都有关系，其中水深的影响最大，非线性项的作用影响与水深成反比；潮汐高潮削弱了风应力作用，从而减弱了风暴潮，低潮时相反[8]。珠江河口天文潮与风耦合作用较为复杂，与理论推导的结果有所差异，非线性增水的极值与高低潮位存在一定的相位差。

4 结语

1) 本文基于SELFE模式和圆形台风风场，考虑了径流、天文潮、台风的耦合作用，建立了粤港澳大湾区风暴潮数值模型。采用实测水文气象资料对“山竹”进行模拟验证，计算最大水位相对误差控制在25%以内，模型计算结果较好地反映了风暴潮增水的过程，可用于粤港澳大湾区风暴潮数值模拟研究。

2) 基于本文建立的风暴潮数值模型，以1822号台风“山竹”为例，讨论了非线性增水特征，结果表明非线性增水对风暴潮增水起到重要的作用，非线性增水极小值出现在高潮位阶段不利于风暴潮增水，极大值出现在低潮位阶段有利于风暴潮增水。