黃毅峰 ，許 婷 ，劉 濤

（1.長沙理工大學水利學院，長沙410076；2.交通運輸部天津水運工程科學研究所工程泥沙交通行業(yè)重點實驗室，天津300456）

近半個世紀以來全球發(fā)生了近10起重大的原油泄露事件。這些原油的泄露造成水質污染，海域環(huán)境惡化，危及海洋生物，造成生態(tài)破壞。人民生活也因此受到影響，企業(yè)因此蒙受巨大的損失。今年4月份美國墨西哥灣鉆井平臺發(fā)生爆炸引發(fā)漏油事件引起了國際社會的廣泛關注。

鑒于溢油事件的頻發(fā)性和嚴重性，國內外有很多學者對溢油事件進行了數值模擬研究。國外的比較著名的有美國的OILMAP系統(tǒng)、英國的OSIS系統(tǒng)、挪威的OSCAR系統(tǒng)、荷蘭的MS4系統(tǒng)、比利時的MU-SLICK系統(tǒng)［1］。國內學者在不同的海區(qū)都進行了溢油數值模擬研究，俞濟清、黃立文等在舟山水域建立了溢油預報模型［2-3］；珠江水域、渤海灣水域、遼東灣水域、膠州灣等水域建立了溢油預報模型［4-7］。文中利用MIKE21 HD模塊建立了甌江口二維水動力模型，并用近幾年最新實測資料驗證模型的可靠性和適用性，模擬結果為甌江口航道二維溢油模型的建立提供水動力基礎數據，再次利用MIKE21 SA模塊建立甌江口航道二維溢油模型，應用“油粒子”模型模擬輸移、風化和熱量遷移等過程，對甌江口航道一期工程溢油泄露事故進行影響預測，選取7種具有代表性情形的事故場景，對不同氣象條件和不同事故發(fā)生地點的工況進行了模擬分析。甌江口以外的海域存在著水產養(yǎng)殖、風景旅游、生態(tài)保護等敏感區(qū)域，一旦發(fā)生溢油事故，將對周圍海洋環(huán)境產生嚴重的影響，應采取溢油應急措施，控制污染影響程度，同時追蹤油膜污染帶，并通知港航監(jiān)督部門。

1 工程海域自然動力條件

1.1 潮汐類型

甌江口是浙江省第二大河流，為強潮河口，據2005年6月水文全潮測驗統(tǒng)計，各潮位站（HK1+H01）/HM2=0.24～0.28，均小于0.5，因此潮汐屬正規(guī)半日潮類型。一晝夜兩個潮，潮高不等現(xiàn)象較為明顯［8］。

1.2 潮流特征

據1999年10月～2005年7月甌江河口及其附近海域水文全潮各站潮流速觀測資料調和分析和計算，各站主要全日分潮流與主要半日分潮流的比值F=（W01+WK1）/WM2均小于0.5，各站淺水分潮流比值G==0.27～0.48，綜合該兩項指標，表明該水域潮流屬正規(guī)淺海半日潮流類型；由于漲落潮流受水域地形限制，基本呈往復流運動。

1.3 風況

本地區(qū)全年常風向為N-NE向，頻率為54.7%，全年平均風速為3.8 m/s；強風向為SSW向，最大風速為32.0 m/s（1975年8月12日）。全年中夏季多為SW向大風，春秋季節(jié)多為偏S向或偏N向大風，以偏N向大風為主，冬季盛行N-NE向大風。根據其多年風資料統(tǒng)計，多年平均≥6級風日數為37 d，多年平均≥7級風日數為8.5 d。

1.4 波浪特征

據洞頭島NE端的甲米礁東岸波浪觀測站（1990年12月～1992年9月）資料統(tǒng)計，該海域常浪向為NNE向，頻率為55%，次常浪向為S-SW向，頻率為20%；涌浪向為ENE向，頻率為26.3%，風浪常浪向主要為N-NE向，頻率55%，次為S-SSW向，頻率為16.5%。最強浪向為NNE-ENE向，年平均H4%波高大于0.6 m，對應年平均周期為0.5～4.2 s，多發(fā)生在8、9月臺風侵襲時期，實測最大波高4.3 m（E向），由9216號臺風造成。

據溫州淺灘圍涂工程1#（北堤處）和2#（溫州淺灘南口水道附近）波浪觀測站資料統(tǒng)計，2001～2004年1#站測得波高2.1～3.0 m出現(xiàn)24次（2001年12月～2002年5月14次、2002年6～7月10次），波高大于3.0 m僅出現(xiàn)2次（2001年12月）；2001～2003年2#站波高2.1～3.0 m出現(xiàn)3次（2001年6月2次、2002年9月1次），兩站波高為0.0～0.5 m的出現(xiàn)年頻率分別為95.0%和91.9%。

