危險化學(xué)品海上漂移擴散數(shù)值模擬研究綜述

郭凱旋

（國家海洋環(huán)境預(yù)報中心，北京100081）

1 引言

隨著我國沿海經(jīng)濟的快速發(fā)展，臨海工業(yè)及港口物流規(guī)模不斷增加，海上危險化學(xué)品運輸需求量越來越大，危險品的種類和數(shù)量也逐年增加。截至2018年12月31日，我國沿海省際運輸化學(xué)品船（含油品和化學(xué)品兩用船，下同）共計288艘、112.9萬載重噸，同比增加16艘、6.72萬載重噸，增幅6.33%[1]；全國萬噸級及以上泊位中，專業(yè)化液體化工泊位217個，比上年增加12個[2]。

自1990年以來我國發(fā)生了多起海上?；沸孤┦鹿剩?997年，Blue Sky No.2散化船在杭州以東200 km洋面沉沒，船上載有的988 t酞酸二辛酯泄漏入海；2001年，韓國籍散化船“大勇”輪在長江口外海域和香港散貨船“大望”輪相撞，造成630 t苯乙烯泄漏，給海洋生態(tài)環(huán)境造成了嚴(yán)重損害；2011年，渤海蓬萊19-3油田發(fā)生溢油事故，造成油田周邊及其西北部面積約6 200 km2的海域海水污染(超第一類海水水質(zhì)標(biāo)準(zhǔn))，其中870 km2海水受到嚴(yán)重污染(超第四類海水水質(zhì)標(biāo)準(zhǔn))；2012年，韓國籍化學(xué)品船“雅典娜”輪在廣東汕尾海域沉沒，船上7 000 t濃硫酸及140 t燃料油發(fā)生泄漏；2018年，巴拿馬籍油船“桑吉”輪與香港散裝貨船“長峰水晶”輪在長江口東約160 n mile處發(fā)生碰撞事故，“桑吉”輪爆燃沉沒，船上的十多萬噸凝析油和大量燃油發(fā)生泄漏，據(jù)國家海洋局的消息，沉海區(qū)域有長約12 km，寬約9 km的油污帶，面積約58 km2，嚴(yán)重影響了東海海洋生態(tài)環(huán)境；2018年，寧波舟山“天桐1號”油輪靠泊泉港東港石化公司碼頭擬接運工業(yè)用裂解碳九，由于操作員違規(guī)操作，造成69.1 t碳九泄漏，雖通過吸油氈回收了約40 t泄漏物，但是剩余的大部分泄漏入海，給當(dāng)?shù)睾Ｑ笊鷳B(tài)環(huán)境及海水網(wǎng)箱養(yǎng)殖造成影響。

海上危化品泄漏事故與其他海洋環(huán)境污染事故相比具有突發(fā)性強、影響范圍廣、影響時間長和泄漏后變化復(fù)雜等特點[3]。突發(fā)性強表現(xiàn)為海上?；沸孤┐蟛糠侄际窃跊]有任何前兆的情況下突然發(fā)生的，瞬間排放大量的污染物，給海洋環(huán)境造成巨大影響。影響范圍廣表現(xiàn)為少量的危化品發(fā)生泄漏，其產(chǎn)生的影響范圍將非常巨大，例如1 t氯的泄漏能夠影響4.8 km2的范圍[3]。影響時間長表現(xiàn)為危化品泄漏對海洋生態(tài)環(huán)境的影響不僅僅發(fā)生在泄漏時，部分?；穼Ｑ笊鷳B(tài)環(huán)境的污染極難消除，例如福島核泄漏事故污染物CS-137的半衰期長達30.07 a[4]。泄漏后變化復(fù)雜表現(xiàn)為?；贩N類繁多且多樣的行為方式造成?；沸孤┤牒：笪锢砗突瘜W(xué)形態(tài)變化十分復(fù)雜。

