WAVEWATCH Ⅲ不同海冰源项的海浪模拟效果对比

2020-10-09苗琪徐福敏俞茂玲

海洋学报 2020年9期

苗琪，徐福敏，俞茂玲

( 1. 中国电建集团华东勘测设计研究院有限公司，浙江杭州 311122；2. 河海大学港口海岸与近海工程学院，江苏南京210098)

1 引言

20 世纪50年代以来，作为北冰洋一部分的波弗特海（Beaufort Sea）近海大陆架被发现具有丰富的石油、天然气储量。随着全球变暖，北极海冰储量大量减少[1–3]，北极航道潜力巨大。北极地区开采所得的矿产利用北极航道经由波弗特海运往世界各地，较传统方式可节约近一半时间，因此近年来随着气候变暖和亚洲天然气需求激增，各国均对该海域的极地研究展现出极大关注。对北极风暴下波弗特海−马更些河（Mackenzie River）河口水域的海浪特性进行分析研究，是油气开采以及北极航道开发中的关键一环。

风暴下短期内海冰密集度、海冰厚度、冰缘线位置均剧烈变化的动态海冰，往往会给北极地区海浪的精确模拟带来巨大困难。国外早期针对波弗特海海域的波浪模拟及海浪特性研究中[4–7]，在海冰的处理上多采用的是未经可靠验证的模型或经验公式，或着眼于夏季无海冰覆盖的近岸区域，对海冰−海浪间的相互作用考虑较为粗糙，且关于北极风暴过程中海冰的短期运动变化、不同风暴下的海浪分布特征的研究依然较少，因此对海冰影响下波弗特海至马更些河河口海域的海浪进行精确数值模拟，以及对如何在海浪数值模拟中正确计量海冰的影响开展研究具有一定的研究意义。

海浪模式是现今较为成熟的海浪计算及预报手段。第三代海浪模式WAVEWATCH Ⅲ在控制方程、程序结构、数值方法等很多方面都作出改进，因此近年来在海浪业务预报上得到了广泛的应用[8–11]。在版本V5.16 中，WAVEWATCH Ⅲ模式已具备多种海冰对海浪的能量耗散源项(IC1−IC4，IS1−IS2)，能在极地地区海浪数值模拟中相对精确地考虑海冰的影响，因此近年来，一些学者开始利用WAVEWATCH Ⅲ模式开展波弗特海海域的海冰−海浪相互作用研究，如Erick 等[12]利用WAVEWATCH Ⅲ模式中的IC3源项对秋季波弗特海的波浪进行后报，并利用后报结果和浮标实测数据分析了不同冰种对海冰的能量耗散能力；Cheng 等[13]利用实测冰情、浪情率定了WAVEWATCH Ⅲ模式中IC3 源项的参数，并对波弗特海海域的海浪进行模拟。这些研究多针对特定海冰源项，而对于WAVEWATCH Ⅲ模式中的多种海冰源项缺乏详细而完整的评述。据此，本文利用WAVEWATCH Ⅲ模式中的不同海冰源项分别对海冰影响下波弗特海夏季开冰期的暴风浪进行数值模拟，以期为极地地区海洋工程建设、海洋防灾减灾等提供一定的参考。

2 WAVEWATCH Ⅲ模式介绍

2.1 控制方程

WAVEWATCH Ⅲ模式中，球坐标系下关于动谱密度的控制方程通常表示为

式中，N为动谱密度；t为时间； ∇x为水平拉普拉斯算子；k为波数矢量；θ为波向；S为源函数；σ为圆频率；Cg为群速度；U为环境流速；d为水深；s为与波向相同的坐标；m为与s垂直的另一坐标。

2.2 海冰源项

WAVEWATCH Ⅲ模式中海冰损耗作用最传统的处理方式是采用经验方法[14]，近似地将海冰密集度大于某临界值（由用户指定）的网格点视作陆地，小于某临界值的视为不受海冰影响的开阔水域，两者之间进行线性插值，并通过修改控制方程使波能在每个海冰网格点处进行百分比传递，该方法记为IC0。WAVEWATCHⅢ V5.16 版本中引入了4 种海冰损耗源项，分别以IC1−IC4 来表示：第一种海冰损耗源项（IC1）对所有的波浪成分都采用相同的耗散率，是一种不依赖于频率的简单损耗过程，在计算时需要输入海冰的密集度以及损耗率；第二种海冰损耗源项（IC2）将海冰视作一个连续的弹性薄盘，能量损耗由冰下边界层的摩擦产生[15–17]，在计算时需要输入海冰密集度、海冰有效剪切模量、黏度；第三种海冰损耗源项（IC3）认为海冰是一个黏滞弹性层[18]，对损耗的计算需要输入海冰的密集度、厚度、有效黏度系数、密度和有效剪切模量；第四种海冰损耗源项（IC4）对能量损耗给出了几组简单的依赖于频率的经验或参数化的公式[19–24]，从而使得源项能够简单灵活却有效地再现冰−浪相互作用过程中海浪能量的衰减。

