浙江地处中国东南沿海, 经济发达、人口稠密, 频繁遭受台风影响。台风会引发巨浪, 并伴随着风暴潮, 不仅直接威胁着海塘的结构安全, 而且因漫堤诱发的海岸洪水还严重影响着低地城镇的防洪安全。因此, 对台风诱发的海浪现象进行准确的后报模拟变得尤为重要, 这不仅有助于对过往灾害事件影响的深入理解, 还可为海塘的工程设计以及未来的防灾减灾管理决策提供科学依据。

数值模拟是研究台风浪的重要手段。目前, WAM (WAve Model)^[1]、WWⅢ(WaveWatch-Ⅲ)^[2]和SWAN^[3]是3种常用的海浪数值模拟模型。其中, SWAN考虑了各种浅水能量源项, 与WAM和WWⅢ相比, 更适用于较小尺度的区域性海洋和近岸浅水区域。已有研究利用SWAN模型, 针对南海^[4]、东海^[5]、黄海^[6]、渤海^[7]等区域性海洋的台风浪进行了模拟, 取得了良好的应用效果。针对近岸水域的研究, 如海南岛周边^[8]、浙江海域^[9]、长江口海域^[10]等, 采用高分辨率的非结构化网格, 以实现对岸线变化特征的精细刻画。

合理可靠的台风场模型是台风浪精确模拟的关键。目前常用的台风场模型有参数化模型和数值模型, 其中参数化模型使用简单、计算高效, 且能合理反映大风区风场的空间变化, 因此在大风、风暴潮和海浪的长期分布、灾害评估、海岸工程设计等相关问题的研究中有着广泛应用^[11-13]。Holland^[14]基于Schloemer的固定剖面气压场模式, 引入剖面参数B, 提出了可变剖面的Holland参数化风场模型。在应用上该模型包含了若干复杂变化的参数, 包括风速剖面参数B、最大风速半径R_max等, 这些参数对于台风强度与尺度的描述至关重要。针对关键参数B和R_max已有多种经验公式, 在建立区域性台风模型时, 掌握台风场的计算结果如何随模型参数变化, 并对不同的参数方案进行比较研究十分必要, 这可以优化出最佳的模型配置, 进而提高台风浪的模拟效果。陈相宇等^[9]利用Holland参数化风场和欧洲中期天气预报中心的再分析数据ERA5构建混合风场, 通过5次台风事件验证了模型在浙江海域对台风浪模拟的合理性, 研究中未考虑不同风场参数对台风浪模拟的影响。马秀玲等^[15]则针对台风“威马逊”, 探讨了B和R_max的不同参数化表达式对南海北部海域台风浪的模拟效果。尽管研究中建议了2个参数的最佳组合方案, 但从台风浪模型应用的区域性差异和验证的广泛性角度考虑, 参数组合方案在浙江海域的应用效果仍有待进一步的研究。

因此, 本文针对浙江沿海台风浪的后报模拟问题, 探讨了Holland模型风速剖面参数B和最大风速半径R_max的不同参数表达式组合方案对模拟结果的影响, 目的是基于较为广泛的资料, 率定并验证这两个参数的最佳组合方案。文中研究的台风浪个例涉及2011年以来影响浙江海域的10场不同路径和不同强度的典型事件, 其中4次事件中有气象站风速和波浪浮标有效波高的观测数据, 这些数据将用于模型参数率定, 全部10次事件中均有卫星观测的大浪区附近的有效波高数据, 它们将用于模型的验证分析。

1 数值模型 1.1 海浪模型

SWAN模式是基于波作用量平衡方程由荷兰Delft理工大学开发的第三代海浪数值模型, 由于充分考虑了浅水作用项, 相较于WAM和WWIII更加适用于近岸海浪生成、传播和耗散的数值模拟, 其在球坐标系下的控制方程如下:

$ \frac{{\partial N}}{{\partial t}} + \frac{{\partial {C_\lambda }N}}{{\partial \lambda }} + \frac{{\partial {C_\varphi }N}}{{\partial \varphi }} + \frac{{\partial {C_\sigma }N}}{{\partial \sigma }} + \frac{{\partial {C_\theta }N}}{{\partial \theta }} = \frac{S}{\sigma } , $

