波浪能是一种可再生的清洁无污染能源, 其开发利用对于缓解不可再生能源危机和减少由化石燃料导致的环境污染具有重要意义(Borthwick, 2016)。波浪能是指海面上波浪的动能和势能, 主要是由风吹过海面引起的(尤再进等, 2017)。波浪能与波高的平方以及周期成正比。波浪能应用范围广泛, 已实现的用途有发电、取水、供热、海水淡化和制氢等(Wang et al, 2011; Ning et al, 2016; 王鑫等, 2020)。

我国的渤海海域受季风影响, 以风浪为主, 波高分布具有明显的季节性变化特征(郑崇伟等, 2013), 波能流密度的分布与有效波高相似, 均是秋冬大, 春夏小(郑崇伟等, 2013; 姜波等, 2017)。以往的研究表明, 渤海中部海域的波能流密度最大, 春季为1.6~ 2.4 kW/m, 夏季为0.7~1.4 kW/m, 秋季为2.4~ 3.4 kW/m, 冬季为4.2~5.4 kW/m (郑崇伟, 2014; 姜波等, 2017)。由此可见, 渤海波浪能最大的季节为冬季(万勇等, 2014), 然而冬季受西伯利亚强寒潮影响时, 一方面会导致风浪较大(郑崇伟, 2014), 另一方面又会导致海水结冰(Zhang et al, 2019, 2020), 海冰的形成会隔绝海-气交换, 导致波浪能的减弱(Zhang et al, 2019, 2020; 张娜等, 2021)。

渤海位于37°~41°N, 是北半球最南端的季节性海冰区(Zhang et al, 2019)。同时, 海冰已成为冬季影响渤海的主要海洋灾害之一(武浩等, 2016)。渤海的初冰日一般为12月上旬, 海冰首先出现在辽东湾近岸, 随后沿西海岸由北向南蔓延到渤海湾, 最后到达莱州湾, 于来年1、2月达到顶峰。其后春季来临, 气温升高, 海冰逐渐向北消融, 终冰日一般为3月中旬(Zhang et al, 2019)。渤海海冰主要集中在辽东湾、莱州湾和渤海湾地区, 海冰的覆盖率随离岸距离的增加而降低(Zhang et al, 2016, 2018)。我国渤海的单层平整冰的最大厚度通常不超过60 cm (白珊等, 1999)。辽东湾作为渤海海域冰情最严重的海域, 在重冰年冰厚一般为25~40 cm (白珊等, 1999)。以往针对渤海波浪能的研究大多没有考虑冬季海冰的影响。根据Zhang等(2020)的研究, 忽略海冰的影响会导致渤海波能流密度的高估, 而这种高估会严重误导电厂的选址。因此, 考虑海冰存在以准确估算渤海冬季波能流密度, 有待进一步研究。

以往的研究表明波流耦合模型能够有效提高近岸水域风浪计算精度与浅水水域波浪模拟精度(夏波等, 2013)。为了模拟更接近实际冰水共存海域动力环境的波浪演化过程, 本课题在Chen等(2018)前期工作的基础上, 采用MCT (model coupling toolkit)耦合器将FVCOM (finite-volume community ocean model)水动力学模型、海冰模型与SWAN (simulating wave nearshore)波浪模型实时耦合, 并在此基础上研究寒潮大风过程中渤海海冰对波浪模拟的影响, 从而实现冰水共存海域下波能流密度的准确评估。

1 数值模型

本文采用非结构化网格有限体积法的海洋模型FVCOM (finite-volume community ocean model)模拟渤海海冰(Chen et al, 2003)。该模型包括水动力学模块和非结构化网格的海冰模块, 其中水动力学模块涉及到潮、流、海温和盐度方程的求解, 海冰模块的求解涉及到海冰热力学、动力学、传输以及本构方程的求解(Chen et al, 2003; Gao et al, 2011; 张庆河等, 2013)。波浪的模拟采用第三代风浪模式SWAN (simulating wave nearshore), 版本号为41.01 (the SWAN Team, 2018)。SWAN模型是基于包含源汇项的动谱平衡方程描述风浪生成及其在近岸区的演化过程, 具体方程描述参见文献(Booij et al, 1999)。

1.1 海冰模型-水动力学模型-波浪模型的耦合原理

本文在Chen等(2018)的波流耦合系统的基础上基于MCT (model coupling toolkit)耦合器建立了海冰模型-水动力学模型-波浪模型的实时耦合模式。MCT耦合器是一套能够使独立的模型实现互相耦合的开源程序工具包, 它采用Fortran90进行编写, 能够支持串行和基于MPI (message passing interface)的并行运算。MCT耦合器能够为每个模型分配系统工作进程, 各子模型通过调用MCT模块实现变量传递。MCT耦合器的具体描述参见文献(Warner et al, 2008)。

