核电站附近海域漂流海藻的漂移路径数值模拟

关春江，王 昆，田新华，徐 鹏，陆宏亮，刘永健，赵行行，孔 芳

(1.国家海洋环境监测中心，辽宁 大连 116023；2.辽宁省海洋水产科学研究院，辽宁 大连 116023；3.辽宁红沿河核电有限公司，辽宁 大连 116001；4.国家海洋局大连海洋环境监测中心站，辽宁 大连 116013；5.葫芦岛市生态环境保护服务中心，辽宁 葫芦岛 125000)

目前，我国的滨海核电站建设[1-3]已迈入了高速发展的时代，位于辽宁省近岸海域的有葫芦岛徐大堡和瓦房店红沿河两座核电站[4-5]。随着越来越多的核电机组投入运行，以往人们大多关注核电站的温排水对附近海域生物的影响，而漂流海藻对核电站正常运行的潜在影响相比之下研究较少。近年来，曾发生过多起生物堵塞核电站冷源取水系统的问题，如国内一些滨海核电站曾发生过因海月水母(Aureliaaurita)、海地瓜(Acaudinamolpadioides)、中国毛虾(Aceteschinensis)堵塞而致使核电系统异常停堆的事件[6-8]。2014年底，广西防城港核电站所在的冷源取水海域钦州湾暴发了球形棕囊藻(Phaeocysticglobosa)赤潮，赤潮藻类聚集严重威胁了核电站的安全，这是国内首次由藻类引发的威胁核电站冷源安全的事件[9]。2018年，海南昌江核电站也曾发生褐藻堵塞冷源取水系统的事件[10]。类似的堵塞风险在国外也有过报道，这些事件的发生都严重威胁着滨海核电站的的正常运行。在此背景下，人们迫切意识到预测漂流海藻堵塞冷源取水系统风险的重要性，也逐渐成为核电系统亟待解决的问题。

截至目前，国内外诸多学者做过这方面的研究，但大部分是基于常规监测和数值模拟相结合的手段进行的。如，Yin 等(2019)利用数值模型研究发现海流运动是控制藻类运移的主要因素[11]；罗晓凡等(2012)利用普林斯顿海洋模型(POM)指出潮汐是藻类漂移路径研究中不可忽略的因素，因此对考虑潮汐运动过程的海洋模式进行质点追踪，进而得出水母漂移路径的结论将更为可靠[12]。Johnson等(2005)利用环流模型对墨西哥湾内的漂流海藻漂移路径进行追踪，发现湾内环流的季节变化对藻类的丰度和分布产生了重要的影响[13]。邱敏志等(2020)采用镭同位素表观年龄模型量化红沿河核电站海域表层水体年龄，继而分析研究区域内水体的主体流向，进一步结合藻类丰度的分布结果，利用镭同位素示踪技术追溯藻类漂移路径[14]。本研究在瓦房店红沿河核电站附近海域，采用浮标观测和数值模拟相结合的方法，构建该海域的二维水动力和粒子追踪数值预测模型，基于近年实际观测的海流和潮位数据，对水动力模型进行验证；进而根据自行完成的浮标现场调查的信息，对该海域同期相同工况下的海藻漂移路径[15-17]进行校验。

最后，在2019年7月大潮期各设计工况下，于核电站两个藻类发源地分别投放示踪粒子，模拟预测漂流海藻对冷源取水口堵塞风险，为管理人员准确预判漂流海藻的运动轨迹，适时采取有效措施规避和减轻漂流海藻对取水口的威胁，提供定量化的数据支撑，为保障核电站冷源系统的安全稳定运行提供技术方法。

1 数据和方法

1.1 数据组成及其来源

经过2015—2018年的实际调查监测，发现核电站南、北两侧的大咀子和江石底藻场为冷源取水的藻类风险源位置，海藻生物量峰值出现在每年的5月份达800吨左右，主要物种为孔石莼(Ulvapertusa)、长石莼(缘管浒苔)(Ulvalinza)、鼠尾藻(Sargassumthunbergii)等。为了调查研究漂流海藻漂移路径，2019年，共投放漂流浮标3次，分别是5月26日于江石底投放1次，6月19日和8月28日于大咀子外海投放2次。通过投放表面漂流浮标对海藻漂移路径进行研究，进而预测不同物理海洋工况下漂流海藻的漂移路径。观测结果显示，2019年5月，在潮流的往复作用和6～9级北风条件下，取水口北侧漂流海藻从江石底漂移到取水口，约需29 h，直线距离约19.8 km。6月份，南风1～4级条件下，从大咀子海藻发源地经过8 h即可漂移到取水口，直线距离约8.9 km。漂流浮标及海上漂浮状态示意图见图1、2。