2 溢油預測模型概述

丹麥水環(huán)境研究所開發(fā)的MIKE 21/3溢油分析（SA）模塊，是基于歐拉-拉格朗日理論體系，通過對油膜在水體中的擴展、傳輸（水流和風場作用）、紊動擴散、分散（夾帶）、蒸發(fā)、乳化、溶解等各種過程的模擬，MIKE 21/3 SA能提供油膜隨時間變化的漂移位置、厚度，以及漂移過程中粘度、油膜表面溫度、傾點等屬性的變化。另外，MIKE 21/3對復雜水體（如各種水工建筑物）的水動力能進行精確模擬，強大的前后處理功能適用于決策系統(tǒng)的開發(fā)。MIKE系統(tǒng)是國際上同類型軟件中最先進的模擬工具之一，在世界各地已有許多成功的應用案例［9］。

2.1 輸移過程

油粒子的輸移包括擴展、漂移、擴散等過程［9］，這些過程是油粒子位置發(fā)生變化的主要原因，而油粒子的組分在這些過程中不發(fā)生變化。

（1）擴展運動。采用修正的Fay重力-粘力公式計算油膜擴展

式中：Aoil為油膜面積，Aoil=πR2oil；Roil為油膜直徑；Ka為系數；t為時間；油膜體積為

初始油膜厚度

（2）漂移運動。油粒子漂移的作用力是水流和風拽力，油粒子總漂移速度由以下權重公式計算

式中：Uw為水面以上10 m處的風速；Us為表面流速；cw為風漂移系數，一般在0.03～0.04。風場數據從氣象部門獲得，而流場從二維水動力模型計算結果獲得。

（3）紊動擴散。假定水平擴散各向同性，一個時間步長內α方向上的可能擴散距離Sα可表示為

式中：［ ］

R1-1為-1～1的隨機數；Dα為α方向上的擴散系數。

2.2 風化過程

油粒子的風化包括蒸發(fā)、溶解和形成乳化物等過程，在這些過程中油粒子的組成發(fā)生改變，但油粒子水平位置沒有變化。

（1）蒸發(fā)。油膜蒸發(fā)受油分、氣溫和水溫、溢油面積、風速、太陽輻射和油膜厚度等因素的影響。假定在油膜內部擴散不受限制（氣溫高于0℃以及油膜厚度低于5～10 cm時基本如此），油膜完全混合，油組分在大氣中的分壓與蒸氣壓相比可忽略不計。

蒸發(fā)率可由下式表示

式中：N為蒸發(fā)率；kei為物質輸移系數；為蒸汽壓；R為氣體常數；T為溫度；M為分子量；ρ為油組分的密度；i為各種油組分。kei由下式估算

式中：k為蒸發(fā)系數；Sci為組分i的蒸氣Schmidts數。

（2）乳化。油向水體中的運動機理包括溶解、擴散、沉淀等。擴散是溢油發(fā)生后最初幾星期內最重要的過程。擴散是一種機械過程，水流的紊動能量將油膜撕裂成油滴，形成水包油的乳化。這些乳化物可以被表面活性劑穩(wěn)定，防止油滴返回到油膜。在惡劣天氣狀況下最主要的擴散作用力是波浪破碎，而在平靜的天氣狀況下最主要的擴散作用力是油膜的伸展壓縮運動。

從油膜擴散到水體中的油分損失量計算

式中：Da為進入到水體的分量；Db為進入到水體后沒有返回的分量

油滴返回油膜的速率為

油中含水率變化可由以下平衡方程表示

式中：R1和R2分別為水的吸收速率和釋出速率。

（3）溶解。溶解率用下式表示

2.3 熱量遷移

蒸氣壓與粘度受溫度影響，而且觀察發(fā)現(xiàn)通常油膜的溫度要高于周圍的大氣和水體。圖1為油膜的熱平衡示意圖。

圖1中1為大氣與油膜之間的傳熱過程；2為大氣與油膜之間熱輻射過程；3為太陽輻射；4為蒸發(fā)熱損失；5為油膜與水體之間的熱量遷移；6為油膜與水體之間散發(fā)和接受的熱輻射。

3 水動力條件模擬

溢油擴散模擬結果是建立在水動力條件模擬的基礎之上，文中采用MIKE21 HD模塊計算水動力條件。

3.1 平面二維水流控制方程

連續(xù)方程

