英曉明，嚴金輝，趙明利

（國家海洋局南海規(guī)劃與環(huán)境研究院，廣東廣州 510300）

0 引言

惠州大亞灣位于粵港澳大灣區(qū)的惠州市，西鄰大鵬灣，東接紅海灣，屬于半封閉海灣（見圖1）。大亞灣灣內(nèi)潮波因受狹長海灣和島嶼地形反射及變形影響，淺水分潮效應較大[1]。根據(jù)惠州海洋站2006—2019 年的長期驗潮資料，該站平均潮差為0.85 m，最大潮差為1.25 m，灣內(nèi)為不正規(guī)半日潮。大亞灣海域年平均氣壓為1 012.0 hPa，多年平均氣溫為22.9 ℃，累年平均降水量為1 734.0 mm，累年平均風速為3.4 m/s，常風向為東南向，強風向為東北向。實測資料表明，2019 年3 月純洲島與灣頂之間海域的垂向平均最大流速為0.17 ～0.18 m/s；大亞灣中部海域垂向平均最大流速為0.18～0.34 m/s。該海域鵝洲站2018 年夏季和冬季浮標觀測資料記錄，夏季有效波高平均值為0.41 m，平均波周期為4.1 s；冬季有效波高平均值為0.23 m，平均波周期為3.5 s。大亞灣內(nèi)島嶼眾多，化工企業(yè)聚集，港口碼頭密集，是粵港澳大灣區(qū)重要的海灣。大亞灣海域海洋生物資源豐富，屬于環(huán)境敏感區(qū)，有大亞灣水產(chǎn)資源自然保護區(qū)、珊瑚、重要產(chǎn)卵場、無居民海島、紅樹林、淡澳溪河口等環(huán)境敏感保護目標。大亞灣是許多名貴魚類的繁殖、產(chǎn)卵和育肥場所[2]。廣東省海洋功能區(qū)劃顯示大亞灣還分布著多個旅游休閑娛樂區(qū)。

圖1 研究區(qū)域和水文觀測站位置示意圖Fig.1 Location of the study area and observation station

隨著中海殼牌、?？松梨诨葜菀蚁┑缺姸嗍椖康穆鋺?，大亞灣過往船舶密度增大，溢油事故風險增加。大亞灣海域一旦發(fā)生船舶溢油事故，將嚴重影響其環(huán)境敏感保護目標并破壞海洋生態(tài)環(huán)境。1970—2018 年登陸粵港澳大灣區(qū)的熱帶氣旋共有40 個，年均為0.8 個[3]，臺風風暴潮發(fā)生頻率較高。2017 年臺風“苗柏”風暴潮期間，大亞灣多艘船舶走錨；2018 年臺風“山竹”風暴潮期間，大亞灣大辣甲西側(cè)海域發(fā)生船舶走錨事故，所幸均沒有發(fā)生環(huán)境污染事故。風暴潮期間因發(fā)生船舶事故導致破壞海洋環(huán)境安全的風險是潛在的，應用數(shù)學模型研究風暴潮期間的船舶溢油事故能夠預測污染物擴散范圍，為溢油事故應急處置提供參考。已有學者對大亞灣溢油和?；沸孤妒鹿曙L險進行了模擬研究[1-5]。然而，已有研究是針對常規(guī)天氣下環(huán)境污染事故風險進行模擬，而極端天氣風暴潮期間溢油事故的模擬研究較少。本文通過數(shù)值模擬方式，對大亞灣海域風暴潮期間的溢油事故進行情景模擬，為溢油事故應急決策提供技術(shù)支撐。

1 數(shù)學模型

1.1 水動力模型