海上?；沸孤┦鹿手幸缬褪鹿首顬轭l繁，因此英、美等發(fā)達國家很早就開展了海上溢油事故相關(guān)理論研究，并通過大量的試驗驗證，建立了眾多預(yù)測模型[5-8]。與溢油相比，其他類型的海上?；沸孤┢茢U散數(shù)值模擬研究相對較少[9-12]。本文將對近年來國內(nèi)外海上危化品泄漏在海洋水體環(huán)境中漂移擴散數(shù)值模擬方面的研究進行整理、分類和介紹。

2 ?；沸孤?shù)值模擬水動力模式

水動力數(shù)值模擬是研究海上?；沸孤┢茢U散的基礎(chǔ)。三維模型由于能夠全面而立體地反映實際河口海岸和近海環(huán)境中海水的運動狀況和基本規(guī)律，因此在海洋水動力數(shù)值模擬研究中得到了廣泛研究和應(yīng)用。伴隨著計算機的飛速發(fā)展以及數(shù)值模擬技術(shù)的不斷完善，目前三維模型已經(jīng)成為水動力數(shù)值模擬的主要發(fā)展方向[13-14]，其中研究應(yīng)用較多的模型有POM（Princeton Ocean Model)、FVCOM（Finite-Volume Coastal Ocean Model）和ROMS（Regional Ocean Modeling System）。

2.1 POM

POM是由美國普林斯頓大學(xué)的Blumberg和Mellor于1977年建立起來的一個三維斜壓原始方程數(shù)值海洋模式[15-17]。POM采用的蛙跳有限差分格式使得模式具有很高的垂向分辨率[18-20]。POM采用時間分裂算法，模式的外模方程是二維的，內(nèi)模方程是三維的，內(nèi)外模分離技術(shù)比完全三維計算節(jié)省很大計算量[18-20]。POM在水平方向采用正交曲線網(wǎng)格，變量空間配置使用“Arakawac C”網(wǎng)格，能更好地擬合岸線邊界，減少“鋸齒”效應(yīng)；在垂向上采用σ坐標(biāo)變換，可體現(xiàn)不規(guī)則海底地形的變化特點；通過干濕網(wǎng)格動邊界技術(shù)，不僅能很好地處理復(fù)雜地形水域的模擬問題，而且可更好地解決三維水動力環(huán)境中大量淺灘的“干出”與“淹沒”等難點問題[15，19-20]。于海亮[18]以POM為基礎(chǔ)并基于Lagrange粒子追蹤方法開發(fā)了模擬海上溢油的三維輸運模型，成功模擬了渤海海域1990年的一次溢油事故。李連峰[19]在POM模式的基礎(chǔ)上建立了三維海上類油型化學(xué)品及沉降型化學(xué)品輸運模型，該模型采用Langevin方程與Fokker-Planck方程相結(jié)合來模擬油粒子的輸運過程，有效地解決了具有隨機性的粒子運動過程。通過在寧波-舟山港海域的應(yīng)用和模擬案例的計算校驗，該模型顯示出良好的效果。

2.2 FVCOM

FVCOM是由陳長勝（美國麻州大學(xué)海洋科技學(xué)院）于2000年成功建立的一套三維、有限體積、自由表面和非結(jié)構(gòu)網(wǎng)格的海洋環(huán)流與生態(tài)數(shù)值模型[20-24]。FVCOM在水平方向上采用非結(jié)構(gòu)化三角形網(wǎng)格，能夠更精確地擬合復(fù)雜曲率的岸界。FVCOM借鑒了許多POM模型的優(yōu)點，例如：垂直方向采用σ坐標(biāo)用于擬合復(fù)雜的海底地形；使用干濕網(wǎng)格判別法處理近岸灘涂演變問題；采用模式分裂法求解方程大大提高計算效率[22，25]。國內(nèi)學(xué)者藏士文[22]、徐國懷[25]在模擬海上危化品泄漏漂移擴散的研究中分別采用FVCOM建立了渤海和寧波-舟山海域的水動力模型，將實際觀測資料與模擬結(jié)果進行對比驗證，發(fā)現(xiàn)潮位、流速和流向的模擬結(jié)果與實際觀測情況吻合性較高、趨勢一致且誤差范圍較小。

2.3 ROMS

