復雜河網水系油粒子模型開發及溢油污染模擬

王 鵬，申 霞，王船海，華祖林，馬乙心，劉曉東

（1.河海大學淺水湖泊綜合治理與資源開發教育部重點實驗室，江蘇南京 210098；2.河海大學環境學院，江蘇南京 210098；3.南京水利科學研究院，江蘇南京 210029；4.河海大學水文水資源學院，江蘇南京 210098）

1 研究背景

隨著全球經濟的發展，世界各國對能源的需求與日俱增，海上石油運輸量迅猛增長，水體溢油事故風險大為提升。為了預測溢油事故發生后的運動軌跡和擴散范圍，科學評估其對生態環境的影響程度，提高溢油事件處置的應急響應能力，通常采用數值模擬方法對溢油自身擴展運動、在水流和風場作用下的漂移和擴散以及各種風化過程進行預測。溢油模型的研究始于1960年代，主要經歷三個階段：首先是油膜擴展模型研究階段，以Fay理論為代表，考慮重力、慣性力、表面張力和黏性力作用，通過不同階段油膜擴展的動力機制建立油膜直徑經驗表達式，取得了開創性成果［1］；第二階段為基于對流擴散方程的數值離散模擬溢油運動，其可能會引入與物理擴散無關的數值擴散，甚至完全掩蓋溢油的實際物理擴散過程，計算結果欠佳；第三階段為Johansen［2］和Elliott［3］等人提出的油粒子模型，該方法將水體中的溢油視為眾多離散油粒子的組合，粒子的平流過程具有拉格朗日性質，而剪切流和紊流引起的擴散過程屬于隨機運動，可用隨機走動法來模擬，即將紊流視為一種隨機流場，每個油粒子在紊流場中的運動類似于流體分子的布朗運動，運動的隨機性使得整個粒子“云團”在水體中擴散，由于該方法對溢油擴散過程的模擬更加接近真實情景，目前應用較為廣泛［4-6］。

在國外，油粒子模型最先開發并應用于海洋、海灣等水體的二維、三維溢油模擬，而在溢油概率較低的河網水系未得到應用。然而，我國河網地區的航運業普遍非常發達，隨著船舶流量和貨運量的快速上升，船舶溢油事故風險驟增［7］。長江、淮河、珠江等流域下游廣泛分布著大量的復雜河網水系，這些地區河道溝渠縱橫交錯，山丘區樹狀河網和平原區網狀河網并存，水系結構高度分異。此外，流域內水閘、泵站等水利工程眾多，水系連通性受人工高度控制，加上海洋潮汐的驅動作用，導致河道水流流向往復不定，水動力特征十分復雜。由于河網水流運動過程通常采用一維水動力模型模擬，僅能輸出沿河道方向的斷面平均流速，無法據此判斷油粒子運動至河網節點后流動方向，導致現有油粒子模型不能直接在河網地區應用。因此，在復雜河網水系采用油粒子模型開展溢油污染模擬，關鍵在于解決油粒子漂移至河網節點后的遷移和分配問題。

本文以一維河網水動力模型為基礎，針對油品在該類水體的擴展、輸移和風化特征，開發適用于河網水系溢油污染模擬的油粒子模型。一方面，使用隨機走動法模擬溢油的慣性擴展、黏性擴展和表面張力擴展等機械擴展過程，解決了油粒子模型不能模擬溢油自身擴展階段的難題［8］。另一方面，通過定義河網節點出流河道流向因子，判斷油粒子運動至節點后的流動去向，實現了油粒子模型與一維河網水動力模型的結合。

2 材料與方法

2.1 模型構建本研究在“數字流域系統”［9-10］基礎上進行油粒子模型開發工作，該系統的模型庫包括流域水文模型、污染負荷模型、水動力模型和水質模型。由于河網地區各類水體的水深均較淺，系統采用一、二維耦合的水動力模型對水流運動過程進行模擬［11］，為油粒子模型提供流速信息。

溢油進入水體后，其行為和歸宿包括兩個過程：空間上的傳輸以及同時發生的油品風化。空間傳輸由自身擴展和輸移兩個階段構成：擴展指油膜由于內部張力而導致的面積增大；輸移過程指在環境動力要素的作用下溢油的遷移運動，包括水平方向的漂移、擴散過程。油品的風化引起溢油組成性質的改變，包含蒸發、乳化、分散、溶解、光氧化、生物降解、吸附沉降、水體混合擴散以及生物體內的代謝作用等。