對于南海地區(qū)臺風風暴潮的數(shù)值模擬，已有學者采用數(shù)學模型MIKE21、高級環(huán)流模型（Advanced Circulation Model，ADCIRC）、Delft 3D 三維水動力-水質(zhì)模型系統(tǒng)和非結(jié)構(gòu)網(wǎng)格海洋環(huán)流與生態(tài)模型（Finite-Volume Coastal Ocean Model，F(xiàn)VCOM）等開展了風暴增水和海浪等研究[6-11]。本文選擇MIKE21 數(shù)值模式進行風暴潮與天文潮耦合情況下的水動力預測模擬，模型網(wǎng)格及模擬范圍見圖2。模型共有35 835個節(jié)點，63 770個網(wǎng)格單元，最小網(wǎng)格邊長約為17 m。

圖2 模型網(wǎng)格和水深地形Fig.2 Model mesh and bathymetry

大亞灣海域具有顯著的潮位“雙峰”現(xiàn)象和漲潮流速“雙峰”現(xiàn)象。研究表明，當潮波由灣口傳播至灣頂時，上述現(xiàn)象主要由M2、M4和M6分潮共同作用造成[12]，M6分潮的異常增大是引起大亞灣內(nèi)部潮位“雙峰”現(xiàn)象的主要原因，2MS6分潮對大亞灣內(nèi)部潮汐的扭曲也具有一定作用[13]。多站位實測潮位調(diào)和分析研究結(jié)果表明，淺水分潮（M4、M6、MS4和2MS6）振幅的增長是造成“雙峰”現(xiàn)象在灣內(nèi)增強的主要原因[14]。模型水位邊界條件采用的M2、S2、N2、K2、K1、O1、P1、Q1、M4、MS4分潮調(diào)和常數(shù)初始值引自全球潮汐模型[15]，此外增加了M6和2MS6兩個分潮，這兩個分潮的調(diào)和常數(shù)初始值取自南海潮流模型結(jié)果[16]。海面風風場及氣壓場采用Holland 模型。臺風氣壓模型如下[9,17]：

式中：Pc為臺風中心氣壓；Pn為臺風外圍氣壓；RMW為最大風速半徑；r為距離臺風中心的距離；Vmax為臺風最大風速；lat 為臺風中心緯度；B為Holland 參數(shù)。臺風對水動力的影響通過海水表面摩擦系數(shù)來體現(xiàn)，摩擦系數(shù)取值范圍為0.000 5～0.001 7。

1.2 溢油模型

溢油事故發(fā)生后，油膜主要在風和海流共同驅(qū)動下發(fā)生平移運動、擴展以及蒸發(fā)、乳化、溶解等風化過程。本文采用MIKE21 HD 模塊構(gòu)建平面二維水動力模型，在此基礎(chǔ)上耦合MIKE21 OS模塊建立溢油模型。關(guān)于模型中的風漂移系數(shù)，一些學者取值0.03～0.04[18-19]，另一些學者取值0.02～0.03[20-22]，也有學者取值0.025～0.044[23]。風漂移系數(shù)和風速有關(guān)，當風速較大時風漂移系數(shù)取大值，因此本文取0.044。

2 模型驗證和流場分析

本文選擇2019年3月20—21日項目附近海域2個潮位觀測站和3個潮流觀測站的觀測值進行常規(guī)情況下的潮位和潮流驗證，選擇2016年7月30日17時（北京時，下同）—8 月2 日23 時赤灣、珠海、深圳和惠州4 個驗潮站的潮位值以及2017 年8 月22 日00 時—24 日05 時赤灣、大萬山、深圳和惠州4 個驗潮站的潮位值進行臺風風暴潮期間潮位和風速、風向的驗證。觀測站位見圖1。風暴潮期間潮位觀測值引自文獻[7]。

2.1 模型驗證和流場特征分析

2019年3 月潮位和潮流模型驗證結(jié)果分別見圖3 和圖4。由圖3 可見，潮位模擬值和觀測值基本一致，T1 站和惠州站的平均潮位誤差分別為10%和-3%，潮位驗證結(jié)果較好，模型能夠較好地模擬大亞灣海域的潮位“雙峰”現(xiàn)象。由圖4 可見，潮流模擬值和觀測值變化趨勢大體一致，V1、V2 和V3 站的平均流速誤差分別為2%、13%和14%，平均流向誤差分別為3%、-6%和44%。潮流驗證結(jié)果總體較好，模型能夠反映大亞灣海域的潮流運動特征。