散射源项有两种，分别以IS1 和IS2 表示：第一种散射源项（IS1）是比较简单的海冰对海浪的守恒作用，仅考虑海冰密集度，海冰只散射海浪能量，但是不耗散海浪的能量；第二种散射源项（IS2）依赖于浮冰尺寸，并且考虑了海浪作用下海冰的破碎从而能够不断更新最大浮冰尺寸[24]。

3 验证资料来源及模型设置

3.1 验证浮标资料

2014年7月末，陆缘海冰带（Marginal Ice Zone，MIZ）项目由美国海军在波弗特海/楚科奇海开展（http://www.apl.washington.edu/project/project.php?id=miz），范围为72°～79°N，137°～170°W。项目过程中共投放3 个配备定位系统的SWIFT（Surface Wave Instrument Floats with Tracking）漂移浮标，并于9月底回收。SWIFT10 和SWIFT11 浮标于7月27 日在阿拉斯加北部约100km 的近岸处从RV“Ukpik”号极地考察船上投放，向西北方向漂移。在9月1 日之前两者漂移路径大致相同，然后SWIFT10 浮标开始进入高密集度海冰区域，直至9月15 日。SWIFT15 浮标在8月5−17 日期间由R/V “Araon”号极地考察船投放并回收，且全程被海冰环绕。各浮标漂移路径见图1。

由于SWIFT15 浮标全程被高密集度海冰包围，几乎处于搁浅状态，大部分时间测得有效波高均为0m，最大有效波高不足0.7m，故本次研究仅选取SWIFT10和SWIFT11 两个浮标用作后续验证。

图1 SWIFT 浮标漂移路径Fig. 1 The track of SWIFT drifting buoys

3.2 模型设置

2014年8−9月有多场风暴袭击波弗特海至马更些河河口水域，且影响范围较广。研究区域过大，通常会导致计算效率过于低下，故研究采用自嵌套的方案，即先用分辨率较低的网格对较大的海域进行计算，然后将计算所得波浪谱作为边界输入到内层分辨率更高、范围更小的模型中进行计算，以提高计算效率。嵌套外层包括整个东西伯利亚海、楚科奇海、波弗特海和部分位于西北太平洋的白令海，内层为波弗特海至马更些河口外水域，据此建立自波弗特海至马更些河口的WAVEWATCH Ⅲ两级嵌套模型，由外层的计算获取传播至波弗特海的所有低频涌浪成份，输入到第二层中，然后在存在大量海冰的波弗特海进行风−冰−浪相互作用模拟研究，并探究WAVEWATCH Ⅲ对海冰影响下的海浪场模拟的准确性。

嵌套计算时，中低频截断频率取为0.04118 Hz，高频截断频率取为0.7904 Hz，以对数分布fn+1=1.1fn划分为32 个频率网格，外层谱方向划分为30 个波向，分辨率为12°，内层谱方向划分为36 个波向，分辨率为10°。

3.2.1 计算区域和水深

本文建立的WAVEWATCH Ⅲ两级自嵌套海浪模型网格采用球面坐标，各层范围为：（1）嵌套Ⅰ区：整个东西伯利亚海、波弗特海和部分位于西北太平洋的白令海，范围为58°～78°N，140°E～110°W，空间分辨率为0.25°×0.25°，网格节点数为441×81；（2）嵌套Ⅱ区：波弗特海至马更些河河口外水域，范围为65°～75°N，125°～175°W，空间分辨率为5′×5′，网格节点数为601×121（图2）。

图2 两级嵌套模型研究区域示意图Fig. 2 Two-level nested computational domain

模型水深数据采用美国国家海洋与大气管理局（NOAA）的ETOPO1 全球地形数据集，其空间分辨率为1′×1′，在各级嵌套区域分别插值至各网格点处。

3.2.2 驱动风场

模型风场数据采用CCMP（Cross-Calibrated Multi-Platform）V2.0 风场，它是CCMP 风场的升级版本，在继承CCMP 风场高时间、空间分辨率的基础上（时间分辨率为6 h，空间分辨率为0.25°×0.25°），又将数据覆盖时间从2011年12月延长至今。对2014年8月1 日00 时至9月30 日18 时的CCMP V2.0 风场进行空间线性插值，作为模型驱动风场。

