基于二维浅水方程和SWAN模型的风暴潮模拟研究

于英霞，张秀姬，宋怀辉

(1.河南科技大学 土木工程学院，河南 洛阳 471023；2.洛阳理工学院 土木工程学院，河南 洛阳 471023)

0 引言

随着全球气候变化以及海洋开发利用的进一步深入，风暴潮灾害的频发程度和危害程度也将呈继续上升趋势[1]。目前，风暴潮灾害已位居中国海洋灾害首位[2]。上海由于其特殊的地理位置，常遭受台风引起的强风暴潮侵袭。因此，建立一个高精度的风暴潮数值预报模型，对于其海岸保护系统的设计、风险评估和应急疏散规划都具有非常重要的意义。

目前，运用最广泛的风暴潮模式主要有英国的海模型(sea model)[3]、美国的飓风暴潮模式(sea,lake and overland surges from hurricanes, SLOSH)[4-5]、有限体积海岸海洋模型(finite-volume coastal ocean model, FVCOM)[6-7]和有限元海洋模式(advanced circulation model,ADCIRC)[8]。中国从20世纪80年代开始，利用引进的风暴潮模型进行大量的数值模拟研究和海洋预报业务。文献[9]利用FVCOM对发生在厦门湾附近海域的“莫兰蒂”台风引起的风暴潮进行模拟，文献[10]利用ADCIRC对上海地区历史风暴潮进行大量数值模拟研究，虽然两者都对增水分布、风生流场等风暴潮特征进行了较深入研究，但未考虑波浪的作用。波浪与风暴潮共存于同一水体，相互作用又相互影响，在具有宽阔大陆架的区域，波浪可使水位增加5%～15%，在具有狭窄大陆架或陡坡的区域，水位可增加35%～40%[11]。文献[12]利用ADCIRC和浅水波浪数值模型(simulating waves nearshore,SWAN)建立耦合模型,对珠江口风暴潮、近岸浪开展数值模拟与预报，文献[13]利用FVCOM和SWAN建立耦合模型,对汕尾沿岸海域影响较大的两次台风风暴潮过程进行模拟，分析浪流耦合作用对汕尾风暴潮增水的影响，虽然模拟精度都有所提高，但受制于引进的风暴潮模型，未能充分考虑波浪与风暴潮的相互作用。近年来，国内学者自主建立的一些风暴潮模型也很少考虑波浪的影响作用[14-16]。另一方面，研究表明不同地区波浪对风暴潮的影响有很大差异[17-19]，而上海沿海地区地处长江口和杭州湾交汇处，风、浪、潮的相互作用复杂，但相关研究较少[20]。

本文基于二维浅水方程建立风暴潮模型，利用高分辨率的三角形非结构网格与SWAN建立耦合模式。该模式充分考虑波浪的辐射应力与表面风应力作用，同时考虑风暴潮过程对波浪的影响。通过对台风温妮(Winnie)期间上海沿海地区复杂的风暴潮增水过程进行模拟分析，讨论了波浪对该地区风暴潮数值模拟的影响，对提高区域风暴潮预报精度以及潮流-波浪耦合研究均具有一定的实际应用价值。

1 台风Winnie简况

1997年第11号台风Winnie于8月10日8时在关岛东北偏东洋面上生成，即北纬15.4°、东经153.8°，当时中心气压995 hPa，近中心最大风力为8级(20 m/s)，朝西北偏西方向移动，14时发展成强热带风暴。11日8时到达北纬16.3°、东经150.1°时发展为台风，当时中心气压975 hPa，近中心最大风力12级(风速33 m/s)，并继续朝西北偏西方向移动。12日20时到达北纬18.1°、东经144.8°时达到鼎盛时期，当时中心气压920 hPa，近中心最大风速达60 m/s，并继续向西北偏西方向移动，向浙江中南部靠近，18日晚21:30在浙江温岭登陆，登陆时中心气压955 hPa，近中心最大风速为40 m/s，台风登陆后中心气压上升，近中心风力减弱，于19日8时降为强热带风暴，14时降为热带风暴，后经浙江北部、安徽、山东、辽宁省于21日2时减弱为低压。

8月18日至8月20日上海地区遭受此台风影响严重。台风Winnie(199711号)影响上海期间，适逢农历七月半天文大潮(一年中天文潮最高时段)、风潮遭遇，沿杭州湾、长江口和黄浦江干流各水文站均出现有记录以来的最高潮位。杭州湾沿岸的各站比原历史记录最高潮位抬高42～64 cm，其中金山嘴站达6.57 m；沿长江口各站抬高13～36 cm，其中外高桥站达5.99 m；沿黄浦江各站抬高24～50 cm，其中吴淞站达5.99 m，黄浦公园站为5.72 m，吴泾站为4.82 m，米市渡站为4.27 m。

