徐圩港區(qū)鄰近海域溢油模擬分析

林偉波，陳曉燕，李婧慧

（1.江蘇省海涂研究中心，江蘇南京210036；2.江蘇省海洋環(huán)境監(jiān)測(cè)預(yù)報(bào)中心，江蘇南京210036）

1 引言

近年來(lái)，隨著我國(guó)近海及海上船舶施工和油氣儲(chǔ)運(yùn)設(shè)施增多，油品泄漏事故增加，不僅給企業(yè)帶來(lái)重大經(jīng)濟(jì)損失，而且對(duì)海洋生態(tài)環(huán)境造成巨大威脅。海上溢油事故發(fā)生后，油膜在海域內(nèi)的運(yùn)動(dòng)及變化是一個(gè)極其復(fù)雜的過(guò)程，它是海洋環(huán)境問(wèn)題的難點(diǎn)之一。國(guó)外自20世紀(jì)60年代開(kāi)始溢油輸運(yùn)擴(kuò)散數(shù)值模擬預(yù)報(bào)方法的研究[1]。Reed等[2]對(duì)20世紀(jì)溢油模型的研究和發(fā)展進(jìn)行了回顧，主要經(jīng)歷了油膜擴(kuò)展理論、對(duì)流擴(kuò)散方法和油粒子法3個(gè)階段?！坝土Ｗ印蹦Ｊ皆谝缬湍M發(fā)展過(guò)程中具有劃時(shí)代的意義，也成為了當(dāng)今主流的溢油模式[3]。國(guó)內(nèi)的學(xué)者在海南島南部[4]、福建東吾洋灣[5]、福寧灣[6]、浙江樂(lè)清灣[7]、長(zhǎng)江口[8]、江蘇濱海[9]、膠州灣及渤海灣等[10-12]不同海域運(yùn)用該方法開(kāi)展了大量的溢油數(shù)值模擬預(yù)報(bào)研究工作。

本文基于丹麥水利研究所（Danish Hydraulic Institute，DHI）開(kāi)發(fā)的MIKE軟件中非結(jié)構(gòu)化網(wǎng)格的水動(dòng)力模型對(duì)徐圩港區(qū)附近海域進(jìn)行潮流數(shù)值模擬，在此基礎(chǔ)上通過(guò)“油粒子”模型對(duì)該區(qū)進(jìn)行溢油模擬，并分析不同情況下溢油油膜的漂移軌跡、掃海面積，為海洋環(huán)境影響評(píng)價(jià)和溢油事故應(yīng)急措施的制定提供科學(xué)依據(jù)。

2 材料與方法

2.1 水動(dòng)力模型控制方程

MIKE21 HD是一個(gè)可用于河口、海灣以及近海潮位及潮流的綜合數(shù)值模型，對(duì)于二維流場(chǎng)模擬已經(jīng)被應(yīng)用于許多工程之中。本文使用該模塊對(duì)徐圩港區(qū)進(jìn)行了水動(dòng)力模擬實(shí)驗(yàn)，潮流數(shù)據(jù)驗(yàn)證良好后輸入溢油模塊，作為溢油水動(dòng)力場(chǎng)的驅(qū)動(dòng)條件。水動(dòng)力模型用連續(xù)方程和動(dòng)量守恒方程來(lái)描述流場(chǎng)和潮位的變化。

連續(xù)方程：

x向動(dòng)量方程：

y向動(dòng)量方程：

式中：x,y為直角坐標(biāo)系坐標(biāo)；t為時(shí)間；ζ為相對(duì)某一基面的水位（單位：m）；h為相對(duì)某一基面的水深（單位：m）；u,v分別為x,y方向上的垂線平均流速分量；Nx為x向水流紊動(dòng)粘性系數(shù)（單位：m2/s）；Ny為y向水流紊動(dòng)粘性系數(shù)（單位：m2/s）；f為科氏力系數(shù)，f=2ωsinφ，φ為緯度；ω為地球自轉(zhuǎn)速度；fb底部摩阻系數(shù)，fb=g/C2，其中C為謝才系數(shù)，g為重力加速度。

2.2 溢油控制方程

（1）隨機(jī)游動(dòng)與擴(kuò)散過(guò)程

在流場(chǎng)計(jì)算的基礎(chǔ)上，將溢油分成有限個(gè)質(zhì)點(diǎn)，每個(gè)質(zhì)點(diǎn)的漂移位置按下式計(jì)算：