图1 海洋表层漂流浮标示意图Fig.1 Schematic diagram of drifting buoy on the sea surface

图2 海洋表层漂流浮标的海上漂移Fig.2 Drifting of the drifting buoy on the sea surface

1.2 研究方法

1.2.1 水动力环境影响模式基本控制方程 水动力场的准确模拟是漂流海藻漂移路径模拟预测的基础与关键，采用二维水动力模型对研究海域的潮流运动情况进行了模拟。方程如下:

①水动力控制方程

(1)

运动方程：

+h·us·S

(2)

+h·vs·S

(3)

式(1)至(3)中：t为时间(h)；η为水位(m)；h=η+d为总水深(m)；d为静水深(m)；u、v为x、y方向的垂向平均流速(m/s)；S为源汇项(kg/s)；pa为大气压力(N/m2)，f为地转科氏力系数(s-1)。

②边界条件

以海域周边岸线和海岛岸线作为闭边界，在闭边界处法向流速为零。

开边界处采用如下的分潮调和形式给定，公式如下：

(4)

式(4)中：ωi是第i个分潮的角速度，fi、u是第i个分潮的交点因子和迟角修正，Hi和gi是调和常数，分别为各分潮的振幅和迟角，V0是分潮的时角。

③初始条件

计算开始时“冷态”启动，水位和流速初始值均为0。

④风场条件

由于本项目取得的风场数据的时间间隔为6 h，并且海上风况变化复杂，呈现短周期波动的季风特征，所以实际计算中要求给定详实的风场条件，在本次模拟中，根据已有的风场资料，分析实施研究期间海域的主要风向，据此对数值模拟过程给以风场条件，足以满足模型插值计算的需要。

⑤水平涡粘系数

采用考虑亚尺度网格效应的Smagorinsky(1963)公式计算水平涡粘系数[18]，公式如下：

(5)

1.2.2 模拟区域及计算参数 网格是离散的基础，采用非结构化三角网格剖分计算域(图略)。为了准确模拟本研究海域的流场分布情况，将地形较为复杂的区域进行局部网格加密，整体计算域由9 815个节点和18 337个三角单元组成，其中最小计算网格面积为60 000 m2，最大计算网格面积为8 400 000 m2。模型计算时间步长根据CFL条件进行调整，确保模型计算收敛，最小计算时间步长取0.1 s。底床糙率通过曼宁数进行控制，曼宁数取为32～40 m1/3/s。

模型运行时考虑漫滩过程，采用限制水深法和保持水量守恒相结合的方法来实现网格点有干出的动边界处理。

2 结果与讨论

2.1 水动力计算结果

在研究海域选取了1个潮位和4个潮流验证点，采用2014年7月大、小潮两个连续历时的实测潮位、潮流数据，作为模型验证资料。实测点位置见图3。验证结果见图4、5所示。

图3 实测站位示意图Fig.3 Locations of measuring stations

图4 潮位验证结果Fig.4 Verification results of tide level

图5 大潮流速、流向验证Fig.5 Verification of flow velocity and direction during spring tide

通过模型验证情况可以看出，模拟计算值与实测值变化趋势基本一致，吻合较好，说明该模型能较好的再现研究海域的实际水动力状况。同时，利用MIKE21后处理模块给出取水口附近海域典型时刻的流场分布，如图6所示。

图6 大潮流场分布图Fig.6 Flow field distribution during spring tide

辽东湾东部海域核电站取水口周围海域潮流为渤海旋转潮流系统的一部分，此处范围较小，海流空间差异较小。取水口周边海域主流向为NE—SW，流速约为1.0 m/s，潮流流速达到涨急和落急最大值时，潮汐达半潮面位置，可以看出该海域潮汐主要呈现往复流运动特征，潮汐呈驻波特征。同时，海域流场受岸线引起的净余流影响较大，海流基本沿着岸线方向运动，在取水口附近的部分小湾内显现典型的涡旋特征，这也许就是漂流海藻大量聚集引起冷源堵塞风险的主要原因。

以上述获得的水动力场为背景，采用粒子追踪模型对漂流海藻的漂移路径进行模拟，与表面漂流浮标采集路径进行对比，分析不同情景模式工况下漂流海藻的漂移路径。粒子追踪通过Langevin方程描述粒子的输运和扩散状态。其方程微分形式如下：

dXt=a(t,Xt)dt+b(t,Xt)ξtdt

(6)

式(6)中：a为漂移项，这里取为流速；b为扩散项，保守型粒子可不考虑。ξ为随机数，取值为[0,1]。

2.2 粒子追踪计算结果

