核電站附近海域漂流海藻的漂移路徑數值模擬

關春江，王 昆，田新華，徐 鵬，陸宏亮，劉永健，趙行行，孔 芳

(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]，使模型能盡量真實客觀的反應自然條件，更加準確的模擬核電站附近海域漂流海藻的漂移運動軌跡，保障核電冷源取水系統安全，保證核電產業的社會安全。同時，本研究方法是將大型藻類的運動當作完全隨水動力的被動漂移擴散。對于有一定主動性的水生生物在海里的運動，應同時考慮生物的自身生活習性的影響，詳細了解抽象為粒子的水生物物種本身的屬性，來盡可能客觀實際的反應物種原有的習性和狀態，此方面在以后的程序設計中需要進一步的補充和完善。