新的拖曳系数参数化方案在飓风“卡特里娜”和“瑞塔”过程海浪数值模拟中的应用研究

李文博，李煊，李庆杰，陈默，王立鹏

（1.山东省海洋生态环境与防灾减灾重点实验室，山东 青岛 266100；2.国家海洋局北海预报中心，山东 青岛 266100；3.海军航空大学，山东 烟台 264000；4.自然资源部烟台海洋中心，山东 烟台 264000；5.国家海洋局北海海洋技术保障中心，山东 青岛 266100）

1 引言

台风过程中伴随的大风、巨浪和风暴潮等海洋灾害会给人类海上活动及海岸设施带来巨大威胁，极端天气条件下海浪的研究及预报具有非常重要的意义。海表面拖曳系数是在高风速下进行海浪数值模拟的关键因子之一，选取合理的拖曳系数将直接决定海浪数值预报或后报的精准度[1-2]。

在海气相互作用中，海面风应力是海浪产生的主要驱动力之一，海气之间的动量通量通常由海表面拖曳系数或者海面空气动力学粗糙长度来表征。最初研究认为海表面拖曳系数是常量，随着研究的深入，人们逐渐发现海表面拖曳系数与风速之间存在一定的数学关系。GARRATT[3]通过研究给出了Charnock无量纲粗糙度表达式中常数α的取值以及海表面拖曳系数与海表面10 m 高度处风速（U10）的参数关系；SMITH[4]用涡相关法分析得出海表面拖曳系数随风速（6～22 m/s）的变化规律；WU[5]提出了海表面拖曳系数的线性参数化方案，并认为在飓风条件下海表面拖曳系数随风速的变化仍遵从线性关系。LARGE 等[6]利用雷诺动量耗散法对海表面动量通量进行了估算，认为海表面拖曳系数在风速为 4～11 m/s 时为常数，大于 11 m/s 后与U10呈线性增长关系。GEERNAERT[7]获取了北海的观测数据，利用涡相关法并考虑波龄的影响进行分析，得出了4～24 m/s 风速条件下海表面拖曳系数的线性表达式。以上研究均认为海表面拖曳系数或动力学粗糙长度随着风速增加呈线性增长。然而近年来观测资料、实验室数据以及数值模拟结果均表明高风速下海表面拖曳系数的增长有所减缓，甚至出现减小的趋势。ALAMARO 等[8]通过实验得出在风速达到25 m/s 后海表面拖曳系数随风速的增大而减小；MAKIN[9]利用Barenblatt提出的有限饱和沫滴悬浮层理论，建立了高风速下海表面动力学粗糙度模型，并指出飞沫水滴组成的飞沫边界层是海表面拖曳系数在高风速下发生衰减的原因；HOLTHUIJSEN等[10]根据观测指出海表面拖曳系数在风速达到40 m/s 后开始随风速的增加而减小，当风速达到80 m/s时海表面拖曳系数将降为零。

在台风浪的模拟过程中，风应力的准确估计及其对海洋上层影响的表征十分重要。海表面拖曳系数是海气动量通量参数化研究的关键，因此本文着眼于高风速下海表面拖曳系数的参数化方案研究，期望得到高风速下更为准确且便于数值模式应用的海表面拖曳系数表达式，以提高台风过程中数值模拟预报的准确性。

2 研究方法

2.1 海浪数值模式

海浪数值模拟和预报工作最早开始于20 世纪40 年代，经过前人的不断探索与技术的发展，海浪数值模拟预报模式已经从第一代发展到第三代。荷兰Delft理工大学（Delft University of Technology）对第二代海浪模式进行了改进和优化，并研制开发出第三代近岸海浪数值模式（Simulation WAves Nearshore，SWAN）。经过多年的改进，SWAN 模式趋于成熟，目前已被广泛使用。SWAN 模式采用全隐式有限差分格式和基于能量守恒原理的平衡方程，它将波浪进入浅水后的传播特点考虑在内的同时保证计算过程中时空网格能够相对稳定地进行，而且计算空间网格和时间步长上不会受到牵制。在平衡方程的各源项中，SWAN 模式不仅考虑了风输入项、四波相互作用、破碎和摩擦项等，还考虑了波浪的深度破碎作用和三波相互作用。该模式包括了当前海浪预报的较新成果[11-13]。

2.2 海表面拖曳系数

