匡翠萍， 謝華浪， 蘇 平， 陳 括

(1. 同濟(jì)大學(xué) 土木工程學(xué)院， 上海 200092； 2. 廣西壯族自治區(qū)交通規(guī)劃勘察設(shè)計(jì)研究院, 廣西 南寧 530029)

?

蓬萊19-3溢油運(yùn)移與歸宿數(shù)值模擬和分析

匡翠萍1， 謝華浪1， 蘇 平2， 陳 括1

(1. 同濟(jì)大學(xué) 土木工程學(xué)院， 上海 200092； 2. 廣西壯族自治區(qū)交通規(guī)劃勘察設(shè)計(jì)研究院, 廣西 南寧 530029)

通過(guò)耦合水動(dòng)力模型建立渤海海上溢油運(yùn)移與歸宿數(shù)學(xué)模型，包含漂移、擴(kuò)散、蒸發(fā)、乳化等過(guò)程，通過(guò)拉格朗日方法模擬漂移過(guò)程、蒙特卡羅方法模擬油粒子擴(kuò)散的隨機(jī)性，采用變化的紊動(dòng)擴(kuò)散系數(shù)使溢油所受紊動(dòng)擴(kuò)散作用隨時(shí)間變化.通過(guò)率定得到蓬萊19-3溢油模擬中風(fēng)拖拽系數(shù)和擴(kuò)散系數(shù)校正參數(shù)b的合理取值范圍分別為1.4%～2.3%和0.38～0.42.對(duì)蓬萊19-3溢油事故進(jìn)行了模擬，通過(guò)溢油附岸位置與時(shí)間驗(yàn)證了模型的準(zhǔn)確性.模擬結(jié)果表明: 由于風(fēng)速較大、方向變化較小，風(fēng)對(duì)漂浮于水面上的溢油產(chǎn)生較大作用，且溢油漂移方向與主風(fēng)向一致；渤海海域內(nèi)的潮流以往復(fù)流、旋轉(zhuǎn)流為主且余流小使其對(duì)溢油運(yùn)移的作用弱于風(fēng)；溢油初期的蒸發(fā)速度較快，后期因揮發(fā)組分的減少以及乳化作用使溢油含水率增大使得蒸發(fā)速度明顯減小；蒸發(fā)和附岸是溢油漂浮組分減少的主要原因，也是溢油的主要?dú)w宿.

渤海； 海上溢油； 運(yùn)移； 數(shù)值模擬

目前，海上石油污染已成為一種全球性海洋污染[1].海上溢油以其突發(fā)性、污染的嚴(yán)重性以及應(yīng)急處理的艱巨性使其在過(guò)去的幾十年里成為海洋生態(tài)環(huán)境最大的破壞因素，且呈現(xiàn)越來(lái)越嚴(yán)重的趨勢(shì)[2].尤其當(dāng)溢油事故發(fā)生于近岸或河口水體環(huán)境時(shí)，被污染的水體及近岸環(huán)境將受到慘重的破壞[3].2011年渤海蓬萊19-3溢油事故造成了極其惡劣的環(huán)境破壞，海水污染面積高達(dá)4 900 km2，造成的海洋生態(tài)損害價(jià)值總計(jì)16.83億元[4].根據(jù)康菲公司公布數(shù)據(jù)，本次事故溢油量為500 t，然而多數(shù)研究學(xué)者認(rèn)為實(shí)際溢油量甚至高達(dá)該數(shù)據(jù)的十倍以上[5].

海上溢油屬于突發(fā)性的海洋污染事故，且傳播速度快、影響范圍大.但由于海洋環(huán)境復(fù)雜多變，溢油事故發(fā)生于不同的時(shí)間和地點(diǎn)都可能產(chǎn)生不同的漂移軌跡和污染范圍[6].因此，對(duì)海上溢油進(jìn)行數(shù)值模擬，有利于預(yù)測(cè)溢油變化趨勢(shì)，制定合理的治理方案，及時(shí)采取應(yīng)對(duì)措施，最大程度地減少溢油事故帶來(lái)的破壞和損失.徐青等[7]通過(guò)GNOME導(dǎo)入潮流和風(fēng)建立了渤海溢油模型，通過(guò)將溢油油膜概化為由多數(shù)圓形斑點(diǎn)組成，并采用拉格朗日方法模擬了渤海蓬萊19-3溢油事故6月11—19日中部分油膜的漂移過(guò)程.鄧增安等[3-6]考慮了風(fēng)、潮流、波浪等因素，研發(fā)了嵌套于中國(guó)數(shù)字海洋原型系統(tǒng)(CDOPS)中的BSOSSM溢油模擬服務(wù)組件，并對(duì)蓬萊19-3溢油事故進(jìn)行了模擬.汪守東[8-9]與王金華等[10]將海上溢油分為表層浮油和懸浮油滴，通過(guò)離散顆粒的拉格朗日算法進(jìn)行計(jì)算，考慮了溢油的漂移、擴(kuò)散、蒸發(fā)和乳化等過(guò)程，建立了渤海海域溢油模型.畢海普等[11]建立了三峽水庫(kù)溢油模型，通過(guò)在模型中設(shè)置一個(gè)衰減參數(shù)來(lái)模擬與分析三峽水庫(kù)溢油事故發(fā)生后各溢油組分在機(jī)械回收和投放分散劑、乳化劑等緊急措施影響下的變化情況.李燕等[12]建立了三維的拉格朗日隨機(jī)步溢油模型，并研究了在不同方案下油粒子垂直擴(kuò)散的變化情況.總體來(lái)說(shuō)，目前各種溢油模型所采用方法與考慮因素各有不同.其中，對(duì)溢油的漂移過(guò)程模擬均采用拉格朗日方法，計(jì)算中主要將風(fēng)和潮流的作用作為溢油漂移的主要影響因素.半封閉的渤海海域中風(fēng)對(duì)溢油漂移的影響是最重要的，一般認(rèn)為風(fēng)拖拽系數(shù)取值范圍在1%～3%之間，但對(duì)渤海區(qū)域風(fēng)拖拽系數(shù)的取值尚未有詳細(xì)的討論和研究.另外，多數(shù)溢油模型由于對(duì)模擬粒子數(shù)限制較大，以至于對(duì)大規(guī)模的溢油模擬未能形成一個(gè)較為直觀的溢油平面密度分布場(chǎng).