2 风暴潮耦合模式的建立

2.1 风暴潮耦合模式

风暴潮过程中潮流、风生流组成的复合流场的运动规律可以用守恒性的浅水方程来描述。本文应用守恒性高、区域适应性强、分辨率高的非结构网格有限体积法离散浅水方程，运用HLLC格式计算界面通量，利用动边界跟踪的模式来处理干、湿交替的边界问题，建立了二维风暴潮水动力数值预报模型[20]。

采用SWAN模拟风浪在向岸边传播过程中各项因素对风浪形态的影响。SWAN 能够系统地模拟波浪产生过程、耗散过程及传播过程，并且源代码公开，被广泛应用于近海岸、湖泊及河口等浅水地区的波浪研究。

图1 风暴潮和波浪耦合数值模式流程图

本文通过运用高分辨率非结构三角形网格将风暴潮模型与SWAN进行双向耦合，能够更好地反映风暴潮与波浪之间的非线性相互作用。耦合计算过程中，每个循环时间间隔内，风暴潮模型首先运行，它将上一个时间间隔中SWAN模型内插波浪参数所得到的表面风应力和辐射应力加入到包括天文潮在内的风暴潮过程的水动力学计算中，得到水位和流速。在这个时间间隔的最后一个时间步长中，风暴潮模型将该间隔内计算得到的风场、水位和流场数据作为计算条件代入到SWAN，然后SWAN利用这些数据计算有潮流影响的台风浪过程，并计算出新的表面风应力和辐射应力，从而传递给风暴潮模型用于下一个时间间隔潮流场的计算并继续反馈给SWAN，实现双向耦合的循环迭代。耦合模式的计算流程如图1所示。

2.2 网格及模式参数

本文研究范围为上海沿海地区，涵盖了小尺度的河道(黄浦江)，中尺度的长江口和杭州湾，以及大尺度的东海部分海域。由于研究区域尺度跨越很大，为提高预报精度，设定了一个相对较小但合理的研究区域。其中，黄浦江和长江上游边界分别取至米市渡和徐六泾，杭州湾边界取至乍浦，外海的北边界取至吕泗以北40 km处，南边界至象山，东边界取至吴淞口以北200 km处。模型的计算区域约为120.9°E～123.6°E，29.5°N ～33.0°N，该区域水深分布图及长江口的局部放大图分别如图2a和图2b所示。

(a) 计算区域

图3 计算区域三角形非结构网格及站点分布

为了更好地适应上海沿海地区河网密布及地形复杂的特点，本文采用三角形非结构网格进行研究区域网格划分。耦合模式的网格图和站点分布见图3，共包含3 162个网格结点，5 634个网格单元。网格分辨率在外海靠近开边界处约为8 km，随着向近岸地区的靠近，网格逐渐加密，到长江口的浅水区域只有几百米，而黄浦江仅为几十米。

模型的初始条件为水位和流速。它们对外界动力的响应十分迅速，因此假设模型计算开始时海面是静止的，即初始时刻所有单元的流速均为0，水位取平均海平面高程。外海开边界利用全球天文潮位预报系统TPXO7.0中的8个主要分潮(M2、S2、K2、N2、K1、O1、P1和Q1)的调和常数来给定。对于黄浦江和长江口的上游边界，分别选用米市渡站和徐六泾站的实测数据；对于杭州湾的开边界，选用乍浦站的实测数据。

3 风暴潮的模拟与分析

利用建立的风暴潮和波浪的耦合模式，对台风Winnie(199711号)在上海近海海域引起的风暴潮增水情况进行了模拟。为了分析风暴潮数值模拟中波浪的作用，分别采用不考虑波浪对风暴潮影响的纯风暴潮模式(记为“非耦合模式”)、考虑波浪和风暴潮相互作用的耦合模式(记为“耦合模式”)进行计算，并将两种模式的计算结果与实测值进行对比分析。

3.1 潮位的模拟与验证

为了分析台风Winnie影响下的潮水位和流速变化过程，计算了1997年8月17日2时至8月20日23时的风暴潮过程。为了获取准确的初始条件，模型启动24 h后才开始预报。

在风暴潮对上海影响严重的69 h模拟中，计算出风暴潮和波浪耦合作用下的潮位过程线，与观测潮位过程的相位一致，水位增幅基本接近。横沙和吴淞两个测站的潮位对比结果分别见图4a和图4b。最大风暴潮增水出现在8月19日0点，计算值与实测值吻合较好，说明耦合模式在风暴潮增水计算方面具有较高精度。

(a) 横沙

