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

(1.國家海洋環(huán)境監(jiān)測中心，遼寧 大連 116023；2.遼寧省海洋水產(chǎn)科學(xué)研究院，遼寧 大連 116023；3.遼寧紅沿河核電有限公司，遼寧 大連 116001；4.國家海洋局大連海洋環(huán)境監(jiān)測中心站，遼寧 大連 116013；5.葫蘆島市生態(tài)環(huán)境保護服務(wù)中心，遼寧 葫蘆島 125000)

目前，我國的濱海核電站建設(shè)[1-3]已邁入了高速發(fā)展的時代，位于遼寧省近岸海域的有葫蘆島徐大堡和瓦房店紅沿河兩座核電站[4-5]。隨著越來越多的核電機組投入運行，以往人們大多關(guān)注核電站的溫排水對附近海域生物的影響，而漂流海藻對核電站正常運行的潛在影響相比之下研究較少。近年來，曾發(fā)生過多起生物堵塞核電站冷源取水系統(tǒng)的問題，如國內(nèi)一些濱海核電站曾發(fā)生過因海月水母(Aureliaaurita)、海地瓜(Acaudinamolpadioides)、中國毛蝦(Aceteschinensis)堵塞而致使核電系統(tǒng)異常停堆的事件[6-8]。2014年底，廣西防城港核電站所在的冷源取水海域欽州灣暴發(fā)了球形棕囊藻(Phaeocysticglobosa)赤潮，赤潮藻類聚集嚴重威脅了核電站的安全，這是國內(nèi)首次由藻類引發(fā)的威脅核電站冷源安全的事件[9]。2018年，海南昌江核電站也曾發(fā)生褐藻堵塞冷源取水系統(tǒng)的事件[10]。類似的堵塞風(fēng)險在國外也有過報道，這些事件的發(fā)生都嚴重威脅著濱海核電站的的正常運行。在此背景下，人們迫切意識到預(yù)測漂流海藻堵塞冷源取水系統(tǒng)風(fēng)險的重要性，也逐漸成為核電系統(tǒng)亟待解決的問題。

截至目前，國內(nèi)外諸多學(xué)者做過這方面的研究，但大部分是基于常規(guī)監(jiān)測和數(shù)值模擬相結(jié)合的手段進行的。如，Yin 等(2019)利用數(shù)值模型研究發(fā)現(xiàn)海流運動是控制藻類運移的主要因素[11]；羅曉凡等(2012)利用普林斯頓海洋模型(POM)指出潮汐是藻類漂移路徑研究中不可忽略的因素，因此對考慮潮汐運動過程的海洋模式進行質(zhì)點追蹤，進而得出水母漂移路徑的結(jié)論將更為可靠[12]。Johnson等(2005)利用環(huán)流模型對墨西哥灣內(nèi)的漂流海藻漂移路徑進行追蹤，發(fā)現(xiàn)灣內(nèi)環(huán)流的季節(jié)變化對藻類的豐度和分布產(chǎn)生了重要的影響[13]。邱敏志等(2020)采用鐳同位素表觀年齡模型量化紅沿河核電站海域表層水體年齡，繼而分析研究區(qū)域內(nèi)水體的主體流向，進一步結(jié)合藻類豐度的分布結(jié)果，利用鐳同位素示蹤技術(shù)追溯藻類漂移路徑[14]。本研究在瓦房店紅沿河核電站附近海域，采用浮標觀測和數(shù)值模擬相結(jié)合的方法，構(gòu)建該海域的二維水動力和粒子追蹤數(shù)值預(yù)測模型，基于近年實際觀測的海流和潮位數(shù)據(jù)，對水動力模型進行驗證；進而根據(jù)自行完成的浮標現(xiàn)場調(diào)查的信息，對該海域同期相同工況下的海藻漂移路徑[15-17]進行校驗。

最后，在2019年7月大潮期各設(shè)計工況下，于核電站兩個藻類發(fā)源地分別投放示蹤粒子，模擬預(yù)測漂流海藻對冷源取水口堵塞風(fēng)險，為管理人員準確預(yù)判漂流海藻的運動軌跡，適時采取有效措施規(guī)避和減輕漂流海藻對取水口的威脅，提供定量化的數(shù)據(jù)支撐，為保障核電站冷源系統(tǒng)的安全穩(wěn)定運行提供技術(shù)方法。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)組成及其來源