耦合模型包括FVCOM水动力学模型、海冰模型和SWAN波浪模型, 由主程序调用各子模型同时独立运行, 当达到用户设定的数据交换时间时, 各子模型通过MCT耦合器进行实时变量交换, 交换结束后, 各子模式继续进行运算, 循环进行直至达到用户设定的结束时间, 结束运算。FVCOM水动力学模型、海冰模型与SWAN波浪模型采用了相同的非结构化三角形网格, 不需进行插值计算。

海冰模型-水动力学模型-波浪模型耦合模式之间的变量交换为: (1) 水动力学模型中的潮位(E_lev)、垂向平均流速分量(U_a、V_a); (2) 海冰模型中的海冰浓度(A_ice); (3) SWAN中的波长(L)、波向(D_ir)、有效波高(H_s)、谱峰波周期(T_p)、平均波周期(T₀₁)、底部波周期(T_bot)以及底部波动轨迹速度(U_bot), 如图 1所示。

1.2 由冰引起的波浪衰减原理

当达到用户设定的交换时间时, 海冰模型将模拟得到的海冰浓度与水动力学模型计算得到的潮位和垂向平均的流速分量传递给MCT耦合器, 由MCT耦合器将这些变量传递给SWAN波浪模型, 以此考虑潮流和海冰对波浪场的影响。同时FVCOM水动力学模型通过MCT耦合器接收SWAN波浪模型计算得到的包含冰影响的波长、波向、有效波高、谱峰波周期、平均波周期、底部波周期以及底部波动轨迹速度, 从而计算出新的潮位和流速传输给海冰模型用于计算新的冰浓度, 而新的冰浓度与新的潮位和流速又用于传输给SWAN波浪模型, 从而实现海冰模型-水动力学模型-波浪模型的耦合。当达到用户设定的结束时间时, 计算终止。

在冬季, 渤海海域寒潮大风频发, 强风是浪的主要驱动条件。在寒潮大风期间, 低气温导致海水结冰, 而海冰的出现会影响大气与海表面之间的动量交换, 从而影响了风浪的形成。根据Connolley等(2004)的结果, 在冰水共存的海域中, 海表面风应力可利用冰浓度的线性修正函数表示。因此本文将海冰模型和水动力学模型耦合模拟得到的冰浓度A_ice作为海表面应力的一部分, 以线性变化加权百分比的方式纳入SWAN模型的海面摩阻风速方程中, 用于模拟冰水共存条件下的波浪演化过程。在冰水共存海域的波浪模型中, 海面摩阻风速U_fric的表达式被改写为:

式中, C_drag为风拖曳力系数(Wu, 1982); A_ice为冰浓度, 由FVCOM模型的海冰模块计算输出(Thorndike et al, 1975); U₁₀为10 m高度处的风速。根据式(1)计算出的海面摩阻风速计算JONSWAP (joint North Sea wave project)谱; 接着再输入谱峰波周期对应的频率计算出当前时刻下的作用密度, 进而得到波能密度的零阶矩, 最终得到考虑冰影响的波浪参数。

1.3 波能流密度的计算公式

单位时间内波浪在传播方向上通过单位波峰宽度上的能量被定义为波能流密度, 用P表示, 单位为kW/m。由定义可得, 波能流密度的计算公式(2)可写为(Brooke, 2003):

式中, ρ为海水密度; g为重力加速度; C_g (f, h)为波群速度; S(f, θ)为波浪谱能量密度; f为波浪频率; h为水深; θ为波的传播方向。

有效波高(Ochi, 1998)通常定义为。因此, 公式(2)可改写为:

式中, 波群速度C_g(f, h)与水深h有关(Barbariol et al, 2013):

k为频率, 对应波数为(Dean et al, 1991)

或者

式中, L为波长; σ为波的角频率; T_e为波能量周期。

结合公式(4)~(6), C_g(f, h)波群速度可以表示为(Liang et al, 2017)

结合公式(3)和公式(7), 波能流密度公式可以进一步表示为

对于深水区(kh > π, 也即时), , 波能流密度公式可以进一步简化为(Liang et al, 2017)

1.4 数据来源及模型网格

FVCOM模型的输入数据包括气象数据(如气温、风、相对湿度、气压、云盖、短波辐射、长波辐射等)、开边界的驱动数据(如潮、海温等)以及初始场数据(海温、盐度)。SWAN模型的输入数据为10 m高度处的风和波浪开边界。