(1)

式中, N是波作用量密度; C_λ和C_φ分别为作用量沿地理空间(经度方向λ和纬度方向φ)的传播速度; C_σ和C_θ分别为作用量沿波谱空间(频率空间σ和方向空间θ)的传播速度; S表示能量源项, 在浅水中可表达为

$ S = {S_{{\text{in}}}} + {S_{{\text{ds,w}}}} + {S_{{\text{ds,b}}}} + {S_{{\text{ds,br}}}} + {S_{{\text{nl4}}}} + {S_{{\text{nl3}}}}, $

(2)

式中, S_in为风能输入; S_{ds, w}为白浪破碎耗散; S_{ds, b}为底摩阻耗散; S_{ds, br}为水深诱导的波浪破碎耗散; S_nl4为四波相互作用; S_nl3为三波相互作用。

1.2 驱动风场

Holland^[14]在Schloemer指数型气压模型的基础上, 引入了一个可调整的剖面参数B, 给出了气压的径向分布表达式:

$ p(r) = {p_{\text{c}}} + ({p_{\text{n}}} - {p_{\text{c}}}){( - \frac{{{R_{\max }}}}{r})^B} , $

(3)

式中, p(r)为距台风中心r处的气压值; r为距台风中心的距离; p_c为台风中心气压; p_n为台风外围环境气压; R_max为最大风速半径。

基于气压场分布, 利用旋恒风原理, 得到梯度风的计算公式为:

$ \begin{aligned} V_{\mathrm{G}} & =\sqrt{\left(p_{\mathrm{n}}-p_{\mathrm{c}}\right) \frac{B}{\rho_{\mathrm{a}}}\left(\frac{R_{\max }}{r}\right)^B \exp \left[-\left(\frac{R_{\max }}{r}\right)^B\right]+\left(\frac{r f}{2}\right)^2} \\ & -\frac{r|f|}{2}, \end{aligned} $

(4)

式中, V_G为梯度风; ρ_a为空气密度; f是科氏力参数(f=2ωsin(ϕ), ω为地球自转角速度)。

通过叠加移行风场以解决台风在移动过程中产生的不对称现象, 台风移行风场采用宫崎正卫模式^[4]:

$ {V_{\text{T}}} = \exp ( - {\text{π}}r/50{\text{ }}000)\left[ {\begin{array}{*{20}{c}} {{V_x}} \\ {{V_y}} \end{array}} \right] , $

(5)

式中, V_T是移行风速, V_x和V_y为台风中心移动速度的正东分量和正北分量。

通过将梯度风场和移行风场叠加得到台风的模型风场:

$ {V_{\text{M}}} = {c_1}{V_{\text{G}}}\left[ {\begin{array}{*{20}{c}} { - \sin (\theta + \beta )} \\ {\cos (\theta + \beta )} \end{array}} \right] + {c_2}{V_{\text{T}}} , $

(6)

式中, V_M为叠加风场的风速; c₁和c₂是风速修正系数, 模拟时取c₁=0.90, c₂=0.80。

鉴于参数化台风模型对台风外围风场反映不足, 为形成更加可靠的驱动风场, 将其与再分析风场相融合。CFSV2 ^[16]是NOAA NCEP于2011年发布的全球尺度的再分析数据, 空间分辨率为0.2°, 时间分辨率为1 h, 在海浪数值模拟研究中有着广泛应用。Shi等^[17]利用CFSR风场建立了中国近海近40年的波浪数据库, 他们的研究显示CFSR对中国近海波浪展现出较好的效果。因此, 本文选取CFSV2作为背景风场, 通过权重系数e, 将Holland参数化风场与CFSV2再分析风场进行融合:

$ {V_{\text{C}}} = \left( {1 + e} \right){V_{\text{M}}} + e{V_{\text{E}}} , $

(7)