2.2.1 浮标投放信息 经过调查发现，取水口附近海域漂流海藻发源地主要为取水口北侧的江石底和南侧的大咀子附近海域，为了调查研究漂流海藻漂移路径，2019年截至目前，于两个发源地共投放漂流浮标3次，具体浮标漂移历时和投放位置见表1、图7所示。

表1 漂流浮标投放时间及位置Tab.1 Launched time and locations of the drifting buoys

图7 漂流浮标投放点Fig.7 Launched locations of the drifting buoys

2.2.2 浮标漂移路径模拟结果与讨论 由于漂流海藻基本上是随水动力被动漂移，所以其在海水中的漂移过程可以看作是质点在表层海流和海面风的共同作用下的漂移输运过程。为了准确模拟表面粒子随表层流和风海流的漂移运动情况，本研究将计算得到的风生流与二维流场叠加，其中风漂移系数的大小与海面上方10 m处风速大小成正比，风漂移系数的取值为0.01～0.06，进而求得风漂移矢量和风生流。漂流海藻漂移运动路径模拟结果如图8至图10所示。

图8 编号318501浮标运动轨迹Fig.8 Trajectory of No.318501 Buoy2019年5月26日11时至2019年6月23日17时。

图9 编号1298929浮标运动轨迹Fig.9 Trajectory of No.1298929 Buoy2019年6月19日10时至2019年8月1日17时。

图10 编号1299908浮标运动轨迹Fig.10 Trajectory of No.1299908 Buoy2019年8月28日21时至2019年9月3日12时。

由编号318501浮标，即图8的模拟结果可知，发源于江石底的漂流海藻在潮流的作用下不断往复运动，受SSW风的影响，浮标总体的趋势朝NNE方向运动，6月1日以后风向发生改变，浮标在NWW向风的驱使下，运动向左偏转90°左右，朝着NWW方向继续运动。分析取得的气象资料发现，在5月26日至6月4日期间，此处海域受气候影响，出现大风天气，最大风速约为20.6 m/s，在模型中加入相应的风场件后，模拟结果与实测运动轨迹基本一致。

由编号1298929浮标的原模拟结果可知(图9)，发源于大咀子的漂流海藻，在S风的作用下，随潮流一直向N向运动，最后到达大凌河入海口附近。模拟结果与浮标实测运动路径大部分吻合，但在7月上旬模拟路径与实测路径偏差较大，通过实测值发现在此期间浮标运动路径向NNE方向发生偏转，到7月中旬，浮标恢复N向运动。分析原因应该是7月上旬研究海域风向发生变化，由S风转为SSW风，造成浮标向NNE方向运动，由于此次研究所叠加的风场资料略显粗糙，无法满足模型精确计算的需求。故而在后来的模型设置中，根据研究海区气候条件，将风场设定为S风-SW风-SE风，最后又变为S向风的非定常条件，发现模拟路径与实测运动路径部分时段稍有偏差，但总体趋势较好(图9)，说明本研究方法能够基本准确反映该浮标整体的运移趋势。

由编号1299908浮标，即图10的模拟结果可知，发源于大咀子的漂流海藻在运动的过程中，受N风的影响，不断向S运动，最终迁移到了长兴岛附近海域，模拟路径与浮标实测运动路径大致吻合，能够反映该浮标的整体运动趋势。根据以上模拟结果可以发现，在叠加合适的风场后，模拟结果与实测路径基本符合，说明浮标在海面漂移的过程中，海面风的影响占据主导作用。在海藻的漂移运动中，由于海藻基本悬浮于表层海水，其单位质量较小，更容易受海面风的影响，因此，在对海藻漂移进行预测时，需要搜集更为详细的风场资料，以确保模型预测的准确性。

2.2.3 取水口海藻堵塞风险模拟结果与讨论 为了分析江石底、大咀子附近海域漂流海藻对取水口造成风险的可能性，在江石底、大咀子附近海域表层释放代表漂流海藻的中性粒子进行预测模拟，基于运动轨迹分析辽东湾东部海域的漂流海藻对取水口的风险威胁程度见图11。

图11 漂流海藻释放点位置Fig.11 Released locations of the drifting sea algae

由于在大潮时，潮流更强，质点运动范围更大，因此选取2019年7月6日—8日的大潮时期进行漂流海藻漂移路径的模拟。由于潮流运动的周期性特征，即使是同一位置，不同潮时释放的中性粒子，其运动轨迹也是不同的。所以选取了8个代表时刻，在不同的风假设条件下释放漂流海藻。其释放时刻及对应潮时见图12。释放的8个时刻分别为：2019年7月7日3时、7时、10时、14时、17时、20时、23时、2019年7月8日1时。并以到达取水口口门中心1.5 km半径范围内定义为威胁区。限于篇幅，本研究只给出两张粒子轨迹图见图13，为了更好地分析不同潮时发源于不同位置的浮游生物对取水口的威胁，对不同潮时释放的中性粒子到达取水口威胁区的时间进行模拟，结果见表2—4。

