危險化學品海上漂移擴散數值模擬研究綜述

郭凱旋

（國家海洋環境預報中心，北京100081）

1 引言

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

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

海上危化品泄漏事故與其他海洋環境污染事故相比具有突發性強、影響范圍廣、影響時間長和泄漏后變化復雜等特點[3]。突發性強表現為海上危化品泄漏大部分都是在沒有任何前兆的情況下突然發生的，瞬間排放大量的污染物，給海洋環境造成巨大影響。影響范圍廣表現為少量的危化品發生泄漏，其產生的影響范圍將非常巨大，例如1 t氯的泄漏能夠影響4.8 km2的范圍[3]。影響時間長表現為危化品泄漏對海洋生態環境的影響不僅僅發生在泄漏時，部分危化品對海洋生態環境的污染極難消除，例如福島核泄漏事故污染物CS-137的半衰期長達30.07 a[4]。泄漏后變化復雜表現為危化品種類繁多且多樣的行為方式造成危化品泄漏入海后物理和化學形態變化十分復雜。

海上危化品泄漏事故中溢油事故最為頻繁，因此英、美等發達國家很早就開展了海上溢油事故相關理論研究，并通過大量的試驗驗證，建立了眾多預測模型[5-8]。與溢油相比，其他類型的海上危化品泄漏漂移擴散數值模擬研究相對較少[9-12]。本文將對近年來國內外海上危化品泄漏在海洋水體環境中漂移擴散數值模擬方面的研究進行整理、分類和介紹。

2 危化品泄漏數值模擬水動力模式

水動力數值模擬是研究海上危化品泄漏漂移擴散的基礎。三維模型由于能夠全面而立體地反映實際河口海岸和近海環境中海水的運動狀況和基本規律，因此在海洋水動力數值模擬研究中得到了廣泛研究和應用。伴隨著計算機的飛速發展以及數值模擬技術的不斷完善，目前三維模型已經成為水動力數值模擬的主要發展方向[13-14]，其中研究應用較多的模型有POM（Princeton Ocean Model)、FVCOM（Finite-Volume Coastal Ocean Model）和ROMS（Regional Ocean Modeling System）。

2.1 POM

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

2.2 FVCOM

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

2.3 ROMS

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

3 海上危化品漂移擴散數值模型

在海洋水動力數值模式基礎上，通過引入危化品運動擴散方程，發展形成了適用于不同危化品類型的海上危化品漂移擴散預測數值模型。海上危險化學品一旦發生泄漏，危化品將在海流和風的作用下進行平流、擴散、揮發、溶解和沉降（或上升）等理化過程[31]。海上危險化學品泄漏后在海中的輸移和擴散模式主要受兩方面影響。一方面是外部環境，如地理環境、水文因素、氣象因素和其他因素[16]。危化品在海流和風的影響下將進行平流擴散和湍流擴散：平流擴散指海水的整體運動引起污染物有規律的平移；湍流擴散指因風和潮流引起的隨機無序的流體運動,運動狀態只能通過概率統計方法加以描述。擴散類型由海流類型決定,并受潮流、風生環流、密度流以及水深和岸線狀況等諸多因素制約[32-33]。在海洋這樣一個大環境中,污染物一旦泄漏，不僅是單一的平流擴散或湍流擴散，而是兩種擴散形式同時作用于污染擴散的整個過程，只不過作用的程度不同[33]。另一方面是化學品的主要物理化學性質，其中最主要的3個因素是溶解度、揮發性和密度。溶解度是決定危化品擴散運動形式的主要因素，是對其他性質進行分析的前提；揮發性決定了危化品入海后是否以蒸氣形式向空氣中擴散，危化品的揮發性用20℃飽和蒸氣壓衡量，通常以3 kPa作為揮發與不揮發的界限；密度（與水的相對密度）則決定危化品入海后是在水面運動還是進入水體甚至沉降至海底[34]。根據危化品泄漏入海后的運動形式，我們將海上危化品分為溶解擴散型、海面漂移型、沉降型和揮發型4種類型。

3.1 溶解擴散型

溶解擴散型危化品溶于海水，并在海水中以三維形式進行擴散。對于溶解型危化品在海水中的遷移擴散研究，目前主要從數值優化求解污染物對流擴散方程這一方向來解決。國內外研究求解對流擴散方程的數值方法大致可分兩類：一類是歐拉法，將研究海域的空間點作為研究對象，利用解析解的方法或者數值求解的方法，求解危化品在該海域的連續性方程、運動方程和輸運方程，從而得出危化品濃度隨時間的變化特征；另一類是拉格朗日法，該方法將研究海域中的危化品概化為具有一定數量和質量的質點顆粒，追蹤每個質點粒子在研究海域流場中的運動軌跡，得到每個時刻和每個質點顆粒所處的空間位置[35]。