式中：X0、Y0為某質(zhì)點(diǎn)的初始坐標(biāo)；U、V為X、Y方向的流速分量；W10為海面10 m高處風(fēng)速；A為風(fēng)向；α為修正系數(shù)；r為隨機(jī)擴(kuò)散系數(shù)，r=R·E，R為0～1之間的隨機(jī)數(shù)，E為擴(kuò)散系數(shù)；B為隨機(jī)擴(kuò)散方向，B=2πR。

（2）延展過(guò)程

延展過(guò)程決定了表面浮油的面積擴(kuò)展，從而進(jìn)一步影響水面油膜的蒸發(fā)、溶解、擴(kuò)散和光氧化作用。延展是湍流擴(kuò)散以及重力、慣性、黏性和表面張力平衡的聯(lián)合作用結(jié)果。采用修正的Fay理論基礎(chǔ)上的重力-粘力公式計(jì)算油膜擴(kuò)展：

式中：Aoil為油膜面積,為油膜半徑；Kα為系數(shù)（率定為0.6）；t為時(shí)間；油膜體積Voil為：

式中：hs為油膜初始厚度。

（3）蒸發(fā)過(guò)程

蒸發(fā)過(guò)程可導(dǎo)致10%～30%的浮油從水面進(jìn)入大氣，具體百分比取決于油種。簡(jiǎn)單時(shí)間型蒸發(fā)計(jì)算式如下：

式中：A為油特征常數(shù)；B為油的溫度特征常數(shù)；T為油溫，單位：℃；t為油齡，單位：min。

（4）水體攜帶過(guò)程

水面浮油暴露在風(fēng)和浪中，浮油會(huì)被攜帶或擴(kuò)散進(jìn)入水體。水體攜帶是一種物理過(guò)程，在破碎浪的作用下，油滴從水面遷移到水體中。水體攜帶速率Qd（kg/（m2·s））和油粒子大小之間的關(guān)系見(jiàn)下式：

式中：C*為與油種和風(fēng)化狀態(tài)相關(guān)的水體攜帶速率經(jīng)驗(yàn)常數(shù)；Dd為單位表面積耗散的破碎波能量，單位：J/m2；S為浮油覆蓋的水面面積分?jǐn)?shù)；F為受破碎浪侵襲的水面面積分?jǐn)?shù)；d為油粒子直徑，單位：m；Δd為油粒子直徑差，單位：m。

（5）乳化過(guò)程

水-油乳化物（或稱(chēng)為乳膠狀物）的形成取決于油的組分和水環(huán)境條件。乳化油可能有80%是以連續(xù)相油存在的微米級(jí)油粒子。一般乳化油的黏度要高于形成乳化油之前的油品黏度。由于水的混入，油/水混合物的體積明顯加大。水混入油相的速率計(jì)算方法見(jiàn)下式：

式中：U為風(fēng)速，單位：m/s；Kem為乳化率常數(shù)；Yw為水分?jǐn)?shù)；Ymax為水在油相中的最大比例。

2.3 水動(dòng)力計(jì)算參數(shù)及網(wǎng)格

模型計(jì)算域從南至北包含了鹽城、連云港、日照所在海域，南北跨度約130 km，東至30 m等深線處。將其劃分成非結(jié)構(gòu)三角形網(wǎng)格，并在徐圩港區(qū)內(nèi)進(jìn)行了加密處理，在港區(qū)內(nèi)部附近分辨率最大達(dá)到30 m，港池內(nèi)水深在1～5 m之間，整個(gè)計(jì)算區(qū)域包括三角形網(wǎng)格39 544個(gè)，網(wǎng)格節(jié)點(diǎn)數(shù)20 554個(gè)（見(jiàn)圖1）。開(kāi)邊界上使用東中國(guó)海模型預(yù)報(bào)得到邊界水位進(jìn)行驅(qū)動(dòng)，閉邊界采用干濕邊界處理。模型使用冷啟動(dòng)，時(shí)間步長(zhǎng)為60 s，每小時(shí)輸出潮流場(chǎng)數(shù)據(jù)，計(jì)算10 d，取后6 d穩(wěn)定流場(chǎng)進(jìn)行驗(yàn)證分析，得到該區(qū)水動(dòng)力場(chǎng)。

2.4 溢油模型參數(shù)設(shè)置

圖1 模型計(jì)算網(wǎng)格（單位：m）