鉴于CCMP V2.0 风场在相关研究中应用暂时较少，将CCMP V2.0 风场的风速与SWIFT 漂移浮标的实测风速进行对比，以验证其可靠性。将CCMP 风速分别插值至每个时刻两个浮标所在位置处，对比结果如图3 所示。可以看出，两条风速曲线趋势和变化几乎完全一致，数据吻合较好，证明研究所用的CCMP V2.0 风场风速与波弗特海的风速吻合较好，能够较好地反映北极风暴过程。

3.2.3 海冰资料

海冰密集度数据资料来源有两种。第一种为NOAA 最优插值海面温度分析数据集(NOAA Optimum Interpolation 1/4 Degree Daily Sea Surface Temperature Analysis, Version 2)，被用于嵌套模型的外层，即从白令海峡至波弗特海的整个海域。该数据集由NOAA/NCDC 提供，空间上覆盖全球，分辨率为0.25°×0.25°，时间上从1981年9月1 日至今，分辨率为1 d。第二种为MASAM2 逐日4km 北极海冰密集度数据集（MASAM2: Daily 4km Arctic Sea Ice Concentration,Version 1），被用于模型的第二层，即波弗特海至马更些河河口。该数据集时间覆盖范围为2012年7月至今，分辨率为1 d，空间覆盖范围为29.08°～90°N，环全球经度，空间分辨率为4km×4km。

海冰厚度数据采用美国国家环境预报中心（NCEP）气候预报系统2.0 版本（CFSV2）的逐时时间序列产品（NCEP Climate Forecast System Version 2（CFSV2）Selected Hourly Time-Series Products），范围覆盖全球，有0.2°×0.2°、0.5°×0.5°、1.0°×1.0°和2.5°×2.5°的空间分辨率可供选择，时间覆盖为2011年1月1 日至今，分辨率为6 h。CFSV2 海冰厚度数据相比实测数据往往偏大[25–26]，故在本次模拟中对其乘以0.6 的折减。

3.2.4 参数设置

计算时，对于海冰的散射作用，由于计算海域中浮冰直径相对波长来说较小（可达2m，但通常小于1m），能带来的散射作用有限，且WAVEWATCH Ⅲ模式中的海冰散射模块尚不能很好地体现海冰对海浪能量的衰减作用，故本次计算不考虑海冰的散射作用。对于海冰对海浪能量的衰减作用，分别采用IC0、IC1、IC2、IC34 种机理进行考虑，其中IC0 额外设置海冰损耗上下限为50%～50%及25%～75%两种方案。其中50%～50%方案为模型的缺省设置方案，而25%～75%方案为一些学者采用过的改进方案。

图3 CCMP V2.0 风场风速验证Fig. 3 Validation of CCMP V2.0 wind speed

由于应用于该大范围海域时，缺乏海冰的剪切模量、黏滞系数等实测属性参数，故除海冰密集度、厚度由卫星资料提供外，其余海冰相关的模式参数均取为模型缺省值。其他源项设置还包括风能输入和白浪耗散项选取Ardhuin 等[27]的WAM4 方案，四波非线性相互作用项选取DIA 方法[28]，三波非线性相互作用项选取LTA 方法[29]，底摩擦项选取JONSWAP 底摩擦公式[30]，水深变浅引起的波浪破碎项选取Battjes-Janssen 的方法[31]，参数均为默认值。

4 模拟结果及对比分析

4.1 模拟结果

由于SWIFT 浮标属于漂移浮标，所处位置随时间在不断变化，需利用WAVEWATCHⅢ模式中的沿轨迹输出功能输出两个SWIFT 浮标的轨迹在对应时刻的波浪要素以进行对比验证。WAVEWATCH Ⅲ模式的沿轨迹输出功能仅能输出一维海浪谱，利用公式编程由离散的海浪谱换算有效波高（其中m0为海浪谱的零阶谱矩），从而与实测数据进行验证比较。图4 为不同海冰方案下两个SWIFT 浮标处的有效波高模拟结果及浮标实测结果。

由图4 可见，在9月1 日之前，由于SWIFT10 浮标和SWIFT11 浮标均处于开阔水域中，海冰对两浮标唯一的影响为当风向背离海冰区时风区长度的不同，但该影响较为微弱，故此时不同的海冰耗散机理对浮标所在位置的海浪模拟并无明显影响，5 种机理结果并无明显区别，部分方案结果几乎完全一致，因此部分方案的结果被遮挡，但各方案结果均较为理想。

9月1−14 日，SWIFT10 浮标实测有效波高开始迅速衰减，从浮标挂载的摄像头录像可以发现，此时SWIFT10 浮标已被海冰所包围（图5），海面上大面积覆盖的海冰对波浪能量产生剧烈的衰减作用。这段时间内，IC0（25%）、IC2、IC3 有效波高模拟结果的趋势与实测数据相同，但数值明显偏高，IC1 趋势与实测相同，且数值相对较为接近。IC2 模拟结果与IC1 非常接近，几乎被完全遮挡（从下文误差分析中也可以看出两方案差距较小，SWIFT10 浮标处IC2 误差略微高于IC1，而在SWIFT11 浮标处略低于IC1）。IC0（50%）在该时段内波要素始终为0。

