1.   引言

随着经济的快速发展，环渤海近岸分布有水电站、核电站、市政供水系统等取水工程，然而冬季取水时经常遭遇巨厚冰体（冰厚大于等于3 m、直径为5 m左右的大圆盘）的威胁^[1-2]。巨厚冰体是由堆积冰冻结而成，在风、潮、流的共同作用下，会导致取水口堵塞，降低冷源取水效率从而影响电站正常运行，因此也称为高风险冰块^[3-4]。在2004−2008年全球取水堵塞事件高达61起，其中绝大多数事件由碎冰、高风险冰块堵塞取水口引起^[5]。除此之外，在极端寒潮天气条件下高风险冰块可能会冲击取水构筑物，引发结构性失效^[1]，从而对冷源取水安全造成影响。高风险冰块已成为冷源取水工程的重要安全隐患，因此建立数值模型提出高风险冰块追踪预测技术可为环渤海及其他寒区海域的冷源取水动态风险预警提供技术支撑。

以往多采用现场观测的手段对定点区域^[1]的冰情进行监测，但无法定量指导预警工作。数值模拟是冰情预警工作的一个必备手段。现有追踪取水口浮冰的工作多是采用离散元方法^[6]针对取水口附近几百米范围内的冰块堆积和碰撞开展研究，忽略了取水口附近的地形水深变化、潮位以及极端天气条件对取水的影响^[7]。然而，高风险冰块在风、潮和流的共同作用下，其运移轨迹模拟是该致险模式下海冰预警的重点^[8-10]。而高风险冰块在24 h内的运移轨迹模拟需要至少10³ km²以上的模型范围，这也就导致现有的浮冰追踪数值模拟方法存在局限性。

因此本文拟基于潮汐模型、海冰模型和拉格朗日粒子追踪模型提出一种高风险冰块追踪预测技术，并以辽东湾某取水口为例探讨高风险冰块运动对冷源取水的影响，为冷源取水安全风险评估提供一种切实可行的研究手段。

2.   研究方法

本文基于非结构化网格有限体积法的海洋模型^[11]（Finite-Volume Community Ocean Model，FVCOM）建立海冰–海洋耦合数值模型，并通过修改拉格朗日粒子追踪模型读取耦合模拟的冰速驱动粒子运动实现对高风险冰块运移轨迹的模拟，其数值模型耦合原理及变量关系传递如图1所示。

图  1  数值模型耦合原理

Figure  1.  Coupling principle of numerical model

下载: 全尺寸图片 幻灯片

潮汐模型是基于FVCOM的动量方程、连续性方程以及状态方程建立的，基本的控制方程及离散方法参考文献[12]。在水平方向使用Smagorinsky公式^[13]进行计算，垂向使用Mellor-Yamada 2.5阶湍流闭合模型^[14]进行计算。通过潮汐模型的计算可以得到潮位、流速等水动力学要素。

将潮汐模型计算的潮位和海表面流速数据传输到海冰模型，循环利用海冰动力学方程、海冰热力学方程和海冰传输方程计算海冰的生消演化，得到海冰密集度、冰厚和冰速数据。其中海冰动力学方程^[15]为

式中，$\sigma $为内部应力张量；m代表单位面积冰的质量之和；$g$代表重力加速度；$f$为科氏力；${u_{\rm{i}}}$和${v_{\rm{i}}}$分别代表$x$方向和$y$方向的冰速；$\zeta $是海表面高程；$ {\sigma _{ij}} $是内部应力张量，下标1、2和$i$、$j$分别代表$x$（下标1）和$y$（下标2）方向上的分量；${A_{\rm{i}}}$为0～1之间的冰密集度；${\tau _{x,{\rm{w}}}}$和${\tau _{x,{\rm{i}}}}$分别为$x$方向上的海表面风应力和冰应力；${\tau _{y,{\rm{w}}}}$和${\tau _{y,{\rm{i}}}}$为$y$方向上的海表面风应力和冰应力^[16]。海表面风应力和冰应力的表达式分别如式（2）和式（3）所示：

式中，${C_{{{\rm{w}}}}}$为风阻力系数，${C_{\rm{w}}} = (0.8 + 0.065 \times \sqrt {{u^2}_{\rm{w}} + {v^2}_{\rm{w}}}) \times {10^{ - 3}}$^[17]；${C_{\rm{i}}}$为冰阻力系数，取值为0.007 5^[18]；${u_{\rm{w}}}$、${v_{\rm{w}}}$分别为$x$、$y$方向风速分量；${\rho _{\rm{a}}} $、${\rho _{\rm{w}}}$分别为大气密度与水密度；$\varphi $为风应力矢量与风向之间的角度；${u_1}$、${v_1}$分别为$x$、$y$方向流速分量；$\omega $为冰流矢量和水之间的角度；$ \Delta u $为冰与流之间的相对流速：