經(jīng)過2015—2018年的實際調(diào)查監(jiān)測，發(fā)現(xiàn)核電站南、北兩側(cè)的大咀子和江石底藻場為冷源取水的藻類風(fēng)險源位置，海藻生物量峰值出現(xiàn)在每年的5月份達800噸左右，主要物種為孔石莼(Ulvapertusa)、長石莼(緣管滸苔)(Ulvalinza)、鼠尾藻(Sargassumthunbergii)等。為了調(diào)查研究漂流海藻漂移路徑，2019年，共投放漂流浮標3次，分別是5月26日于江石底投放1次，6月19日和8月28日于大咀子外海投放2次。通過投放表面漂流浮標對海藻漂移路徑進行研究，進而預(yù)測不同物理海洋工況下漂流海藻的漂移路徑。觀測結(jié)果顯示，2019年5月，在潮流的往復(fù)作用和6～9級北風(fēng)條件下，取水口北側(cè)漂流海藻從江石底漂移到取水口，約需29 h，直線距離約19.8 km。6月份，南風(fēng)1～4級條件下，從大咀子海藻發(fā)源地經(jīng)過8 h即可漂移到取水口，直線距離約8.9 km。漂流浮標及海上漂浮狀態(tài)示意圖見圖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 水動力環(huán)境影響模式基本控制方程 水動力場的準確模擬是漂流海藻漂移路徑模擬預(yù)測的基礎(chǔ)與關(guān)鍵，采用二維水動力模型對研究海域的潮流運動情況進行了模擬。方程如下:

①水動力控制方程