式中, V_C为融合风场; V_M为Holland模型风场; V_E为CFSV2背景风场; e是权重系数($ e={C}^{4}/（1+{C}^{4}） $), C是考虑台风影响范围的系数($ C = r/n{r_{\max }} $), 模拟时系数n取为9。

1.3 关键风场参数及组合方案

在参数化台风场模型中, 剖面参数B反映了风速空间分布的剖面形状, 最大风速半径R_max刻画了大风区的控制范围。在台风移动过程中两个参数随时间和空间不断变化, 其参数化表达直接影响着台风风场的计算效果。

关于风速剖面参数B, Holland基于澳大利亚海域台风数据, 建议了B的合理取值范围为1.0~2.5。Hubbert等^[18]提出了适用于澳大利亚海域的参数B计算式:

$ B = 1.5 + (980 - {P_{\text{c}}})/120. $

(8)

Willoughby等^[19]利用美国国家海洋和大气管理局(NOAA)飓风研究部(HRD)发布的飞行层探测资料, 对参数B进行了最小二乘法拟合, 给出的计算式如下:

$ \begin{aligned} B= & 1.0036+0.0173 V_{\max } \\ & -0.0312 \ln R_{\max }+0.0087 \varphi . \end{aligned} $

(9)

基于NOAA-HRD飞行层探测资料, Powell等^[20]进一步筛选了数据, 给出了参数B的表达式为

$ B = 1.881 - 0.005\;57{R_{\max }} - 0.012\;95\varphi . $

(10)

关于最大风速半径参数R_max, Willoughby等^[19]给出的计算式如下:

$ {R_{\max }} = 51.6\exp ( - 0.022\;3{V_{\max }} + 0.028\;1\varphi ) . $

(11)

Vickery等^[21]结合NOAA-HRD飞行层探测数据和H*Wind资料, 拟合给出了参数R_max的计算式为

$ {R_{\max }} = \exp (3.015 - 6.291 \times {10^{ - 5}}\Delta {p^2} + 0.033\;7\varphi ) . $

(12)

考虑到不同B和R_max的参数化组合对台风场的刻画存在较大差异, 所以选取上述3个风速剖面参数B表达式和2个最大风速半径参数R_max表达式, 二者交叉组合构成了不同的计算方案, 具体见表 1。

1.4 数值模型参数设置

模拟的空间范围为18°N~34°N, 114°E~134°E, 向北抵达长江口, 向南包括台湾岛, 向西包含琉球群岛。模型水深(图 1)在外海采用美国国家地理数据中心(U.SNGDC)发布的分辨率为1′×1′的ETOPO1全球地形数据, 在浙江沿海采用航保部海图水深进行细化。计算网格采用非结构化网格, 以合理刻画复杂岸线的变化。模型网格共计39 589个节点, 77 053个三角形单元, 由东南部西北太平洋向东北部中国东南沿海逐渐加密: 外海网格尺度为10 km, 近岸网格尺度为1 km, 到浙江沿岸网格尺度加密到500 m。模型计算采用冷启动的方式, 为消除冷启动对模拟结果的影响, 模型起算时间从台风进入研究区域开始提前3 d, 结束时间为气旋消失时间, 模拟的时间步长取为2 min。模型输出包括有效波高、平均波周期、平均波向、谱峰波向等, 时间分辨率取为1 h。在物理过程的考虑上, 风能输入和白浪损耗采用Komen方案, 海底摩擦损耗使用Collins公式, 考虑近岸三波相互作用和浅水波浪破碎的影响, 相关参数设置均为模型建议的默认值^[22]。

2 台风浪研究个例、观测资料与误差评估指标 2.1 研究个例

自2011年起浙江沿海发生多场热带气旋, 筛选出其中强度达到台风及以上且被Jason卫星观测到的10场台风作为台风浪研究个例, 具体包括: 1109号超强台风“梅花”、1211号强台风“海葵”、1323号强台风“菲特”、1509号超强台风“灿鸿”、1616号强台风“马勒头”、1618号超强台风“暹芭”、1718号台风“泰利”、1909号超强台风“利奇马”、1918号超强台风“米娜”和2106号强台风“烟花”。在这10场台风中: 超强台风有5场, 强台风有3场, 台风有1场, 按照强度不同将其路径分别绘制于图 2(a)—2(c)中, 可见有4场台风的移动路径为西北向, 西向、北向和转向型台风各有2场。