式中，${u_{\rm{i}}}$和${v_{\rm{i}}}$分别代表$x$方向和$y$方向的冰速；${u_1}$、${v_1}$分别为$x$、$y$方向流速分量。

将潮汐模型计算出的流速分量${u_1}$、${v_1}$与海冰模型计算得到的海冰密集度${A_{\rm{i}}}$用于耦合求解冰–水共存海域的海表面应力方程^[16]：

式中，海表面应力分量${\tau _{{\rm{s}},x}}$、${\tau _{{\rm{s}},y}}$和海表面高程（$\zeta $）共同更新下一个时刻潮汐模型中的流速分量，从而实现海冰和潮汐模型的耦合。将耦合模拟的冰速替换掉式（6）原始读取的流速变量，采用冰速作为动力条件驱动高风险冰块粒子运动，垂向速度设置为0，修正后的求解方程如下：

式中，$\vec x$为第$t$时刻的粒子位置；${\rm{d}}\vec x /{\rm{d}}t$为粒子位置随时间的变化率；${\vec u _{\rm{i}}}\left[ {\vec {x} \left( t \right),t} \right]$为潮汐模型与海冰模型耦合模拟的冰速矢量场。

由于冰是浮在水面上的，其上表面和下表面分别受到风和流的作用，采用耦合模拟的冰速场代替拉格朗日粒子追踪模型原始读取的流速场来驱动粒子运动，不仅可以合理模拟风和流驱动下的粒子运移轨迹，而且可以描述多个粒子之间的堆叠和碰撞特性。最终得到高风险冰块在风、潮和流作用下的运移轨迹，从而实现高风险冰块的追踪模拟。

3.   FVCOM冰–海耦合模型验证

本文是在渤海模型的基础上，通过对工程区进行局部加密从而实现取水口附近高风险冰块的追踪数值模拟，最小网格为5 m，取水口位置见图2。渤海模型的模型范围参见文献[19]。模型验证包括潮位、流速、冰层、冰速4个部分的验证。冰层、冰面积验证资料通过决策树方法反演MODIS卫星遥感图像获得^[20]。冰速验证资料来源于JZ20-2平台^[21]。潮位和流速验证资料来源于现场观测数据（图2）。

图  2  渤海水深（单位：m）与测点图

P1（40.30°N，122.08°E）为潮位测站，P2（40.75°N，121.20°E）为流速测站，P3（40.47°N，121.33°E）为JZ20-2平台冰速测站

Figure  2.  The water depth (unit: m) and measured stations in the Bohai Sea

P1 (40.30°N, 122.08°E) is the tidal level station, P2 (40.75°N, 121.20°E) is the tidal current station, and P3 (40.47°N, 121.33°E) is the ice velocity station of JZ20-2 platform

下载: 全尺寸图片 幻灯片

3.1   潮位验证

潮位的数据验证时间为2010年1月1日至31日。根据图3，P1测点模拟与实测的最大值相对误差为5.70%，平均值相对误差为9.40%，验证结果表明实测潮位与模拟潮位值吻合较好。

图  3  测点P1（40.30°N，122.08°E）潮位验证图

Figure  3.  Tidal level verification of Station P1 (40.30°N, 122.08°E)

下载: 全尺寸图片 幻灯片

3.2   流速验证

流速验证的时间取自2009年11月12日12时至13日12时。

根据图4，P2测点的实测最大流速为0.39 m/s，模拟值与实测值的相对误差为0.60%；实测平均流速为0.26 m/s，模拟值与实测值的相对误差为11.20%。流向的模拟值与实测值均以50°（NE）和220°（SW）为主。验证结果表明，模拟的流速与实测结果吻合较好。

图  4  测点P2（40.75°N，121.20°E）流速模拟值与实测值对比

Figure  4.  The simulated and measured current velocity of Station P2 (40.75°N, 121.20°E)

下载: 全尺寸图片 幻灯片

3.3   海冰分布验证

海冰分布的验证时间为2010年1月13日，模拟的海冰面积为30.15 ×10³ km²，占整个渤海面积的30.90%，比MODIS反演的海冰面积偏小3.98%。其中辽东湾海冰面积约为18.01×10³ km²，比MODIS反演的海冰面积偏小6.23%。渤海湾海冰面积约为6.53×10³ km²，比MODIS反演的海冰面积偏小2.33%。莱州湾海冰面积约为4.41×10³ km²，比MODIS反演的海冰面积偏小0.75%。

验证结果表明，辽东湾、渤海湾和莱州湾数值模拟的海冰分布与卫星遥感图像基本一致，如图5所示。