利用Delft3D通過(guò)對(duì)水動(dòng)力模型進(jìn)行耦合建立了渤海海上溢油模型，模型計(jì)算包括漂移、擴(kuò)散、蒸發(fā)、乳化等變化過(guò)程.本模型將溢油作為離散化的油粒子，通過(guò)拉格朗日方法模擬漂移過(guò)程、蒙特卡羅方法模擬湍流擴(kuò)散過(guò)程，且采用隨時(shí)間變化的紊動(dòng)擴(kuò)散系數(shù)使溢油所受湍流作用亦隨時(shí)間變化，使得大規(guī)模的渦流和環(huán)流對(duì)溢油的擴(kuò)散作用得以體現(xiàn)出來(lái)，并采用大量的油粒子進(jìn)行模擬使溢油分布密度具有較高分辨率，對(duì)于判斷不同區(qū)域的污染程度差異更為清晰直觀.通過(guò)模型率定得到了風(fēng)拖拽系數(shù)和擴(kuò)散系數(shù)校正參數(shù)b的合理取值范圍，在此基礎(chǔ)上對(duì)2011年蓬萊19-3溢油事故進(jìn)行了驗(yàn)證模擬，研究了渤海海域風(fēng)和潮流對(duì)溢油運(yùn)移的影響，并分析了溢油歸宿的變化特征.

1 數(shù)學(xué)模型

1.1 Delft3D-FLOW水動(dòng)力模型

FLOW模型的基本方程為二維或三維的非線性淺水方程，采用淺水假定和布西奈斯克(Boussinesq)近似由不可壓縮流體的三維Navier-Stokes方程推求得出.由于垂向加速度遠(yuǎn)遠(yuǎn)小于重力加速度，故在垂直方向上的動(dòng)量方程中忽略了垂直加速度的影響，并由此推導(dǎo)出在靜壓假定下的水流方程[13].平面二維水動(dòng)力模型的基本方程為

(1) 連續(xù)方程

(1)

式中：ζ為水位；d為基準(zhǔn)面下的深度，總水深H=d+ζ；u和ν分別表示水平方向上ζ和η向的流速分量；Gξξ和Gηη分別為正交曲線坐標(biāo)系和Cartesian坐標(biāo)系之間的轉(zhuǎn)換系數(shù)；Q為單位面積的源匯項(xiàng).

(2) 動(dòng)量方程

ξ方向的動(dòng)量方程

(2)

η方向的動(dòng)量方程

(3)式中：ξ和η為平面曲線坐標(biāo)；f=2Ωsinφ，為科氏力參數(shù)，由地球自轉(zhuǎn)角速度和緯度值決定；ρ0為水的密度；Pξ和Pη為壓力梯度；Fξ和Fη為水平雷諾應(yīng)力的不平衡性；Mξ和Mη為外來(lái)源匯項(xiàng)引起的附加動(dòng)量.

1.2 Delft3D-PART模塊

Delft3D-PART模型是一個(gè)包含隨機(jī)運(yùn)動(dòng)的粒子追蹤模型，它將溶解或懸浮于水中的物質(zhì)通過(guò)一定數(shù)量離散的油粒子來(lái)表示，每個(gè)油粒子代表一定質(zhì)量的溢油.油粒子的運(yùn)動(dòng)包括漂移與擴(kuò)散.漂移主要由流和風(fēng)的作用引起；擴(kuò)散則主要由湍流引起[14].本溢油模型計(jì)算還包括蒸發(fā)和乳化等變化過(guò)程.

1.2.1 漂移

漂移包括平流與風(fēng)的拖拽產(chǎn)生的移動(dòng).每一個(gè)油粒子的平流速度通過(guò)所在計(jì)算網(wǎng)格節(jié)點(diǎn)的速度(由潮流模型獲得)線性內(nèi)插得到.風(fēng)的拖拽計(jì)算中，風(fēng)的影響僅僅作用于溢油，而對(duì)流場(chǎng)不再產(chǎn)生影響.風(fēng)的拖曳作用，即將溢油因風(fēng)產(chǎn)生的速度增量表達(dá)為風(fēng)速的百分比進(jìn)行計(jì)算，其關(guān)系見(jiàn)式(4).漂移過(guò)程通過(guò)拉格朗日方法進(jìn)行模擬，其表達(dá)式如式(5)所示.