運動方程

式中：x、y為直角坐標系坐標；t為時間；h為水深（基準面到床面的距離）；ζ為潮位（基準面到自由水面的距離）；u、v分別為x、y方向的垂線平均流速分量；f為科氏系數；g為重力加速度；Ex、Ey分別為x、y方向的水平紊動粘性系數；τx、τy分別為波流共同作用下床面剪切力在x、y方向的分量；Sxx、Sxy、Syy分別為各方向的波浪輻射應力。

3.2 計算域的確定及網格剖分

該海域島嶼眾多，潮波運動受地形控制顯著、徑流影響明顯的特點，要復演該水域的潮流運動，必須模擬甌江口整體海域的潮波傳播過程，建立足夠大的計算區(qū)域。因此，此次數學模型的計算區(qū)域西邊界取至120°30′經度線，東邊界至121°45′經度線，北至樂清灣，南至南麂島南側27°10′緯度線。東西向寬約120 km，南北向為150 km。模型采用正四邊形網格，計算網格步長為Dx=Dy=50 m，并就復雜的陸地和島嶼岸線進行二次精細光滑處理，以保證計算精度。

3.3 模型驗證

為了驗證模型的合理性，采用2006年10月甌江口水文全潮測驗（每次測驗均為大、中、小3個潮型）中的大潮資料進行了驗證。水文全潮共布置了13個潮位站、14條水文垂線（圖2）。其中1#～7# 垂線施測大、中、小潮，8#～14# 垂線僅施測大潮，測驗內容為潮位、流速、流向等。該次水文全潮規(guī)模宏大、內容豐富、測點涵蓋了甌江口大范圍海域。部分測站大潮潮位及流速流向驗證曲線見圖3～圖4。

通過潮位、流速、流向的驗證，模型中無論計算的位相、量值還是過程線，均與原型吻合良好，符合《海岸與河口潮流泥沙模擬技術規(guī)程》［10］的要求，可以用來進行溢油擴散數值模擬計算研究。

3.4 工程海域潮流場模擬

采用已驗證的潮流數學模型對工程海域潮流場進行模擬計算。從工程海域的流場看，甌江口及口外的海域漲落潮流主要受岸線和地形控制及工程方案布置的影響，總體呈現(xiàn)往復流運動；由于江口外的島嶼眾多，導致附近區(qū)域的流態(tài)比較復雜。漲潮的時候潮流從外海區(qū)域向甌江口門涌入，受小門島、大門島和洞頭列島的狀元岙島和霓嶼島影響分成兩股主要的漲潮流，其中大門島的南北兩側形成了甌江口潮流的主流和口門的分流區(qū)，漲潮流在流經小門島被分為一南一北兩支，北支流在流向沿岸時受岸線約束又分為南北2個流向，向南漲潮流跟大門島和小門島之間的支流匯合成口門分流。落潮時基本呈現(xiàn)出與漲潮相反的運動過程。當落潮時，外海潮位降低，大小門島的西側及南北側的落潮流向外海運動，落潮流經小門島分為南北兩支流，匯合小門島東側海域，而后跟大門島南側的主流匯合向外海方向運動，至低平潮沿岸形成一定范圍的露灘。工程海域漲、落急流場分布見圖5。

4 溢油影響預測

在水動力條件模擬的基礎上，利用MIKE21 SA模塊建立甌江口航道二維溢油模型，應用“油粒子”模型模擬輸移、風化和熱量遷移等過程，對甌江口航道一期工程溢油泄露事故進行影響預測，預測時刻選擇大潮的4個潮時（高潮時、落急時、低潮時、漲急時），考慮甌江口地區(qū)全年常風向為N－NE向，頻率為54.7%，全年平均風速為3.8 m/s；強風向為SSW向，最大風速為32.0 m/s。溢油計算的結果受水動力、風向、溢油事故發(fā)生位置、發(fā)生時間等影響。綜合考慮以上各影響因素，共計36個計算工況。限于篇幅，這里選取7種代表工況進行結果展示，其掃海范圍及漂移軌跡見圖6～圖12，圖中所示為時間總長2個潮周期（即24 h），時間分別取2 h、6 h、12 h和24 h的結果。不同時間段過后其溢油掃海面積見表1。

表1 事故溢油掃海面積預測結果Tab.1 Sweeping area of oil spill accident prediction