海气相互作用的一个重要形式是海气界面的动量交换过程，风速脉动会产生动量通量，海表面风应力是驱动上层海洋环流和风浪的主要因素，因此，如何表征风应力对海气界面的影响至关重要。在海气相互作用的研究中，通常开展海气动量通量的参数化研究，并用块体公式来计算海面风应力。在这一研究过程中，海气动量通量一般由海表面拖曳系数或者海面空气动力学粗糙长度进行参数化。公式如下：

式中：Cd,z为海表面Z高度处拖曳系数；Uz为海表面Z高度处平均风速；U0为海表面处平均风速，通常取值为0，故上式可写为：

拖曳系数描述的是动量通量在海-气界面的传输程度。在流体动力学中，阻力系数即表面拖曳系数Cd是用于量化流体环境（例如空气或水）中物体阻力的无量纲量，较低的阻力系数表示物体具有较少的空气动力学或流体动力学阻力，阻力系数总是与特定表面积相关联。Cd公式如下：

式中：Fd是阻力；ρ是流体的质量密度；u是物体相对于流体的运动速度；A是作用区域。任何对象的阻力系数均受表面摩擦阻力和寄生阻力两个基本因素影响，其中寄生阻力是在拖曳物体穿过大气中流体介质过程中由气态介质的空气动力阻力引起。Cd并不是常数，而是作为流速、流动方向、物体位置、物体尺寸、流体密度和流体粘度的函数。流速、流体粘度和物体尺寸也被结合到被称为雷诺数（Re，惯性力和粘性力的比值）的无量纲量中，以便量化给定条件的流动，Cd是Re的相关量。在可压缩流体中，Cd同时也是马赫数（Ma，表示超过边界的流速与声音的局部速度的比率的无量纲量）的函数。对于某些形状，Cd只取决于Re、Ma和流动方向，当Ma很低时，阻力系数与Ma无关。

基于Cd与U10具有相关性这一前提，依据有限的实测资料和不同的参数化方法，前人给出了低风速条件下Cd随风速线性变化的关系。图1为不同研究者获得的Cd与U10的线性关系图。

图1 Cd与U10线性关系图Fig.1 The linear relationship between Cd and U10

低风速下Cd线性化参数方案基本成型，但高风速下Cd与U10的关系仍是需要继续探讨的问题。近年来基于外海观测数据、实验数据和数值模型数据的研究都显示出高风速下Cd随风速的增长趋于平缓，并且在风速的持续增长下Cd将出现减小的趋势。究其原因，可能是由于白帽破碎和飞沫造成的。

在SWAN 模式中，默认使用WU[5]得出的关于风速变化的Cd表达式：

2.3 海表面拖曳系数的参数化方案

COARE 3.0[14]模型是国际上广泛采用的海气动量通量计算模型，其中海表面动力粗糙度长度有3种参数化方案，分别由 YELLAND 等[15]、TAYLOR等[16]和OOST 等[17]提出，本文分别称为 YT96 方案、TY01方案和O02方案。

（1）YT96方案仅考虑了风速的影响：

式中：α是一个通用常数或者无量纲化粗糙长度；v是运动粘度；Z0为海表面动力粗糙度长度;u*为摩擦速度。

（2）TY01 方案在风速的基础上同时考虑了波浪的影响：

式 中 ：HS为 有 效 波 高 ，0.015u10N)；LP为主波波长，tw= 0.729u10N;HS/LP近似为主波波陡；v为运动粘度；u10N为海面10 m风速。

（3）O02方案也考虑了波浪的影响，但是选取了与TY01方案不同的波浪参数：

式中：LP为主波波长；CP为主波向速度为波龄；v为运动粘度。

Cd的参数化方法通常有3种：涡相关法、廓线法和惯性耗散法。廓线法内容为：在高风速条件下达到稳定的中性层结大气，大气边界层的垂直风廓线可以由以下对数关系式表达：

式中：u(z)表示海表面Z高度处的风速；u*为摩擦速度；k为卡曼系数，一般取值为0.4；z0为海面粗糙度。在流体力学中，摩擦速度通常用来与真实速度进行比较，摩擦速度可以用来重写剪切应力。根据摩擦速度的定义：

式中：τ 表示流体的任意层中的剪切应力；ρ表示流体的密度。结合式（2），拖曳系数可以写为：

再结合式（8）可以得到海表面拖曳系数与海面粗糙度的关系为：

