長江干流銅陵段突發(fā)水污染情景模擬分析

齊于順，劉仁志*，張啟月，崔保山，郭忠，王飛，方麗，方陵

1.北京師范大學(xué)環(huán)境學(xué)院

2.銅陵市生態(tài)環(huán)境局

突發(fā)環(huán)境事件是指由于污染物排放或自然災(zāi)害、生產(chǎn)安全事故等因素，導(dǎo)致有毒有害物質(zhì)大量進(jìn)入環(huán)境介質(zhì)，突然造成或可能造成環(huán)境質(zhì)量下降，危及公眾身體健康和財產(chǎn)安全，需要采取緊急措施予以應(yīng)對的事件[1]。與其他普通污染事件不同的是，突發(fā)環(huán)境事件具有隨機(jī)性、復(fù)雜性、高強度、高風(fēng)險的特點[2]，會在短時間內(nèi)造成周邊環(huán)境污染物超標(biāo)，甚至?xí)斐砷L期影響[3]。其中，突發(fā)水環(huán)境事件最易發(fā)生，李旭等[4]對2011—2017年突發(fā)環(huán)境污染事件的受體類型進(jìn)行統(tǒng)計發(fā)現(xiàn)，水是主要污染受體，占比高達(dá)85.59%。長江沿岸城市“化工圍江”、上下游排污口與取水口交錯分布等布局性風(fēng)險突出，導(dǎo)致長江水環(huán)境風(fēng)險形勢不容樂觀[5-6]。此外，由于長江水體流動性強的特點，突發(fā)水環(huán)境事件的影響范圍更大，響應(yīng)時間更短，如果不能及時妥善處置，短時間內(nèi)可能導(dǎo)致飲用水水源地污染、水域生態(tài)環(huán)境破壞，甚至造成跨界污染，產(chǎn)生極大的不良社會影響[7-8]。目前突發(fā)水環(huán)境事件仍是影響供水和人民安全的重要威脅因素[9-10]。因此，有必要通過多種方式提前做好突發(fā)水環(huán)境事件的預(yù)測預(yù)警與應(yīng)急準(zhǔn)備，有效防范突發(fā)水污染事故，降低突發(fā)水環(huán)境風(fēng)險[11]。

水環(huán)境模型可對水環(huán)境系統(tǒng)及其內(nèi)部發(fā)生的復(fù)雜過程進(jìn)行定量化描述[12]，模擬污染物在一定時間和空間范圍內(nèi)的遷移轉(zhuǎn)化過程，計算污染物到達(dá)下游環(huán)境敏感目標(biāo)的時間、污染團(tuán)的歷經(jīng)時間和污染物的濃度分布等[13]，有助于幫助應(yīng)急管理人員了解事故可能的發(fā)展趨勢并做出合理的處置。常用的河流水質(zhì)模型主要有EFDC、WASP、QUAL2E、QUAL2K、MIKE、BASINS、SWAT 等[13-17]。其中，丹麥水力研究所開發(fā)的平面二維模型MIKE21，因用戶界面友好，具有多種計算網(wǎng)格、模塊供用戶選擇，且有強大的前、后處理功能，被廣泛用于河流[18-19]、湖泊[20-21]、水庫[22]、海灣[23]等水質(zhì)模擬。如袁玥等[18]應(yīng)用MIKE21模型對長江蘄春段非正常排污情況進(jìn)行了模擬，舒長莉等[24]應(yīng)用MIKE21模型對贛江南昌段突發(fā)污染事故進(jìn)行了模擬，均實現(xiàn)了對污染影響后果的量化表達(dá)。為貫徹落實習(xí)近平總書記“不搞大開發(fā)，共抓大保護(hù)”的重要指示精神，著力解決長江水環(huán)境、水生態(tài)、水風(fēng)險突出問題[25]，為長江經(jīng)濟(jì)帶高質(zhì)量發(fā)展提供保證[26]，筆者以長江生態(tài)環(huán)境保護(hù)修復(fù)駐點跟蹤研究工作為基礎(chǔ)，利用MIKE21構(gòu)建長江干流銅陵段突發(fā)水污染事故水動力-水質(zhì)耦合模型，考慮了多種可能的泄漏情景，并對不同水文期進(jìn)行模擬分析，實現(xiàn)對突發(fā)水污染事故的動態(tài)模擬，以期為長江干流銅陵段水域突發(fā)水污染事故風(fēng)險防范、預(yù)測預(yù)警和應(yīng)急處置提供參考。