3.2 潮流场的变化特征分析

非耦合模式和耦合模式计算所得台风Winnie最高潮位时刻长江口潮流场分别见图5a和图5b，杭州湾潮流场分别见图6a和图6b。非耦合模式下长江口与杭州湾交汇处产生了3.16 m/s的最大流速，耦合模式下该区域最大流速为3.28 m/s，比未考虑波浪作用的情况增速0.12 m/s，且流速超过3 m/s的区域略微增大。同时，长江口入海处流速分布也发生一些改变，北槽中西侧和佘山附近的流速都增大了0.2 m/s。非耦合模式下杭州湾开边界处有高达3 m/s往东北方向的流速，且区域较大。而耦合模式下该处最大流速仅有2.2 m/s。附近区域流速在迎面而来的潮流影响下迅速转向，且速度不到1 m/s。但上部区域在波浪的影响下流速由不到0.8 m/s增速至1.2 m/s左右。流速等值线的分布也发生了较大改变，流速大于2.0 m/s的区域虽然减少，但1.0～1.4 m/s流速区域增大很多。

(a) 非耦合模式

(a) 非耦合模式

从流场等值线的分布与流速极值发生的位置和大小来看，波浪作用在长江口区域对流场的影响较小，而在杭州湾其作用却十分显著。最高潮位时刻的风场是由东往西的方向，表面风应力随风拖曳力系数的增大而增强，致使流向长江口和杭州湾的潮流流速增大。然而,长江口水深较浅，波浪在向长江口传播时，水深变浅，再加上床面摩阻力的影响，在浅水中产生变形和破碎，波浪破碎以后，造成波高的减小，波浪作用也相应减小，所以该区域潮流场大部分区域变化较小。而杭州湾的水深较长江口相对较大，岸线较平缓，波浪在向岸传播过程中，引起波高在传播过程中的空间变化较长江口小，波浪的空间变化也相对较小，因此传播至杭州湾的波浪能量损失较小，在其影响下流入杭州湾的潮流流速增大。但最高潮位时刻该潮流正巧遇上杭州湾的退潮，这就造成杭州湾开边界处退潮的最大流速大大降低。

3.3 最高潮位变化特征分析

非耦合模式和耦合模式计算所得最高潮位时刻长江口潮位分布分别见图7a和图7b，杭州湾潮位分布分别见图8a和图8b。最高潮位时刻长江口潮位升高明显，北支支口涌入的潮量在其最狭窄段的下游形成6 m以上的高潮位，且5 m以上的高潮位范围也从北上支口漫延至横沙。杭州湾潮位变化不大，仅4.4～5.2 m的潮位范围有所缩小。这是由于耦合模式潮流场的影响，使得潮水在北支支口大量汇集，而杭州湾的潮水则往长江口退去。在受波浪影响的表面风应力作用下，从外向杭州湾推进的潮流增强，使得杭州湾高潮位等值线向开边界处推移。

(a) 非耦合模式

(a) 非耦合模式

选取吴淞和横沙两个潮位站作为验证点，将两种模式的计算结果与实测值进行对比，结果分别如图9a和图9b所示。与非耦合模式相比，耦合模式两个站点的计算潮位都与实测潮位更吻合，特别是最高潮位，两个站点的最高潮位计算值都更接近实测值。

(a) 横沙

吴淞、横沙、北槽中和长兴4个潮位站两种模式最高潮位误差对比结果见表1。由表1可知：非耦合模式4个站点最高潮位的相对误差为8.8%～10.1%，耦合模式预报精度明显有提高，相对误差为6.2%～7.9%，平均相对误差为7.2%，较非耦合模式提高2.3%。这是因为波浪造成的水面小尺度的波动改变了水面的表观粗糙度，同时波浪破碎会消耗大量的波能，辐射应力也随之而变，所以耦合模式能更好表现台风对大尺度潮流场的驱动力，从而改善了潮位增水。

表1 耦合模式和非耦合模式最高潮位误差对比

4 结论

(1)构建的风暴潮与波浪耦合模式，充分考虑了风暴潮与波浪的相互作用，并能有效模拟出台风Winnie引起的风暴潮潮位变化过程，经验证4个站点最高潮位平均相对误差仅7.2%，相对于未考虑波浪作用的纯风暴潮模型模拟精度提高了2.3%。

(2)波浪对上海沿海地区风暴潮潮流场作用受地形影响较大，最高潮位时刻，长江口潮流场受波浪影响较小，而杭州湾则影响较大。

(3)台风Winnie影响期间，波浪对上海地区不同水域风暴潮最高潮位影响也有区别，由于波浪和天文潮的综合作用，潮水向长江口汇集，造成长江口潮位明显增加。