林偉波，陳曉燕，李婧慧

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

1 引言

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

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

2 材料與方法

2.1 水動力模型控制方程

MIKE21 HD是一個可用于河口、海灣以及近海潮位及潮流的綜合數(shù)值模型，對于二維流場模擬已經(jī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高處風速；A為風向；α為修正系數(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ā)過程可導致10%～30%的浮油從水面進入大氣，具體百分比取決于油種。簡單時間型蒸發(fā)計算式如下：

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

（4）水體攜帶過程

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

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

（5）乳化過程

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

式中：U為風速，單位：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）。開邊界上使用東中國海模型預報得到邊界水位進行驅(qū)動，閉邊界采用干濕邊界處理。模型使用冷啟動，時間步長為60 s，每小時輸出潮流場數(shù)據(jù)，計算10 d，取后6 d穩(wěn)定流場進行驗證分析，得到該區(qū)水動力場。

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

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

表1 溢油風況及潮時組合

工程實施過程中的施工船舶燃料油泄漏會造成海域泄漏事故，不同船型、不同事故情況的溢油量都不同，具有較大的隨機性。此次模擬的溢油油品設(shè)定為燃料油，密度為920 kg/m3，最大含水率為0.85，以10 t燃料油作為泄漏源強。溢油類型為固定點源瞬時溢油，油粒子數(shù)為1 000。油膜漂移過程中抵達岸邊視為被攔截，不再計算溢油量。在本工程實施時，溢油最大風險可能發(fā)生位置位于四區(qū)導堤龍口南側(cè)區(qū)。根據(jù)連云港海洋站1974—2005年風速、風頻率資料，夏季常風向為E，風速為4.8 m/s，冬季常風向為N，風速為4.5 m/s，最不利風向為SE，風速為10.7 m/s。本文模擬了夏季常風向、冬季常風向及最不利風向下漲潮和落潮時刻分別發(fā)生溢油的共計5種情景，具體溢油風況及潮時組合見表1。利用水動力結(jié)果作為溢油模型驅(qū)動，考慮風況及溢油潮時組合，溢油模型時間步長為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向退出海州灣；潮流的流向與等深線或岸線的交角較大，即潮流的沿岸運動趨勢較小，而以離岸、向岸的往復運動為主。

3.3 溢油模擬結(jié)果

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

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

表2 溢油風險預測分析

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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