2.1.1 對流過程 對流過程考慮水流流速和風漂流流速對油粒子漂移的綜合影響。其對流過程的基本方程如下［12］：

式中：Si為油粒子位置；Uc為水流表面流速，m/s；Uw為風漂流流速，m/s；α為水流運動對油膜漂移的影響因子，通常取1.1～1.2。

（1）一維河網模型。一維河網模型水流流速采用一維河網水動力模型計算得到，風引起的油膜漂移速度大小表示為：

式中：β為風漂流系數，通常取3%～ 4%；w10為水面上10 m處風速，m/s；φ為風向與水流流向的夾角。

河網是由縱橫交錯的河道構成的水系，從物理結構上看，河道交匯處形成了若干節點。根據河道水流流向與河網節點的關系，可以將與節點相連的河道分為入流河道和出流河道。其中入流河道指的是水流流入節點的河道，水流流出節點的河道稱為出流河道。由于河網地區地勢低平，河道比降小，加上受潮汐和水利工程調控的影響，造成河道水流流向不定，因此與節點相連的河道是入流河道，還是出流河道，需要根據河網水動力模型的計算結果進行動態判斷。

對于出、入流河道交匯形成的節點，在河網數值模擬中通常將其作為零維模型對象處理，由于零維模型無法模擬水動力過程，導致油粒子隨水流運動至河網節點后，立即隨出流河道的水流繼續運動，因此節點大小和屬性不會對油粒子漂移和歸宿產生影響。同時，由于出流河道數量可能不止一條，為了預測油粒子通過節點后的運動軌跡，除了利用一維河網水動力模型預測流速之外，還需要判斷油粒子運動至河網節點后的去向。

由于河網水系的河流寬度普遍較窄，溢油可以在較短的時間內在寬度方向上擴散均勻，加之絕大部分溢油漂浮于水面，則油粒子在節點處的流向概率與出流河道的寬度及其流速成正比。此外，考慮到這類地區的河道水深通常相差不大，因此，模型假設油粒子流入某條河道的概率與該條河道的出流流量大小成正比，即以與節點相連的各出流河道流量為權重判斷油粒子流向。具體計算方法如下：

首先，將與某一節點相連的所有出流河道的流量從小到大排列，定義各條出流河道的流向因子為dfi，并按下式計算：

式中：dfi為第i條出流河道的流向因子；qi為第i條出流河道的流量；n為出流河道數量。

其次，定義各條出流河道的流向因子判斷區間。其中，第1 條出流河道的流向因子判斷區間為［0，df1］，第i條出流河道的流向因子判斷區間為

最后，對于運動至河網節點處的每一個油粒子，生成［0，1］的均勻隨機數Rd，如果Rd位于某條出流河道的流向因子判斷區間，則油粒子將隨水流漂移至該條出流河道。

（2）平面二維模型。對于平面二維模型，在風力作用下，油粒子漂移方向與風向成0°～40°夾角，此時風引起的油膜漂移速度可表示為［12］：

式中：D為考慮風向偏角的轉換矩陣；θ為地轉科氏力引起的風向偏角。

2.1.2 機械擴展和紊動擴散 溢油機械擴展及剪切流和紊流引起的粒子紊動擴散過程，均采用隨機走動方法模擬。通過在流場中追蹤各質點的運動軌跡，得到每一時刻各個油粒子所處的空間位置，統計各時刻油粒子的位置可得到各時刻溢油的空間分布。

根據分子紊動擴散理論，在一維空間情況下，隨機走動方差與擴散系數之間的關系可以表示為［12］：

式中：x為隨機走動距離，m；K為擴散系數，m2/s；Rn為均值為0，標準差為1.0的正態分布隨機數；Δt為時間步長，s。

因此，對于一維河網模型，油粒子每一個時間步長的隨機走動速度可采用下式計算：

式中：Vt為紊流擴散的隨機走動速度，m/s；De為油膜機械擴展系數，m2/s；DT為紊動擴散系數，m2/s。

同理，對于平面二維模型，每一個時間步長的隨機走動速度可表示為［8］：

式中δ為［0，π］之間的均勻分布隨機角。