(4)

(5)

1.2.2 擴(kuò)散

在釋放粒子后的初始階段，粒子團(tuán)相對(duì)較小，這時(shí)粒子的擴(kuò)散作用只是由小規(guī)模的湍流作用引起.當(dāng)“粒子云團(tuán)”充分?jǐn)U散開來(lái)，大規(guī)模的渦流將對(duì)粒子的擴(kuò)散產(chǎn)生明顯的影響.故本模型采用隨時(shí)間變化紊動(dòng)擴(kuò)散系數(shù)D，使得模擬溢油所受湍流過(guò)程亦隨時(shí)間變化，以更好地模擬大規(guī)模渦流和環(huán)流對(duì)后期溢油的擴(kuò)散作用.紊動(dòng)擴(kuò)散系數(shù)D的冪函數(shù)表達(dá)式為

D=atb

(6)

式中：t為時(shí)間，以粒子釋放時(shí)間為初始時(shí)刻(t=0)；a和b為紊動(dòng)擴(kuò)散系數(shù)的校正參數(shù)，且0

湍流擴(kuò)散是在水動(dòng)力速度場(chǎng)的基礎(chǔ)上使油粒子產(chǎn)生一個(gè)附加位移.這個(gè)附加位移的方向是隨機(jī)選取的，通過(guò)蒙特卡羅方法來(lái)模擬.蒙特卡羅方法即利用計(jì)算機(jī)通過(guò)構(gòu)造符合一定規(guī)則的隨機(jī)數(shù)以模擬隨機(jī)現(xiàn)象進(jìn)行近似計(jì)算的方法.在紊動(dòng)擴(kuò)散的過(guò)程中，其隨機(jī)性的模擬基于一個(gè)白噪聲過(guò)程，其隨機(jī)數(shù)的取值范圍為-1～+1，且均值為0.每個(gè)油粒子在一個(gè)時(shí)間步長(zhǎng)Δt內(nèi)的擴(kuò)散位移為

(7)

式中：ΔS為單位時(shí)間步長(zhǎng)的擴(kuò)散位移；R為隨機(jī)數(shù)，取值范圍為[-1,1].

1.2.3 蒸發(fā)

PART溢油模型將溢油分為揮發(fā)組分與不揮發(fā)組分，不揮發(fā)組分不發(fā)生蒸發(fā)過(guò)程，揮發(fā)組分將產(chǎn)生一階衰減的蒸發(fā)過(guò)程.表達(dá)如下:

(8)

式中：Fvol為揮發(fā)組分分?jǐn)?shù);Fv為已經(jīng)蒸發(fā)掉的組分分?jǐn)?shù);k為蒸發(fā)率.

1.2.4 乳化

乳化過(guò)程本身是相對(duì)快速的，在實(shí)驗(yàn)室中乳化過(guò)程從發(fā)生到乳化完成大約需要0.1～3 h[15].但形成油包水乳化物的最長(zhǎng)時(shí)間一般在10～100 h之間，這主要是因?yàn)橐缬桶l(fā)生后往往要經(jīng)過(guò)一段時(shí)間乳化作用才逐漸發(fā)生.乳化幾乎是一個(gè)不可逆轉(zhuǎn)的過(guò)程，它使溢油變?yōu)槊芏容^大、具有高粘度的半固體物質(zhì).乳化過(guò)程根據(jù)Mackay等人[16]所提出的算法進(jìn)行計(jì)算，其吸水率為

(9)

乳化過(guò)程不僅起到改變粘度的作用，它還將影響到溢油的蒸發(fā)過(guò)程和水下分散過(guò)程.根據(jù)Fingas[17]的研究，乳化過(guò)程能夠使溢油粘度增大2～3個(gè)數(shù)量級(jí)，溢油擴(kuò)散到水體中的過(guò)程減慢，而蒸發(fā)過(guò)程將近停止.為了使乳化過(guò)程和蒸發(fā)過(guò)程之間形成關(guān)聯(lián)，計(jì)算中乳化物的含水率的增大將減小蒸發(fā)率，模型通過(guò)減小溢油中蒸發(fā)組分的大小來(lái)實(shí)現(xiàn)，相應(yīng)的計(jì)算公式如下：

(10)

式中：Fvol為乳化發(fā)生前溢油的蒸發(fā)組分；Few為乳化發(fā)生后溢油的蒸發(fā)組分.當(dāng)乳化物達(dá)到最大含水率時(shí)，F(xiàn)ew取值為0，蒸發(fā)過(guò)程終止.

1.2.5 附岸

附岸過(guò)程即溢油粘附于岸邊并不再進(jìn)入水體的過(guò)程.當(dāng)溢油進(jìn)入近岸海域后將發(fā)生附岸過(guò)程，并對(duì)海岸生態(tài)環(huán)境造成一定的破壞.本模型通過(guò)設(shè)定一個(gè)粘附概率進(jìn)行模擬.當(dāng)油粒子進(jìn)入岸邊或沉至床底時(shí)會(huì)獲得一個(gè)介于0～1之間的一個(gè)隨機(jī)數(shù)，當(dāng)隨機(jī)數(shù)小于粘附概率時(shí)則油粒子粘附于岸上并不再參與計(jì)算.

