平潭竹嶼灣水交換能力和溢油風險數值模擬

張振偉，劉必勁*，曾銀東，張世民

(1. 廈門理工學院土木工程與建筑學院，福建 廈門 361024; 2. 福建省海洋預報臺，福建 福州 350503；3. 國家海洋局廈門海洋預報臺，福建 廈門 361008)

平潭竹嶼灣(圖1)位于平潭綜合實驗區中部，北靠霞嶼村鵝頭山(臨近城市核心區)，東接平潭舊城區，南與敖東鎮及南寨山接壤，西南接竹嶼出海口。國家相關部門支持在平潭竹嶼灣開展生態整治和修復工作，工作內容主要包括海灣深度清淤疏浚、海綿污染防控系統建設、竹嶼灣水環境監控能力建設等。生態整治和修復工作實施后，可以有效提高平潭防災減災能力、改善竹嶼灣水環境等。本次數模主要針對竹嶼灣水交換能力和溢油風險影響范圍進行研究，為竹嶼灣的水環境監控能力建設及管理提供支撐。

有限體積海岸海洋模型(Finite Volume Coast Ocean Model, FVCOM)模型是應用有限體積法對守恒形式的三維水動力方程進行離散而構建的海洋動力計算模型[1-2]。模型基于非結構三角形網格進行離散，這種網格可以較好的擬合不規則水域邊界。許多學者開展了海域水體交換和溢油等方面的研究。趙亮等(2002)應用水動力模型(Estuary, Coast Ocean Model, ECOM)基于水質點運移的物理交換，針對膠州灣的水交換能力開展了計算研究[3]。李真(2009)應用FVCOM模型對圍墾引起的羅源灣的水動力環境進行了模擬研究[4]。王興剛等(2015)開展了連云港港主體港區水交換能力的研究[5]。遲萬清等(2018)研究了不同風向作用下膠州灣水體交換規律[6]。王雪等(2017)應用FVCOM模型研究了膠州灣各個子區域之間的水交換情況[7]。劉曉東等(2017)應用FVCOM模型開展了廈門圍頭灣海域溢油風險計算分析[8]。

1 材料與方法

1.1 FVCOM模型

為了提高淺水區域的垂向分辨率，FVCOM模型垂向坐標應用了σ坐標變換。FVCOM模型中σ坐標有多種形式可以選擇，本次模擬中我們采用了垂向均勻分布的σ坐標。FVCOM模型的動量方程和連續方程表達為：

(1)

(2)

(3)

水平耗散項表達為：

(4)

(5)

(6)

式(1)至(6)中：x、y、σ為三維笛卡爾直角坐標系的正東(m)、正北(m)和垂直方向坐標；u、v、ω為x、y、σ方向上的速度分量(m/s)；ζ為自由水面(m)；D為總水深(m)；H為靜水水深(m)；q為湍動能(m/s)；l為湍流宏觀尺度(m)；Fx、Fy分別為動量x方向、動量y方向的水平擴散項(m/s2)；ρ為水體密度(kg/m3)；f為科氏力系數(s-1)； g為重力加速度(m/s2)；Km為垂向渦粘系數(m2/s)；Kh為垂向熱擴散系數(m2/s)；Am為水平渦粘系數(m2/s)，Ah為水平熱擴散系數(m2/s)；它們可以由修正的MY-2.5湍流閉合子模型計算[9-10]。

圖1 平潭竹嶼灣及其鄰近海域水深地形和計算網格Fig. 1 Bathymetry of Zhuyu Bay and its adjacent waters and grids of the model

1.2 污染擴散模型

FVCOM模型中污染擴散方程表達為：

(7)

式(7)中：C為物質濃度，KH為垂向渦粘系數；Fc為水平擴散項；C0為點源濃度。

(8)

式(8)中：ts、te分別為物質釋放的開始和結束時間；i為三角網格節點編號，N為物質釋放節點的總數。

水交換能力反映了水體的自凈能力，與水環境、水生態密切相關。評價水體交換能力的指標有很多，如平均滯留時間[11]、半交換時間[12]、沖洗時間等。平均滯留時間是指水質點離開研究區所需的平均時間。半交換時間是指保守物質濃度降至初始濃度一半時所需時間。沖洗時間是指研究區內保守物質濃度降為零時所需時間。假定t0時研究區域內某保守物質的濃度為C0，而后某時刻t(t>t0)存留于研究區域內的物質濃度變為C(t)，保守物質平均滯留時間(θ)定義為：

(9)

式(9)中：當C(t)=50%C0時，Δt=t-t0即為水體的半交換時間；當C(t)=0時，Δt=t-t0即為水體的沖洗時間。

1.3 粒子追蹤模型

拉格朗日粒子追蹤模型，可以模擬粒子在風、波浪、海流等作用下的漂移運動軌跡，假設粒子不會發生碰撞，也不會發生混合，則粒子擴散方程寫為：

(10)

u=uc+κ·uw+ur，v=vc+κ·vw+vr