機械擴展系數和紊動擴散系數的計算采用Sayed［8］提出的方法。對于機械擴展系數，將溢油的機械擴展過程分為慣性、黏性、表面張力三個階段，首先計算油粒子實際厚度及油粒子厚度的臨界值，再根據兩者之間的大小關系確定擴展階段，最后計算不同擴展階段對應的機械擴展系數De。紊動擴散系數考慮摩阻速度和水深的影響。

2.1.3 風化過程 溢油事故發生后，除了伴隨著擴展、對流、擴散等動力過程外，油品還經歷如蒸發、乳化、分散、溶解、光氧化及生物降解等風化過程，使油膜質量、油膜物理化學性質等發生一系列變化，這些變化主要和油品自身性質以及海況條件如風、波浪、水流、氣溫以及生物活動等有關。本研究建立的溢油預測模型主要用于突發溢油污染應急預測，分析溢油發生后短期內的油膜遷移擴散規律，為短期模型，而光氧化、生物降解等過程相比其它風化過程是一個非常緩慢的過程。此外，分散會影響溢油的垂向分布特征，但是不會對溢油在水體中的整體質量變化產生影響。因此，本研究重點對溢油初期風化作用較突出的蒸發和乳化兩個過程進行模擬。

（1）蒸發過程。蒸發是油品中的石油烴輕組分從液態變為氣態的過程，是溢油初期階段油品與大氣物質交換的重要過程，是使油品殘留量大幅減少的重要途徑。蒸發過程與油品組分、油膜厚度以及環境狀況等因素相關。對于原油泄露，模型采用目前廣泛使用的Stiver等［13］提出的計算模式。對于柴油或汽油等燃料油泄露導致的蒸發損失，采用Fingas 模型［14］計算，該模型將石油及其產品的蒸發表示為時間、溫度和180℃時質量蒸餾比的函數，認為蒸發過程滿足平方根函數或對數函數形式。

（2）乳化過程。溢油乳化是另一項重要的風化過程，指油類吸收水而形成油包水乳化液的過程。乳化使油滴體積增加2～5倍，主要受風速、波浪、油膜厚度、環境溫度、油風化程度等因素的影響。乳化作用通常在油膜拓展較大、厚度較薄時發生。風浪能量打碎油膜，水滴分散到油中，形成油包水的乳化液，呈黑褐色黏性泡沫乳油狀漂浮于水面。乳化物含水率采用Mackay 等提出的乳化作用方程［15］，該方程也被其它大多數溢油模型所采用。

2.1.4 河岸吸附 河網地區的河道寬度普遍較小，因此溢油在自身擴展階段就可能接觸到水陸邊界。當溢油接觸到河岸時，有可能被河岸本身或灘涂、坑洼和植物吸附。研究表明［16］，河岸對溢油的吸附作用有三種情況：完全吸收、完全反射和部分吸收。然而，由于河岸在材質、形狀、水力條件、植被等方面存在較大差異，對其吸收能力只能進行粗略分析。為了簡化計算，本研究對于河岸吸附作用模擬只考慮完全吸收和完全反射兩種情況，通過假設油粒子吸附概率，判斷其吸附狀態［17］。

2.2 輸入與輸出根據以上對溢油模型對流、擴散和風化過程的分析，模型輸入數據包括以下幾個方面：（1）溢油事件時空特征：溢油發生時間、溢油點坐標等；（2）溢油質量及油粒子數量；（3）油品初始理化性質，具體包括油品的初始溫度、密度、蒸發參數；（4）氣象數據，包括風速、風向等。

為了反映溢油對水環境的影響范圍和程度，模型輸出部分包括各計算單元的油粒子數量、質量、體積。其中油粒子體積考慮蒸發和乳化過程的影響，其剩余體積可表示為：

式中：V0為油粒子的初始體積，m3；Vi為蒸發和乳化后油粒子的剩余體積，m3；Fvi為第i個油粒子的蒸發率；Ywi為第i個油粒子乳化物的含水率。

初始密度由用戶根據油品種類設定，考慮溫度、蒸發和乳化作用的影響。