2 數(shù)學(xué)模型建立與驗(yàn)證

本模型計(jì)算范圍為以大連—煙臺(tái)為開邊界的整個(gè)渤海海域，采用正交曲線網(wǎng)格，網(wǎng)格尺度約為1 000 m.計(jì)算網(wǎng)格與邊界如圖1所示.

圖1 計(jì)算區(qū)域與網(wǎng)格

2.1 水動(dòng)力模型及驗(yàn)證

水動(dòng)力模型通過(guò)環(huán)渤海多個(gè)測(cè)站的潮位潮流資料進(jìn)行了驗(yàn)證.以下僅列出2011年5月20—21日秦皇島與山海關(guān)的潮位驗(yàn)證如圖2和圖3所示，以及2011年5月20日11點(diǎn)～21日11點(diǎn)V01和V10測(cè)站的潮流驗(yàn)證(圖4).

圖2 秦皇島潮位驗(yàn)證(2011-05-20—21)

圖3 山海關(guān)潮位驗(yàn)證(2011-05-20—21)

圖4 V01和V10潮流流速流向驗(yàn)證 (2011-05-20 11:00—5.21 11:00)

采用Willmott[18]的統(tǒng)計(jì)學(xué)方法對(duì)水動(dòng)力模型計(jì)算結(jié)果進(jìn)行效率評(píng)價(jià)，計(jì)算公式如下：

(11)

圖5與圖6分別為2011年6月4日漲急時(shí)刻與落急時(shí)刻渤海流場(chǎng)圖.渤海灣、遼東灣和萊州灣潮流以往復(fù)流為主，渤海灣漲落潮流以W—E向?yàn)橹?，遼東灣漲落潮流以NE—SW向?yàn)橹鳎R州灣漲落潮流以SW—NE向?yàn)橹?，近岸海域由于岸線影響復(fù)雜多變.渤海中北部以旋轉(zhuǎn)流為主，潮流流向呈順時(shí)針旋轉(zhuǎn)，漲落潮流向以NW—SE向?yàn)橹?秦皇島海域的潮流具有明顯的往復(fù)流特征，漲、落潮方向與岸線平行，漲潮方向?yàn)槲髂戏较颍涑狈较驗(yàn)闁|北方向，漲、落潮平均流速均較小，并且差別不大.空間分布上，不同海域漲落急時(shí)刻并不同步，在三大海灣內(nèi)漲落潮流速變幅較大，而渤海中部區(qū)域?qū)掗煟魉僮兎^小，秦皇島海域處于無(wú)潮點(diǎn)附近，潮流流速變幅更小.總體上一個(gè)潮周期內(nèi)渤海的余流較小.

圖5 2011年6月4日漲急時(shí)刻流場(chǎng)

圖6 2011年6月4日落急時(shí)刻流場(chǎng)

2.2 溢油模型

溢油模型將溢油概化為大量的油粒子，每個(gè)油粒子代表一定質(zhì)量的溢油, 在拉格朗日方法的基礎(chǔ)上增加了蒙特卡羅法計(jì)算溢油擴(kuò)散隨機(jī)性.根據(jù)國(guó)家海洋局調(diào)查報(bào)告[4]，蓬萊19-3(圖1)溢油事故發(fā)生于2011年6月4日，6月21日基本控制溢油，之后較長(zhǎng)時(shí)間內(nèi)仍有少量溢油溢出.2011年7月中下旬在東戴河岸灘發(fā)現(xiàn)長(zhǎng)約4 km，寬約0.5 m呈不均勻帶狀分布的油污帶；在昌黎黃金海岸岸灘高潮線附近發(fā)現(xiàn)長(zhǎng)約1.2 km零星分布的油污帶；在河北唐山淺水灣岸灘高潮線附近發(fā)現(xiàn)長(zhǎng)約500 m，寬約1～1.5 m呈帶狀分布的油污帶，低潮線附近有長(zhǎng)約300 m，寬約1.5～2 m的油污帶.即蓬萊19-3溢油最終在7月中下旬進(jìn)入東戴河至淺水灣沿岸一帶.根據(jù)溢油最終附岸的時(shí)間與位置對(duì)模擬結(jié)果進(jìn)行驗(yàn)證.

2.2.1 基本參數(shù)設(shè)置

模擬對(duì)象為2011年蓬萊19-3溢油事故，作為連續(xù)溢油進(jìn)行考慮，溢油釋放時(shí)間為6月4日—7月31日，溢油總量為500 t，總粒子數(shù)為50 000個(gè).溢油密度參數(shù)采用了蓬萊19-3油田中未熟-低熟油密度值(935 kg·m-3)[18].具體參數(shù)設(shè)置如表1所示，其中蒸發(fā)率為每天的蒸發(fā)部分所占比例.

表1 模型參數(shù)設(shè)置

2.2.2 風(fēng)拖拽系數(shù)、擴(kuò)散系數(shù)校正參數(shù)b率定