图  5  2010年1月13日数值模拟的渤海海冰分布（橙色线围成的区域）与MODIS卫星遥感图像对比

Figure  5.  Comparison between simulated ice cover surrounded by orange lines of the Bohai Sea and MODIS satellite remote sensing image on January 13, 2010

下载: 全尺寸图片 幻灯片

3.4   冰速验证

冰速验证时间为2004年1月22日23:00至25日00:00，如图6所示。JZ20-2平台附近U方向实测最大冰速为0.78 m/s，模拟值与实测值的相对误差为0.50%；V方向实测最大冰速为0.61 m/s，模拟值与实测值的相对误差为–7%。

图  6  JZ20-2（P3）平台附近海域u_i和v_i方向冰速模拟值与实测值对比

Figure  6.  The simulated and measured values of ice velocity in the direction of u_i and v_i in the sea area near JZ20-2 (P3) platform

下载: 全尺寸图片 幻灯片

模拟冰速的平均值为0.49 m/s，与实测值的相对误差为1.50%。冰向模拟值以16°和330°为主，与实测基本一致。验证结果表明，冰速的模拟结果与实测吻合较好。

4.   高风险冰块的运移轨迹

在辽东湾某冷源取水口附近结冰海域随机释放了1 200个拉格朗日示踪粒子，作为高风险冰块，模拟了2010年1月13日0时至16日0时寒潮期间高风险冰块24 h、48 h和72 h的运移轨迹及进入取水口的概率。本文通过将随机释放的粒子分为8个区域，并用不同颜色标识区分从而实现对高风险冰块的追踪与溯源，当粒子进入取水口分界线时，粒子将被统计为进入到取水口内，见图7。

图  7  取水口位置判定线及初始时刻粒子分布

Figure  7.  Determination line of water intake inlet and particle distribution at initial time

下载: 全尺寸图片 幻灯片

4.1   取水口外粒子的运移轨迹

根据图8，取水口外随机释放的1200个高风险粒子在0～24 h（红线）的运移范围是47.86 km²，24～48 h（蓝线）的运移范围是125.13 km²，48～72 h（黑线）的运移范围是138.68 km²。在此期间，美国国家环境预报中心（National Centers for Environmental Prediction，NCEP）观测^[22]风速最大值为11.70 m/s，平均风速为5.88 m/s，风向以NNE向为主；最大海表面流速为1.27 m/s，平均海表面流速为0.51 m/s；模拟冰速最大值为1.27 m/s，平均冰速为0.53 m/s，冰向和海表面流向基本一致。高风险冰块的运移轨迹基本呈现NE–SW走向，与海表面流速、流向基本一致，如图9所示。

图  8  0～24 h（a）、24～48 h（b）和48～72 h（c）所有释放粒子的移动轨迹

Figure  8.  Trajectories of all released particles in 0−24 h (a), 24−48 h (b) and 48−72 h (c)

下载: 全尺寸图片 幻灯片

图  9  2010年1月13日0时至16日0时取水口外的风、流、冰情况

Figure  9.  Wind, current and ice velocity outside the water intake from 0:00 on January 13 to 0:00 on January 16, 2010

下载: 全尺寸图片 幻灯片

4.2   取水口内粒子的运移轨迹

根据表1，当不考虑取水流量时，取水口内的流速较小，0～72 h内海表面流速最大值和平均值分别为0.09 m/s和0.04 m/s，冰速最大值和平均值分别为0.14 m/s和0.01 m/s。当考虑取水流量时，因取水泵房上沿位于海表面以下5.70 m，取水口内的海表面流速和冰速的最大值不超过0.15 m/s，平均值均不超过0.05 m/s，与不考虑取水流量时相差不大。从图10中粒子在取水口内0～72 h的移动轨迹可以发现，第0～24 h粒子在NNE向风和落潮流的作用下从取水口外进入到取水口内，第24～48 h在取水口内15.39×10³ m²范围内移动，粒子在第48～72 h几乎不再移动。

表  1  模拟时段取水口外、内流速和冰速情况（单位： m/s）

Table  1.  Current and ice velocity outside and inside the water intake during the period of simulation (unit: m/s)

类型	海表面流速 最大值	海表面流速 平均值	冰速 最大值	冰速 平均值
取水口外	1.27	0.51	1.27	0.53
取水口内	0.09	0.04	0.14	0.01
取水口内（考虑取水流量）	0.15	0.05	0.15	0.01

下载: 导出CSV 

| 显示表格

图  10  0~24 h（红线）、24~48 h（蓝线）和48~72 h（黑线）取水口内追踪粒子的移动轨迹

Figure  10.  Trajectories of particles in the water intake in 0–24 h (red line), 24–48 h (blue line) and 48–72 h (black line)