圖3 2019年3月20日10時—21日12時潮位驗證結(jié)果Fig.3 Tidal level verification results from 10:00 on March 20 to 12:00 on March 21,2019

圖4 2019年3月20日10時—21日12時潮流驗證結(jié)果Fig.4 Tidal current verification results from 10:00 on March 20 to 12:00 on March 21,2019

大亞灣海域漲急時刻和落急時刻的流場分布見圖5。無論是漲急還是落急時刻，大亞灣中部海域和范和灣灣口的流速較大。漲急時刻的平均流速約為0.23 m/s；落急時刻平均流速約為0.18 m/s。漲急時刻，潮流主要由東南向進入大亞灣，灣內(nèi)中央主流向北；落急時刻，灣內(nèi)中央主流向南，潮流主要由東南向流出大亞灣。

圖5 大亞灣海域漲急和落急時刻流場分布Fig.5 Distribution of tide current velocity field in Daya Bay at flood tide and ebb tide

2.2 風暴潮驗證

本文采用1713號臺風“天鴿”風暴潮的風速和風向過程對模式進行風的參數(shù)驗證，驗證結(jié)果見圖6，由圖可見，模擬的風速和風向與實測值趨勢相同，整體吻合較好。1604 號臺風“妮妲”和1713 號臺風“天鴿”風暴潮期間的潮位過程驗證結(jié)果見圖7，由圖可見，在兩次臺風風暴潮期間，潮位模擬值和觀測值基本接近，模型能夠反映臺風風暴潮的潮位變化過程。

圖6 1713號臺風“天鴿”風暴潮期間風速過程比較結(jié)果Fig.6 Comparison results of wind speed during storm surge of 1713 Typhoon"Hato"

圖7 1604號臺風“妮妲”和1713號臺風“天鴿”風暴潮期間潮位過程比較結(jié)果Fig.7 Comparison results of tidal level during storm surge of 1604 Typhoon"Nida"and 1713 Typhoon"Hato"

3 結(jié)果分析

由于近年來大亞灣海域沒有發(fā)生典型等級以上的涉海溢油事故案例，這里僅進行情景事故模擬。假設(shè)臺風風暴潮期間發(fā)生特別重大的船舶污染事故，參照《防治船舶污染海洋環(huán)境管理條例》，船舶發(fā)生溢油情景事故的溢油量為1 000 t。漏油時間為1 h，在模型中溢油由240 000個“油粒子”表示。溢油方式選擇固定點源連續(xù)性溢油，油品選擇柴油，密度取0.89×103kg/m3。

3.1 1604號臺風“妮妲”風暴潮溢油情景模擬結(jié)果

1604 號臺風“妮妲”風暴潮模擬時間為2016 年7 月29 日20 時—8 月3 日05 時。在臺風行進過程中，假設(shè)在8月2日凌晨（臺風中心臨近大亞灣前3 h）大辣甲西側(cè)海域發(fā)生溢油事故。風暴潮過程期間溢油掃海面積、殘油量和油膜中心位置見表1，風暴潮中止后油膜掃海面積見圖8。

表1 1604號臺風“妮妲”風暴潮過程期間溢油掃海面積、殘油量和油膜中心位置Tab.1 Oil spill sweeping area,residual oil volume and oil film center position during the storm surge process of 1604 Typhoon "Nida"

圖8 1604號臺風“妮妲”風暴潮過程期間溢油29 h后油膜掃海范圍Fig.8 Scope of oil film sweeping after 29 h oil spill during the storm surge of 1604 Typhoon"Nida"