(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為地轉(zhuǎn)科氏力系數(shù)(s-1)。

②邊界條件

以海域周邊岸線和海島岸線作為閉邊界，在閉邊界處法向流速為零。

開邊界處采用如下的分潮調(diào)和形式給定，公式如下：

(4)

式(4)中：ωi是第i個分潮的角速度，fi、u是第i個分潮的交點因子和遲角修正，Hi和gi是調(diào)和常數(shù)，分別為各分潮的振幅和遲角，V0是分潮的時角。

③初始條件

計算開始時“冷態(tài)”啟動，水位和流速初始值均為0。

④風(fēng)場條件

由于本項目取得的風(fēng)場數(shù)據(jù)的時間間隔為6 h，并且海上風(fēng)況變化復(fù)雜，呈現(xiàn)短周期波動的季風(fēng)特征，所以實際計算中要求給定詳實的風(fēng)場條件，在本次模擬中，根據(jù)已有的風(fēng)場資料，分析實施研究期間海域的主要風(fēng)向，據(jù)此對數(shù)值模擬過程給以風(fēng)場條件，足以滿足模型插值計算的需要。

⑤水平渦粘系數(shù)

采用考慮亞尺度網(wǎng)格效應(yīng)的Smagorinsky(1963)公式計算水平渦粘系數(shù)[18]，公式如下：

(5)

1.2.2 模擬區(qū)域及計算參數(shù) 網(wǎng)格是離散的基礎(chǔ)，采用非結(jié)構(gòu)化三角網(wǎng)格剖分計算域(圖略)。為了準確模擬本研究海域的流場分布情況，將地形較為復(fù)雜的區(qū)域進行局部網(wǎng)格加密，整體計算域由9 815個節(jié)點和18 337個三角單元組成，其中最小計算網(wǎng)格面積為60 000 m2，最大計算網(wǎng)格面積為8 400 000 m2。模型計算時間步長根據(jù)CFL條件進行調(diào)整，確保模型計算收斂，最小計算時間步長取0.1 s。底床糙率通過曼寧數(shù)進行控制，曼寧數(shù)取為32～40 m1/3/s。

模型運行時考慮漫灘過程，采用限制水深法和保持水量守恒相結(jié)合的方法來實現(xiàn)網(wǎng)格點有干出的動邊界處理。

2 結(jié)果與討論

2.1 水動力計算結(jié)果

在研究海域選取了1個潮位和4個潮流驗證點，采用2014年7月大、小潮兩個連續(xù)歷時的實測潮位、潮流數(shù)據(jù)，作為模型驗證資料。實測點位置見圖3。驗證結(jié)果見圖4、5所示。

圖3 實測站位示意圖Fig.3 Locations of measuring stations

圖4 潮位驗證結(jié)果Fig.4 Verification results of tide level

圖5 大潮流速、流向驗證Fig.5 Verification of flow velocity and direction during spring tide

通過模型驗證情況可以看出，模擬計算值與實測值變化趨勢基本一致，吻合較好，說明該模型能較好的再現(xiàn)研究海域的實際水動力狀況。同時，利用MIKE21后處理模塊給出取水口附近海域典型時刻的流場分布，如圖6所示。

圖6 大潮流場分布圖Fig.6 Flow field distribution during spring tide

遼東灣東部海域核電站取水口周圍海域潮流為渤海旋轉(zhuǎn)潮流系統(tǒng)的一部分，此處范圍較小，海流空間差異較小。取水口周邊海域主流向為NE—SW，流速約為1.0 m/s，潮流流速達到漲急和落急最大值時，潮汐達半潮面位置，可以看出該海域潮汐主要呈現(xiàn)往復(fù)流運動特征，潮汐呈駐波特征。同時，海域流場受岸線引起的凈余流影響較大，海流基本沿著岸線方向運動，在取水口附近的部分小灣內(nèi)顯現(xiàn)典型的渦旋特征，這也許就是漂流海藻大量聚集引起冷源堵塞風(fēng)險的主要原因。

以上述獲得的水動力場為背景，采用粒子追蹤模型對漂流海藻的漂移路徑進行模擬，與表面漂流浮標采集路徑進行對比，分析不同情景模式工況下漂流海藻的漂移路徑。粒子追蹤通過Langevin方程描述粒子的輸運和擴散狀態(tài)。其方程微分形式如下：

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

(6)

式(6)中：a為漂移項，這里取為流速；b為擴散項，保守型粒子可不考慮。ξ為隨機數(shù)，取值為[0,1]。

2.2 粒子追蹤計算結(jié)果

2.2.1 浮標投放信息 經(jīng)過調(diào)查發(fā)現(xiàn)，取水口附近海域漂流海藻發(fā)源地主要為取水口北側(cè)的江石底和南側(cè)的大咀子附近海域，為了調(diào)查研究漂流海藻漂移路徑，2019年截至目前，于兩個發(fā)源地共投放漂流浮標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 浮標漂移路徑模擬結(jié)果與討論 由于漂流海藻基本上是隨水動力被動漂移，所以其在海水中的漂移過程可以看作是質(zhì)點在表層海流和海面風(fēng)的共同作用下的漂移輸運過程。為了準確模擬表面粒子隨表層流和風(fēng)海流的漂移運動情況，本研究將計算得到的風(fēng)生流與二維流場疊加，其中風(fēng)漂移系數(shù)的大小與海面上方10 m處風(fēng)速大小成正比，風(fēng)漂移系數(shù)的取值為0.01～0.06，進而求得風(fēng)漂移矢量和風(fēng)生流。漂流海藻漂移運動路徑模擬結(jié)果如圖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的模擬結(jié)果可知，發(fā)源于江石底的漂流海藻在潮流的作用下不斷往復(fù)運動，受SSW風(fēng)的影響，浮標總體的趨勢朝NNE方向運動，6月1日以后風(fēng)向發(fā)生改變，浮標在NWW向風(fēng)的驅(qū)使下，運動向左偏轉(zhuǎn)90°左右，朝著NWW方向繼續(xù)運動。分析取得的氣象資料發(fā)現(xiàn)，在5月26日至6月4日期間，此處海域受氣候影響，出現(xiàn)大風(fēng)天氣，最大風(fēng)速約為20.6 m/s，在模型中加入相應(yīng)的風(fēng)場件后，模擬結(jié)果與實測運動軌跡基本一致。

由編號1298929浮標的原模擬結(jié)果可知(圖9)，發(fā)源于大咀子的漂流海藻，在S風(fēng)的作用下，隨潮流一直向N向運動，最后到達大凌河入?？诟浇ＤM結(jié)果與浮標實測運動路徑大部分吻合，但在7月上旬模擬路徑與實測路徑偏差較大，通過實測值發(fā)現(xiàn)在此期間浮標運動路徑向NNE方向發(fā)生偏轉(zhuǎn)，到7月中旬，浮標恢復(fù)N向運動。分析原因應(yīng)該是7月上旬研究海域風(fēng)向發(fā)生變化，由S風(fēng)轉(zhuǎn)為SSW風(fēng)，造成浮標向NNE方向運動，由于此次研究所疊加的風(fēng)場資料略顯粗糙，無法滿足模型精確計算的需求。故而在后來的模型設(shè)置中，根據(jù)研究海區(qū)氣候條件，將風(fēng)場設(shè)定為S風(fēng)-SW風(fēng)-SE風(fēng)，最后又變?yōu)镾向風(fēng)的非定常條件，發(fā)現(xiàn)模擬路徑與實測運動路徑部分時段稍有偏差，但總體趨勢較好(圖9)，說明本研究方法能夠基本準確反映該浮標整體的運移趨勢。

