長江宜昌段溢油事故影響數學模型分析研究

胡　旭，辛小康，尹　煒，陳和春，王繼保，李光浩

（1.三峽大學　水利與環(huán)境學院，湖北　宜昌　443002；2.長江水資源保護科學研究所，武漢　430051）

長江是我國主要的河流運輸通道，客貨運輸量大。交通運輸部發(fā)布《2016年交通運輸行業(yè)發(fā)展統(tǒng)計公報》顯示[1]，內河運輸完成貨物量達35.72億噸，占水路總運輸貨物量的55.97%；內河港口貨物吞吐量完成47.76億噸，增長3.1%。隨著國家依托黃金水道推動長江經濟帶發(fā)展戰(zhàn)略的實施，內陸自貿實驗區(qū)相繼開放，貨運量在今后仍將繼續(xù)增加[2]。在同條件及預警技術管理水平不變前提下，這無疑增大了溢油事故風險發(fā)生的概率。眾所周知，溢油事故產生的石油類污染物會給整個河流水體造成巨大水環(huán)境污染，嚴重破壞水環(huán)境平衡，危害水生生物。相關研究表明，中華鱘在含油濃度為10 mg/L時，死亡時間僅為2天[3]。

為了更好地應對內河溢油事故風險，有必要對內河溢油事故進行模擬研究。目前，工程界和學術界對溢油事故的模擬研究多采用數值模擬方法[4-5]。現有研究對象多數為入?？诩昂Ｑ蟮貐^(qū)。對于環(huán)境差異較大、水文條件特殊、河道情況復雜多變的內河，溢油事故研究較少[6-9]。本文采用MIKE21軟件，建立平面二維溢油模型，模擬葛洲壩至下游枝江江段豐水期和枯水期突發(fā)溢油事故，并分析不同溢油量條件下油膜的擴展、輸移范圍及其對下游造成的影響。

1　模型計算原理

1.1　水動力模塊

MIKE21是由丹麥水力研究院研究開發(fā)的應用于近海海域、海灣、港區(qū)及河道等水環(huán)境模擬的系統(tǒng)。水動力模塊（HD）是其核心，可以模擬由于各種力的作用而產生的水位及水流變化。其控制方程是基于不可壓縮流體和Reynolds值均布的Navier-Stokes方程，并遵從Boussinesq假設和靜水壓力假設。二維非恒定淺水方程組如下：

式中，t為時間，η為水面高程，h為總水深，g為重力加速度，ρ為水的密度，ρ0為（淡）水的參考密度，f=2Ωsin?為科氏力系數（Ω為旋轉角速率，?為地理緯度），Pa為大氣壓強，Sij為輻射應力張量，S和(us,vs)分別為點源的排放量和速度，uˉ和vˉ分別為沿x，y方向上的流速在深度上的平均值；(τsx,τsy)和(τbx,τby)為水面風應力張量和河床床面應力張量，Tij為側向應力。

該模塊利用隱式交替方向ADI（Alternating di?rection implicit）技術對上述連續(xù)方程和動量方程進行離散，所得矩陣用追趕法求解，各項微分和主要系數均采用中心差分格式，Taylor級數展開的截斷誤差可達二至三階精度。

1.2　溢油模塊

MIKE21的Oil Spill模塊（OS模塊）是在HD模塊基礎之上建立的溢油模型，用于模擬溢油擴散運動規(guī)律和歸宿。溢油模塊采用“油粒子”對溢油量進行概化，采用拉格朗日方法描述油粒子的輸移軌跡，包括擴展、漂移、擴散等過程[9]，其輸移方程見公式（4）-（7）。同時，在輸移過程中油粒子也會發(fā)生如蒸發(fā)、乳化、溶解等風化過程，該軟件通過計算油粒子質量損失來體現。

（1）擴展：油膜擴展運動采用修正的Fay重力-粘力公式計算。公式如下：

式中：Ao為油膜面積；t為時間；Kα為擴散系數；Vo為油膜體積；Ro為油膜半徑，h0為初始油膜厚度，取10 cm。

（2）漂移：影響油粒子漂移速度的主要因素是水流和風作用力。漂移速度為:

式中：UW為水面上10 m處的風速，Us為表面流速，cw為漂移系數，一般取值0.02～0.03。

（3）擴散：由于單個粒子不能被分成幾片，因此擴散的過程被解釋為在隨機方向上的運動。對于二維的情況，可以將隨機走動的距離形式表示為一個時間步長α方向上的擴散距離。計算公式如下：

式中：Sα為在α方向上的一個時間步長內可能擴散走動的距離，Dα為α方向上的擴散系數，R為-1到1的隨機數。Δtp為計算時間步長。

2　研究區(qū)域概況

本文模擬范圍為長江葛洲壩至下游枝江江段。上游起始斷面為葛洲壩大壩，下游終止斷面為枝江昌門溪，范圍全長80 km。根據《全國重要江河湖泊水功能區(qū)劃（2011-2030年）》，該江段包含長江宜昌飲用水源、工業(yè)用水區(qū)和長江宜昌中華鱘保護區(qū)兩個水功能區(qū)。根據水功能區(qū)管理要求，長江宜昌飲用水源、工業(yè)用水區(qū)水質目標為Ⅱ類，水質目標要求較高。模擬范圍如圖1。

圖1　模擬范圍內各水功能區(qū)、監(jiān)測斷面、概化溢油事故點相對位置示意圖Fig.1　Relative position of each water function area,monitoring section,and generalized point of oil spill accident

3　模型參數設定及計算方案

3.1　網格劃分

本文采用2015年實測水下地形數據。葛洲壩至枝江區(qū)間的長江干流，河道形狀相對規(guī)則，變化梯度不大。為了保證更高計算精度和收斂性，模型網格劃分時采用四邊形網格，網格尺寸為200 m×40 m；支流清江入匯口段，地形復雜，區(qū)域較小，采用三角形網格劃分，三角形網格控制面積小于1 000 m2。計算區(qū)域共剖分計算網格單元16 093個，并根據該實測水下地形進行插值，如圖2和圖3所示。

圖2　葛洲壩-枝江區(qū)間水質數學模型計算網格Fig.2　Computing grid for water quality mathematical model between Gezhouba and Zhijiang section

圖3　葛洲壩-枝江區(qū)間計算網格地形插值Fig.3　Terrain interpolation for computing grid between Gezhouba and Zhijiang section

3.2　水動力模塊邊界條件及參數設定

葛洲壩是三峽的反調節(jié)水庫，為日調節(jié)水庫，設計水文條件可參照三峽出庫流量。三峽水庫2003年才開始初期蓄水，2010年才開始175 m試驗性蓄水，水文情勢尚不穩(wěn)定。于是本文便采用三峽水庫近10年最枯、最豐月的月均流量作為設計水文條件。本模型枯水期上游邊界設計流量為三峽2月份平均出庫流量5 280 m3/s，參照枝江水位站水位—流量關系曲線，得出模型下游邊界設計水位為33.42 m。豐水期上游邊界設計流量為三峽8月份平均出庫流量21 357 m3/s，對應下游邊界設計水位為36.15 m。參照多年平均水溫，枯水期設定水溫10℃，汛期設定水溫20℃。

設定濕水深度為0.100 m，淹沒深度為0.050 m，干水深度為0.005 m，根據辛小康[10]等對長江宜昌江段水動力模型的率定結果，曼寧系數為0.031，渦粘性系數的Smagorinsky系數為0.28，擴散系數D設定為1 m2/s。

3.3　溢油模塊參數設定

本文首先計算各個油粒子的位置變化和組成變化，然后統(tǒng)計各個網格上的油粒子個數和組分含量，模擬出油膜的濃度時空分布和組分變化。MIKE21中溢油模型將單個油粒子質量定義為揮發(fā)性輕組分的質量、不揮發(fā)性重組分的質量、蠟狀物的質量以及瀝青的質量。本文只考慮輕組分和重組分，兩者比例為3∶2，即1個1 000 kg油粒子，包含輕組分600 kg，重組分400 kg。由于石油類污染物易發(fā)生蒸發(fā)、溶解、乳化、沉淀、生物降解、光氧化等特點，本文參考姜衛(wèi)星[11]2007年對狹長型感潮河流溢油模型的研究，采用縱向擴散系數為0.25 m2/s、橫向為0.10 m2/s。風速和風向，可根據氣象資料獲得。