事故發(fā)生后，溢油將立即影響到附近珊瑚區(qū)、大亞灣水產(chǎn)資源自然保護區(qū)中部核心區(qū)、中部緩沖區(qū)和大辣甲無居民海島等4 個敏感目標，2 h 后影響到大亞灣水產(chǎn)資源自然保護區(qū)南部緩沖區(qū)，3 h后影響到大亞灣水產(chǎn)資源自然保護區(qū)南部核心區(qū)，8 h后影響到小辣甲附近的珊瑚，10 h 后影響到許洲無居民海島和重要產(chǎn)卵場，11 h 后影響到大亞灣水產(chǎn)資源自然保護區(qū)北部實驗區(qū)和沙魚洲無居民海島，12 h后影響到沙魚洲東側(cè)的珊瑚，13 h 后影響到大亞灣灣頂岸線，14 h 后影響到鵝洲和赤洲附近的珊瑚，16 h后影響到馬鞭洲附近海域，18 h后影響到寶塔洲無居民海島。臺風“妮妲”風暴潮中止后（溢油后29 h）油膜掃海范圍約為256.38 km2，殘油量為640.10 t。

3.2 1713號臺風“天鴿”風暴潮溢油情景模擬結(jié)果

1713 號臺風“天鴿”風暴潮模擬時間為2017 年8 月21 日02 時—24 日17 時。在臺風行進過程中，假設(shè)在8 月23 日06 時（臺風中心臨近大亞灣南側(cè)）平海電廠碼頭附近海域發(fā)生溢油事故。風暴潮過程期間溢油掃海面積、殘油量和油膜中心位置見表2，風暴潮中止后油膜掃海面積見圖9。

表2 1713號臺風“天鴿”風暴潮過程期間溢油掃海面積、殘油量和油膜中心位置Tab.2 Oil spill sweeping area,residual oil volume and oil film center position during the storm surge process of 1713 Typhoon"Hato"

圖9 1713號臺風“天鴿”風暴潮過程期間溢油35 h后油膜掃海范圍Fig.9 Scope of oil film sweeping after 35 h oil spill during the storm surge of 1713 Typhoon"Hato"

事故發(fā)生后，溢油將在1 h 后影響到大亞灣水產(chǎn)資源自然保護區(qū)北部實驗區(qū)，2 h后影響到大亞灣水產(chǎn)資源自然保護區(qū)中部緩沖區(qū)和桑洲無居民海島，3 h后影響到大亞灣水產(chǎn)資源自然保護區(qū)中部核心區(qū)和大辣甲無居民海島，4 h后影響到大辣甲附近的珊瑚，5 h后影響到大亞灣水產(chǎn)資源自然保護區(qū)中部緩沖區(qū)，7 h 后影響到大亞灣核電站北側(cè)岸線，8 h后影響到小辣甲附近的珊瑚，9 h后影響到小辣甲無居民海島，13 h后影響到重要產(chǎn)卵場、赤洲附近的無居民海島和珊瑚，14 h 后再次影響到大亞灣水產(chǎn)資源自然保護區(qū)北部實驗區(qū)，15 h 后影響到沙魚洲和許洲無居民海島，17 h后影響到馬鞭洲北側(cè)的珊瑚，18 h 后影響到淡澳溪河口和鵝洲附近的珊瑚，20 h后影響到大亞灣灣頂岸線。風暴潮中止后（溢油后35 h）油膜掃海范圍約為115.62 km2，殘油量為626.94 t。

3.3 流速變化分析

本節(jié)通過分析臺風風暴潮期間和無臺風時的流速變化探討臺風風暴潮的影響。由圖10 可見，1604 號臺風“妮妲”風暴潮導致大亞灣內(nèi)大部分區(qū)域流速增加，尤其在灣口區(qū)域、馬鞭洲西側(cè)區(qū)域和大亞灣灣頂區(qū)域，灣口區(qū)域流速最大增加值超過0.9 m/s，馬鞭洲西側(cè)區(qū)域流速最大增加值超過0.8 m/s，大亞灣灣頂區(qū)域流速最大增加值超過0.5 m/s。1713 號臺風“天鴿”風暴潮導致大亞灣內(nèi)大部分區(qū)域流速增加，尤其在灣口區(qū)域、大亞灣核電站附近區(qū)域和大亞灣灣頂區(qū)域，灣口區(qū)域流速最大增加值超過1.0 m/s，大亞灣核電站附近區(qū)域超過0.8 m/s，大亞灣灣頂區(qū)域超過0.5 m/s。臺風風暴潮期間流速的增大加速了油膜的漂移和擴散。