從圖6～圖12和表1中均可以看出，溢油事故發(fā)生后其掃海范圍及遷移路徑與溢油發(fā)生地點、發(fā)生時間、風向情況密切相關，不同情況下，其掃海范圍及遷移路徑差別較大。根據施工期和營運期溢油風險預測可得：（1）不同風向的影響時。風對溢油油膜漂移軌跡的影響比較大。靜風條件下，油膜的漂移路徑主要受控制于潮流動力，在NE向風的作用下，油膜會發(fā)生明顯向西北移動現(xiàn)象，而在SSW風向下，油膜同樣會發(fā)生明顯向東北至東南該區(qū)間方向偏轉；（2）不同潮時的影響。溢油的漂移軌跡受溢油時刻的影響極大，在潮周期的不同階段發(fā)生溢油事故時，溢油油膜漂移軌跡可能完全不同。可以通過掃海范圍及油膜漂移軌跡圖明顯看出不同潮時所帶來的不同溢油事故后果。若發(fā)生溢油事故的地點不同，不同潮時所帶來的影響差別大小也不盡相同。

5 結論

利用MIKE21 HD模塊，建立甌江口水域水動力數學模型模擬潮流特性。在此基礎上利用MIKE 21/3溢油分析（SA）模塊，對甌江口航道二期工程溢油泄露事故進行影響預測。計算結果表明：（1）不同風向的影響。風對溢油油膜漂移軌跡的影響比較大。靜風條件下，油膜的漂移路徑主要受控制于潮流動力，在NE向風的作用下，油膜會發(fā)生明顯向西北移動現(xiàn)象，而在SSW風向下，油膜同樣會發(fā)生明顯向東北至東南區(qū)間方向偏轉。（2）不同潮時的影響。溢油的漂移軌跡受溢油時刻的影響極大，在潮周期的不同階段發(fā)生溢油事故時，溢油油膜漂移軌跡可能完全不同?？梢酝ㄟ^掃海范圍及油膜漂移軌跡圖明顯看出不同潮時所帶來的不同溢油事故后果。若發(fā)生溢油事故的地點不同，不同潮時所帶來的影響差別大小也不盡相同。

［1］徐艷東.海上溢油風化過程及其預測模型研究［D］.青島：中國海洋大學，2006.

［2］俞濟清，魏敏，黃立文，等.中國舟山港溢油模擬信息系統(tǒng)研制［J］.武漢理工大學學報，2002，26（4）：463-465.YU J Q，WEI M，HUANG L W，et al.The spilled oil simulation information system designed for Zhoushan harbor［J］.Journal of Wuhan University of Technology，2002，26（4）：463-465.

［3］黃立文，宋文銅，孫俊，等.實時流場預報及在海面溢油軌跡預測中的應用［J］.中國航海，2000（2）：36-41.HUANG L W，SONG W T，SUN J，et al.Real time current field forecast and its application to marine oil-spill trajectory prediction［J］.Navigation of China，2000（2）：36-41.

［4］熊德琪，楊建立，嚴世強.珠江口區(qū)域海上溢油應急預報信息系統(tǒng)的開發(fā)研究［J］.海洋環(huán)境科學，2005（2）：63-66.XIONG D Q，YANG J L，YAN S Q.R&D on oil spill contingency forecast information system for Pearl River estuary［J］.Marine Environmental Science，2005（2）：63-66.

［5］張波，吳冠.渤海溢油應急預報系統(tǒng)［J］.海洋環(huán)境科學，1998（3）：72-75.ZHANG B，WU G.Contingent predicting system of oil spill for Bohai Sea［J］.Marine Environment Science，1998（3）：72-75.

［6］竇振興，宮成，張存智，等.海上溢油動態(tài)預測應急系統(tǒng)［J］.中國海洋平臺，1995（5）：31-33.DOU Z X，GONG C，ZHANG C Z，et al.The Dynamic Predication and Emergency System of Oil Spill at Sea［J］.China Offshore Platform，1995（5）：31-33.

［7］婁安剛，王學昌，孫長青，等.膠州灣海面溢油軌跡的數值模擬［J］.黃渤海海洋，2001（1）：1-8.LOU A G，WANG X C，SUN C Q，et al.A simulation of Oil spill Trajectory on Sea Surface of Jiaozhou Bay［J］.Journal of Oceanograpgy of Huanghai&Bohai Seas，2001（1）：1-8.

［8］趙洪波，許婷.甌江口航道二期治理潛堤工程潮流泥沙鹽度數學模型研究［R］.天津：交通運輸部天津水運工程科學研究所，2009.

［9］姜衛(wèi)星.黃浦江溢油事故的數值模擬研究［D］.上海：同濟大學，2007.

［10］JTJ/T233-98，海岸與河口潮流泥沙模擬技術規(guī)程［S］.