1 研究區(qū)概況

銅陵市位于安徽省中南部、長江下游（117°04’E～118°09’E，30°38’N～31°09’N），是長江經(jīng)濟(jì)帶重要節(jié)點城市。長江干流銅陵段全長142.6 km，過境的長江水資源豐富，是銅陵市主要供水來源，最小日平均流量為10 700 m3/s，最大日平均流量為68 400 m3/s，年平均流量為29 250 m3/s（2019年），水量年內(nèi)變化顯著。研究河段內(nèi)環(huán)境敏感目標(biāo)包括銅陵市淡水豚國家級自然保護(hù)區(qū)和3處集中式生活飲用水水源地（銅陵市水廠、新三水廠和第五水廠）。淡水豚國家級自然保護(hù)區(qū)范圍上始銅陵市樅陽縣老洲鎮(zhèn)，下至銅陵市義安區(qū)金牛渡，全長58 km，其中核心區(qū)、緩沖區(qū)和實驗區(qū)面積分別為9 534、6 360、15 624 hm2[27]。橫港掃把溝地區(qū)緊臨長江，是銅陵市最早的工業(yè)聚集區(qū)之一，由于歷史原因，該區(qū)域功能布局不合理，工業(yè)用地與居住用地混雜，多家企業(yè)距長江直線距離不足1 km。目前該區(qū)域仍存留多家風(fēng)險企業(yè)，也沒有成立工業(yè)園區(qū)管委會，這給下游集中式飲用水水源地、淡水豚國家級自然保護(hù)區(qū)以及長江水體水質(zhì)安全帶來重大威脅。研究區(qū)概況如圖1所示。

圖1 研究區(qū)概況Fig.1 Schematic map of the study area

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

研究區(qū)水文數(shù)據(jù)采用2019年大通水文站實測逐日流量數(shù)據(jù)和橫港、荻港2個水位站實測逐日水位數(shù)據(jù)；地形文件為1∶25 000實測水下地形數(shù)據(jù)；風(fēng)險企業(yè)、飲用水水源保護(hù)區(qū)和淡水豚保護(hù)區(qū)等相關(guān)資料由銅陵市生態(tài)環(huán)境局提供。

2.2 MIKE21模型基本原理

長江銅陵段江面寬度遠(yuǎn)大于深度，且不存在明顯分層現(xiàn)象，因此采用二維水動力模型可以滿足研究需要[28]。選取MIKE21 Flow Model FM模型進(jìn)行模擬，二維水動力模型基于三向不可壓縮和Reynolds值均布的Navier-Stokes方程，并服從于Boussinesq假定和靜水壓力的假定，水動力控制方程組（淺水方程）為[24,29-31]：

式中：t為時間；x、y為笛卡爾坐標(biāo)； η為水位，m；d為靜止水深，m；h為 總水深，h=d+η，m；u、v分別為x、y方向上的速度分量；f為哥氏力系數(shù)，f=2ωsinφ，其中 ω 為地球自轉(zhuǎn)角速度，φ為當(dāng)?shù)鼐暥?；g為重力加速度，取9.81 m /s2；ρ和 ρ0分別為水體和參考水的密度， k g/m3；Pa為大氣壓力，Pa；sxx、sxy、syx、syy分 別為x軸的法線方向、y軸的法線方向、xy表面上的切線方向和yx表面上的切線方向的輻射應(yīng)力分量；Ax、Ay為應(yīng)力分項；Txx、Txy、Tyx、Tyy為黏滯切應(yīng)力分量；S為源匯項；us、vs為源項水流流速，m/s；uˉ、vˉ為沿水深平均的流速；τsx、τsy為水面風(fēng)應(yīng)力在x、y方向上的分量;τbx、τby為河床 ( 海底) 底部應(yīng)力的分量。