2.3 研究區選擇太湖流域位于長江下游地區，流域面積36 900 km2，流域內河網密布，湖泊眾多，水面率達17%，平原地區河道密度達3.2 km/km2，屬于典型的平原河網水系。為了提升流域在水資源保護、水環境治理等方面的管理水平和決策能力，經過二十多年的研制和不斷完善，以基于雙對象共享結構的數字流域系統為支撐［9］，開發了“太湖流域水量水質決策支持系統”。系統模型庫由水文模型、污染負荷模型、水動力模型、水質模型組成，可對流域產匯流過程、污染物在陸域和水體的遷移和轉化過程進行模擬。系統分別采用零維、一維和二維模型對全流域的97 個小型湖泊、952 條河道及太湖湖體進行概化，并包含188 座閘泵工程。其中一維模型共概化河道斷面4306 個，河道總長7958 km，二維模型生成網格2339個。太湖流域概化河網及太湖計算網格如圖1所示。

圖1 太湖流域概化河網及太湖計算網格

2.4 模型率定和驗證以“太湖流域水資源綜合規劃數模研制項目”為依托［18］，采用沿長江13座水利工程、沿杭州灣3 座水利工程、望亭水利樞紐及太浦閘全年引排水量資料作為“太湖流域河網水動力模型”的邊界條件，對2000年太湖及12個水位代表站、4個流量代表站的全年水位和流量過程進行了率定，同時利用環太湖流量巡測資料和流域供、用、耗、排水量調查資料對出入太湖水量和流域水量進行了平衡分析；此外，還利用1998、1999年水量實測和調查資料，對模型計算結果進行了驗證。率定和驗證結果均表明，太湖及大部分代表站的特征水位、全年水位和流量過程線預測值與實測值相差較小，各控制線水量平衡分析成果能夠反映太湖流域的水流實際運動狀況［19］。

由于溢油事故發生具有突發性和不可預見性，尤其對于水流流向多變的河網水系，很難在這類地區獲得水位、流量、溢油濃度的高頻率同步觀測資料，因此，基于觀測數據對溢油模型預測結果進行率定驗證的難度較高。本文通過設定溢油情景的方法開展數值試驗，重點對油粒子在河網水系遷移和分配過程的合理性進行分析。

2.5 溢油情景設置將太湖流域河網水系作為溢油污染數值試驗場景，耦合油粒子模型與“太湖流域河網水動力模型”，通過設置假想溢油事件，開展全流域溢油漂移和歸宿模擬。采用流向因子判斷法模擬油粒子在河網水系的遷移和分配過程，分析不同典型時刻油粒子空間分布特征，探討溢油預測結果的合理性及預測方法的科學性，評估自然風化過程和人為打撈措施對油粒子漂移和擴散的影響。

太湖流域自2002年開始實施“引江濟太”工程，以提升流域水資源與水環境承載能力。太浦河作為“引江濟太”工程的主要輸水通道、上海和嘉興等地的重要水源地及蕪申線（三級航道）重要組成部分，同時具有泄洪、排澇、供水、航運等多種功能，一旦發生溢油事故，將對水源地水質及周邊地區供水安全產生不利影響。因此，假定在太浦河太浦閘下游約5 km 處發生貨船碰撞溢油事故，泄漏油品為柴油，泄漏量為2 t，泄漏持續時間10 min，油粒子個數5000 個。模型計算時間為2000年1月1日8∶00至2000年1月20日8∶00，假設溢油事件發生在1月1日10∶00，時間步長5 min。

3 結果與討論

按上述設定的溢油情景開展河網突發溢油污染模擬。首先，以某一節點為例，分析油粒子運動至河網節點后的流動去向，驗證出流河道油粒子數量與出流流量大小的關系。圖2為2000年1月5日油粒子漂移至某河網節點處，并隨水流繼續向下游漂移的過程。如圖2所示，油粒子在1月5日00∶00 運動至太浦河與4條支流交匯形成的節點，隨后一部分油粒子沿太浦河繼續向下游運動，另一部分油粒子沿支流2運動。此時，支流1、支流3和支流4均為入流河道，太浦河及支流2為出流河道，兩條河道的油粒子數量分別為1683 個和114 個，流量分別為297.3 m3/s 和20.3 m3/s。由此可見，出流河道油粒子數量與其流量大小成正比，說明流向因子判斷法能夠以出流河道流量為權重，正確模擬油粒子在河網水系的遷移和分配過程。

圖2 油粒子在河網節點的遷移和分配