溢油運(yùn)動(dòng)包括漂移與擴(kuò)散兩部分：漂移主要由風(fēng)應(yīng)力和潮流引起，其中風(fēng)應(yīng)力起著非常關(guān)鍵的作用；擴(kuò)散主要為紊動(dòng)擴(kuò)散，在初始階段由小規(guī)模的湍流作用引起，當(dāng)粒子“云團(tuán)”充分?jǐn)U散開來(lái)后，大規(guī)模的渦流將對(duì)粒子的擴(kuò)散產(chǎn)生明顯的影響.紊動(dòng)擴(kuò)散所引起的附加位移方向是隨機(jī)的，可用隨機(jī)擴(kuò)散作用來(lái)模擬.在溢油數(shù)值模擬中，風(fēng)拖拽系數(shù)和擴(kuò)散系數(shù)對(duì)溢油的運(yùn)移與擴(kuò)散起著關(guān)鍵性的作用.不同的石油具有不同的組分和性質(zhì)，同一環(huán)境條件下性質(zhì)不同的溢油在漂移、擴(kuò)散等變化過(guò)程中也會(huì)呈現(xiàn)出不同的結(jié)果.對(duì)于不同性質(zhì)的溢油，其風(fēng)拖拽系數(shù)和擴(kuò)散系數(shù)是時(shí)空變化的，但是在一定的范圍內(nèi)變化需要通過(guò)率定確定.

由于風(fēng)在溢油的漂移當(dāng)中起著最重要的作用，風(fēng)拖拽系數(shù)的大小對(duì)于溢油模擬的準(zhǔn)確性也非常關(guān)鍵.在以往經(jīng)驗(yàn)當(dāng)中，風(fēng)拖拽系數(shù)一般在1～3%之間進(jìn)行取值，根據(jù)不同的溢油性質(zhì)進(jìn)行選擇和校準(zhǔn)[14].為了探索渤海風(fēng)拖拽系數(shù)的合理范圍，本次研究設(shè)置了不同的風(fēng)拖拽系數(shù)，根據(jù)溢油的附岸時(shí)間與實(shí)際情況進(jìn)行比較來(lái)確定.從圖7可以看出，在6種風(fēng)拖拽系數(shù)下的模擬結(jié)果均與實(shí)際情況較為接近，當(dāng)Cωd=1.4%～2.3%時(shí)，溢油在7月中下旬進(jìn)入東戴河至淺水灣沿岸一帶，與實(shí)際情況較為相符；而當(dāng)Cωd=1.1%、Cωd=2.6%時(shí)，溢油略微偏早或偏晚進(jìn)入近岸海域.在適當(dāng)范圍內(nèi)對(duì)風(fēng)拖拽系數(shù)進(jìn)行取值對(duì)于溢油的漂移路徑?jīng)]有太大影響；而由于較大的拖拽系數(shù)下溢油的漂移速度較快，在漂移相同的距離下向兩邊擴(kuò)散的范圍較小，故溢油的分布形狀較為狹長(zhǎng).根據(jù)模擬結(jié)果與實(shí)際情況的對(duì)比，對(duì)密度大、粘度高的蓬萊19-3溢油進(jìn)行溢油模擬時(shí)，風(fēng)拖拽系數(shù)應(yīng)取1.4%～2.3%.以下蓬萊19-3溢油模擬計(jì)算中風(fēng)拖拽系數(shù)取值為Cωd=2%.

b Cωd=1.4%

c Cωd=1.7%

d Cωd=2.0%

e Cωd=2.3%

f Cωd=2.6%

紊動(dòng)擴(kuò)散系數(shù)式(6)中含兩個(gè)可變參數(shù)a和b.根據(jù)式(6)可知：參數(shù)a的影響較小，可取默認(rèn)值1.0；對(duì)于時(shí)間跨度較長(zhǎng)的溢油模擬，參數(shù)b起主要作用.本文主要針對(duì)參數(shù)b的影響進(jìn)行探討，通過(guò)設(shè)置不同的參數(shù)b進(jìn)行模擬與對(duì)比，根據(jù)溢油的附岸范圍和實(shí)際情況進(jìn)行比較以確定合適的參數(shù)范圍.模擬結(jié)果如圖8所示，結(jié)果表明：參數(shù)b的改變對(duì)溢油的漂移速度和路徑影響不大，對(duì)溢油的擴(kuò)散范圍影響較大.當(dāng)b=0.38～0.42時(shí)，溢油于東戴河至淺水灣沿岸一帶附岸，與實(shí)際情況相符；當(dāng)b=0.3,0.35時(shí)，淺水灣沿岸未出現(xiàn)溢油附岸，溢油擴(kuò)散范圍偏??；當(dāng)b=0.45時(shí)，溢油附岸范圍兩端分別略微超出了東戴河與淺水灣，擴(kuò)散范圍偏大.因此，在利用本模型對(duì)密度大、粘度高的蓬19-3溢油進(jìn)行模擬時(shí)，參數(shù)b合理的取值范圍為0.38～0.42.以下蓬萊19-3溢油模擬取b=0.4.經(jīng)計(jì)算，時(shí)間步長(zhǎng)設(shè)為1 h，則1 h后紊動(dòng)擴(kuò)散系數(shù)D=26.5，則最大擴(kuò)散位移ΔS=756 m；24 h后D=94.3，則ΔS=1 427 m；10 d后D=236.9，ΔS=2 262 m；50 d后D=451.0，ΔS=3 122 m.模擬期間紊動(dòng)擴(kuò)散系數(shù)取值大概在26.5～500范圍內(nèi).由于擴(kuò)散位移的方向是隨機(jī)選取的，故擴(kuò)散位移主要使擴(kuò)散范圍逐漸增大，對(duì)溢油質(zhì)心位移的影響較小.可見(jiàn)，溢油后期的紊動(dòng)擴(kuò)散系數(shù)增大，油粒子的擴(kuò)散位移也明顯增大.這樣大渦促使溢油后期迅速擴(kuò)散的作用便得以體現(xiàn)出來(lái).