uˉ、vˉ計算公式如下：

水質(zhì)模塊建立于二維空間下的輸移擴(kuò)散方程為：

式中：C為各典型污染物濃度；Dx、Dy為各典型污染物在x、y上的擴(kuò)散系數(shù)；P為各典型污染物降解項。

MIKE21模型采用有限體積法對計算區(qū)域進(jìn)行空間離散，將該連續(xù)統(tǒng)一體細(xì)分為若干個不重疊的三角形或四邊形單元。模型計算的時間和精度取決于計算數(shù)值方法所使用的求解格式精度，淺水方程的時間積分和輸移擴(kuò)散方程基于半隱格式求解，相應(yīng)平流項采用顯式格式求解，垂直對流項采用全隱格式求解。受顯式格式穩(wěn)定性的限制，為保持模型計算的穩(wěn)定性，模型中時間步長的設(shè)定必須保證CFL（Courant-Friedrich Levy)數(shù)小于1。淺水方程和輸移擴(kuò)散方程在笛卡爾坐標(biāo)上的CFL分別定義為：

式中： Δx和 Δy為x、y方向的特征長度； Δt為時間間距。 Δx和 Δy近似于三角形網(wǎng)格的最小邊長，水深和流速值為三角形網(wǎng)格中心的取值。

2.3 模型建立與驗證

利用MIKE21建立水動力-水質(zhì)耦合模型，模型構(gòu)建流程如圖2所示?；阢~陵市水利工程圖和實測水下地形數(shù)據(jù)，利用ArcMap與AutoCAD軟件提取生成長江銅陵段水陸邊界線（Land.xyz文件），并讀取相應(yīng)水深散點數(shù)據(jù)（Water.xyz文件）。在Mesh Generator 中導(dǎo)入水陸邊界線（Land.xyz文件）并進(jìn)行光滑處理，采用非結(jié)構(gòu)三角網(wǎng)格對河道地形進(jìn)行處理，然后生成網(wǎng)格；插入地形水深散點數(shù)據(jù)（Water.xyz文件）并進(jìn)行地形插值后導(dǎo)出mesh文件（圖3），計算區(qū)域共有10 896個三角網(wǎng)格，6 886個網(wǎng)格節(jié)點。模型上邊界采用2019年大通水文站逐日流量時間序列，下邊界采用2019年荻港水位站逐日水位時間序列。選擇Flow Model FM建立水動力-水質(zhì)耦合模型，經(jīng)過多次調(diào)試后選擇合適參數(shù)；將模擬時長設(shè)為24～30 h，時間步長選擇60 s；為保證模型穩(wěn)定運行，設(shè)置CFL為0.8，最大時間步長為60 s；干水深(hdry)為0.005m，淹沒水深（hflood）為0.05m，濕水深(hwet)為0.1m；底床摩擦力采用曼寧系數(shù)，取值32 m1/3/s；初始水位根據(jù)不同水文期選擇模擬對應(yīng)開始日期水位值；其他參數(shù)采用推薦或默認(rèn)值。在對流擴(kuò)散模塊，添加需要模擬的風(fēng)險物質(zhì)硫酸，從更安全的角度考慮預(yù)測結(jié)果，將硫酸降解系數(shù)設(shè)為0[32]；并在源匯項中添加風(fēng)險物質(zhì)初始濃度和泄漏流量等相關(guān)數(shù)據(jù)。

圖2 模型構(gòu)建流程Fig.2 Model construction flowchart