對于復雜河網水系，由于其河流寬度普遍比較狹窄，加之油膜在溢油初始階段的機械擴展速度較快，可以在溢油后的較短時間內沿河寬方向擴散均勻，因此，油粒子沿出流河道的漂移數量與出流河道寬度成正比，即出流河道寬度越大，向該河道漂移的油粒子數量越多。此外，出流河道流速大小也是影響油粒子漂移方向的關鍵因素，即出流河道流速越大，油粒子向該河道漂移的概率越高。同時，河網地區各河道水深通常相差不大。因此，綜合考慮油粒子流向概率與河寬、流速及水深的關系，以出流河道流量為權重判斷油粒子流向的方法是合理的。

目前，油粒子模型在河流溢油污染模擬中開展了一些應用研究［20］，取得了較好的預測效果，但是在油粒子漂移方向判斷上，其與現有應用于海洋、海岸等開敞水域的油粒子模型并無本質區別，也是基于二維或三維水動力模型的流場預測結果，按照對流方程計算其漂移軌跡，仍然無法應用于眾多河流交匯形成的復雜河網水系，使油粒子模型的應用場景受到一定限制。本研究通過定義出流河道流向因子及其判斷區間，在流向概率與流量大小間建立聯系，從實際預測效果來看，較好地解決了一維河網水動力模型預測結果無法決定油粒子漂移方向的問題。

為了進一步分析自然風化過程對溢油漂移和歸宿預測的影響，將不同典型時刻油粒子的空間分布繪于圖3。如圖所示，在溢油發生的初始階段，由于各支流流向大部分為流入太浦河，因此，從1月1日10∶00 到1月5日00∶00，油粒子總體隨水流沿太浦河向下游移動，部分溢油短時間進入支流，最終又匯入太浦河，因此形成了兩條溢油帶，總長度約為900 m（圖2（a）），在蒸發和河岸吸附等風化因素的作用下，油粒子總質量減小為854 g，相對于初始質量減少了57%。隨后，油粒子開始向周邊河網水系擴散，進入太浦河支流，影響范圍逐漸擴大（圖2（b））。1月8日11∶00 溢油進入攔路港，此時距事故發生169 h，油粒子總質量為528 g。1月9日13∶00 油粒子進入黃浦江，總質量進一步減至398 g。最后，受上游徑流和下游潮汐的雙重影響，溢油沿黃浦江往復震蕩并逐漸向下游入江口漂移。

圖3 典型時刻油粒子空間分布

除了自然風化過程，人為打撈措施也是影響油粒子漂移和擴散的重要因素，模型通過設置打撈范圍和打撈率，模擬圍油欄、吸油氈、消油劑等溢油回收和消除措施。假設在溢油事故發生24 h 后采取打撈措施，打撈率設置為50%。圖4為溢油打撈情景下典型時刻油粒子空間分布圖，由圖可見，油粒子進入攔路港和黃浦江的時間與未采取打撈措施的方案幾乎沒有差別，進一步說明對流過程是控制溢油在河網水系輸移的主導因素。從油粒子空間分布來看，打撈措施對減少局部水體的油粒子數量有較大影響，并且該影響隨著溢油輸移范圍的擴大而增強。

圖4 溢油打撈情景下典型時刻油粒子空間分布

4 結論

近年來，油粒子模型在溢油輸運和歸宿模擬方面獲得了巨大成功，已被廣泛應用于海洋、海灣等開敞水域的二維、三維溢油數值模擬。然而，由于河網地區的水流運動模擬通常采用一維水動力模型，導致現有的油粒子模型不能與其耦合，無法在河網水系開展溢油污染模擬。本研究通過定義河道流向因子及其判斷區間，確定油粒子運動至河網節點處的流動去向，開發了一種能夠與一維河網水動力模型耦合的油粒子模型。通過在太湖流域開展溢油數值試驗，分析了自然風化過程和人為打撈措施對油粒子漂移和擴散的影響。結果表明，流向因子判斷法能夠正確模擬油粒子在河網水系的遷移和分配過程，基于該方法開發的油粒子模型適用于水流條件多變的復雜河網水系，可對溢油在該類地區的擴散范圍及其水質影響進行科學評估。目前，該研究成果已成功集成于《太湖流域水環境綜合治理總體方案》重點項目“太湖流域水資源監控與保護預警系統”，具備了在該流域開展溢油污染模擬的應用條件，為提升復雜河網水系突發水污染預警與決策的科學水平提供了有力的技術支撐。