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

林偉波，陳曉燕，李婧慧

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

1 引言

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

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

2 材料與方法

2.1 水動力模型控制方程

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

連續(xù)方程：

x向動量方程：

y向動量方程：

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

2.2 溢油控制方程

（1）隨機游動與擴散過程

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

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

（2）延展過程

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

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

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

（3）蒸發(fā)過程

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

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

（4）水體攜帶過程

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

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

（5）乳化過程

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

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

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

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

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

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

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

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

3 結(jié)果與討論

3.1 模型驗證

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

圖2 潮位、潮流驗證測站位置

圖3 贛榆港區(qū)H1潮位站潮位驗證過程線

圖4 西連島H2潮位站潮位驗證過程線

圖5 開山島H3潮位站潮位驗證過程線

圖6 V1測站流速與流向計算值與實測值的比較

圖7 V2測站流速與流向計算值與實測值的比較

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

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

3.3 溢油模擬結(jié)果

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

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

表2 溢油風(fēng)險預(yù)測分析

圖8 大潮流場分布（單位：m/s）

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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