由編號1299908浮標，即圖10的模擬結(jié)果可知，發(fā)源于大咀子的漂流海藻在運動的過程中，受N風(fēng)的影響，不斷向S運動，最終遷移到了長興島附近海域，模擬路徑與浮標實測運動路徑大致吻合，能夠反映該浮標的整體運動趨勢。根據(jù)以上模擬結(jié)果可以發(fā)現(xiàn)，在疊加合適的風(fēng)場后，模擬結(jié)果與實測路徑基本符合，說明浮標在海面漂移的過程中，海面風(fēng)的影響占據(jù)主導(dǎo)作用。在海藻的漂移運動中，由于海藻基本懸浮于表層海水，其單位質(zhì)量較小，更容易受海面風(fēng)的影響，因此，在對海藻漂移進行預(yù)測時，需要搜集更為詳細的風(fēng)場資料，以確保模型預(yù)測的準確性。

2.2.3 取水口海藻堵塞風(fēng)險模擬結(jié)果與討論 為了分析江石底、大咀子附近海域漂流海藻對取水口造成風(fēng)險的可能性，在江石底、大咀子附近海域表層釋放代表漂流海藻的中性粒子進行預(yù)測模擬，基于運動軌跡分析遼東灣東部海域的漂流海藻對取水口的風(fēng)險威脅程度見圖11。

圖11 漂流海藻釋放點位置Fig.11 Released locations of the drifting sea algae

由于在大潮時，潮流更強，質(zhì)點運動范圍更大，因此選取2019年7月6日—8日的大潮時期進行漂流海藻漂移路徑的模擬。由于潮流運動的周期性特征，即使是同一位置，不同潮時釋放的中性粒子，其運動軌跡也是不同的。所以選取了8個代表時刻，在不同的風(fēng)假設(shè)條件下釋放漂流海藻。其釋放時刻及對應(yīng)潮時見圖12。釋放的8個時刻分別為：2019年7月7日3時、7時、10時、14時、17時、20時、23時、2019年7月8日1時。并以到達取水口口門中心1.5 km半徑范圍內(nèi)定義為威脅區(qū)。限于篇幅，本研究只給出兩張粒子軌跡圖見圖13，為了更好地分析不同潮時發(fā)源于不同位置的浮游生物對取水口的威脅，對不同潮時釋放的中性粒子到達取水口威脅區(qū)的時間進行模擬，結(jié)果見表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 靜風(fēng)狀態(tài)下不同時刻釋放中性粒子的運動軌跡Fig.13 Trajectory diagrams of neutral particles released at different times under static wind condition

表2 靜風(fēng)狀態(tài)下中性粒子從釋放位置到達取水口威脅區(qū)所需時間Tab.2 Drifting time of neutral particles from the released position to the warning area near water intake under static wind condition

續(xù)表