ROMS是一個三維自由表面非線性原始方程近海區(qū)域模式，是在垂向靜壓近似和Boussinesq假定下求解自由表面下Reynolds平均的Navier-Stokes方程[14，26-28]。ROMS水平方向采用正交曲線網(wǎng)格,垂直方向采用σ坐標(biāo)系統(tǒng)[14，27]。為了提高計算效率節(jié)約計算時間，ROMS使用經(jīng)過正壓（快）和斜壓（慢）模式之間的特殊處理和耦合的分離顯式時間步方案求解動量流體靜力原始方程[28-29]。ROMS包含非線性計算內(nèi)核算法等多種算法，并有MY-2.5混合參數(shù)方案等多種垂直混合參數(shù)化方案供選擇[28-30]。由于ROMS具有集合預(yù)報功能，可與多種海洋及氣象模式進行耦合使用，因此，近年來被廣泛應(yīng)用于近海海洋環(huán)流、海洋生態(tài)環(huán)境以及海冰等領(lǐng)域的模擬研究[27，29]。崔可夫[14]基于ROMS模式建立了渤海灣三維潮流模型，數(shù)值模擬結(jié)果在水位、流速和流向方面均與實際觀測資料吻合良好，為渤海灣水體交換和污染物遷移擴散研究提供了可靠的水動力參數(shù)。陳超[27]運用ROMS模型，對大連灣危化品輸運過程進行精細化模擬，模擬結(jié)果和實測值吻合較好。

3 海上?；菲茢U散數(shù)值模型

在海洋水動力數(shù)值模式基礎(chǔ)上，通過引入?；愤\動擴散方程，發(fā)展形成了適用于不同?；奉愋偷暮Ｉ衔；菲茢U散預(yù)測數(shù)值模型。海上危險化學(xué)品一旦發(fā)生泄漏，?；穼⒃诤Ａ骱惋L(fēng)的作用下進行平流、擴散、揮發(fā)、溶解和沉降（或上升）等理化過程[31]。海上危險化學(xué)品泄漏后在海中的輸移和擴散模式主要受兩方面影響。一方面是外部環(huán)境，如地理環(huán)境、水文因素、氣象因素和其他因素[16]。?；吩诤Ａ骱惋L(fēng)的影響下將進行平流擴散和湍流擴散：平流擴散指海水的整體運動引起污染物有規(guī)律的平移；湍流擴散指因風(fēng)和潮流引起的隨機無序的流體運動,運動狀態(tài)只能通過概率統(tǒng)計方法加以描述。擴散類型由海流類型決定,并受潮流、風(fēng)生環(huán)流、密度流以及水深和岸線狀況等諸多因素制約[32-33]。在海洋這樣一個大環(huán)境中,污染物一旦泄漏，不僅是單一的平流擴散或湍流擴散，而是兩種擴散形式同時作用于污染擴散的整個過程，只不過作用的程度不同[33]。另一方面是化學(xué)品的主要物理化學(xué)性質(zhì)，其中最主要的3個因素是溶解度、揮發(fā)性和密度。溶解度是決定?；窋U散運動形式的主要因素，是對其他性質(zhì)進行分析的前提；揮發(fā)性決定了?；啡牒：笫欠褚哉魵庑问较蚩諝庵袛U散，?；返膿]發(fā)性用20℃飽和蒸氣壓衡量，通常以3 kPa作為揮發(fā)與不揮發(fā)的界限；密度（與水的相對密度）則決定?；啡牒：笫窃谒孢\動還是進入水體甚至沉降至海底[34]。根據(jù)?；沸孤┤牒：蟮倪\動形式，我們將海上?；贩譃槿芙鈹U散型、海面漂移型、沉降型和揮發(fā)型4種類型。

3.1 溶解擴散型