气象数据来自美国国家环境预报中心(National Centers for Environmental Prediction, NCEP), 其时间和空间分辨率分别为1 h和0.2° (Saha et al, 2010)。FVCOM模型的开边界潮位数据来自Chinatide (李孟国等, 2007), 开边界和初始场的海温和盐度数据来自HYCOM (HYbrid coordinate ocean model) (Chassignet et al, 2007)。SWAN模型的开边界波浪数据由第五代欧洲中期天气预报中心ECMWF (European Centre for Medium-range Weather Forecasts)再分析数据提供, 时间分辨率为1 h, 空间分辨率为0.5° (Hersbach et al, 2016)。以往的研究表明, 以上数据与实测值吻合较好, 能够为数值模拟提供合理的模型驱动条件(Zhang et al, 2020)。

地形水深来自渤海1︰15万海图地理信息系统数据(1999年), 结合近岸较高分辨率的海图数据以及实测地形数据加密, 最小空间分辨率为800 m, 如图 2所示。

2 模型验证

模型验证包括海冰和波浪两部分。海冰验证资料采用MODIS (moderate-resolution imaging spectroradiometer)卫星遥感图像, 波浪验证基于全国海洋立体观测网中的浮标观测波浪资料。

2.1 海冰模型验证

2012年2月8日, 渤海海冰覆盖范围最大, 模拟的渤海海冰面积约为29.79×10³ km², 其中辽东湾占比71.1%, 渤海湾占比22.3%, 莱州湾占比6.6%。模拟的渤海海冰面积占整个渤海面积的30%, 比MODIS的海冰面积偏小1.2%。渤海湾、辽东湾和莱州湾数值模拟的海冰分布及最大离岸距离与卫星遥感图像基本一致, 如图 3所示。

2.2 波浪模拟结果的验证

波浪模拟结果的验证取自2012年2月5~8日的一个寒潮大风过程, 受低气温影响, 此时渤海海冰达到年度最大值。如图 4所示, 2月6~8日期间, 实测日最大风速在14.4~15.0 m/s之间, 风向为NE和NW向。波浪模拟结果验证包括以下2种。(1) 不考虑海冰影响的波浪模拟结果; (2) 通过海冰模型-水动力学模型-波浪模型实时耦合模拟得到的波浪结果。从图 4中可以看出, 第2种实时耦合模式模拟得到的有效波高值与实测值最为接近, 尤其是最大值几乎与实测值吻合。第1种无冰影响的有效波高模拟值在大风期间要远远大于实测值, 最大有效波高比实测值偏大15%。大风期间不考虑海冰影响的有效波高模拟值与实测值的平均相对误差是实时耦合模式的1.7倍。

3 海冰对波能流密度的影响

2012年2月5日至8日, 渤海最低气温为−24℃, 平均气温为−20 ℃。渤海三个湾均被海冰覆盖, 分布范围达到最大值。从图 5可知, 2012年2月5~8日, 辽东湾近岸的海冰出现概率最大高达95%, 渤海湾近岸的海冰出现概率最大为50%, 而莱州湾的海冰出现概率最大为40%, 结冰区域面积占到了渤海总面积的三分之一。

3.1 平均波能流密度及其受海冰影响的空间分布

根据图 6, 2012年2月5日0时至2012年2月8日0时, 在考虑海冰影响的情况下, 辽东湾、渤海湾和莱州湾平均波能流密度的最大值分别为5.9, 6.1和5.6 kW/m, 分别位于各自的湾口位置处, 比不考虑海冰影响时分别减小了10.6%、1.6%和1.8%。辽东湾、渤海湾和莱州湾近岸区域的平均波能流密度较不考虑海冰影响时最多分别减少了100%、60%和50%, 见表 1。

3.2 平均波能流密度沿中轴线受海冰影响的程度

为了深入研究海水结冰对渤海三个湾波能流密度的影响, 图 7给出了渤海各湾中轴线处平均波能流密度受冰影响的程度, 中轴线位置见图 2a。2012年2月5日0时至2012年2月8日0时, 距辽东湾湾顶129.2 km是海冰沿中轴线处出现的最大离岸距离。其中沿辽东湾中轴线方向距湾顶78.2 km内, 海冰出现概率为95%, 其平均波能流密度受海冰影响最为严重, 减少了100%。而在距辽东湾顶78.2~129.2 km范围内, 随着海冰出现概率逐渐减小, 平均波能流密度受海冰影响的程度也越来越小, 从100%减小到32%。而在距辽东湾顶129.2~255.1 km范围内, 虽然没有海冰覆盖, 但平均波能流密度仍然受海冰影响而减小, 减小程度从32%到10%。这主要是由于海冰的出现阻隔了海-气之间的热量和动量交换, 相当于减小了渤海的有效风区, 从而使得无冰海域的平均波能流密度也会减小。