2.2 观测资料

在上述10场台风个例中, 收集到了4场台风影响下的气象站风速(A01—A05, 5站)和波浪浮标有效波高(B01—B04, 4站)的观测资料, 4场台风路径、5个气象站和4个浮标站的空间分布及其与台风路径的位置关系分别见图 3(a)和图 3(b), 其中: 气象站点A01—A03位于大洋, A04和A05位于沿岸; 浮标站B01—B04位于浙江外海深水域, 这些风/浪资料将用于台风参数的率定分析, 以确定1组最佳参数组合方案。上述10次台风活动均被Jason卫星(包括Jason2和Jason3)所捕捉, 具有高度计有效波高观测数据, 其轨道分布见图 3(c), 这些卫星高度计记录的有效波高数据将用于进一步验证最佳组合方案下台风浪模拟的空间分布效果。

2.3 误差评估指标

采用误差指标定量分析台风风场模型及其驱动海浪的模拟效果, 具体包括平均偏差(BIAS)、均方根偏差(RMSB)和相关系数(R), 表达式如下:

$ {\text{BIAS}} = \frac{1}{N}\sum\nolimits_{i = 1}^N {\left( {{S_i} - {O_i}} \right)} , $

(13)

$ {\text{RMSB}} = \sqrt {\frac{1}{N}\sum\nolimits_{i = 1}^N {{{\left( {{S_i} - {O_i}} \right)}^2}} } , $

(14)

$ R = \frac{{\sum\nolimits_{i = 1}^N {\left( {{S_i} - \bar S} \right)\left( {{S_i} - \bar O} \right)} }}{{\sqrt {\sum\nolimits_{i = 1}^N {{{\left( {{S_i} - \bar S} \right)}^2}{{\left( {{S_i} - \bar O} \right)}^2}} } }} , $

(15)

式中, N为样本量; S为模拟值; O为实测值; “–”为求平均运算符。

3 台风浪模拟结果分析 3.1 基于气象站实测风速的模型率定分析

首先, 从风速角度对台风场模型进行率定。针对1323号“菲特”、1509号“灿鸿”、1909号“利奇马”和1918号“米娜”4场台风, 使用不同的参数组合方案对风场进行模拟计算。为节省篇幅, 对于每场台风, 在大洋和沿海气象站中分别选取1个台风路径附近的气象站数据, 与模拟风速进行对比, 时间过程见图 4(a)~图 4(h)。可以看出: 6种参数组合方案下台风风场模型的模拟结果与气象站资料总体吻合良好, 但是在台风风速峰值附近的大风区, 各组次模拟结果存在着明显差异。R_max的不同造成了选用R2的3个组次(B1R2、B2R2、B3R2)对于大风区风速的模拟普遍高于选用R1的3个组次(B1R1、B2R1、B3R1)。选用相同R_max的组次时间过程曲线较为相似, 这在“威马逊”台风浪的模拟中(马秀玲等^[15])也有相似的结果, 这表明B参数对台风风速的模拟影响相对较小, B3R1方案下的台风风场模型对风速变化刻画得最好, 与各气象站数据均较为贴合, 适合浙江沿海台风风场模型的模拟计算。

为直观对比各参数组合方案下风场模拟效果的差异, 将所有气象站风速模拟值与实测值进行对比, 并绘制散点密度图如图 5。由此可见, 对于10 m/s以下风速, 各参数组合方案下的模拟值与实测值之间的差异相对较小, 但对10 m/s以上(尤其是20 m/s以上)大风速的刻画差异较大, 这说明了对台风场参数组合率定的必要性。较最大风速半径选用R1的3个组次(B1R1、B2R1、B3R1), 选用R2的3个组次(B1R2、B2R2和B3R2)出现了显著的对20 m/s以上大风速的高估, 且在大风区的散点分布更显分散。对比B1R1、B2R1、B3R1可以发现, B3R1在大风区的表现最优, 散点更加集中在y=x直线附近, 由于风速剖面参数B的组合不同, 相关性也高于B1R1和B2R1。