枯水期葛洲壩下游江段石油類背景濃度采用宜昌黃陵廟斷面2010-2014年2月份實測平均值0.03 mg/L。豐水期采用該斷面2010-2014年8月份實測平均值0.01 mg/L。

3.4　計算方案

本次模擬設計溢油點位置為該區(qū)域航運交通要塞三江引航道和大江引航道的匯合點，坐標經緯度為（E111°16′4.8″，N30°42′10.8″），垂向水深為零。參考我國《船舶油污染事故等級》（JT/458-2001）[12]，考慮最不利情形，即油船造成的水域油污染事故。根據油船大事故和重大事故的劃分依據，擬取設計總溢油量分別為5 t和10 t，并分別模擬枯水期和豐水期兩種不同水文條件下的溢油事故情形，共4種模擬工況，提取每種工況1 h，3 h，10 h后的溢油模擬計算結果進行對比分析。

4　計算結果與分析

4.1　枯水期

枯水期瞬時溢油為10 t時，事故發(fā)生后1 h，3 h和10 h以后油膜位置和面積如圖4所示?？梢钥闯觯斠缬?0 t后1 h，油膜尚未完全展開，油膜厚度≥0.3μm的范圍為600 m×300 m；油膜厚度大部分為100μm～0.1 mm的范圍為400 m×200 m。溢油后3 h，油膜厚度≥0.3μm的范圍為3600 m×300 m；油膜厚度為100μm～0.1 mm的范圍為3 200 m×280 m。溢油10 h以后，油類物質經過擴散，污染帶范圍進一步增加，但由于乳化、蒸發(fā)、風化、光解等綜合作用，油膜厚度有所減小，油膜厚度≥0.3μm的范圍可達到14 000 m×400 m；油膜厚度為100μm～0.1 mm的范圍為600 m×200 m。

圖4　枯水期溢油10 t后1 h，3 h，10 h油膜厚度分布圖Fig.4　Oil film thickness distribution at 1 h,3 h and 10 h after 10 t spillage during low water period

枯水期瞬時溢油5 t后1 h，3 h和10 h油膜演進情況見圖5?？梢钥闯觯斠缬秃? h，油膜厚度≥0.3μm的范圍為400 m×240 m；油膜厚度大部分為100μm～0.1 mm，其范圍為300 m×200 m。溢油后3 h，油膜厚度≥0.3μm的范圍為2 200 m×300 m；油膜厚度為100μm～0.1 mm的范圍為2 000 m×260m。溢油10 h以后，油類物質經過擴散，并有少量附著在河岸，污染帶范圍進一步增加，油膜厚度有所減小，油膜厚度≥0.3μm的范圍可達到10 000 m×400 m，而油膜厚度在100μm～0.1 mm幾乎消失。

圖5　枯水期溢油5 t后1 h，3 h，10 h油膜厚度分布圖Fig.5　Oil film thickness distribution at 1 h,3 h and 10 h after 5 t spillage during low water period

4.2　豐水期

豐水期瞬時溢油10 t時，事故發(fā)生后1 h，3 h和10 h以后油膜位置和面積如圖6所示?？梢钥闯?，當溢油10 t后1 h，油膜厚度為100μm～1 mm，范圍為800 m×100 m。溢油后3 h，油膜面積達到4 400 m×200 m，油膜厚度超過100μm。溢油10 h以后，油類物質經過擴散，污染帶范圍進一步增加，油膜范圍達到10 000 m×400 m，同樣由于乳化、蒸發(fā)、風化、光解等綜合作用，油膜厚度減小，基本小于100μm。

圖6　豐水期溢油10 t后1 h，3 h，10 h油膜厚度分布圖Fig.6　Oil film thickness distribution at 1 h,3 h and 10 h after 10 t spillage during high water period

豐水期瞬時溢油5 t時，事故發(fā)生后1 h，3 h和10 h以后油膜位置和面積如圖7所示?？梢钥闯?，當溢油5 t后1 h，油膜尚未完全展開，油膜厚度為100μm～1 mm，范圍800 m×100 m。溢油后3 h，油膜面積達2 600 m×200 m，油膜厚度超過100μm。溢油10 h以后，油類物質經過擴散，飄移，污染帶范圍進一步增加，油膜范圍達到8 000 m×300 m，而油膜厚度減小，基本小于100μm。