圖3 水動力模型網(wǎng)格地形高程Fig.3 Grid terrain elevation of hydrodynamic model

因長江銅陵段水文情勢變化顯著，有明顯的豐、平、枯水文情勢節(jié)律性變化，且枯水期部分時間段水位站處于不工作狀態(tài)。故對豐、平、枯水文時期部分時間段進(jìn)行模擬，并選取2019年橫港水位站逐日水位數(shù)據(jù)進(jìn)行驗證，將實測水位數(shù)據(jù)與模擬水位數(shù)據(jù)進(jìn)行對比。經(jīng)對照發(fā)現(xiàn)（圖4），模擬水位與實際水位數(shù)值相差不大，平均相對誤差為2.4%，曲線趨勢基本相似，波峰與波峰相對，波谷與波谷相對。模擬值與實測值吻合性良好，可作為水質(zhì)模塊的水動力基礎(chǔ)。

圖4 橫港水位站2019年模擬水位與實測水位Fig.4 Simulated and measured water level of Henggang water-level station in 2019

2.4 泄漏情景設(shè)定與水質(zhì)模擬

橫港掃把溝地區(qū)多家風(fēng)險企業(yè)的儲罐內(nèi)存在大量風(fēng)險物質(zhì)，給下游集中式飲用水水源地、淡水豚國家級自然保護(hù)區(qū)以及長江水體水質(zhì)帶來巨大安全隱患。針對某化工企業(yè)硫酸儲罐，設(shè)定3種泄漏情景，分別為罐體20%管徑破裂、罐體100%管徑破裂和儲罐完全泄漏，泄漏時間均設(shè)為10 min。在豐水期、平水期和枯水期對3種泄漏情景分別進(jìn)行模擬，具體泄漏時間分別設(shè)定為7月10日、4月16日和10月12日的06:00。所設(shè)定的泄漏情景如表1所示。

表1 3種設(shè)定泄漏情景Table 1 Three supposed leakage situations

硫酸罐體20%管徑和100%管徑破裂時泄漏速率用伯努利方程[33]計算（液體在泄漏口沒有急驟蒸發(fā)），公式如下：

式中：QL為液體泄漏速率， k g/s ；P1為容器內(nèi)介質(zhì)壓力，Pa；P0為環(huán)境壓力，Pa；ρ為泄漏液體密度， k g/m3；h為 裂口之上液位高度， m ；Cd為液體泄漏系數(shù)，取0.62；A為裂口面積， m2。

根據(jù)實地調(diào)研，判斷突發(fā)泄漏情況下的可能入江位置如圖1所示，將不同泄漏情景下計算得到的泄漏流量-時間序列等相關(guān)數(shù)據(jù)輸入構(gòu)建的水動力-水質(zhì)耦合模型，可得到不同水文期的風(fēng)險物質(zhì)遷移擴(kuò)散分布結(jié)果。

3 結(jié)果與討論

采用MIKE21水動力-水質(zhì)耦合模型對3種泄漏情景在豐水期、平水期和枯水期的風(fēng)險物質(zhì)遷移擴(kuò)散情況分別進(jìn)行了模擬，定量化表達(dá)出風(fēng)險物質(zhì)到達(dá)下游環(huán)境敏感目標(biāo)的時間、最大污染峰團(tuán)濃度和污染持續(xù)時間等。泄漏硫酸對長江水質(zhì)的影響可通過對水體pH的影響來反映[32]，GB 3838—2002《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》規(guī)定地表水Ⅰ類～Ⅴ類水域水體pH為6～9，當(dāng)pH為6時，通過計算得到對應(yīng)硫酸的濃度為0.049 mg/L，即硫酸濃度大于0.049 mg/L時會對下游各風(fēng)險受體和長江水質(zhì)造成影響。通過水質(zhì)模型模擬結(jié)果可得到儲罐泄漏后各時間節(jié)點的濃度，進(jìn)而預(yù)測和判斷影響程度。下面針對同一水文期不同泄漏情景和同一泄漏情景在不同水文期的水質(zhì)模擬結(jié)果進(jìn)行說明。