表 2对6个参数组合方案下风速模拟值与实测值之间的误差进行了计算。由此可见, 在大洋气象站(A01、A02、A03)的模拟相对优于近岸气象站(A04、A05)。其中, A01和A02站的模拟效果最佳, 各组次在A01有着最高的相关系数(R)和相对较低的均方根偏差(RMSB)。B3R1组次从6个组次中脱颖而出, 在各气象站风速模拟值的平均偏差、均方根偏差和相关性分别为0.7 m/s、3.7 m/s和0.85。

3.2 基于波浪浮标实测有效波高的模型率定分析

然后, 从有效波高模拟角度对台风浪模型进行率定。针对4场台风浪过程, 分别使用6种参数组合方案下的风场驱动SWAN模型, 对台风浪进行过程模拟, 将有效波高模拟结果分别与B01—B04四个浮标站实测值进行对比分析。

图 6给出了台风浪模型有效波高模拟值与浮标观测值的时间过程对比。可以看出, 波高的对比结果与风速的对比结果相一致: 最大风速半径参数R_max对模拟结果影响较大, 选用Willoughby表达式的3个组次和选用Powell表达式的3个组次模拟结果差异明显, 选用相同R_max表达式的组次之间较为相似。B1R1、B2R1、B3R1对于4场台风浪有效波高的模拟总体上优于B1R2、B2R2、B3R2, 更加贴合实测曲线。B1R2、B2R2、B3R2这3个组次普遍存在低海况下对有效波高的低估和高海况下对有效波高的高估。相较于B1R1和B2R1, B3R1组与实测值变化趋势最为吻合, 模拟效果最佳。

图 7给出了各组次有效波高模拟值与实测值的散点密度分布。可以看出, 相较R1组, R2组对于有效波高＜6.0 m的模拟出现了总体低估的现象, 而在有效波高＞6.0 m时, 则有整体上的高估现象。R1组模拟值与实测值的相关系数较R2组更高, 散点的分散程度也更小。B3R1对6.0 m以上大浪的模拟比B1R1和B2R1更好, 散点的分布更为集中在对角线附近, 相关性也最高, 线性回归的决定系数R²为0.846, 相关系数为0.92。

表 3进一步给出了6个参数组合方案下有效波高模拟值与浮标实测值的误差指标, 由此可见: R1组的模拟效果略优于R2组, 在R1的3个组次中, B3R1在每个浮标处的表现均为最优。

综上分析可见, 无论是从风速角度还是从有效波高角度来看, B3R1参数组合方案下的模拟效果均为最佳, 适用于浙江沿海台风浪的模拟后报。

3.3 基于Jason卫星有效波高观测结果的模型验证分析

进一步利用卫星轨道有效波高观测数据, 对B3R1方案下的波浪场进行验证分析。图 8给出了10场台风影响期间卫星观测的有效波高沿轨分布及其与对应时刻有效波高模拟结果的比较, 其中, 台风“梅花”和“利奇马”分别有2条符合的卫星轨道数据, 其余台风均有1条; 为风浪场分析需要, 台风路径、卫星轨道生成时刻的模拟结果(包括风向和平均波向)的空间分布也一并绘于图中。

由图 8可知, 台风路径右手侧风向和平均浪向大致相同, 而在左手侧差异较大。由图 8(d)、(e)、(g)、(h)、(j)、(k)可以看出台风风场和浪场呈现出非对称性: 台风路径右手侧波高大于左手侧, 造成这种不对称分布的原因是台风移动方向右手侧风速强于左手侧, 且右手侧风、浪同向, 海浪始终被锁定在风场中, 持续不断获得能量^[23-24]。在台风中心附近的大浪区, B3R1方案下各台风浪有效波高模拟值空间分布同卫星波高空间分布吻合良好, 对于部分场次的模拟, 在台风外围出现模拟值偏低的现象, 这可能是因为在模拟时未考虑边界入射波浪所致。