9月14 日以后，SWIFT10 浮标重新回到开阔水域，IC0（50%）方案仍围绕观测值出现剧烈上下波动，而其余方案均较为理想。与SWIFT10 浮标不同，SWIFT11 浮标全程位于冰密集度较低的开阔海域内，从图4b 可以看出该浮标的有效波高模拟结果趋势与9月14 日后的SWIFT10 浮标结果类似，即IC0（50%）方案出现异常波动，而其余方案均拟合较好。

IC0（50%）的参数设置，根据其原理，是一种不连续的处理方式，会对海冰密集度大于50%区域的海浪能量产生剧烈衰减，而对海冰密集度小于50%的区域衰减较微弱。根据卫星及遥感数据，浮标所在的波弗特海域在计算时间段内海冰密集度多为40%～70%，这就导致在该密集度范围下，IC0（50%）方案对密集度变化过于敏感，结果极易产生突变及波动。在9月1−14 日，海冰密集度大于50%，故该时段内IC0（50%）方案模拟结果波高与周期均为0。

图4 不同海冰方案下两个SWIFT 浮标有效波高对比Fig. 4 Comparison of simulated significant wave height by two SWIFT bouys of different ice source terms

IC0（25%）方案虽然也是一种不连续的处理方式，会对海冰密集度大于75%区域的海浪能量产生剧烈衰减，而对海冰密集度小于25%的区域衰减较微弱，但其能在两密集度范围之间均匀过渡。因而虽然机理相同，但并未出现类似IC0（50%）方案的波高突变。

4.2 不同海冰模式下模拟结果误差分析

为更直观地定量比较不同机理对有效波高的模拟水平，确定该海域最佳的海冰机理，首先分别对不同机理的模拟结果和实测数据进行显著性水平检验，假定模拟结果与实测数据之间相关性为0，计算得p值均小于0.05，证明存在相关性。其次，选取均方根误差（Root Mean Square Error，RMSE）、相关系数（Correlation Coefficient，CC）和绝对平均误差（Mean Absolute Error，MAE）3 个参数来统计分析实测值和计算值之间的误差，分别定义为

图5 9月1−14 日环绕SWIFT10 的海冰Fig. 5 Sea ice pictures taken around SWIFT10 from September 1st to 14th

表1 和表2 分别为SWIFT10 和SWIFT11 浮标有效波高模拟值与浮标值的相关系数、均方根误差及绝对平均误差。

由表1 和表2 可以看出，各方案下SWIFT11 浮标处的有效波高与实测数据之间的均方根误差和绝对平均误差普遍小于SWIFT10，而相关系数高于SWIFT10，证明WAVEWATCHⅢ模式对有效波高的模拟效果在远离冰缘、水域相对开阔时比在被海冰控制下更好。此外，从SWIFT10 浮标处的有效波高结果可以看出，在离海冰较近，受其影响较大的水域，IC1 方案相较其他方案而言，模拟结果误差更小，相关系数更高。

表1 SWIFT10 有效波高模拟值与浮标值的均方根误差（RMSE）、相关系数（CC）和绝对平均误差（MAE）Table 1 RMSE, CC and MAE of observed and simulated significant wave heights in SWIFT10

表2 SWIFT11 有效波高模拟值与浮标值的的均方根误差（RMSE）、相关系数（CC）和绝对平均误差（MAE）Table 2 RMSE, CC and MAE of observed and simulated significant wave heights in SWIFT11

尽管已有相关研究[25]显示，在有较为完备的海冰属性参数时，IC3 源项表现出的对海浪的能量衰减特征更加符合实际情况，但在本次研究中，由于研究区域范围较大，并无能力获得除海冰密集度、厚度以外的其他海冰数据，如剪切模量、黏性系数等，而这些时空均变化的海冰数据恰恰是IC3 源项必不可少的输入参量，均取为缺省值显然不合理。相反，IC1 源项对海冰数据要求较低，因此在应用于大范围水域、资料匮乏的前提下，计算结果反而与实测数据更为接近。

本文所采用参数的模拟结果并非最优解，其他海冰−波浪源项在海冰数据充足、参数调试合理的情况下可能会得到更优的结果。但是在以下限定条件下：（1）应用于北极波佛特海这一特定海域的有效波高模拟，（2）缺乏海冰剪切模量、黏性系数等实测参数，（3）模型参数均采用默认值，使用IC1 海冰源项来计量海冰对海浪的衰减作用最为合理，能够为后续浅水地区的进一步海浪嵌套模拟提供更加精准的边界条件。