为研究高风速下Cd与风速的关系，本文使用美国国家浮标数据中心（National Data Buoy Center，NDBC）的高风速实测数据作为研究资料。本文选用了9 个飓风过程期间8 个浮标的12 份浮标数据。所选浮标的有效观测数据中，最大风速达到40.9 m/s，风速在25～41 m/s 的有效数据超过25 个，这表明高风速区域海表面拖曳系数的研究范围有了很大拓展，数据的准确度也有一定保证。根据TY01和O02方案并结合式（8）和式（10），可分别计算出两种方案下的Cd，在此基础上对Cd与U10的散点关系进行二项式拟合，即可得出两种方案下Cd与U10的参数化公式。由于YT96 方案本身即为与风速相关的关系式，故不再对其进行拟合，将YT96 方案作为新的参数化关系的对比方案之一。图2 和图3 为TY01 和O02 方案计算所得的Cd及其拟合函数的曲线图。

图2 TY01方案所得Cd及其拟合曲线Fig.2 Cd and its fitting curve obtained from the TY01 scheme

图3 O02方案所得Cd及其拟合曲线Fig.3 Cd and its fitting curve obtained from the O02 scheme

所得拟合函数曲线中，O02 方案中高风速区的Cd随风速增大而增大的趋势有所减缓，而TY01 方案在风速增大到一定程度后Cd的增长更为显著。根据所得散点图进行分段二项式拟合，得到两个Cd与U10的参数化公式：

TY01方案为：

对比两种方案所得的Cd与U10的散点图及参数化公式，统计参数包括误差平方和（Sum of Squares for Error，SSE）、相关系数（R2）、调整后的R2以及均方根误差（Root Mean Square Error，RMSE），结果见表1。通过综合比较分析可知，TY01 方案下由U10计算Cd存在两个问题：Cd在高风速下的计算结果数值过大，有失准确；Cd与U10的变化趋势在高风速处陡增，与实际相反。故本文淘汰TY01 方案所得的参数化公式，选择O02 方案的拟合结果（以下称新O02 方案）进行下一步的实测对比验证。实测数据选自中国南海海洋观测平台，站点位置距离海岸约6.5 km，测量数据范围全长67 m，水下14 m。

表1 TY01和O02两个方案拟合函数对比统计结果Tab.1 Comparison of the statistical results of the fitting functions obtained by the TY01 and O02 schemes

将新O02 方案与YT96 方案以及根据实测u*和U10得出的Cd进行对比，结果见图4。从图中可以看出两个函数曲线与实测结果的变化趋势一致，数量级相同，故新O02 方案能够合理地反映Cd与U10的变化关系。

图4 YT96、新O02方案与实测数据对比图Fig.4 Comparison chart between YT96 and new O02 schemes and measured data

新 O02 方案的R2为 0.997 0，SSE 和 RMSE 分别为2.273 3 × 10-5和8.250 1 × 10-5，数量级均为10-5，表明该方案下的参数化拟合函数与实测数据的拟合效果良好。

3 参数化拖曳系数在台风浪数值模拟中的验证

3.1 飓风“卡特里娜”和“瑞塔”介绍

风暴“卡特里娜”于 2005 年 8 月 23 日（北京时，下同）在巴哈马群岛附近生成，24 日增强为飓风后，于佛罗里达州南部以1 级飓风强度登陆，数小时内其强度迅速增强为5 级，近中心最高持续风速为150 n mile/h，29日在密西西比河口登陆后强度减弱为3级，之后系统加速向东北移动，31日在俄亥俄州转化为温带气旋。风暴“瑞塔”于2005 年9 月18 日在特克斯与凯科斯群岛形成，20 日早上升级为1 级飓风，21 日进一步升级为5 级飓风，24 日在德克萨斯州萨宾帕斯与路易斯安那州卡梅伦之间的海岸地区登陆。虽然飓风“瑞塔”登陆时已减弱为3级飓风，但其最大风速超过190 km/h，给沿海一些中小城镇造成了严重的水灾和不同程度的财产损失。

3.2 海浪模型的设置和数据来源

本文选定的计算模拟区域为79°～98°W，18°～32°N。水深采用全球水深地形模型ETOPO1 数据，利用非结构网格，网格边界分辨率是1/20°。采用美国国家海洋和大气管理局（National Oceanic and Atmospheric Administration，NOAA）提供的浮标数据作为实测数据进行对比验证。根据飓风路径及NDBC 所提供的浮标位置找到6 个有效浮标数据，浮标经纬度位置及水深见表2。

表2 浮标信息表Tab.2 Buoys information table