圖7　豐水期溢油5 t后1 h，3 h，10 h油膜厚度分布圖Fig.7　Oil film thickness distribution at1 h,3 h and 10 h after 5 t spillage during high water period

4.3　結果分析

根據模擬結果可知，對于所有工況，漏油5 t以上均會產生明顯的油膜厚度≥100μm的區(qū)域，且持續(xù)時間都超過3 h。溢油1 h后油膜處于或剛完成溢油擴散的初始階段，油膜展開范圍不大，厚度均大于100μm。溢油3 h后，油膜范圍繼續(xù)擴大，油膜影響范圍均大于2 000 m×260 m，最大可達3 600 m×300 m。溢油10 h以后，油類物質經過擴散，污染帶范圍進一步增加，無論是枯水期還是豐水期，兩種溢油量油膜影響范圍均大于8 000 m×300 m，最大可達10 000 m×400 m，但由于物理、化學等綜合作用，油膜厚度減小較明顯，普遍小于100μm。在同等水力條件下，溢油量越多油膜相對更厚，但油膜擴散范圍的差異較小，這可能與河道較狹窄有關。

對比枯水期和豐水期可知，枯水期由于水量偏小，流速小，油污范圍較小，持續(xù)的時間較長，而豐水期水量大，流速大，油膜掃過的范圍相對更大更遠，油膜厚度減小速率更大，持續(xù)的時間也相對較短?？菟谟湍哌^的范圍可達到猇亭一帶，豐水期可達到枝城一帶。無論是枯水期還是豐水期，當瞬時溢油量為5 t及以上時，均會對中華鱘自然保護區(qū)產生明顯影響。

5　結論

我國內河流域環(huán)境普遍復雜多變，與海洋有較大差別，對于內河溢油事故模擬，需要考慮的因素眾多。本文根據實測地形資料借用MIKE21軟件對宜昌江段進行溢油事故模擬預測，考慮了不同水文條件，不同溢油量的情形，而模擬結果也表明這對油膜厚度、飄移范圍、持續(xù)時間的差別有較為顯著的影響，這對內河溢油應急對策具有指導意義，也為內河溢油事故對水環(huán)境影響評估提供了很好的參考價值。

[1]中華人民共和國交通運輸部.交通運輸行業(yè)發(fā)展統(tǒng)計公報[R].北京：綜合規(guī)劃司，2016.

[2]袁群.淺談長江航運油污染現狀、原因及其對策[J].水運管理，2004（5）：22-23.

[3]倪朝輝，翟良安.石油對魚類等水生生物的毒性[J].淡水漁業(yè)，1997，27（6）：38-40.

[4]陳家興，杜娟，付金鋒.基于GIS的三峽庫區(qū)溢油模型研究[J].中國水運，2010，10（9）：96-98.

[5]高龍驤.內河船舶溢油擴散的建模與仿真[D].上海：復旦大學，2013.

[6]趙琰鑫，王永桂，張萬順，等.河道溢油污染事故二維數值模型研究[J].人民長江，2012，43（15）：81-84.

[7]GARCIA M R，TOVAR F H.Computer modeling of oil spill trajectories with a high accuracy method[J].Spill Science&Technology Bulletin，1999，5（5/6）：323-330.

[8]巫麗俊.感潮江段事故溢油二維數值模擬[D].南京：河海大學，2006.

[9]張帆，黃立文，鄧建，等.重慶主城區(qū)江段溢油模型及數值試驗研究[J].武漢理工大學學報（交通科學與工程版），2011，35（1）：87-90.

[10]辛小康，葉閩，尹煒.長江宜昌江段水污染事故的水庫調度措施研究[J].水電能源科學，2011，29（06）：46-48，95.

[11]姜衛(wèi)星.黃浦江溢油事故的數值模擬研究[D].上海：同濟大學，2007.

[12]中華人民共和國交通部.船舶油污染事故等級：JT/458-2001[S].北京：中國標準出版社，2001.