(11)

式(10)、(11)中：原坐標為(x0，y0)的粒子經時間Δt=t-t0后，漂移到坐標(x,y)；u和v分別是粒子運動的東、北分量，它由流速uc、風速uw、粒子隨機運動速度ur組成，κ為風對粒子拖曳系數，取0.025。通過跟蹤各粒子坐標(x,y)的各位置，確定運移范圍，統計分析粒子掃海面積。

2 結果與討論

2.1 潮汐潮流計算結果

本研究的計算區域范圍為24.82°—25.90°N，119.10°—120.23°E。為保證模擬精度，海域水深根據中華人民共和國海事局出版的最新海圖及竹嶼灣附近海域實測水深插值得到。模型的水深地形和網格如圖1所示。在開邊界處模型計算網格分辨率約為4 km，在近海岸處網格進行了加密，竹嶼灣海域網格分辨率約為30 m左右。模型垂直方向采用均勻分布的符號坐標，為與實測層數一致垂向分為7個σ層，摩擦系數選用ORIG選項進行計算。

模型開邊界的潮位時間序列從OTIS全球潮汐模型 (OSU Tidal Inversion Software)數據庫中提取，在生成潮位時間序列時選擇了M2、S2、K2、N2、K1、O1、P1、Q1八個主要分潮。

模型模擬了時間自9月1日—14日的潮汐潮流運動情況，該時間范圍包含了大潮和小潮的觀測期限，小潮觀測時間(UTC)9月4日03時—5日03時。大潮觀測時間(UTC)自9月12日03時—13日03時，水文調查站位見圖2。圖3給出了平潭海洋站大潮測流期間潮位模擬與實測值的對比，由圖可知潮位時間序列計算結果與實測結果符合較好。

圖2 平潭竹嶼灣及其鄰近海域水文調查站點位Fig. 2 Locations of tidal current observation stations in Zhuyu Bay and its adjacent waters, Pingtan

圖3 平潭竹嶼灣2018年9月潮位觀測結果與模擬結果對比Fig. 3 Comparison of observed and modeled water elevations of Zhuyu Bay in September 2018

圖4 大潮和小潮流速、流向驗證結果Fig. 4 Modeled current validation against observations during spring and neap tides

圖4給出了大潮流速、流向和小潮流速、流向驗證結果，由圖中結果可以看出，計算流速、流向總體上與實測值吻合較好。海壇海峽內潮波屬于前進波，潮位越高流速越大。高潮位時海峽內海流自北向南流，流速大；落潮半潮位時海流在分流尾嶼區域形成南北分流，分流尾嶼北側海流向北流，而南側海流則向南流，海峽中部流速較小；低潮位時海峽內海流自南向北流，水道上流速較大。因潮差大，海峽兩岸露灘面積大。自低潮再經過約3 h，至漲潮半潮位時，海峽北東口海流自北向南流與海峽北上漲潮海流在分流尾嶼區域匯合，海峽中部流速較小[13]。潮位潮流驗證結果表明本研究建立的數值模型，能夠很好的模擬研究區內潮位、潮流運動特征，可以應用模型開展區域水環境特征的計算研究。

圖5給出了大潮高潮竹嶼灣內潮流場的矢量圖。由圖可以看出，竹嶼灣內的潮流流速在大潮高潮落潮半潮位和高潮漲潮半潮位時的流速最大，流速超過1.0 m/s，說明潮流動力很強。在大潮高潮漲平和落平時刻流速小，最小流速量值在0.3 m/s左右。因竹嶼灣呈狹長狀態，分析竹嶼灣中特征點的流速流向過程線可以看出潮流呈往復流特征，潮汐呈駐波特征。

圖5 平潭竹嶼灣及其鄰近海域大潮水深平均潮流場Fig. 5 Tidal current velocity vector and elevation in Zhuyu Bay and its adjacent waters

2.2 水體交換計算結果

利用建立的竹嶼灣水動力模型，構建竹嶼灣對流擴散模型，進行竹嶼灣水交換能力指標的計算。計算時以竹嶼灣口為界，在灣內釋放可溶于水的保守物質，灣外物質濃度設為0.0 mg/L，灣內物質濃度設為1.0 mg/L，釋放時間選為小潮高潮時刻(UTC時間9月3日08時)。計算水半交換時間是對于計算域內每個網格，當網格節點的濃度第一次小于0.5 mg/L時，記下該時刻，這一時刻與保守物質釋放時刻的差即為水半交換時間。對于平均滯留時間，需要按照式(9)進行積分計算，對于每個網格節點，記錄其保守物質第一次變為0.0 mg/L的時刻，積分的時間為保守物質釋放時刻至保守物質變為0.0 mg/L的時刻，這樣可以得到每個網格節點的平均滯留時刻。對于沖洗時間，每個網格節點也是不一樣的，這里沒有針對每一個節點開展計算，而是取竹嶼灣交換能力最差的位置處的沖洗時間作為整個海灣的沖洗時間。圖6、7給出了竹嶼灣水體半交換時間和平均滯留時間的分布情況，由圖可知，竹嶼灣大部分水域水體半交換時間小于1.0 d，平均滯留時間3.0 d左右。竹嶼灣全部水體的沖洗時間為15.0 d。總體上，竹嶼灣水體交換能力較強。