本文SWAN 模式中的风场采用模型风场与CCMP（Cross-Calibrated，Multi-Platform）风场作为背景风场构造的合成风场，模型风场由圆形对称风场Jelesnianski-II 模型表达式得出[18-19]，背景风场采用经过线性插值后的CCMP 风场，时间步长为1 h，空间分辨率为5′×5′。采用随风速变化的权重系数e构造合成风场，构造公式为：

式中：Vm为模型风场；Ven为背景风场；c是一个与台风影响范围有关的参数；系数n一般取9 或10；r表示计算点到台风中心的距离；Rmax为最大风速半径。合成风场综合了模型风场与背景风场的优势，在飓风中心附近模型风场中起决定性作用，而在台风外围背景风场中起主导作用。

误差分析表明，飓风“卡特里娜”合成风场与实测风场的平均绝对误差（Mean Absolute Error，MAE）为 1.93 m/s,平均相对误差（Mean Relative Error，MRE）为20.4%；飓风“瑞塔”合成风场与实测风场的 MAE 为 1.78 m/s，MRE 为 18.8%。上述数据说明合成风场基本能够反映真实风场的结构和变化。

3.3 数值实验结果及分析

利用新O02 方案分别对飓风“卡特里娜”与“瑞塔”的台风浪过程进行数值模拟，并与YT96 方案和SWAN 模式默认方案进行对比。图5 为飓风“卡特里娜”过程有效波高实测数据与模拟值的对比，图6为飓风“瑞塔”过程有效波高实测数据与模拟值的对比。根据图像可知：

图5 飓风“卡特里娜”过程有效波高实测数据与模拟值的对比Fig.5 Comparison of measured and simulated values of significant wave height during hurricane"Katrina"

图6 飓风“瑞塔”过程有效波高实测数据与模拟值的对比Fig.6 Comparison of measured and simulated values of significant wave height during hurricane"Rita"

（1）对于飓风“卡特里娜”台风浪过程，YT96 方案在浮标42001 和42002 处的模拟结果更好，而在浮标 42019、42035、42039 和 42040 处，新 O02 方案的参数化关系对台风浪的模拟更接近实测值。

（2）对于飓风“瑞塔”台风浪过程，YT96 方案在浮标42002 和42019 处的模拟结果更好，而在浮标42001、42035、42039和42040处，新O02方案的参数化关系对台风浪的模拟更接近实测值。

（3）综合来说，新O02 方案参数化关系的模拟效果优于SWAN 模式默认方案，而与YT96 方案的模拟结果各有千秋。

通过飓风“卡特里娜”和“瑞塔”海浪数值模拟有效波高与实测数据的对比发现，新O02 方案中海浪数值模拟结果与实测浮标数据的趋势较为一致，且与实测结果符合较好，说明新的拖曳系数参数化方案能够满足海浪数值模拟的应用要求。有效波高的对比检验利用MAE 和MRE 参数进行，检验结果见表3和表4。

表3 飓风“卡特里娜”和“瑞塔”有效波高模拟值与实测数据平均绝对误差Tab.3 MAE of significant wave height between simulated values and measured data during hurricanes"Rita"and"Katrina"

表4 飓风“卡特里娜”和“瑞塔”有效波高模拟值与实测数据平均相对误差Tab.4 MRE of significant wave height between simulated values and measured data during hurricanes"Rita"and"Katrina"

经过计算，在YT96、SWAN和新O02方案中，有效波高模拟值和实测数据对比统计的MAE 分别为0.40 m、0.50 m 和0.38 m，MRE 分别为9.4%、12.2%和9.0%。

4 结论

本文利用9 次飓风过程中的12 份浮标实测资料，采用 COARE 3.0 算法得到基于 TY01 和 O02 两种考虑波浪方案下的海表面拖曳系数，通过二项式拟合得到两种系数关于风速的参数化关系式，其中TY01 方案的系数结果在高风速下未表现出“饱和”特征，O02方案的结果符合已有的观测规律，并且其与南海观测平台实测数据符合良好，因此选用基于O02的新参数化方案用于数值模拟研究。

为进一步验证高风速下海表面拖曳系数参数化表达式的准确性，本文选取了飓风“卡特里娜”和“瑞塔”两个过程进行海浪数值模拟。从模拟结果与实测数据的对比中可以发现，新O02 参数化海表面拖曳系数的模拟结果与NDBC 浮标实测数据符合更好。新的参数化关系方案的MRE 与SWAN 模式默认方案相比提高了3.2%，而与YT96 方案相比提高了0.4%，所以新的参数化关系方案能更准确地描述海表面拖曳系数在高风速下的变化，可以更好地用来进行高风速条件下的海浪数值模拟。