随着我国交通路网从内陆向岛屿延伸，近岸海域建设了越来越多跨海桥梁，例如杭州湾大桥、东海大桥、港珠澳大桥等。这些桥梁的建设为实现我国东部沿海，以及“21世纪海上丝绸之路”沿线经济区域的互联互通起着举足轻重的作用。不同于内陆区域，这些桥梁所处的近海环境复杂，常常遭受海洋动力灾害的侵袭。强台风作用下，风、波浪、风暴潮形成彼此关联的极端环境荷载，严重威胁桥梁结构安全[1-2]。飓风Ivan、Katrina和Ike形成的极端波浪和风暴潮先后造成美国多座跨海桥梁被毁[3-4]。2018年台风“山竹”也曾让万千国人为港珠澳大桥捏了一把汗。
基于历史观测数据的概率分析是研究风、波浪、风暴潮荷载参数的常用手段[5-6]。由于设备有限、数据保密等原因，很多近岸海域无法获取到充足的风、浪、风暴潮实测数据。因此，掌握强台风作用下近岸海域风、浪、风暴潮耦合的数值模拟技术，不仅可以通过数值分析手段获取任意海域在台风下的风、浪、风暴潮数据，更可为桥梁结构设计、施工及运营管理提供重要的环境荷载参数。第三代近岸风浪模式和水动力模式的提出，实现了对台风作用下风暴潮和波浪场等发展过程的数值模拟[7]。为提高波浪和风暴潮数值预报水平，第三代风浪模式作为最常用的波浪模式，常被耦合到不同的海洋环流模型中[8]。Dietrich等[9]开发了采用非结构性网格进行飓风期间波浪-风暴潮模拟的数值计算模式，对Katrina和Rita台风引发的风暴潮进行了后报。众多研究表明[9-10]，采用海洋-波浪双向耦合模型预报风浪、风暴潮能达到较高的精度，可用于研究台风下的波浪和风暴潮环境参数。本文基于开源的SWAN+ADCIRC波流耦合模式，在经验风场的驱动下，建立了台风作用下波浪-风暴潮耦合数值模型，模拟强台风“玛莉亚”下台湾海峡及近岸海域的波浪-风暴潮耦合发展过程，重点研究了台风下波浪-风暴潮耦合模拟的驱动风场的构建方法和近岸海浪、风暴潮的分布特点，基于实测数据验证了模拟结果的可靠性。
1 波浪-风暴潮耦合数值模拟方法
第三代近岸风浪模式SWAN和水动力模式ADCIRC是开展波浪-风暴潮耦合数值模拟的常见工具。SWAN模式采用引入源项的动谱平衡建立控制方程来模拟风浪生成及在近岸的传播过程。其中源项包括风浪成长；三波和四波相互作用；白帽耗散、底部摩擦、水深变浅引起的波浪破碎、植被耗散、风浪遇到水下障碍物的绕射和反射、波浪绕射等过程。ADCIRC模式是一种连续伽辽金、有限元、浅水模型，可用于求解台风影响下的水位。因为二者皆为开源软件，在国内外相关设计、研究工作中得到了广泛应用。
1.1 SWAN风浪模式
SWAN模式的控制方程[11]为：
式中：N为波作用密度谱；σ为相对波频；θ为波向角；cx和cy为波浪传播速度的x和y向分量；cσ和cθ为σ、θ空间的波浪传播速度；S为谱密度表示的源项。
在本模式参数计算中，频率谱范围取0.03 Hz ~1.42 Hz；方向谱空间分为36段，即10°的分辨率；风浪成长线性增长部分采用Cavaleri与Malanotte-Rizzoli模式；指数增长部分采用Komen模式；同时考虑白帽耗散、三波和四波相互作用、波浪破碎、底部摩擦过程。其中：底部摩擦采用JONSWAP模式，摩擦系数取0.067；其余耗散项均取默认值。
1.2 ADCIRC水动力模式
ADCIRC模式根据广义波动连续性方程和动量守恒方程求解水位和流速，并在计算中考虑风、浪、潮汐和径流的影响。其控制方程[9]为：
式中：ζ为从海平面起算的水位高度；U、V为两个方向的平均流速；H为总水深；ps为水面大气压强；ρ0为水的密度；g为重力加速度；η为牛顿平衡潮汐势能；α为地球弹性影响因子；τbx、τby为底部切应力项；τsx、τsy为风应力项；τsx、τsy为波浪辐射应力；Mx、My为侧应力梯度项；Dx、Dy为扩散项在笛卡尔坐标下x、y方向的分量。
1.3 SWAN+ADCIRC耦合模式
SWAN+ADCIRC耦合模式是将SWAN模式和ADCIRC模式在相同的非结构化网格下进行耦合计算[9]。当考虑台风风场和潮汐的共同驱动时，首先由ADCIRC计算水位和潮流数据，并将计算得到的数据实时传递给SWAN；SWAN根据上述数据、驱动风场和外海波浪场，通过式(5)~式(7)所示谱方程计算得到波浪辐射应力张量，进而得到波浪信息，并将波浪辐射应力传递给ADCIRC；ADCIRC根据式(8)~式(9)计算波浪辐射应力，并代入下一个时间步的台风风场和潮汐信息开始下一步的水位和潮流计算。
式中：n为群速度与相速度的比值；E为能量谱密度。
上述耦合过程可以充分考虑波浪对风暴潮的影响以及风暴潮与潮汐之间的非线性关系。本文建立的数值模型主要考虑台风风场和天文潮汐驱动，同时ADCIRC计算步长设置为2 s，SWAN的计算步长为600 s。
2 台风波浪-风暴潮耦合数值模型
2.1 台风资料
以图1所示2018年8号强台风“玛莉亚”下的台湾海峡及近岸海域作为研究对象。玛莉亚台风于北京时间2018年7月4日20时在西太平洋生成，7月11日9时以强台风等级在我国福建沿海登陆，登陆时最大持续风速为42 m/s，中心最低气压为960 hPa。本文模拟了2018年7月7日0时~2018年7月12日0时，共5天台风“玛莉亚”下上述海域的波浪-风暴潮发展过程(本文时间基准都采用世界时间，世界世界=北京时间-8 h)。计算所需的台风轨迹、中心气压、及台风最大持续风速等数据来自中国台风网(http://www.typhoon.org.cn/)。为验证数值模拟的准确性，选取福建省海洋预报台(http://www.fjmf.gov.cn/)的风、浪实测数据作为数值模拟对照依据。图1给出了玛莉亚台风的轨迹以及风、浪浮标的位置，其中#1浮标位于近岸海域，#2浮标位于远海海域。
图1 台风轨迹及浮标位置图
Fig.1 Typhoon track and the position of buoy
2.2 计算海域网格划分
选取图2(a)所示位于112.6°E~131.6°E、14.3°N~31.9°N之间的区域作为计算海域。根据全球地形数据集ETOPO1 (https://maps.ngdc.noaa.gov)，建立地形模型。考虑到近岸水深对波浪传播影响较大，本研究重点关注的平潭海峡公铁两用大桥所在的平潭海峡附近海域采用国家海道测量局提供的海图进行岸线和水深的精细化修正。采用网格划分软件SMS对上述区域进行非结构化网格离散，如图2所示。其中最大网格尺寸为40 km，而最小网格尺寸为60 m。计算模型节点总数为69822个，单元总数134168个。
图2 模型区域及网格示意图
Fig.2 Model domain with unstructured triangular grids
2.3 驱动风场模型
风场是波浪和风暴潮形成的主要驱动力。因此，风场的准确描述，是开展波浪-风暴潮数值模拟的前提。气旋形成的风场主要由两部分构成，一部分为台风移动引起的移行风场，一部分为气旋平衡得到的梯度风场[12]，即：
式中：vm为台风移动风场；vg为台风梯度风场。本文vm由宫崎正卫模型[13]确定，计算公式如下：
式中，vmc为台风中心移动风速矢量。梯度风场vg可通过求解动量方程得到。当忽略海表面的粗糙度，考虑气压和科氏力平衡得到的vg为：
式中：ρa为空气密度(取1.29kg/m3)；f为科氏力系数；r为计算点到台风中心的距离；P为计算位置处气压。目前，国内常用的气压模型主要有：Takahashi模型、Fujita模型、Jelesnianski模型、Holland模型等。由于Holland模型，在最大风速半径Rmax的基础上，引入了形状系数B，具有更普遍的适用性，所以本文选择Holland模型确定梯度风场。Holland模型[14]表达式为：
式中：P0为台风中心气压；Pn为台风外围气压(取1010 hPa)；ΔP=Pn-P0为外围气压与中心气压差。
相比其他气压模式，Holland模型需要确定2个关键参数：最大风速半径Rmax和形状系数B。它们与台风中心气压，台风位置等有关[15]。基于实测数据，前人提出了许多最大风速半径、形状系数关于中心气压、地理位置等的经验公式。为了获得较准确的模型参数，Holland模型最大风速半径由比较适用于西北太平洋不同的最大风速半径经验公式得到，形状系数由中国台风网提供的最大持续风速，通过式(14)[14]计算得到。
式中：vgmax为最大梯度风速，近似取最大持续风速；e为自然底数。通过查阅文献，适用于西北太平洋的最大风速半径经验公式主要有：
李瑞龙[16]：
卢安平[17]：
Zhou等[18]：
将式(11)、式(12)代入式(10)得到台风风场经验模型为：
式中：Vx、Vy为笛卡尔坐标系内模型风场风速在x方向(东)和y(北)方向的分量；Vm为移行风场风速；Vg为梯度风场风速；θ为计算点与台风中心的连线与x轴的夹角(逆时针)；ϕ为台风移动方向与x轴的夹角(逆时针)；β为流入角，取20°，C1、C2为修正系数，本文取1.0和0.71，C2取0.71是由大气边界层理论，梯度风速和10 m风速之间转换关系确定[19]，并得到了下文实测数据的验证。
有研究[12]表明，台风参数模型在台风中心附近具有较高的精度，但随着距台风中心距离的增加，精度下降；相比台风参数模型，风场再分析产品在远离台风中心处，具有较高精度，在台风中心附近，精度较差。因此，为了获得更精确的风场，本文将经验风场与ERA-Interim风场叠加，获得最终波浪-风暴潮耦合模型的驱动风场，叠加过程如图3所示。ERA-Interim风场是欧洲中期天气预报中心(ECMWF：https://www.ecmwf.int/)的海表面风场再分析数据。叠加风场的原则为：在台风中心以经验风场为主，在台风外围以ERA-Interim风场为主，叠加方式如下[20]：
式中：vsup为叠加风场风速矢量；vERA为ERA-Interim风场风速矢量；n为常数，一般取9或10，本文取10。
最大风速半径分别采用式(15)~式(17)确定的模型风场叠加ERA-Interim风场获得图1所示的风速浮标的模拟风速与实测风速对比情况如图4所示。为定量分析数值模拟效果，时间段7月10日~7月12日内的模拟结果采用相关系数R和均方根误差RMSE对模拟风速和实测风速的一致性进行评价：
图3 叠加风场的计算过程示意图(世界时间2018-07-11 00:00)
Fig.3 Modeling process of driven wind field at 00:00 on July 11, 2018 (UTC)
式中：Mi为实测风速；Si为模拟风速；
为实测风速均值；
为模拟风速均值；n为实测风速的总样品量。
由图4知，最大风速半径采用李瑞龙[16]、卢安平[17]和Zhou等[18]的经验公式计算的玛莉亚台风期间的模拟风速与实测风速在#1和#2浮标吻合都较好。这主要原因是各最大风速半径经验公式在台风中心气压较小时，差距较小。由#2浮标风速验证知：当台风靠近#2浮标时，数值风速低估了实际风速，说明此时经验公式计算的最大风速半径偏小；当台风远离#2浮标后，数值风速高估了实际风速，说明此时经验公式计算的最大风速半径偏大，这间接地表明了台风结构的不对称性，台风中心前侧与后侧最大风速半径可采用不同的经验公式描述。由于#1浮标距离台风路径较远，ERA-Interim风场起到了很好的修正作用，所有在台风靠近1#浮标前数值风速与实测风速吻合较好。结合#1、#2浮标风速时程验证和相关系数、均方根误差分析，可以看出，最大风速半径采用Zhou等[18]的经验公式得到的风场与实测吻合情况最好，实测风速与模拟风速在#1浮标相关系数和均方根误差分别为0.92、2.21 m/s；在#2浮标分别为0.93、3.20 m/s，所以下文玛莉亚台风风场的最大风速半径由Zhou 等[18]的经验公式确定。
图4 实测风速与模拟风速比较图
Fig.4 The comparison of calculated wind field with measurement
2.4 天文潮驱动模型
天文潮造成水位的涨落，对近岸海域波浪的模拟有重要影响。本文选用美国TPXO模式[21]对天文潮进行预测，具体为：采用TPXO9_atlas模式提供的M2、S2、N2、K2、K1、O1、P1、Q1分潮的调和常数来预测图2网格边界上的水位，作为波浪-风暴潮耦合模式的天文潮汐驱动。
图5给出了采用TPXO模式提供边界水位的ADCIRC模式预测的#1浮标水位与TPXO模式及国家海洋信息中心预测的#1浮标水位的对比情况。由图5知，采用ADCIRC模式预测的#1浮标水位与TPXO模式预测水位几乎重合，与国家海洋信息中心预测水位吻合较好，整体趋势一致。
图5 天文潮验证
Fig.5 Variation of astronomical tide
3 数值模拟结果与讨论
本节采用SWAN+ADCIRC波流耦合模式，在叠加风场和天文潮汐的驱动下，模拟了强台风“玛莉亚”下近岸海域的波浪和风暴潮发展过程。
3.1 波浪计算结果与讨论
图6对比了模拟与实测得到的有效波高，并采用相关系数R和均方根误差RMSE对模拟效果进行了评价。由图6知，#1、#2浮标实测有效波高与模拟有效波高的相关系数分别为0.80、0.97，均方根误差分别为0.50 m、1.12 m。可见本数值模拟的#1、#2浮标有效波高与实测有效波高相关性较高，均方根误差较小，说明数值模拟结果与实测吻合较好。结合图4风速验证，可以看出，有效波高发展规律与风速发展规律相似，#2浮标有效波高存在与风速类似的双峰结构，同时近海数值模拟有效波高和实测有效波高，在小风速时差别小，大风速时差别大；而远海波高在风速小时差别大，风速大时差别小。
图7分别给出了10日18时、11日0时、11日6时，算例海域的有效波高分布。由图7可知，当台风向近岸靠近时，有效波高增加，由于平潭海域位于台风轨迹左侧，风速与海浪的传播方向相反，同时受近岸水深浅化的影响，平潭海域有效波高在11日达到最大后衰减，并在靠近台风轨迹位置波高衰减较缓。同时，岛屿对波浪的传播有明显的遮蔽效应，所以模拟近岸海域的台风浪时，应注意近岸岛屿水深和岸线的处理，可能会显著影响近岸海域波浪模拟精度。
图6 实测有效波高与模拟有效比较图
Fig.6 The comparison of calculated significant wave height with measurement
图7 不同时刻近岸海浪分布图
Fig.7 Wind field and significant wave height distribution of typhoon Maria at different times
3.2 风暴潮计算结果与讨论
台风期间，向岸的大风驱动海水向岸流动，受岸线阻挡，使得海平面上升。因此通过数值模拟研究风暴潮，对于近岸结构灾害预警十分重要。图8给出了#1、#2浮标数值模拟得到的风暴增水情况。
图8 不同浮标风暴潮增水随时间变化图
Fig.8 Time history of storm surge at different buoys
由图8可知，#1、#2浮标位置处的风暴潮发展规律相似，于10日0时开始出现明显的风暴增水，并在11日0时台风穿过#2号浮标，增水达到最大，此时#1号浮标增水为0.77 m，#2增水为0.93 m。从图8还可以发现在10日和11日之间存在一个增长较缓的时间段，这是由于#1、#2浮标距台风中心距离减小，但是台风强度减低所导致。
图9 不同时刻近岸风暴潮分布图
Fig.9 Wind field and surge distribution of typhoon Maria at different time
图9分别给出了不同时刻近岸风暴潮的分布。由图9可知，风暴潮是强风导致的海面异常升高现象，受水深的影响不显著，与台风强度及距离台风中心的距离关系密切，台风轨迹左侧的风暴潮在靠近台风中心处和风向朝向岸线的岸线附近增水较大。由图9还可以明确，台风过程中风暴潮与风速、风向关系较为复杂，同一时刻的风速、风向与风暴潮没有直接对应关系。
4 结论
文中采用Holland模型风场叠加宫崎正卫模型风场和ERA-Interim风场模拟台风风场，通过SWAN+ADCIRC波流耦合模式，模拟了强台风“玛莉亚”下台湾海峡及近岸海域的波浪-风暴潮发展过程，结合风、浪实测资料对数值模拟结果进行了验证，探讨了李瑞龙[16]、卢安平[17]和Zhou等[18]提出的最大风速半径经验公式对台风“玛莉亚”风场的适用性，并研究了近岸海浪、风暴潮的分布规律，文章得出的主要结论如下：
(1) 最大风速半径采用Zhou等[18]提出的经验公式确定Holland模型风场，最终得到的叠加风场与台风“玛莉亚”实测风场的吻合度最好。
(2) 叠加风场和TPXO潮汐模式驱动SWAN+ADCIRC耦合模式得到的有效波高与实测值吻合较好，该数值模拟方法可以较好地模拟近岸海域波浪的生成与发展过程。
(3) 岛屿对波浪的传播有明显的遮蔽效应，模拟近岸海域的台风浪时，应注意近岸岛屿水深和岸线的处理，可能对近岸海域波浪预测带来较大的影响。
(4) 台风期间，强风会导致海面的异常升高形成风暴潮；风暴潮受水深的影响不显著，与台风强度及距离台风中心的距离关系密切，台风轨迹左侧的风暴潮在靠近台风中心处和风向朝向岸线的岸线附近增水较大。
综上，本文模型模拟精度较好，能有效地用于确定台风作用下近岸海域的极端风、浪、风暴潮环境荷载，可为开展风-浪-风暴潮环境荷载概率分析，提供有效手段。但是实际台风结构往往是不对称的，采用对称的Holland模型风场，对不对称的台风结构进行描述，本身就是一种近似。同时，本文限于篇幅仅针对一次台风过程进行了数值模拟，本文结论是否能够推广到其他台风和其他海域还有待进一步验证。
[1]Ti Z L, Wei K, Qin S Q, et al.Assessment of random wave pressure on the construction cofferdam for sea-crossing bridges under tropical cyclone [J].Ocean Engineering, 2018, 160: 335―345.
[2]张家瑞, 魏凯, 秦顺全.基于贝叶斯更新的深水桥墩波浪动力响应概率模型[J].工程力学, 2018, 35(8):138―171.Zhang Jiarui, Wei Kai, Qin Shunquan.An Bayesian updating based probabilistic model for the dynamic response of deep-water bridge piers under wave loading[J].Engineering Mechanics, 2018, 35(8): 138―171.(in Chinese)
[3]Xu G J, Chen Q, Chen J H.Prediction of solitary wave forces on coastal bridge decks using artificial neural networks [J].Journal of Bridge Engineering, 2018, 23(5):04018023-1―04018023-13
[4]Qeshta I M I, Hashemi M J, Gravina R, et al.Review of resilience assessment of coastal bridges to extreme wave-induced loads [J].Engineering Structures, 2019,185: 332―352.
[5]姚博, 全涌, 顾明, 等.混合气候地区极值风速分析方法研究[J].工程力学, 2018, 35(5): 86―92.Yao Bo, Quan Yong, Gu Ming, et al.Study on the analysis method of extreme wind speed in mixed climate areas [J].Engineering Mechanics, 2018, 35(5): 86―92.(in Chinese)
[6]孙富学, 许向楠, 史文海, 等.温州滨海平坦地貌近地台风特性实测研究[J].工程力学, 2018, 35(9): 73-80.Shu Fuxue, Xu Xiangnan, Shi Wenhai, et al.Field measurements of typhoon characteristics near ground in Wenzhou coastal flat terrian [J].Engineering Mechanics,2018, 35(9): 73―80.(in Chinese)
[7]Ti Z L, Wei K, Qin S Q, et al.Numerical simulation of wave conditions in nearshore island area for sea-crossing bridge using spectral wave model [J].Advances in Structural Engineering, 2018, 21(5): 756―768.
[8]Chen T, Zhang Q, Wu Y, et al.Development of a wave-current model through coupling of FVCOM and SWAN [J].Ocean Engineering, 2018, 164: 443―454.
[9]Dietrich J C, Tanaka S, Westerink J J, et al.Performance of the unstructured-mesh, SWAN+ADCIRC model in computing hurricane waves and surge [J].Journal of Scientific Computing, 2012, 52(2): 468―497.
[10]Wang Y, Mao X, Jiang W.Long-term hazard analysis of destructive storm surges using the ADCIRC-SWAN model: A case study of Bohai Sea, China [J].International Journal of Applied Earth Observation and Geoinformation, 2018, 73: 52―62.
[11]Booij N, Ris R C, Holthuijsen L H.A third-generation wave model for coastal regions: 1.Model description and vaildation [J].Journal of Geophysical Research: Oceans,1999, 104(C4): 7649―7666.
[12]Pan Y, Chen Y P, Li J X, et al.Improvement of wind field hindcasts for tropical cyclones [J].Water Science and Engineering, 2016, 9(1): 58―66.
[13]Miyazaki M, Ueno T, Unoki S.Theoretical investigations of typhoon surges along the Japanese coast [J].Oceanographic Magazine, 1962, 13(2) : 103―117.
[14]Holland G J.An analytic model of the wind and pressure profiles in hurricanes [J].Monthly Weather Review,1980, 108(8): 1212―1218.
[15]方根深, 赵林, 梁旭东, 等.基于强台风“黑格比”的台风工程模型场参数在中国南部沿海适用性研究[J].建筑结构学报, 2018, 39(2): 106―113.Fang Genshen, Zhao Lin, Liang Xudong, et al.Applicability analysis of typhoon field parameters in engineering model for south coastal region of China based on strong typhoon Hagupit 0814 [J].Journal of Building Structures, 2018, 39(2): 106―113.(in Chinese)
[16]李瑞龙.基于改进的台风关键参数的台风极值风速预测[D].哈尔滨: 哈尔滨工业大学, 2007.Li Ruilong.Prediction of typhoon extreme wind speeds based on improved typhoon key parameters [D].Harbin:Harbin Institute of Technology, 2007.(in Chinese)
[17]卢安平.登陆台风近地风压场实测重构及其对工程场地极值风速影响分析[D].上海: 同济大学, 2012.Lu Anping.Near ground air pressure of landing typhoon measured reconstruction and its impact on the extreme value wind velocity evaluations of engineering fields[D].Shanghai: Tongji University, 2012.(in Chinese)
[18]Zhou T Y, Tan Y, Chu A, et al.Integrated model for astronomic tide and storm surge induced by typhoon for Ningbo coast [C].Sapporo Japan: International Society of Offshore and Polar Engineers, 2018: 1124―1129.
[19]Wei K, Arwade S R, Myers A T, et al.Effect of wind and wave directionality on the structural performance of non-operational offshore wind turbines supported by jackets during hurricanes [J].Wind Energy, 2017, 20:289―303.
[20]李健, 侯一筠, 孙瑞.台风模型风场建立及其模式验证[J].海洋科学, 2013, 37(11): 95―102.Li Jian, Hou Yijun, Sun Rui.Surge model caused by 0814 typhoon and mold wind field established [J].Marine Sciences, 2013, 37(11): 95―102.(in Chinese)
[21]Egbert G D, Erofeeva S Y.Efficient inverse modeling of barotropic ocean tides [J].Journal of Atmospheric &Oceanic Technology, 2002, 19(2): 183―204.