甌江口航道海域溢油擴(kuò)散數(shù)值模擬

黃毅峰 ，許 婷 ，劉 濤

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

近半個(gè)世紀(jì)以來全球發(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)化和熱量遷移等過程，對(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ì)周圍海洋環(huán)境產(chǎn)生嚴(yán)重的影響，應(yīng)采取溢油應(yīng)急措施，控制污染影響程度，同時(shí)追蹤油膜污染帶，并通知港航監(jiān)督部門。

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

1.1 潮汐類型

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

1.2 潮流特征

據(jù)1999年10月～2005年7月甌江河口及其附近海域水文全潮各站潮流速觀測(cè)資料調(diào)和分析和計(jì)算，各站主要全日分潮流與主要半日分潮流的比值F=（W01+WK1）/WM2均小于0.5，各站淺水分潮流比值G==0.27～0.48，綜合該兩項(xiàng)指標(biāo)，表明該水域潮流屬正規(guī)淺海半日潮流類型；由于漲落潮流受水域地形限制，基本呈往復(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端的甲米礁東岸波浪觀測(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#（溫州淺灘南口水道附近）波浪觀測(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)境研究所開發(fā)的MIKE 21/3溢油分析（SA）模塊，是基于歐拉-拉格朗日理論體系，通過對(duì)油膜在水體中的擴(kuò)展、傳輸（水流和風(fēng)場(chǎng)作用）、紊動(dòng)擴(kuò)散、分散（夾帶）、蒸發(fā)、乳化、溶解等各種過程的模擬，MIKE 21/3 SA能提供油膜隨時(shí)間變化的漂移位置、厚度，以及漂移過程中粘度、油膜表面溫度、傾點(diǎn)等屬性的變化。另外，MIKE 21/3對(duì)復(fù)雜水體（如各種水工建筑物）的水動(dòng)力能進(jìn)行精確模擬，強(qiáng)大的前后處理功能適用于決策系統(tǒng)的開發(fā)。MIKE系統(tǒng)是國(guó)際上同類型軟件中最先進(jìn)的模擬工具之一，在世界各地已有許多成功的應(yīng)用案例［9］。

2.1 輸移過程

油粒子的輸移包括擴(kuò)展、漂移、擴(kuò)散等過程［9］，這些過程是油粒子位置發(fā)生變化的主要原因，而油粒子的組分在這些過程中不發(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ù)從氣象部門獲得，而流場(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)化過程

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

（1）蒸發(fā)。油膜蒸發(fā)受油分、氣溫和水溫、溢油面積、風(fē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)最重要的過程。擴(kuò)散是一種機(jī)械過程，水流的紊動(dòng)能量將油膜撕裂成油滴，形成水包油的乳化。這些乳化物可以被表面活性劑穩(wěn)定，防止油滴返回到油膜。在惡劣天氣狀況下最主要的擴(kuò)散作用力是波浪破碎，而在平靜的天氣狀況下最主要的擴(kuò)散作用力是油膜的伸展壓縮運(yùn)動(dòng)。

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

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

油滴返回油膜的速率為

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

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

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

2.3 熱量遷移

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

圖1中1為大氣與油膜之間的傳熱過程；2為大氣與油膜之間熱輻射過程；3為太陽輻射；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方向的垂線平均流速分量；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)，必須模擬甌江口整體海域的潮波傳播過程，建立足夠大的計(jì)算區(qū)域。因此，此次數(shù)學(xué)模型的計(jì)算區(qū)域西邊界取至120°30′經(jīng)度線，東邊界至121°45′經(jīng)度線，北至樂清灣，南至南麂島南側(cè)27°10′緯度線。東西向?qū)捈s120 km，南北向?yàn)?50 km。模型采用正四邊形網(wǎng)格，計(jì)算網(wǎng)格步長(zhǎng)為Dx=Dy=50 m，并就復(fù)雜的陸地和島嶼岸線進(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條水文垂線（圖2）。其中1#～7# 垂線施測(cè)大、中、小潮，8#～14# 垂線僅施測(cè)大潮，測(cè)驗(yàn)內(nèi)容為潮位、流速、流向等。該次水文全潮規(guī)模宏大、內(nèi)容豐富、測(cè)點(diǎn)涵蓋了甌江口大范圍海域。部分測(cè)站大潮潮位及流速流向驗(yàn)證曲線見圖3～圖4。

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

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

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

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

在水動(dòng)力條件模擬的基礎(chǔ)上，利用MIKE21 SA模塊建立甌江口航道二維溢油模型，應(yīng)用“油粒子”模型模擬輸移、風(fēng)化和熱量遷移等過程，對(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é)果展示，其掃海范圍及漂移軌跡見圖6～圖12，圖中所示為時(shí)間總長(zhǎng)2個(gè)潮周期（即24 h），時(shí)間分別取2 h、6 h、12 h和24 h的結(jié)果。不同時(shí)間段過后其溢油掃海面積見表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í)，溢油油膜漂移軌跡可能完全不同?？梢酝ㄟ^掃海范圍及油膜漂移軌跡圖明顯看出不同潮時(shí)所帶來的不同溢油事故后果。若發(fā)生溢油事故的地點(diǎn)不同，不同潮時(shí)所帶來的影響差別大小也不盡相同。

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í)，溢油油膜漂移軌跡可能完全不同?？梢酝ㄟ^掃海范圍及油膜漂移軌跡圖明顯看出不同潮時(shí)所帶來的不同溢油事故后果。若發(fā)生溢油事故的地點(diǎn)不同，不同潮時(shí)所帶來的影響差別大小也不盡相同。

［1］徐艷東.海上溢油風(fēng)化過程及其預(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)的開發(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］.