3.1 不同泄漏情景下水質(zhì)模擬結(jié)果

以3種不同泄漏情景在豐水期模擬結(jié)果為例進(jìn)行說明，模擬結(jié)果如圖5～圖7所示。3種泄漏情景下最大污染峰團(tuán)到達(dá)下游敏感受體所需時間相同，到達(dá)三水廠取水口、市水廠取水口、淡水豚核心區(qū)、五水廠取水口和市出境斷面的時間分別為89、93、332、677、751 min，這是由于同一時間段的水動力結(jié)果相同。但在儲罐完全泄漏情景下，風(fēng)險物質(zhì)到達(dá)下游敏感受體時間稍早于其他2種泄漏情景，風(fēng)險物質(zhì)完全流過敏感受體時間稍晚于其他2種泄漏情景，即影響時間更長。此外，該情景下最大污染峰團(tuán)濃度也明顯高于其他2種泄漏情景下2～3個數(shù)量級。在三水廠取水口處，20%管徑破裂時最大污染峰團(tuán)濃度為1.891 mg/L，100%管徑破裂時最大污染峰團(tuán)濃度為47.273 mg/L，儲罐完全泄漏時最大污染峰團(tuán)濃度為1 381.39 mg/L，這是由于不同泄漏情景下風(fēng)險物質(zhì)泄漏總量不同。隨著時間的延長，風(fēng)險物質(zhì)被稀釋，到達(dá)下游各敏感受體時的最大污染峰團(tuán)濃度會逐漸降低。在3種不同泄漏情景下，風(fēng)險物質(zhì)到達(dá)市出境斷面時的最大污染峰團(tuán)濃度分別為0.090 5、2.262 4和66.123 4 mg/L，均超過0.049 mg/L的標(biāo)準(zhǔn)要求，會造成不同程度的跨界污染。淡水豚核心保護(hù)區(qū)范圍最大，其受影響持續(xù)時間也最長，分別為 251、313、343 min。

圖5 豐水期罐體20%管徑破裂時水質(zhì)模擬結(jié)果Fig.5 Simulation results of water quality when 20% pipe diameter of tank breaks in wet period

圖6 豐水期罐體100%管徑破裂時水質(zhì)模擬結(jié)果Fig.6 Simulation results of water quality when 100% pipe diameter of tank breaks in wet period

3.2 不同水文期泄漏水質(zhì)模擬結(jié)果

以不同水文期儲罐完全泄漏情景的模擬結(jié)果為例進(jìn)行說明，模擬結(jié)果如圖7～圖9所示。根據(jù)模擬結(jié)果可知，豐水期風(fēng)險物質(zhì)到達(dá)下游敏感受體時間最短，其中風(fēng)險物質(zhì)到達(dá)三水廠取水口僅用75 min，影響持續(xù)時間為43 min。相較于豐水期，平水期和枯水期的流速更小，風(fēng)險物質(zhì)需更長時間到達(dá)下游敏感受體，如平水期、枯水期風(fēng)險物質(zhì)到達(dá)三水廠取水口的時間分別為103和111 min。但在平水期和枯水期，風(fēng)險物質(zhì)對敏感受體的影響持續(xù)時間會更長，其對三水廠取水口的影響持續(xù)時間分別為65和72 min。同時因流量更小，風(fēng)險物質(zhì)在市水廠取水口、三水廠取水口和淡水豚核心區(qū)的最大污染峰團(tuán)濃度會更高，枯水期的濃度分別為1 883、1 736和426 mg/L。而到達(dá)下游更遠(yuǎn)的五水廠取水口時，豐水期、平水期、枯水期濃度依次降低，其主要原因可能是停留時間較長，風(fēng)險物質(zhì)被稀釋[24]。平水期和枯水期，在五水廠取水口、市出境斷面處會先后出現(xiàn)2次污染峰團(tuán)，這是由于江心洲兩側(cè)水體流速相差較大，部分風(fēng)險物質(zhì)隨右側(cè)支流先行到達(dá)敏感受體，隨后左側(cè)支流中的風(fēng)險物質(zhì)再次到達(dá)敏感受體。