a b=0.3

b b=0.35

c b=0.38

d b=0.4

e b=0.42

f b=0.45

3 模擬結(jié)果與分析

模擬采用風(fēng)場(chǎng)由美國(guó)NOAA提供的9個(gè)渤海及臨近區(qū)域站點(diǎn)的風(fēng)場(chǎng)數(shù)據(jù)，其中渤海中部站點(diǎn)(120°E,39.047°N)的風(fēng)速風(fēng)向如圖9所示，平均風(fēng)速為4.10 m·s-1，最大風(fēng)速為12.04 m·s-1；常風(fēng)向?yàn)槟舷颍舷蝻L(fēng)占總數(shù)的25.4%，平均風(fēng)速為5.02 m·s-1；南偏東、南、南偏西風(fēng)向占總數(shù)的57.3%，平均風(fēng)速為4.61 m·s-1.

在前面參數(shù)率定(即風(fēng)拖拽系數(shù)Cωd=2%、擴(kuò)散系數(shù)校正參數(shù)b=0.4)的基礎(chǔ)上，對(duì)蓬萊19-3溢油事故進(jìn)行模擬，結(jié)果如圖10.溢油發(fā)生后擴(kuò)散范圍逐漸增大，前端擴(kuò)散面積較大，平面分布密度較小，約為0.000 02 kg·m-2.溢油平臺(tái)周圍及其北偏西方向溢油分布密度較大，大約為0.000 05～0.000 5 kg·m-2.結(jié)合圖9與圖10進(jìn)行分析可知：6月22日以前以南風(fēng)為主，溢油運(yùn)移方向與風(fēng)向大致相同為正北方向；6月22日—7月4日，風(fēng)向以東南偏東向和南向?yàn)橹鳎缬椭饾u轉(zhuǎn)為向西北方向運(yùn)移，并靠近西北沿岸一帶；7月4日—7月16日，風(fēng)向以南向?yàn)橹?，溢油往北運(yùn)移進(jìn)入東戴河至黃金海岸沿岸一帶；7月16日—7月22日，風(fēng)向以南偏東向?yàn)橹?，溢油往西北方向運(yùn)移，附岸范圍向南擴(kuò)展至淺水灣沿岸；隨后風(fēng)向以南向?yàn)橹?，溢油擴(kuò)散范圍略微往北偏移.總體上，溢油油污于7月中下旬進(jìn)入淺水灣至東戴河沿岸一帶，其中東戴河與黃金海岸沿岸一帶受到的溢油污染最為嚴(yán)重，淺水灣受到的污染較輕，模擬結(jié)果與實(shí)際情況相符.