圖10 1604號臺風“妮妲”和1713號臺風“天鴿”風暴潮過程期間和無臺風的流場和流速變化Fig.10 Comparison of current fields and current velocity change among the storm surge process of 1604 Typhoon"Nida",1713 Typhoon"Hato"and the scenario without typhoon

4 結(jié)論和建議

大亞灣海域島嶼眾多，潮汐具有獨特的“雙峰”現(xiàn)象。本文通過MIKE21 數(shù)值模式，建立了高分辨率的惠州大亞灣附近海域的風暴潮和溢油數(shù)學模型，通過引入M4、MS4、M6和2MS6淺水分潮，精確模擬了大亞灣潮位的“雙峰”特征。通過Holland 模型對1713號臺風“天鴿”風暴潮的風場進行了模擬，驗證結(jié)果較好，運用模型模擬分析了在1604 號臺風“妮妲”和1713 號臺風“天鴿”兩場典型臺風風暴潮期間，大亞灣海域發(fā)生1 000 t 特大溢油情景事故對敏感目標的影響。

在1604號臺風“妮妲”風暴潮期間，假設(shè)大辣甲西側(cè)海域發(fā)生溢油事故，油膜最大擴散范圍約為256.38 km2，溢油發(fā)生后，3 h 內(nèi)將影響附近珊瑚區(qū)、大亞灣水產(chǎn)資源自然保護區(qū)中部核心區(qū)、大辣甲無居民海島、中部緩沖區(qū)、大亞灣水產(chǎn)資源自然保護區(qū)南部緩沖區(qū)和大亞灣水產(chǎn)資源自然保護區(qū)南部核心區(qū)等敏感目標，臺風風暴潮中止后，殘油為溢油量的64.0%；在1713 號臺風“天鴿”風暴潮期間，假設(shè)平海電廠碼頭附近海域發(fā)生溢油事故，油膜最大擴散范圍約為115.62 km2，溢油發(fā)生后，3 h 內(nèi)將影響附近大亞灣水產(chǎn)資源自然保護區(qū)北部實驗區(qū)、大亞灣水產(chǎn)資源自然保護區(qū)中部緩沖區(qū)、桑洲無居民海島、大亞灣水產(chǎn)資源自然保護區(qū)中部核心區(qū)和大辣甲無居民海島等敏感目標，臺風風暴潮中止后，殘油為溢油量的62.7%。本文建立的溢油模型可以應用于常規(guī)天氣和臺風風暴潮情況下的溢油擴散模擬，能夠為溢油事故應急決策提供技術(shù)支撐。

本文所采用的模型沒有考慮大洋環(huán)流對近岸區(qū)域的影響，同時由于大亞灣溢油事故較少，缺少驗證資料，因此沒有深入分析溶解、擴散等參數(shù)對溢油的影響，這是本文存在的不足，后續(xù)將根據(jù)資料情況進一步開展研究。

對于溢油事故風險防控，應做好應急預案，充分準備溢油防控應急物資，平時做好應急演練；加強船舶管控，在臺風風暴潮來臨前做好防臺準備；在臺風風暴潮期間若發(fā)生溢油事故，除借助衛(wèi)星遙感監(jiān)測技術(shù)外，還可以通過數(shù)值模擬及時預測溢油擴散范圍和到達敏感目標的時間；在臺風風暴潮過后，及時安排船舶、人員利用除油物資開展油膜打撈等應急處置。