表1 溢油風(fēng)況及潮時(shí)組合

工程實(shí)施過(guò)程中的施工船舶燃料油泄漏會(huì)造成海域泄漏事故，不同船型、不同事故情況的溢油量都不同，具有較大的隨機(jī)性。此次模擬的溢油油品設(shè)定為燃料油，密度為920 kg/m3，最大含水率為0.85，以10 t燃料油作為泄漏源強(qiáng)。溢油類(lèi)型為固定點(diǎn)源瞬時(shí)溢油，油粒子數(shù)為1 000。油膜漂移過(guò)程中抵達(dá)岸邊視為被攔截，不再計(jì)算溢油量。在本工程實(shí)施時(shí)，溢油最大風(fēng)險(xiǎn)可能發(fā)生位置位于四區(qū)導(dǎo)堤龍口南側(cè)區(qū)。根據(jù)連云港海洋站1974—2005年風(fēng)速、風(fēng)頻率資料，夏季常風(fēng)向?yàn)镋，風(fēng)速為4.8 m/s，冬季常風(fēng)向?yàn)镹，風(fēng)速為4.5 m/s，最不利風(fēng)向?yàn)镾E，風(fēng)速為10.7 m/s。本文模擬了夏季常風(fēng)向、冬季常風(fēng)向及最不利風(fēng)向下漲潮和落潮時(shí)刻分別發(fā)生溢油的共計(jì)5種情景，具體溢油風(fēng)況及潮時(shí)組合見(jiàn)表1。利用水動(dòng)力結(jié)果作為溢油模型驅(qū)動(dòng)，考慮風(fēng)況及溢油潮時(shí)組合，溢油模型時(shí)間步長(zhǎng)為3 600 s，每小時(shí)輸出溢油結(jié)果，模擬6個(gè)潮周期（72 h）內(nèi)油膜漂移運(yùn)動(dòng)。溢油模型考慮水平擴(kuò)散、蒸發(fā)和乳化過(guò)程，水平擴(kuò)散系數(shù)取0.01 m2/s，蒸發(fā)常數(shù)A和B取2.67和0.06，乳化率常數(shù)取10-6，最大水分?jǐn)?shù)取0.75。溢油蒸發(fā)比例隨時(shí)間變化在10%～26%之間。

3 結(jié)果與討論

3.1 模型驗(yàn)證

潮位驗(yàn)證數(shù)據(jù)潮位采用2015年1月23日13時(shí)—1月29日15時(shí)（北京時(shí)，下同）的現(xiàn)場(chǎng)實(shí)測(cè)資料，共設(shè)3個(gè)潮位站，分別為贛榆港區(qū)H1、西連島H2和開(kāi)山島H3站，流速及流向驗(yàn)證數(shù)據(jù)為徐圩港區(qū)V1及灌河口外部V2測(cè)站2015年1月23日13時(shí)—1月24日16時(shí)的大潮實(shí)測(cè)資料，2015年1月28日11時(shí)—1月29日15時(shí)的小潮期實(shí)測(cè)資料，潮流、潮位測(cè)站分布見(jiàn)圖2。圖中流速的正負(fù)分別代表潮流的漲落，落潮流速為正值，漲潮流速為負(fù)值。根據(jù)驗(yàn)證結(jié)果可知，該模型模擬所得潮位及流速流向與實(shí)測(cè)數(shù)據(jù)驗(yàn)證較為吻合，可較好地反映當(dāng)?shù)厮畡?dòng)力形態(tài)，計(jì)算所得水動(dòng)力場(chǎng)可用于溢油模型（見(jiàn)圖3—7）。

圖2 潮位、潮流驗(yàn)證測(cè)站位置

圖3 贛榆港區(qū)H1潮位站潮位驗(yàn)證過(guò)程線

圖4 西連島H2潮位站潮位驗(yàn)證過(guò)程線

圖5 開(kāi)山島H3潮位站潮位驗(yàn)證過(guò)程線

圖6 V1測(cè)站流速與流向計(jì)算值與實(shí)測(cè)值的比較

圖7 V2測(cè)站流速與流向計(jì)算值與實(shí)測(cè)值的比較

3.2 水動(dòng)力模擬結(jié)果