渤海海域夏季以南風(fēng)和東南風(fēng)為主，風(fēng)速較大且風(fēng)向較為穩(wěn)定，溢油的漂移方向與主風(fēng)向大致相同，表明風(fēng)對(duì)漂浮于水面上的溢油的運(yùn)移起著很明顯的作用.然而，渤海海域的潮流流速不大，尤其是秦皇島海域處于無(wú)潮點(diǎn)區(qū)域，潮流流速更小.渤海潮流以往復(fù)流和旋轉(zhuǎn)流為主，漲潮潮流與落潮潮流流向大致相反，且流速差異不明顯.一個(gè)潮周期內(nèi)溢油在漲落潮流的作用下會(huì)來(lái)回?cái)[動(dòng)，余流小、凈位移小.因此，在渤海海域內(nèi)潮流對(duì)溢油的漂移不及風(fēng)的作用.

圖9 風(fēng)速風(fēng)向過(guò)程 (2011-06-04—7.31)

a 6月14日

b 6月24日

c 7月4日

d 7月14日

e 7月24日

f 7月31日

溢油模型計(jì)算過(guò)程中，溢油主要分為漂浮組分、蒸發(fā)組分和附岸組分.圖11為整個(gè)溢油模擬過(guò)程中各組分的質(zhì)量變化圖.其中溢油總量設(shè)定為線性增長(zhǎng)，16 d(6月19日)基本控制溢油，溢油總量達(dá)290 t；后期依然持續(xù)有少量溢油溢出.在溢油前期，溢油漂浮組分和蒸發(fā)組分占溢油總量的絕大部分.漂浮組分在初始階段隨著溢油的發(fā)生呈線性增長(zhǎng)，36 d(7月10日)達(dá)到最大值204 t，約占溢油總量的50.2%.之后隨著漂浮于水面的溢油逐漸發(fā)生蒸發(fā)與附岸使其逐漸下降，第56 d(7月31日)，漂浮組分下降至約155.0 t，占溢油總量30.6%.溢油的蒸發(fā)主要發(fā)生于溢油初期，由于前25 d溢油量較大，故蒸發(fā)量增長(zhǎng)較快，于第25 d便達(dá)到約160 t，約占溢油總量48.2%；之后隨著溢油揮發(fā)組分減少，且乳化過(guò)程使溢油含水率增大，使得蒸發(fā)在溢油的后期階段速度減緩，第57 d(7月31日)蒸發(fā)量達(dá)到最大值約為223 t，占溢油總量的44.0%.其中，第22 d(6月26日)出現(xiàn)了12.04 m·s-1的最大風(fēng)速，溢油的蒸發(fā)出現(xiàn)了短暫的陡增，漂浮組分也隨之出現(xiàn)驟減，可見(jiàn)風(fēng)對(duì)促進(jìn)溢油蒸發(fā)也有一定作用.至第35 d(7月9日)之后，溢油進(jìn)入沿岸海域并逐漸發(fā)生附岸，于7月31日達(dá)126 t左右，占溢油總量的24.9%.

圖11 溢油組分質(zhì)量變化圖

4 結(jié)論

通過(guò)Delft3D中的FLOW模塊建立了渤海水動(dòng)力模型，在水動(dòng)力模型的基礎(chǔ)上利用PART模塊耦合建立了渤海海上溢油模型.率定了該溢油模型中風(fēng)拖拽系數(shù)、擴(kuò)散系數(shù)校正參數(shù)b的合理范圍，對(duì)2011年蓬萊19-3溢油事故進(jìn)行了模擬與驗(yàn)證，并對(duì)溢油運(yùn)移與歸宿進(jìn)行了相關(guān)探討.得出以下結(jié)論：

(1) 渤海海域潮流以往復(fù)流和旋轉(zhuǎn)流為主，其中渤海中北部以旋轉(zhuǎn)流為主，呈順時(shí)針旋轉(zhuǎn)，漲落潮流向以NW-SE向?yàn)橹?渤海海域潮流流速較小，且漲落潮流速相差小，變幅也小.其中，秦皇島海域處于無(wú)潮點(diǎn)附近，潮流流速更小.

(2) 對(duì)于密度大粘度高的蓬萊19-3溢油風(fēng)拖拽系數(shù)合理范圍為1.4～2.3%；擴(kuò)散系數(shù)公式中校正參數(shù)b對(duì)溢油的擴(kuò)散范圍影響較大，其合理取值范圍為0.38～0.42.

(3) 渤海夏季風(fēng)速較大、風(fēng)向較為穩(wěn)定的南風(fēng)和東南風(fēng)對(duì)漂浮于水面上的溢油作用較大，且溢油漂移方向與主風(fēng)向大致相同.而渤海潮流以往復(fù)流和旋轉(zhuǎn)流為主，溢油在漲落潮流的作用下會(huì)出現(xiàn)來(lái)回?cái)[動(dòng)，凈位移小，潮流對(duì)溢油的漂移不及風(fēng)的作用.

(4) 溢油初期的蒸發(fā)速度較快；由于揮發(fā)組分的減少以及乳化作用使溢油含水率增大，溢油后期的蒸發(fā)速度明顯減小.進(jìn)入沿岸海域后溢油容易發(fā)生附岸，附岸和蒸發(fā)使得溢油漂浮組分逐漸減少.最終，漂浮組分、蒸發(fā)組分和附岸組分分別占溢油總量的30.6%、44.0%和24.9%，蒸發(fā)和附岸是溢油的主要?dú)w宿.

[1] 張麗萍,王輝.不同海面狀況海洋石油污染處理方法優(yōu)化配置[J].海洋技術(shù),2006,25(3):1.

ZHANG Liping, WANG Hui. Optimization design of countermeasures to cleanup the marine oil pollution in different sea surface conditions.[J].Ocean Technology,2006,25(3):1.

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

XU Yandong. A study on theweathering and predicting model of spilled oils at sea[D]. Qingdao: Ocean University of China,2006.

[3] DENG Zenggan, YU Ting, JIANG Xiaoyi,etal. Bohai Sea oil spill model: a numerical case study[J]. Marine Geophysical Research.2013,34(2):115.

[4] 國(guó)家海洋局.蓬萊19-3油田溢油事故聯(lián)合調(diào)查組關(guān)于事故調(diào)查處理報(bào)告[R].北京:國(guó)家海洋局,2011[2015-02-03]. http://www.soa.gov.cn/xw/hyyw_90/201211/t20121109_884.html.

State Oceanic Administration People’s Republic of China. The joint investigation group report on the investigation and handling of the Penglai 19-3 oil spill accident[R].Beijing: State Oceanic Administration People’s Republic of China, 2011[2015-02-03]. http://www.soa.gov.cn/xw/hyyw_90/201211/t20121109_884.html.

