黃毅峰 ，許 婷 ，劉 濤

（1.長(zhǎng)沙理工大學(xué)水利學(xué)院，長(zhǎng)沙410076；2.交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所工程泥沙交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，天津300456）

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

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

1 工程海域自然動(dòng)力條件

1.1 潮汐類(lèi)型

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

1.2 潮流特征

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

1.3 風(fēng)況

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

1.4 波浪特征

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

據(jù)溫州淺灘圍涂工程1#（北堤處）和2#（溫州淺灘南口水道附近）波浪觀(guān)測(cè)站資料統(tǒng)計(jì)，2001～2004年1#站測(cè)得波高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 溢油預(yù)測(cè)模型概述

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

2.1 輸移過(guò)程

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

（1）擴(kuò)展運(yùn)動(dòng)。采用修正的Fay重力-粘力公式計(jì)算油膜擴(kuò)展

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

初始油膜厚度

（2）漂移運(yùn)動(dòng)。油粒子漂移的作用力是水流和風(fēng)拽力，油粒子總漂移速度由以下權(quán)重公式計(jì)算

式中：Uw為水面以上10 m處的風(fēng)速；Us為表面流速；cw為風(fēng)漂移系數(shù)，一般在0.03～0.04。風(fēng)場(chǎng)數(shù)據(jù)從氣象部門(mén)獲得，而流場(chǎng)從二維水動(dòng)力模型計(jì)算結(jié)果獲得。

（3）紊動(dòng)擴(kuò)散。假定水平擴(kuò)散各向同性，一個(gè)時(shí)間步長(zhǎng)內(nèi)α方向上的可能擴(kuò)散距離Sα可表示為

式中：［ ］

R1-1為-1～1的隨機(jī)數(shù)；Dα為α方向上的擴(kuò)散系數(shù)。

2.2 風(fēng)化過(guò)程

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

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

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

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

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

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

從油膜擴(kuò)散到水體中的油分損失量計(jì)算

式中：Da為進(jìn)入到水體的分量；Db為進(jìn)入到水體后沒(méi)有返回的分量

油滴返回油膜的速率為

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

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

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

2.3 熱量遷移

蒸氣壓與粘度受溫度影響，而且觀(guān)察發(fā)現(xiàn)通常油膜的溫度要高于周?chē)拇髿夂退w。圖1為油膜的熱平衡示意圖。

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

3 水動(dòng)力條件模擬

溢油擴(kuò)散模擬結(jié)果是建立在水動(dòng)力條件模擬的基礎(chǔ)之上，文中采用MIKE21 HD模塊計(jì)算水動(dòng)力條件。

3.1 平面二維水流控制方程

連續(xù)方程

運(yùn)動(dòng)方程

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

3.2 計(jì)算域的確定及網(wǎng)格剖分

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

3.3 模型驗(yàn)證

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

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

3.4 工程海域潮流場(chǎng)模擬

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

4 溢油影響預(yù)測(cè)

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

表1 事故溢油掃海面積預(yù)測(cè)結(jié)果Tab.1 Sweeping area of oil spill accident prediction

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

5 結(jié)論

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

［1］徐艷東.海上溢油風(fēng)化過(guò)程及其預(yù)測(cè)模型研究［D］.青島：中國(guó)海洋大學(xué)，2006.

［2］俞濟(jì)清，魏敏，黃立文，等.中國(guó)舟山港溢油模擬信息系統(tǒng)研制［J］.武漢理工大學(xué)學(xué)報(bào)，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］黃立文，宋文銅，孫俊，等.實(shí)時(shí)流場(chǎng)預(yù)報(bào)及在海面溢油軌跡預(yù)測(cè)中的應(yīng)用［J］.中國(guó)航海，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］熊德琪，楊建立，嚴(yán)世強(qiáng).珠江口區(qū)域海上溢油應(yīng)急預(yù)報(bào)信息系統(tǒng)的開(kāi)發(fā)研究［J］.海洋環(huán)境科學(xué)，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］張波，吳冠.渤海溢油應(yīng)急預(yù)報(bào)系統(tǒng)［J］.海洋環(huán)境科學(xué)，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］竇振興，宮成，張存智，等.海上溢油動(dòng)態(tài)預(yù)測(cè)應(yīng)急系統(tǒng)［J］.中國(guó)海洋平臺(tái)，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］婁安剛，王學(xué)昌，孫長(zhǎng)青，等.膠州灣海面溢油軌跡的數(shù)值模擬［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］趙洪波，許婷.甌江口航道二期治理潛堤工程潮流泥沙鹽度數(shù)學(xué)模型研究［R］.天津：交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所，2009.

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

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