图12 漂流海藻释放时刻示意图Fig.12 Schematic diagram of released time for the drifting sea algae ①至⑧分别代表漂流海藻释放的8个典型时刻，①为一个潮汐历程的低低潮时(7月7日3时)，②为涨急时(7月7日7时)，③为高高潮时(7月7日10时)，④为落急时(7月7日14时)，⑤为高低潮时(7月7日17时)，⑥为一个太阴日的另一个涨急时(7月7日20时)，⑦为低高潮时(7月7日23时)，⑧为另一个落急时(7月8日1时)。

图13 静风状态下不同时刻释放中性粒子的运动轨迹Fig.13 Trajectory diagrams of neutral particles released at different times under static wind condition

表2 静风状态下中性粒子从释放位置到达取水口威胁区所需时间Tab.2 Drifting time of neutral particles from the released position to the warning area near water intake under static wind condition

续表

表3 SSE四级风况下中性粒子从释放位置到达取水口威胁区所需时间Tab.3 Drifting time of neutral particles from the released position to the warning area near water intake under SSE level 4 wind conditions

表4 NNE四级风况下中性粒子从释放位置到达取水口威胁区所需时间Tab.4 Drifting time of neutral particles from the released position to the warning area near water intake under NNE level 4 wind conditions

上述表中的模拟结果是基于中性粒子追踪的方式探讨了不同风况下、不同时刻的漂流海藻25 h内对取水口的风险大小，对发源于大咀子和江石底海域的漂流海藻在静风、SSE四级风况、NNE四级风况下的运动轨迹模拟，分析可知：①发源于大咀子或江石底海域的漂流海藻在潮流和风的影响下，基本上在10 h以内能够到达距取水口口门1.5 km的海域内，需要加强关注，并适时采取有效措施减轻漂流海藻对取水口的威胁。②在SSE四级风况和NNE四级风况下，漂流海藻不仅受表层流作用，还受到海面风直接对漂流海藻施加的作用力，以及海面风通过改变表层海水的运动间接对漂流海藻产生的作用力。通过数值模拟结果，可以看到，在施加风场后，漂流海藻在做往复运动的基础上，均有朝着外海运动的趋势，方向同风向。③由于漂流海藻发源地及取水口都位于岸边，通过模拟结果发现，代表漂移海藻的中性粒子释放后，由于近岸流的影响，已模拟的不同风况下，漂流海藻初次到达取水口风险区所需的时间大致相同。但施加NNE、SSE风场后，漂流海藻在取水口风险区停留的时间较静风状态下有所减少，基本上一个涨落潮周期内，在潮流和风的共同作用下离开取水风险区，朝着外海运动。

3 结论

本研究主要是聚焦于红沿河核电站海域漂流海藻漂移运动路径的研究。首先建立了辽东湾海域的二维水动力模型和粒子追踪模型，基于近年观测的水动力实测数据和浮标信息对模型进行调参和验证，获得了符合该区域客观状况的数值模型。研究结果显示，2019年5月，在6～9级北风和水动力条件下，漂流海藻从冷源取水口北侧的江石底发源地漂流到取水口，约需29 h，直线距离约19.8 km；6月份，在南风1～4级条件下，从大咀子海藻发源地经过8 h即可漂流到取水口，直线距离约8.9 km。然后，在大潮条件下，利用该模型探讨了不同工况下，漂流海藻25 h内对取水口堵塞风险的大小，结果发现，发源于大咀子或江石底海域的漂流海藻，基本上10 h以内能够到达距核电站取水口口门半径1.5 km的威胁区海域范围内，最短时间仅为2 h，应该引起足够重视。今后的研究中，可在模型中给入更为详实的平面风场数据，或建立风场模型[19-24]，使模型能尽量真实客观的反应自然条件，更加准确的模拟核电站附近海域漂流海藻的漂移运动轨迹，保障核电冷源取水系统安全，保证核电产业的社会安全。同时，本研究方法是将大型藻类的运动当作完全随水动力的被动漂移扩散。对于有一定主动性的水生生物在海里的运动，应同时考虑生物的自身生活习性的影响，详细了解抽象为粒子的水生物物种本身的属性，来尽可能客观实际的反应物种原有的习性和状态，此方面在以后的程序设计中需要进一步的补充和完善。