對流擴散方程為：

式中：C為危化品濃度；ux、uy和uz分別為X、Y和Z方向的流速；Dx、Dy和Dz分別為X、Y和Z方向的擴散系數；S0為危化品的源和匯。

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

3.2 海面漂移型（類油型）

海面漂移型危化品不溶于海水，在海面主要以二維形式漂移擴散。人們對于這類危化品的研究主要集中在海上溢油，因此海面漂移型又可簡稱“類油型”。海上溢油的漂移擴展過程主要分慣性力擴展、粘性力擴展和表面張力擴展3個階段進行研究，早期具有代表性的是Fay公式[15-16，36-37]。

慣性力擴展階段：

粘性力擴展階段：

表面張力擴展階段：

式中：D為油膜擴散直徑（單位：m）；g為重力加速度（單位：m/s2）；V為溢油總體積（單位：m3）；t為從泄漏開始計算的時間（單位：s）；β=1-ρ0/ρw，ρ0和ρw分別為危化品密度和海水密度（單位：103kg/m3）；Vw為海水運動粘性系數（單位：m2/s）；δ為擴展系數；K1、K2和K3分別為3個階段的比例系數。

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

Fay理論及其修正模型對于在開闊海域瞬時泄漏情況下的類油型污染物模擬取得了很好的結果，適合類油型危化品在泄漏初期階段的模擬預測[38]。連續泄漏油膜在距離泄漏源一定距離后，側向擴散起主導作用，而對于水下泄漏，浮力的作用至關重要，由于Fay模式忽略這些作用，所以它并不適用于連續泄漏和水下泄漏的模擬[39]。

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

3.3 沉降型

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

海洋水動力要素是影響危化品在海水中擴散和漂移過程的主要因素。擴散是指危化品泄漏入海后由于濃度差和湍流等物理因素而形成污染物向整個水體遷移的混合過程；漂移指泄漏入海的危化品在風、表層和次表層流以及波浪作用下的拉格朗日漂移過程，漂移運動主要由平流條件決定[16，48]。陳協明[48]在研究化學品的溢漏擴散模式時指出，沉降型化學品在海面上的漂移主要受風的切應力、表層/次表層流和余流（波生余流和潮余流）控制。

沉降過程一般可以分為3個階段：沉降射底、觸底和底浪的傳播[16]。危化品入水后發生破碎并迅速下降，形成了充分成長的高速射流。在這一過程中，周圍海水與危化品發生混合，沉降射流的體積不斷擴大；沉降射流觸底時將形成以碰撞點為中心向四周擴展的高密度底浪；開始時底浪的傳播速度很快，隨著傳播距離的增加，底浪的厚度和傳播速度迅速減少，底浪中攜帶的危化品迅速下降并沉積[16，48]。

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

3.4 易揮發型

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

計算機仿真模擬是研究易揮發型危險化學品泄漏擴散的又一熱點領域。計算機仿真模型考慮的因素較多，模擬時間較長，因此預測精確度很高。此類模型最具代表性的就是CFD（Computational Fluid Dynamics）和ALOHA（Areal Locations Of Hazardous Atmospheres）[55]。CFD具有可視化能強和計算方法成熟等特點[43],在實際中被廣泛應用。Ikealumba等[56]在應用CFD研究海上液化天然氣的泄漏擴散規律時發現，海洋波動帶來的不穩定性比風場更容易促進液化天然氣的擴散。ALOHA模型包含兩種不同的擴散模型：高斯模型和重氣云模型，能根據危化品理化性質和泄漏量等參數自動選擇擴散模型。ALOHA不僅能夠模擬揮發型危化品的擴散過程，還能模擬危化品源強的變化規律及擴散區域的濃度變化規律[5]。王志憲[5]在應用ALOHA模型分析液氯泄漏災害的影響因素時指出，ALOHA對氣體危化品泄漏后的擴散范圍和危險濃度閾值的估算預測具有較高的準確度。

4 總結與展望

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

將危化品劃分為4種類型是一種十分理想的狀態。數值模型中只考慮了最主要的擴散形式，忽略了次要擴散過程，因此可以將復雜的運動過程簡化為單一的運動形式，雖然大大減小了計算工作量，但也會嚴重影響模擬結果的準確性。現實生活中危化品種類繁多，泄漏到海洋環境中的擴散過程可能會同時發生揮發、溶解、沉降和漂浮等復雜的運動過程。為了更真實地模擬海上危化品泄漏擴散的實際情況，未來應將危化品的泄漏擴散過程作為一個整體，同時考慮沉降、溶解和揮發等運動過程，構建相互耦合的數值模型。