2012年2月5日0时至2012年2月8日0时, 渤海湾沿着中轴线距离湾顶的海冰最大覆盖范围为35.3 km, 其中距离湾顶21.6 km范围内的平均波能流密度受海冰影响最为严重, 最大减少了52%。在距渤海湾顶21.6~35.3 km范围内, 随着海冰出现概率的逐渐减小, 平均波能流密度受冰影响的减小程度从52%到21%。而在距渤海湾顶35.3~117.7 km范围内, 虽然没有海冰覆盖, 但平均波能流密度仍然受海冰影响而减小, 减小程度从21%到5%。

莱州湾与渤海湾情况类似。2012年2月5日0时至2012年2月8日0时, 莱州湾沿着中轴线距离湾顶的海冰最大覆盖范围为25.8 km, 在距湾顶8.0 km范围内, 平均波能流密度受海冰影响最为严重, 最大减少了50%。在距莱州湾顶8.0~25.8 km范围内, 随着海冰出现概率的逐渐减小, 平均波能流密度受冰影响的减小程度从50%到8%。而在距莱州湾顶25.8~77.1 km范围内, 虽然没有海冰覆盖, 但平均波能流密度仍然受海冰影响而减小, 减小程度从8%到3%。

3.3 无冰覆盖区域波能流密度随时间的演化

渤海的海冰主要位于渤海湾, 辽东湾和莱州湾的近岸区域, 且辽东湾的冰情最重, 如图 3所示。海冰的存在阻隔了海-气之间的热量和动量交换, 减小了渤海的有效风区, 可能导致无冰海域的波能流密度减小。为了探讨无冰覆盖海域的波能流密度在考虑和不考虑海冰影响时随时间演化的规律, 挑选了2012年2月5日0时至2012年2月8日0时渤海海冰年度最大期间渤海湾Q1点与老铁山水道Q2点进行对比分析。这两个点在这段时间内均无海冰覆盖, 因此我们可以通过对比考虑和不考虑海冰影响两种情况下波能流密度随时间的演化进而判断海冰是否能通过减小有效风区(也即开阔水域)的占比而使得无冰区域的波能流密度减小。

在考虑海冰影响的情况下, 位于渤海湾Q1点的波能流密度在2012年2月5日0时至2012年2月8日0时的平均值为2.5 kW/m, 比不考虑海冰影响时减小了24%, 最大值为8.4 kW/m, 比不考虑海冰影响时减小了18%, 如图 8a所示。老铁山水道Q2点在考虑海冰影响下的波能流密度平均值为4.3 kW/m, 比不考虑海冰影响时减小了19%; 最大值为15.4 kW/m, 比不考虑海冰影响时减小了14%, 如图 8b所示。其中, 从2月5日0时至2月6日7时期间平均风速约为4.0 m/s, 海冰对波能流密度的影响几乎可以忽略。在2月6日8时至2月8日0时期间, 平均风速为10.6 m/s, 在考虑海冰影响的情况下渤海湾Q1点的平均波能流密度为4.0 kW/m, 比不考虑海冰影响时减小了25%。老铁山水道Q2点的平均波能流密度为7.8 kW/m, 比不考虑海冰影响时减小了20%。由此可见在强风作用下海冰的消浪效果是十分显著的。

与图 8a的渤海湾Q1点相比, 老铁山水道的波能流密度几乎是它的两倍, 即使是在海冰分布年度最大的时刻, 老铁山水道的波能流密度仍然更适合用于波能发电。

4 结语

根据模型验证结果可知, 本文建立的海冰模型-水动力学模型-波浪模型的实时耦合模式能够较好地模拟渤海冬季寒潮大风期间的海冰及波能流密度的分布。根据模拟结果可知, 寒潮大风期间渤海约1/3的海域被海冰覆盖, 然而其平均波能流密度受冰影响的区域占整个渤海的76%, 辽东湾、渤海湾和莱州湾的平均波能流密度受冰影响最多分别减小了100%, 60%和50%。

老铁山水道虽然冬季不结冰, 但在寒潮大风期间, 波能流密度的平均值仍然比不考虑海冰影响时减小了19%。这说明海水结冰不仅导致结冰区域的有效波高和波能流密度会受到海冰影响而减小, 也会导致不结冰海域的有效风区(开阔水域面积)减小, 从而造成不结冰海域的波能流密度减小。此外, 老铁山水道的波能流密度几乎是渤海湾口的两倍, 即使是在海冰严重的时期, 其波能流密度仍然符合发电需求。