对上述10场台风浪有效波高的模拟值与卫星观测值进行散点密度统计分析, 如图 9(a)所示。可见, 有效波高模拟值与观测值整体吻合良好, 散点集中分布在对角线附近, 线性回归的决定系数R²为0.811, 相关系数达到0.90; 模型对8.0 m以下有效波高的模拟效果良好, 对8.0 m以上有效波高的模拟出现了略偏小的现象, 这可能与高海况下卫星反演有效波高精度不足有关。为考察有效波高模拟的空间分布效果, 计算了卫星轨道上各数据点的模拟值较观测值的偏差以及该数据点与台风中心的距离, 以数据点与台风中心距离为横坐标, 绘制了有效波高偏差随距离的空间变化, 如图 9(b)。可以看出, 样本点出现在距台风中心约1 000 km范围内, 其中200 km以内和900 km以上的样本数量相对较少, 200~900 km分布区内的样本数量最多。考虑到样本数量空间分布不均, 表 4进一步以200 km为间隔, 统计给出了各区间内的样本数量、有效波高模拟值的平均偏差(BIAS)与相对平均偏差(RBIAS, 定义为平均偏差与实测波高均值比)、均方根偏差(RMSB)与相对均方根偏差(RRMSB, 定义为均方根偏差与实测波高均值比)以及相关系数(R)。定量的误差分析表明: 有效波高平均偏差为−0.30 m, 相对平均偏差为−7.0%, 各区间内的相对平均偏差最大为−8.1%, 发生在距台风中心600~800 km内; 均方根偏差为0.74 m, 相对均方根偏差为17.6%, 相对均方根偏差最大为21.5%, 出现在距台风中心400~600 km内; 相关系数为0.90, 各区间上相关系数的最小值为0.78, 发生在距中心400~600 km内。应该指出: 在近台风中心200 km内的大风浪区, 有效波高模拟值具有更低的相对均方根误差(13.6%)和相对较好的相关系数(0.88), 这说明了模型在大风浪区具有较高的精度; 在距中心400~600 km内, 模拟效果略差, 这可能与模型风场与背景风场的融合方法有关, 是后续研究需要关注的问题。

4 结论

融合Holand参数化模型风和CFSR背景风构建台风风场, 以此作为驱动条件, 使用非结构化网格的SWAN模型, 对浙江海域的台风浪进行了数值模拟。为验证模型的准确性和可靠性, 选择了2011年以来发生的10次典型路径和强度的台风浪事件作为个例研究对象, 通过对模拟结果与观测资料(包括气象站风速、波浪浮标以及卫星观测资料)进行广泛的对比分析, 探讨了参数化模型风速剖面参数B和最大风速半径参数R_max组合方案对台风浪模拟的影响, 并确定了B和R_max两个参数组合的最佳方案。

(1) 最大风速半径分别采用了Willoughby和Powell表达式, 结果显示两者之间存在较大差异, 选用Willoughby表达式的3个组次(B1R1、B2R1、B3R1)在模拟效果上表现更佳。

(2) 风速剖面参数分别采用了Hubbert、Willoughby和Vickery表达式, 结果表明, Willoughby的R_max表达式搭配Vickery的B表达式(即B3R1)在浙江海域台风浪的模拟中表现最佳, 其与气象站风速和浮标有效波高实测值在时间过程上最为贴合, 在散点分布上最为集中, 在误差分析中也有最小的均方根偏差和最高的相关性系数。

(3) 基于Jason卫星观测资料对B3R1方案的验证分析显示: 在近台风中心1 000 km范围内, 模型对有效波高空间分布的模拟效果较好, 其与观测值的相关系数可达0.90, 相对平均偏差和相对均方根误差分别为−7.0%和17.6%; 特别值得注意的是, 相较近台风中心400~600 km范围内的风浪影响区, 模型对近中心200 km内的大风浪区有效波高的模拟表现更佳。