表3 SSE四級風(fēng)況下中性粒子從釋放位置到達取水口威脅區(qū)所需時間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四級風(fēng)況下中性粒子從釋放位置到達取水口威脅區(qū)所需時間Tab.4 Drifting time of neutral particles from the released position to the warning area near water intake under NNE level 4 wind conditions

上述表中的模擬結(jié)果是基于中性粒子追蹤的方式探討了不同風(fēng)況下、不同時刻的漂流海藻25 h內(nèi)對取水口的風(fēng)險大小，對發(fā)源于大咀子和江石底海域的漂流海藻在靜風(fēng)、SSE四級風(fēng)況、NNE四級風(fēng)況下的運動軌跡模擬，分析可知：①發(fā)源于大咀子或江石底海域的漂流海藻在潮流和風(fēng)的影響下，基本上在10 h以內(nèi)能夠到達距取水口口門1.5 km的海域內(nèi)，需要加強關(guān)注，并適時采取有效措施減輕漂流海藻對取水口的威脅。②在SSE四級風(fēng)況和NNE四級風(fēng)況下，漂流海藻不僅受表層流作用，還受到海面風(fēng)直接對漂流海藻施加的作用力，以及海面風(fēng)通過改變表層海水的運動間接對漂流海藻產(chǎn)生的作用力。通過數(shù)值模擬結(jié)果，可以看到，在施加風(fēng)場后，漂流海藻在做往復(fù)運動的基礎(chǔ)上，均有朝著外海運動的趨勢，方向同風(fēng)向。③由于漂流海藻發(fā)源地及取水口都位于岸邊，通過模擬結(jié)果發(fā)現(xiàn)，代表漂移海藻的中性粒子釋放后，由于近岸流的影響，已模擬的不同風(fēng)況下，漂流海藻初次到達取水口風(fēng)險區(qū)所需的時間大致相同。但施加NNE、SSE風(fēng)場后，漂流海藻在取水口風(fēng)險區(qū)停留的時間較靜風(fēng)狀態(tài)下有所減少，基本上一個漲落潮周期內(nèi)，在潮流和風(fēng)的共同作用下離開取水風(fēng)險區(qū)，朝著外海運動。

3 結(jié)論

本研究主要是聚焦于紅沿河核電站海域漂流海藻漂移運動路徑的研究。首先建立了遼東灣海域的二維水動力模型和粒子追蹤模型，基于近年觀測的水動力實測數(shù)據(jù)和浮標信息對模型進行調(diào)參和驗證，獲得了符合該區(qū)域客觀狀況的數(shù)值模型。研究結(jié)果顯示，2019年5月，在6～9級北風(fēng)和水動力條件下，漂流海藻從冷源取水口北側(cè)的江石底發(fā)源地漂流到取水口，約需29 h，直線距離約19.8 km；6月份，在南風(fēng)1～4級條件下，從大咀子海藻發(fā)源地經(jīng)過8 h即可漂流到取水口，直線距離約8.9 km。然后，在大潮條件下，利用該模型探討了不同工況下，漂流海藻25 h內(nèi)對取水口堵塞風(fēng)險的大小，結(jié)果發(fā)現(xiàn)，發(fā)源于大咀子或江石底海域的漂流海藻，基本上10 h以內(nèi)能夠到達距核電站取水口口門半徑1.5 km的威脅區(qū)海域范圍內(nèi)，最短時間僅為2 h，應(yīng)該引起足夠重視。今后的研究中，可在模型中給入更為詳實的平面風(fēng)場數(shù)據(jù)，或建立風(fēng)場模型[19-24]，使模型能盡量真實客觀的反應(yīng)自然條件，更加準確的模擬核電站附近海域漂流海藻的漂移運動軌跡，保障核電冷源取水系統(tǒng)安全，保證核電產(chǎn)業(yè)的社會安全。同時，本研究方法是將大型藻類的運動當(dāng)作完全隨水動力的被動漂移擴散。對于有一定主動性的水生生物在海里的運動，應(yīng)同時考慮生物的自身生活習(xí)性的影響，詳細了解抽象為粒子的水生物物種本身的屬性，來盡可能客觀實際的反應(yīng)物種原有的習(xí)性和狀態(tài)，此方面在以后的程序設(shè)計中需要進一步的補充和完善。