下载: 全尺寸图片 幻灯片

5.   高风险冰块进入取水口的概率

5.1   1 200个初始释放粒子进入取水口的百分比

统计2010年1月13日至16日进入取水口内的粒子百分比（图11）可以发现，在粒子释放的第1个小时内，取水口内的粒子百分比增长很快，从初始时刻的1.80%增长到了粒子总数的5.60%。

图  11  释放不同数量粒子进入取水口内的百分比

Figure  11.  Percentage of particles entering the water intake released in different quantities at the initial time

下载: 全尺寸图片 幻灯片

粒子释放的第10个小时，在自东北向西南方向的落潮流和WNW向风作用下，取水口内形成一逆时针涡旋（图12）。有一部分粒子被涡旋流场带出了取水口，由于流向与风向在这段时间内不同向，冰向在风和流的相互作用下不断变化，因此粒子在30 min内又重新回到取水口内。

图  12  取水口附近粒子释放第10个小时的流场图

Figure  12.  The current field near the water intake at the 10th hour after the release of particles

下载: 全尺寸图片 幻灯片

在粒子释放17 h 55 min时，取水口外的风向与流向、冰向基本一致为WSW–SW向之间，大量粒子在自西南向东北的涨潮流作用下突破了取水口分界线进入到取水口内，此时的取水口内粒子占比为14.62%，达到了模拟时刻的最大值。仅仅5 min后，这些粒子又在风和流的作用下被带出取水口，此时取水口内粒子数量占初始释放总数的5.60%，在此之后取水口内粒子总数未发生变化。观察高风险冰块粒子来源可以发现进入取水口内的粒子仅来自于最靠近取水口的3个区域：区域5、区域7和区域8（图7）。

当考虑取水流量时，进入取水口内的粒子百分比为5.93%，较没考虑取水流量时仅仅增加了5.89%，如图11所示。

5.2   释放粒子数量对结果的影响

为了探究进入取水口内的粒子百分比是否受初始释放粒子数量的影响，在取水口附近释放了6 000个粒子，粒子初始时刻位置及分区与图7一致，每个分区内的粒子数量为图7的5倍，并与1 200个释放粒子的工况进行对比。初始时刻进入取水口内的粒子占比均为1.80%，最终时刻1 200个释放粒子与6 000个释放粒子在取水口内的粒子占比分别为5.60%和5.68%。

如图11所示，与1 200个释放粒子工况相似，在粒子释放的第1个小时内增长很快，从初始时刻的1.80%迅速增长到了粒子总数的5.68%。在粒子释放的第9 h 55 min至11 h 5 min，位于取水口北部的一部分粒子被风和流带出取水口，在取水口分界线附近不断移动，取水口内的粒子占比两次降至0.90%之后又回到了5.68%，此时取水口内的粒子主要来源于区域5、区域7和区域8。在粒子释放24 h后，取水口内的粒子不再发生移动，并且取水口外的粒子也没有进入，结束时刻的粒子位置见图13。

图  13  6 000个粒子结束时刻位置图

不同颜色代表图7所示的粒子初始位置

Figure  13.  Position of 6 000 particles at end time

Different colors represent the initial position of particles as shown in Fig. 7

下载: 全尺寸图片 幻灯片

通过1 200个粒子和6 000个粒子两种工况可以发现：1 200个释放粒子数量的结果是趋于稳定的，两种工况进入取水口内的粒子数量占比最终相差仅为0.08%。此外，两种工况下最后时刻进入取水口的粒子均来源于区域5、区域7和区域8。

6.   结语

本文通过修改拉格朗日粒子追踪模型读取潮汐模型和海冰模型耦合模拟的冰速，作为动力条件驱动高风险冰块粒子运动，模拟并分析了寒潮天气下高风险冰块在24 h、48 h和72 h内的运移轨迹及其进入取水口的概率，主要结论如下：

（1）由于取水口的走向与流向和风向几乎垂直，因此进入取水口内的高风险冰块粒子数量不超过6%。

（2）进入取水口内的粒子均来自距取水口分界线300 m范围内，在此范围内可采取一定措施进行防范，保障取水安全。

（3）进入取水口内的流速大大减小，且出现最大海表面流速小于0.10 m/s的涡旋，导致动力条件大大减弱，使得进入取水口内的高风险冰块在24 h后移动速度大大减小，进一步降低了对取水安全的影响。

（4）采用冰–海耦合模拟的冰速驱动粒子运动不仅能够合理描述冰的运移轨迹，而且能够描述冰块之间的堆叠与碰撞特性，而这是无法通过原拉格朗日粒子追踪模型描述的。本研究能够实现高风险冰块粒子的追踪与溯源，为冷源取水安全风险评估提供了一种切实可行的技术手段。