溶解擴散型?；啡苡诤Ｋ?，并在海水中以三維形式進行擴散。對于溶解型?；吩诤Ｋ械倪w移擴散研究，目前主要從數(shù)值優(yōu)化求解污染物對流擴散方程這一方向來解決。國內(nèi)外研究求解對流擴散方程的數(shù)值方法大致可分兩類：一類是歐拉法，將研究海域的空間點作為研究對象，利用解析解的方法或者數(shù)值求解的方法，求解危化品在該海域的連續(xù)性方程、運動方程和輸運方程，從而得出?；窛舛入S時間的變化特征；另一類是拉格朗日法，該方法將研究海域中的?；犯呕癁榫哂幸欢〝?shù)量和質(zhì)量的質(zhì)點顆粒，追蹤每個質(zhì)點粒子在研究海域流場中的運動軌跡，得到每個時刻和每個質(zhì)點顆粒所處的空間位置[35]。

對流擴散方程為：

式中：C為?；窛舛龋籾x、uy和uz分別為X、Y和Z方向的流速；Dx、Dy和Dz分別為X、Y和Z方向的擴散系數(shù)；S0為?；返脑春蛥R。

在水平運動尺度遠遠大于垂向尺度和?；反瓜蚧旌媳容^均勻時，為求解方便往往忽略垂向濃度變化，將對流擴散方程簡化為二維模型，以此提高計算效率。杜海濤[33]在對大連灣溶解、輕質(zhì)和保守液體化學(xué)品溢漏后歸宿行為的研究中使用了簡化后的二維模型，得出的模擬結(jié)果符合大連灣實際情況。伴隨著計算機性能的提高，越來越多的學(xué)者使用三維模型進行研究，于海亮[15]在研究海上液體?；沸孤┤S污染擴散的過程中建立了σ坐標(biāo)系下的溶解型?；啡S輸運模型，并對模型中的對流項分別采用中心差分和一階Smolarkiewicz迎風(fēng)格式離散進行優(yōu)化求解，結(jié)果表明一階Smolarkiewicz迎風(fēng)格式和改進的二維有限體積法相比中心差分格式具有一定的優(yōu)勢。

3.2 海面漂移型（類油型）

海面漂移型?；凡蝗苡诤Ｋ诤Ｃ嬷饕远S形式漂移擴散。人們對于這類?；返难芯恐饕性诤Ｉ弦缬停虼撕Ｃ嫫菩陀挚珊喎Q“類油型”。海上溢油的漂移擴展過程主要分慣性力擴展、粘性力擴展和表面張力擴展3個階段進行研究，早期具有代表性的是Fay公式[15-16，36-37]。

慣性力擴展階段：

粘性力擴展階段：

表面張力擴展階段：

式中：D為油膜擴散直徑（單位：m）；g為重力加速度（單位：m/s2）；V為溢油總體積（單位：m3）；t為從泄漏開始計算的時間（單位：s）；β=1-ρ0/ρw，ρ0和ρw分別為?；访芏群秃Ｋ芏龋▎挝唬?03kg/m3）；Vw為海水運動粘性系數(shù)（單位：m2/s）；δ為擴展系數(shù)；K1、K2和K3分別為3個階段的比例系數(shù)。

泄漏后的油品在海水中形成油膜并在海水運動的作用下產(chǎn)生漂移，同時油膜擴散范圍逐漸擴大,油膜漂移距離的長短通常用油膜等效圓中心位移來判斷[15]。其位移s為：

Fay理論及其修正模型對于在開闊海域瞬時泄漏情況下的類油型污染物模擬取得了很好的結(jié)果，適合類油型?；吩谛孤┏跗陔A段的模擬預(yù)測[38]。連續(xù)泄漏油膜在距離泄漏源一定距離后，側(cè)向擴散起主導(dǎo)作用，而對于水下泄漏，浮力的作用至關(guān)重要，由于Fay模式忽略這些作用，所以它并不適用于連續(xù)泄漏和水下泄漏的模擬[39]。