Liu 等(2004)針對不同釋放時間對平均滯留時間造成的誤差做了分析[14]，其分析結果證明平均滯留時間與初始濃度無關，但保守物質的釋放時刻不同，對平均滯留時間的計算結果會造成影響，其影響范圍在5%～10%之間。

2.3 溢油擴散計算結果

竹嶼灣溢油事故風險主要來自于海域的運輸船只，假定運輸船只發生溢油事故，3 h內漏油5 t，模擬泄露油品為燃料油。溢油點選在竹嶼灣灣口航道區域。采用瞬時排放的方式，計算時不考慮粘結，也不考慮揮發，模擬溢油點位于竹嶼灣口水道，坐標為(25.515°N，119.697°E)，油粒子位置每隔10 min輸出一次。

溢油數模旨在了解一旦發生溢油事故，溢油的擴散分布趨勢。本次溢油模擬主要模擬靜風和不利風向(NNE風向，風速9.2 m/s)條件下，高平潮、低平潮、漲急、落急等典型潮時發生溢油時，溢油泄漏48 h內各典型時刻的影響范圍及到達附近環境敏感點的時間。海壇海峽及附近區域敏感目標包括海壇島北部保留區、福清灣農漁業區、山岐澳中華鱟保護區、塘嶼列島海洋保護區等。

圖6 竹嶼灣水體半交換時間分布Fig. 6 Distribution of half-life time in Zhuyu Bay

圖7 竹嶼灣平均滯留時間分布Fig. 7 Distribution of average residence time in Zhuyu Bay

表1給出了靜風和不利風時刻，油粒子掃海面積的統計結果，油粒子在某一單元內出現，則該單元的面積作為油粒子掃海面積的貢獻值。圖8給出了不利風條件下，油粒子48 h擴散范圍圖。

表1 平潭竹嶼灣溢油掃海面積Tab. 1 Area of oil film after oil spill over 48 h in Zhuyu Bay

圖8 不利風向下竹嶼灣及其鄰近海域4個不同潮時刻溢油48 h軌跡分布Fig. 8 Distribution of oil particles trajectory after oil spill over 48h at 4 tide moments under the averse wind conditions in Zhuyu Bay and its adjacent waters

計算分析表明靜風條件下，油粒子總體上與潮流運動趨勢一致。漲急、高平潮和低潮情況下，油粒子主要呈現出在海壇海峽內的往復運動，總體上油粒子在靠近海壇海峽平潭一側運動。落急時刻，油粒子一旦從竹嶼灣口水道運動至海壇海峽水道，油粒子會迅速向海壇海峽南口區域運動，其擴散范圍達到了平潭壇南灣、平潭塘嶼島、平潭草嶼島、高山灣等區域，其面積可達40.800 km2。

由圖8可知，不利風對油粒子擴散影響顯著。不利風條件下，油粒子在NNE風向的作用下，顯現出向海壇海峽福清一側運動的趨勢。油粒子的擴散范圍可以到達平潭壇南灣、平潭塘嶼島、平潭草嶼島、高山灣等區域。在漲急、高平潮和落急時刻，油粒子的掃海面積均超過了20 km2，落急時刻掃海面積最大，達到了43.100 km2。低平潮時刻釋放的油粒子擴散范圍最小，其掃海面積為6.880 km2。因此，不利風對油粒子的漂移擴散影響顯著，風的作用擴大了油粒子的影響范圍。

3 結論

本研究應用FVCOM海洋動力模型開展了平潭竹嶼灣水體交換和溢油擴散計算分析。通過在小潮高潮時刻釋放示蹤劑，對竹嶼灣水體交換能力開展了研究，計算結果表明竹嶼灣大部分水域的水半交換時間小于1.0 d，最大3.5 d；平均滯留時間大部分區域小于3.0 d，最大4.0 d；水體沖洗時間為15.0 d。

溢油擴散計算結果表明，靜風條件下落急時刻發生溢油時48 h油粒子掃海面積最大，達到了40.800 km2。油粒子的擴散范圍到達了平潭壇南灣、塘嶼島、草嶼島及高山灣等區域，對周圍海域影響大。其他時刻發生溢油時油粒子主要在海壇海峽平潭一側往復活動。不利風(NNE向)作用下，溢油粒子掃海面積明顯比靜風條件大，不同時刻發生溢油時油粒子均可以擴散到平潭壇南灣、塘嶼島、草嶼島及高山灣附近海域，低平潮掃海面積最小6.880 km2，落急時刻發生溢時油粒子掃海面積最大，為43.100 km2。其他時刻發生溢油時，油粒子掃海面積均大于24.200 km2。