圖7 豐水期儲罐完全泄露時水質(zhì)模擬結(jié)果Fig.7 Water quality simulation results when the tank breaks completely in wet period

圖9 枯水期儲罐完全泄露時水質(zhì)模擬結(jié)果Fig.9 Water quality simulation results when the tank breaks completely in dry period

圖8 平水期儲罐完全泄露時水質(zhì)模擬結(jié)果Fig.8 Water quality simulation results when the tank breaks completely in normal period

綜上，風(fēng)險物質(zhì)對下游風(fēng)險受體的影響程度與泄漏總量和水文期密切相關(guān)，所得結(jié)論與長江其他段類似研究[29,31]相同。同一水文期，風(fēng)險物質(zhì)泄漏總量越大，對下游敏感受體的影響時間越長，污染團(tuán)濃度也越高。泄漏風(fēng)險物質(zhì)在平水期和枯水期時，對受體的影響時間更長；在豐水期時，到達(dá)受體的時間更短，即應(yīng)急反應(yīng)時間更少。

4 結(jié)論

（1）利用MIKE21水動力-水質(zhì)耦合模型，對長江銅陵段某化工企業(yè)儲罐內(nèi)風(fēng)險物質(zhì)在不同水文期、3種不同泄漏情景下的突發(fā)水污染事故進(jìn)行了模擬，所建立的水動力模型模擬水位結(jié)果與實測結(jié)果吻合性良好，曲線趨勢基本相似，表明模型選取參數(shù)較為合理，可作為水質(zhì)模塊的水動力基礎(chǔ)。

（2）在設(shè)定的3種不同泄漏情景下，儲罐完全泄漏時因其泄漏總量大，導(dǎo)致風(fēng)險物質(zhì)濃度最高，影響時間最長。在不同水文期的模擬結(jié)果顯示，豐水期風(fēng)險物質(zhì)到達(dá)下游敏感受體的時間最短，留給應(yīng)急人員處置時間也最少；平水期和枯水期風(fēng)險物質(zhì)的最大污染峰團(tuán)濃度會更高，影響時間會更長，且在五水廠取水口和市出境斷面處先后有2次污染峰團(tuán)到達(dá)。

（3）本研究基于MIKE21構(gòu)建了適用于長江干流銅陵段的水動力-水質(zhì)模型，通過對風(fēng)險企業(yè)內(nèi)多種風(fēng)險物質(zhì)在不同水文情勢和泄漏情景下的模擬分析，可建立當(dāng)?shù)氐哪M預(yù)測預(yù)警資料庫，為突發(fā)水污染事故的預(yù)測預(yù)警、應(yīng)急處置提供決策參考。結(jié)合本次模擬結(jié)果，因風(fēng)險物質(zhì)在豐水期到達(dá)受體時間更短，平水期和枯水期對受體影響時間更長，當(dāng)接到突發(fā)水環(huán)境事件報告時，在豐水期應(yīng)立即采取措施保護(hù)風(fēng)險受體，在平水期和枯水期則首先對泄漏風(fēng)險源進(jìn)行控制，從而減少影響時長。此外，有必要通過加強日常環(huán)境安全檢查等措施，及時發(fā)現(xiàn)并整改問題，從而有效降低對水環(huán)境的影響程度。后續(xù)研究還需采集更為豐富的水文數(shù)據(jù)、水質(zhì)監(jiān)測數(shù)據(jù)等，不斷進(jìn)行完善以進(jìn)一步提高模型模擬的準(zhǔn)確性。