Johanseen等[39-43]提出的“油粒子”模型是當(dāng)今世界主流的溢油模式。該模型把溢油分散成有限個油粒子來模擬實際溢油的漂移擴散過程。“油粒子”模型不僅解決了溢油在重力擴展停止后的擴散現(xiàn)象，還實現(xiàn)了海上油膜破碎分離現(xiàn)象的模擬，對溢油在海洋環(huán)境中受到的海洋動力因素的影響模擬得更加準(zhǔn)確[42-44]。國內(nèi)學(xué)者Yang等[45]基于國家海洋環(huán)境預(yù)報中心自主研發(fā)的溢油預(yù)報模型，采用Lagrange追蹤法模擬油粒子的三維時空運動，建立了渤海灣三維溢油預(yù)報模型，較好地重現(xiàn)了2011年6—8月渤海蓬萊19-3溢油事故的發(fā)展過程。Li等[46]基于“油粒子”模型重現(xiàn)了“桑吉”輪事故溢油的漂移擴散過程，模擬結(jié)果與實際觀測相吻合。雖然“油粒子”模型代表了當(dāng)今溢油模擬的方向，但是其缺點也很明顯：所需的模擬粒子數(shù)過于龐大，導(dǎo)致計算量非常大；忽略了溢油自身的重力擴展作用。黃娟等[47]和劉偉峰等[43]研究發(fā)現(xiàn)在溢油大規(guī)模泄漏初期，由于油膜面積在短時間內(nèi)急劇擴大，溢油自身的擴展效應(yīng)顯著大于湍流擴散效應(yīng)，如果僅通過水平擴散的方式來模擬油膜擴展過程，將導(dǎo)致“油粒子”擴散速度響應(yīng)過慢，不能真實反映溢油量對擴散面積的影響。

3.3 沉降型

沉降型化學(xué)品主要指密度比海水大且不溶或微溶的化學(xué)品，入水后分散成小液滴在浪、潮和流的作用下隨水體運動，其運動形式基本可分為擴散和漂移、沉降、沉積和再懸浮[16，48-49]。

海洋水動力要素是影響?；吩诤Ｋ袛U散和漂移過程的主要因素。擴散是指?；沸孤┤牒：笥捎跐舛炔詈屯牧鞯任锢硪蛩囟纬晌廴疚锵蛘麄€水體遷移的混合過程；漂移指泄漏入海的?；吩陲L(fēng)、表層和次表層流以及波浪作用下的拉格朗日漂移過程，漂移運動主要由平流條件決定[16，48]。陳協(xié)明[48]在研究化學(xué)品的溢漏擴散模式時指出，沉降型化學(xué)品在海面上的漂移主要受風(fēng)的切應(yīng)力、表層/次表層流和余流（波生余流和潮余流）控制。

沉降過程一般可以分為3個階段：沉降射底、觸底和底浪的傳播[16]。?；啡胨蟀l(fā)生破碎并迅速下降，形成了充分成長的高速射流。在這一過程中，周圍海水與?；钒l(fā)生混合，沉降射流的體積不斷擴大；沉降射流觸底時將形成以碰撞點為中心向四周擴展的高密度底浪；開始時底浪的傳播速度很快，隨著傳播距離的增加，底浪的厚度和傳播速度迅速減少，底浪中攜帶的?；费杆傧陆挡⒊练e[16，48]。

危化品沉降到海底后會形成沉積，在流與浪產(chǎn)生的剪切力的共同作用下，沉積的?；窌l(fā)生再懸浮，并回到水體中重新進行擴散、漂移和沉降過程[16，48-49]。沉積物的再懸浮過程十分復(fù)雜，目前還沒有成熟的預(yù)測方法，國外學(xué)者McNeil等[50-52]使用水槽裝置模擬水體底泥的再懸浮過程，國內(nèi)學(xué)者陳協(xié)明[48]和李連峰[19]使用經(jīng)驗參數(shù)（沉積物的5%）計算再懸浮通量。

3.4 易揮發(fā)型