[5] KUANG Cuiping, PAN Yi, HU Ying,etal. Primary analysis of the 2011ConocoPhillips oil spill in Bohai Bay, China [C]// Advances in Biomedical Engineering——2012 Asia Pacific Conference on Environmental Science and Technology (APEST 2012). Kuala Lumpur: Information Engineering Research Institute,2012:1-5.

[6] DENG Zenggan, YU Ting, SHI Suixiang,etal. Numerical study of the oil spill trajectory in Bohai Sea, China[J]. Marine Geodesy,2013,36(4): 351.

[7] XU Qing, LI Xiaofeng, WEI Yongliang,etal. Satellite observations and modeling of oil spill trajectories in the Bohai Sea[J]. Marine Pollution Bulletin,2013,71(1/2):107.

[8] WANG Shoudong, SHENYongming, ZHENG Y.H. Two-dimensional numerical simulation for transportand fate of oil spills in seas[J]. Ocean Engineering,2005,32(13):1556.

[9] WANG Shoudong, SHEN Yongming, GUO Yakun,etal. Three-dimensional numerical simulation for transport of oil spills in seas[J]. Ocean Engineering,2008,35(5): 503.

[10] WANG Jinhua, SHEN Yongming. Development of an integrated model system to simulate transport and fate of oil spills in seas[J].Sciences China-Technological Sciences,2010,53(9): 2423.

[11] BI Haipu, SI Hu. Numerical simulation of oil spill for the Three Gorges Reservoir in China[J]. Water and Environment Journal,2014,28(2):183.

[12] LI Yan, ZHU Jiang, WANG Hui. The impact of different vertical diffusion schemes in a three-dimensional oil spill model in the BohaiSea[J]. Advances in Atmospheric Sciences,2013, 30(6): 1569.

[13] WL|DelftHydraulics. User manual Delft3D-FLOW[M].Delft: WL|Delft Hydraulics, 2013.

[14] WL|DelftHydraulics. User manual D-WaqPART[M].Delft: WL|Delft Hydraulics, 2013.

[15] Fingas M, Filedhouse B, Mullin J. Water-in-oil emulsions results of formation studies and applicability to oil spill modeling[J].Spill Science & Technology，1999，5(1): 81.

[16] Mackay D, Paterson S, Trudel K. A mathematical model of oil spill behavior on water with natural and chemical dispersion. [R]. Ottawa: Environment Canada, 1977.

[17] Fingas M F. Chemistry of oil and modelling of spills[J].Journal of Advances in Marine Technology Conference，1994, 11: 41.

[18] Willmott C J. On the validation of models[J].Physical Geography,1981, 2(2):184.

[19] 王廣源,周心懷,王昕,等.蓬萊19-3/25-6油田未熟-低熟油特征與成因[J].石油實(shí)驗(yàn)地質(zhì), 2014,36(2):231.

WANG Guangyuan, ZHOU Xinhuai, WANG Xin,etal. Characteristics and genesis of immature and low-matire oils in Penglai 19-3/25-6 oilfields[J]. Petroleum Geology & Experiment, 2014, 36(2):231.

Numerical Simulation and Analysis of Transport and Fate of Penglai 19-3 Oil Spill

KUANG Cuiping1, XIE Hualang1, SU Ping2, CHEN KUO1

(1. College of Civil Engineering, Tongji University, Shanghai 200092, China； 2. Guangxi Communications Planning Surveying and Designing Institute, Nanning 530029, China)

The numerical model of oil spill transport and fate in the Bohai Sea is set up by coupling with the hydrodynamic model. The change processes of spilled oil included in the model are drift, dispersion, evaporation, emulsification and so on. The drift process is simulated by the Lagrange method, while the turbulent dispersion process is considered by using the ‘Monte Carlo method’, turbulent dispersion of spilled oil changes over time by using the time-varying turbulent dispersion coefficient. The reasonable ranges of wind drag coefficient and the dispersion correction parameterbfor Penglai 19-3 oil spill model obtained by calibration are 1.4～2.3% and 0.38～0.42, respectively. Then, the Penglai 19-3 oil spill accident in the Bohai Sea is used to verify the oil spill model by comparing the time and location of spilled oil arrival at shore, and the simulation results are consistent with the actual situation. The numerical results show that the wind, because of its large velocity and small change in direction, plays a major role on the drift of spilled oil floating on the water surface, and the drift direction of spilled oil is roughly the same with the main direction of wind. The effect of the tidal current, which is dominated by reversing current and rotary current with a very weak residual current, to the drift of spilled oil in Bohai Sea is less than the effect of the wind. The evaporation rate of spilled oil is quite large in the initial stage after oil spill, but significantly reduced in the later stage as a result of the reduction of the volatile fraction and the increase of the water content caused by emulsification. The reduction of floating oil is mainly caused by evaporation and sticking that are also the main fates of spilled oil.

Bohai Sea； oil spill； transport and dispersion； numerical simulation

2015-02-06

國(guó)家海洋公益項(xiàng)目(201305003)

匡翠萍(1966—)，女，教授，博士生導(dǎo)師，工學(xué)博士，主要研究方向?yàn)楹涌诤０豆こ毯铜h(huán)境.E-mail:cpkuang@#edu.cn

X55; P76

A