油類(lèi)在海面漂移，表層流場(chǎng)為其主要的驅(qū)動(dòng)力，圖8為大潮期漲落潮表層流場(chǎng)分布圖，漲潮時(shí)，外海潮流基本以NE～SW方向進(jìn)入海州灣；落潮時(shí)，潮流則基本以SW～NE向退出海州灣；潮流的流向與等深線或岸線的交角較大，即潮流的沿岸運(yùn)動(dòng)趨勢(shì)較小，而以離岸、向岸的往復(fù)運(yùn)動(dòng)為主。

3.3 溢油模擬結(jié)果

表2和圖9—11顯示了當(dāng)船舶分別在漲潮階段及落潮階段發(fā)生泄漏時(shí)，油膜粒子經(jīng)72 h后的漂移軌跡及掃海范圍。

通過(guò)預(yù)測(cè)可以看出，在夏季常風(fēng)條件下，當(dāng)泄漏發(fā)生在漲潮階段，油膜向北側(cè)圍堤內(nèi)漂移，72 h后油膜掃海面積0.6 km2；當(dāng)泄漏發(fā)生在落潮階段，油膜在落潮流和風(fēng)作用下，油膜先向西北漂移出防波堤口門(mén)后，再往西漂移，72 h后油膜掃海面積約15.8 km2。

表2 溢油風(fēng)險(xiǎn)預(yù)測(cè)分析

圖8 大潮流場(chǎng)分布（單位：m/s）

圖9 工程局部大潮流場(chǎng)分布（單位：m/s）

圖10 夏季常風(fēng)條件下1～72 h掃海面積和漂移軌跡

在冬季常風(fēng)條件下，當(dāng)泄漏發(fā)生在漲潮階段，油膜向北側(cè)圍堤內(nèi)漂移，72 h后油膜掃海面積2.1 km2；當(dāng)泄漏發(fā)生在落潮階段，油膜在落潮流和風(fēng)作用下，油膜先向西北側(cè)漂移，出口門(mén)后，在N風(fēng)的作用下，往口門(mén)內(nèi)側(cè)移動(dòng)，72 h后油膜掃海面積約11.8 km2。

從夏季和冬季風(fēng)常風(fēng)向的實(shí)驗(yàn)中可以看出漲潮時(shí)溢油僅在圍堤內(nèi)漂移，因此本文對(duì)漲潮時(shí)最不利SE風(fēng)向的情況不再計(jì)算。

在最不利風(fēng)向條件下，當(dāng)泄漏發(fā)生在落潮階段，油膜在落潮流和風(fēng)作用下，油膜先向西北側(cè)漂移，在3.5 h后抵達(dá)東西防波堤口門(mén)，之后在SE風(fēng)和潮流的共同作用下，往西北方向移動(dòng)，在11 h后到達(dá)田灣核電站取水工程，此時(shí)油膜距離泄漏位置漂移路程約為25.8 km，72 h后油膜掃海面積約28.2 km2。

圖11 冬季常風(fēng)條件下1～72 h掃海面積和漂移軌跡

圖12 落潮期最不利風(fēng)條件下1～72 h掃海面積和漂移軌跡

4 結(jié)論

利用MIKE21 HD模塊建立了連云港海域水動(dòng)力模型，較好地模擬了該海域的潮流情況。在此基礎(chǔ)上，利用拉格朗日“油粒子”模型對(duì)徐圩港區(qū)附近溢油情景進(jìn)行了模擬，由數(shù)值模擬結(jié)果可得到以下結(jié)論：

（1）油類(lèi)在水中運(yùn)動(dòng)受多種因素共同影響，不同潮時(shí)不同風(fēng)況下溢油的漂移軌跡、漂移路程及掃海面積都存在較大差異。

（2）徐圩港區(qū)附近溢油漂移擴(kuò)散主要受潮流和風(fēng)場(chǎng)影響，72 h內(nèi)油膜最大掃海面積及漂移路程均出現(xiàn)在落潮期最不利風(fēng)時(shí)溢油，分別為28.2 km2及25.8 km。徐圩港區(qū)內(nèi)落潮時(shí)發(fā)生溢油，在夏季風(fēng)E風(fēng)和最不利風(fēng)SE風(fēng)的條件下，油膜會(huì)經(jīng)過(guò)口門(mén)飄向西北側(cè)海域，對(duì)該區(qū)生態(tài)環(huán)境會(huì)造成一定影響。

（3）不同風(fēng)況下溢油，油膜受風(fēng)及潮流共同驅(qū)使，油膜漂移方向與表面海流和風(fēng)所引起的流速之矢量和的方向大致相同。