根據(jù)揮發(fā)性?；沸孤┖蠡旌蠚庠婆c環(huán)境大氣的密度差異，將混合氣云分為浮性氣云（密度比空氣?。?、中性氣云（密度與空氣相近）和重質(zhì)氣云（密度比空氣大）3類。研究學(xué)者往往將浮性氣云和中性氣云歸為一類進行研究[53]，簡稱為“非重氣云”。非重氣云團擴散模型中最具代表性的是高斯（Gauss）模型，其理論基礎(chǔ)是湍流擴散梯度輸送理論，忽略重力對氣體擴散的影響，認為擴散主要由空氣的湍流決定，污染物在擴散截面的濃度呈正太分布，擴散系數(shù)K為常數(shù)[5，54]。高斯模型由于所需參數(shù)少、運算量小和計算快捷等優(yōu)點，被廣泛應(yīng)用于小尺度的污染物模擬。重氣云團的擴散比中性和浮性氣云復(fù)雜的多，主要經(jīng)歷閃蒸和空氣夾帶、密度差作用下的重量沉降、重力沉降的地表作用和被動擴散4個階段[54]。重氣云團的擴散研究不僅要考慮氣體重力因素還要考慮重氣云流動擴散阻力等因素的影響。孫莉等[6]通過收集大量重氣云擴散實驗結(jié)果加以分析整理，解析出能較好反應(yīng)重氣云瞬時和連續(xù)釋放規(guī)律的方程式，并成功建立了唯象模型（BM模型）。BM模型是最簡單的根據(jù)經(jīng)驗判斷的重氣云擴散模型，可以根據(jù)重氣云擴散曲線圖和方程式快速得出相對應(yīng)的擴散數(shù)據(jù)，因此被廣泛應(yīng)用。

計算機仿真模擬是研究易揮發(fā)型危險化學(xué)品泄漏擴散的又一熱點領(lǐng)域。計算機仿真模型考慮的因素較多，模擬時間較長，因此預(yù)測精確度很高。此類模型最具代表性的就是CFD（Computational Fluid Dynamics）和ALOHA（Areal Locations Of Hazardous Atmospheres）[55]。CFD具有可視化能強和計算方法成熟等特點[43],在實際中被廣泛應(yīng)用。Ikealumba等[56]在應(yīng)用CFD研究海上液化天然氣的泄漏擴散規(guī)律時發(fā)現(xiàn)，海洋波動帶來的不穩(wěn)定性比風(fēng)場更容易促進液化天然氣的擴散。ALOHA模型包含兩種不同的擴散模型：高斯模型和重氣云模型，能根據(jù)?；防砘再|(zhì)和泄漏量等參數(shù)自動選擇擴散模型。ALOHA不僅能夠模擬揮發(fā)型?；返臄U散過程，還能模擬?；吩磸姷淖兓?guī)律及擴散區(qū)域的濃度變化規(guī)律[5]。王志憲[5]在應(yīng)用ALOHA模型分析液氯泄漏災(zāi)害的影響因素時指出，ALOHA對氣體危化品泄漏后的擴散范圍和危險濃度閾值的估算預(yù)測具有較高的準(zhǔn)確度。

4 總結(jié)與展望

本文簡要概述了前人關(guān)于海上危險化學(xué)品泄漏擴散數(shù)值模擬方面的研究成果，并將危化品分為漂浮型、溶解型、沉降型和揮發(fā)型4大類，分別介紹其泄漏后在海洋環(huán)境中的擴散形式及歸宿特點；通過最廣泛的數(shù)值模型闡述其理論方法，并對其優(yōu)缺點和特性進行討論。

將?；穭澐譃?種類型是一種十分理想的狀態(tài)。數(shù)值模型中只考慮了最主要的擴散形式，忽略了次要擴散過程，因此可以將復(fù)雜的運動過程簡化為單一的運動形式，雖然大大減小了計算工作量，但也會嚴(yán)重影響模擬結(jié)果的準(zhǔn)確性?，F(xiàn)實生活中?；贩N類繁多，泄漏到海洋環(huán)境中的擴散過程可能會同時發(fā)生揮發(fā)、溶解、沉降和漂浮等復(fù)雜的運動過程。為了更真實地模擬海上?；沸孤U散的實際情況，未來應(yīng)將危化品的泄漏擴散過程作為一個整體，同時考慮沉降、溶解和揮發(fā)等運動過程，構(gòu)建相互耦合的數(shù)值模型。