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

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

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

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

1　模型計算原理

1.1　水動力模塊

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

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

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

1.2　溢油模塊

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

（1）擴(kuò)展：油膜擴(kuò)展運(yùn)動采用修正的Fay重力-粘力公式計算。公式如下：

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

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

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

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

式中：Sα為在α方向上的一個時間步長內(nèi)可能擴(kuò)散走動的距離，Dα為α方向上的擴(kuò)散系數(shù)，R為-1到1的隨機(jī)數(shù)。Δtp為計算時間步長。

2　研究區(qū)域概況

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

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

3　模型參數(shù)設(shè)定及計算方案

3.1　網(wǎng)格劃分

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

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

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

3.2　水動力模塊邊界條件及參數(shù)設(shè)定

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

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

3.3　溢油模塊參數(shù)設(shè)定

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

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

3.4　計算方案

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

4　計算結(jié)果與分析

4.1　枯水期

枯水期瞬時溢油為10 t時，事故發(fā)生后1 h，3 h和10 h以后油膜位置和面積如圖4所示。可以看出，當(dāng)溢油10 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以后，油類物質(zhì)經(jīng)過擴(kuò)散，污染帶范圍進(jìn)一步增加，但由于乳化、蒸發(fā)、風(fēng)化、光解等綜合作用，油膜厚度有所減小，油膜厚度≥0.3μm的范圍可達(dá)到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油膜演進(jìn)情況見圖5。可以看出，當(dāng)溢油后1 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以后，油類物質(zhì)經(jīng)過擴(kuò)散，并有少量附著在河岸，污染帶范圍進(jìn)一步增加，油膜厚度有所減小，油膜厚度≥0.3μm的范圍可達(dá)到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所示?？梢钥闯?，當(dāng)溢油10 t后1 h，油膜厚度為100μm～1 mm，范圍為800 m×100 m。溢油后3 h，油膜面積達(dá)到4 400 m×200 m，油膜厚度超過100μm。溢油10 h以后，油類物質(zhì)經(jīng)過擴(kuò)散，污染帶范圍進(jìn)一步增加，油膜范圍達(dá)到10 000 m×400 m，同樣由于乳化、蒸發(fā)、風(fēng)化、光解等綜合作用，油膜厚度減小，基本小于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所示?？梢钥闯觯?dāng)溢油5 t后1 h，油膜尚未完全展開，油膜厚度為100μm～1 mm，范圍800 m×100 m。溢油后3 h，油膜面積達(dá)2 600 m×200 m，油膜厚度超過100μm。溢油10 h以后，油類物質(zhì)經(jīng)過擴(kuò)散，飄移，污染帶范圍進(jìn)一步增加，油膜范圍達(dá)到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　結(jié)果分析

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

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

5　結(jié)論

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

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

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

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

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

[5]高龍驤.內(nèi)河船舶溢油擴(kuò)散的建模與仿真[D].上海：復(fù)旦大學(xué)，2013.

[6]趙琰鑫，王永桂，張萬順，等.河道溢油污染事故二維數(shù)值模型研究[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]巫麗俊.感潮江段事故溢油二維數(shù)值模擬[D].南京：河海大學(xué)，2006.

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

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

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

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