基于MIKE21模型的漢江中下游水環(huán)境影響風(fēng)險(xiǎn)研究

趙曉磊

(上?？睖y(cè)設(shè)計(jì)研究院有限公司，上海 200335)

1 引言

由于水資源過度開發(fā)和缺乏保護(hù)，我國(guó)江河湖海正遭受嚴(yán)重的水環(huán)境問題，迫切需要加強(qiáng)治水工作[1]。遵循新時(shí)代治水要求，基于“安全、清潔、健康”的理念，我國(guó)圍繞水資源管理、水環(huán)境治理、水生態(tài)保護(hù)、水景觀建設(shè)和水管理措施等方面，提出控源截污、內(nèi)源治理和生態(tài)修復(fù)等治理工程來提升水環(huán)境質(zhì)量和加強(qiáng)水資源保護(hù)[2]。湖北省結(jié)合本省實(shí)際，于2017年針對(duì)漢江流域正式執(zhí)行《湖北襄陽(yáng)市漢江流域水環(huán)境保護(hù)條例》，并充分運(yùn)用水質(zhì)同步監(jiān)測(cè)、遙感技術(shù)和地表水環(huán)境數(shù)學(xué)模型等手段，搭建漢江流域預(yù)警預(yù)報(bào)平臺(tái)[3～5]，以此達(dá)到精準(zhǔn)預(yù)測(cè)、科學(xué)分析、嚴(yán)格評(píng)判的作用，打通漢江流域預(yù)警通道“最后一公里”。

漢江是長(zhǎng)江中游最長(zhǎng)的一級(jí)支流，全長(zhǎng)1577 km，近年來漢江中下游水環(huán)境事故頻發(fā)[6, 7]，如硅藻水華事件、水體富營(yíng)養(yǎng)化污染等，導(dǎo)致流域水質(zhì)變差，嚴(yán)重影響沿江用水安全。本研究基于漢江中下游多年水質(zhì)數(shù)據(jù)，應(yīng)用MIKE21軟件[8, 9]，建立二維水動(dòng)力水質(zhì)耦合模型，對(duì)漢江中下游水動(dòng)力和水質(zhì)變化過程進(jìn)行模擬，分析污染物擴(kuò)散及其影響變化，以期增強(qiáng)流域預(yù)警預(yù)報(bào)能力，為流域水生態(tài)環(huán)境風(fēng)險(xiǎn)管控提供理論支撐。

2 漢江中下游水環(huán)境現(xiàn)狀評(píng)價(jià)

根據(jù)2017年正式執(zhí)行的《湖北襄陽(yáng)市漢江流域水環(huán)境保護(hù)條例》要求，漢江干流要求滿足《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》(GB3838-2002)Ⅱ類水質(zhì)標(biāo)準(zhǔn)，支流則滿足Ⅲ類水質(zhì)標(biāo)準(zhǔn)。

2.1 干流水質(zhì)現(xiàn)狀分析

基于收集到的漢江中下游干流丹江口壩下、襄陽(yáng)、余家湖及仙桃等斷面2008～2014年水質(zhì)監(jiān)測(cè)數(shù)據(jù)，選取CODMn、氨氮、TP作為代表指標(biāo)，對(duì)各斷面進(jìn)行水質(zhì)分析。

監(jiān)測(cè)結(jié)果(圖1)顯示，丹江口壩下斷面CODMn濃度在1.6～4 mg/L之間，氨氮濃度在0.062～0.546 mg/L之間，TP濃度在0.01～0.06 mg/L之間，除2011年7月份氨氮超標(biāo)外，其余月份水質(zhì)均滿足Ⅱ類水；襄陽(yáng)斷面CODMn濃度在1～3 mg/L之間，氨氮濃度在0.057～0.456 mg/L之間，TP濃度在0.01～0.08 mg/L之間，水質(zhì)較好，均能達(dá)到Ⅱ類水質(zhì)；余家湖斷面CODMn濃度在1.6～4 mg/L之間，氨氮濃度在0.068～0.91 mg/L之間，TP濃度在0.02～0.28 mg/L之間，可見余家湖斷面大部分時(shí)間段能夠達(dá)到Ⅱ類水質(zhì)，但2008及2009年水質(zhì)較差，部分月份出現(xiàn)Ⅲ類甚至Ⅳ類水，超標(biāo)因子為氨氮和TP；仙桃斷面CODMn濃度在1.4～4.2 mg/L之間，氨氮濃度在0.025～0.279 mg/L之間，TP濃度在0.001～0.16 mg/L之間，除個(gè)別月份CODMn和TP超標(biāo)外，其他時(shí)段水質(zhì)均能滿足Ⅱ類標(biāo)準(zhǔn)。

總體而言，從CODMn、氨氮、TP3個(gè)指標(biāo)來講，漢江中下游干流各斷面水質(zhì)基本較好，除2008～2009年余家湖斷面水質(zhì)較差外，其它斷面基本能夠滿足Ⅱ類水質(zhì)要求。

2.2 入?yún)R支流水質(zhì)現(xiàn)狀分析

根據(jù)收集到的漢江中下游入?yún)R支流小清河、唐白河斷面2008～2014年水質(zhì)監(jiān)測(cè)數(shù)據(jù)，選取CODMn、氨氮、TP作為代表指標(biāo)，對(duì)入?yún)R支流各斷面進(jìn)行水質(zhì)分析。

監(jiān)測(cè)結(jié)果(圖2)顯示，小清河斷面CODMn濃度在2～21 mg/L之間，氨氮濃度在0.09～10.326 mg/L之間，TP濃度在0.03～0.86 mg/L之間，小清河斷面2008～2012年水質(zhì)較好，而2013～2014年水質(zhì)較差，未能達(dá)到Ⅲ類水質(zhì)要求；唐白河斷面CODMn濃度在1.9～14.4 mg/L之間，氨氮濃度在0.106～3.67 mg/L之間，總磷濃度在0.03～0.79 mg/L之間，CODMn在大部分時(shí)段均能達(dá)到Ⅲ類標(biāo)準(zhǔn)，而氨氮、TP大部分監(jiān)測(cè)時(shí)段濃度較高，水質(zhì)未能達(dá)到Ⅲ類水質(zhì)要求。

總體而言，從CODMn、氨氮、TP 3個(gè)指標(biāo)來講，漢江中下游入?yún)R支流各斷面水質(zhì)大部分監(jiān)測(cè)時(shí)段均無法保證滿足Ⅲ類水質(zhì)，部分時(shí)段水質(zhì)污染較為嚴(yán)重，尤其是氮污染比較嚴(yán)重，需要特別重視。

3 模型構(gòu)建

3.1 研究方案

本研究范圍為漢江中下游干流丹江口壩下斷面至漢川斷面，沿途經(jīng)襄陽(yáng)市、宜城市、鐘祥市、潛江市、仙桃市等縣市，跨東經(jīng)111°29′～113°57′，北緯30°22′～32°34′，呈西北東南方向分布。

收集重要水位站(襄陽(yáng)、余家湖、宜城、皇莊、沙洋、仙桃)逐日水位數(shù)據(jù)、流量站(襄陽(yáng)、余家湖、皇莊、仙桃)逐日流量數(shù)據(jù)以及重要水質(zhì)站(襄陽(yáng)、余家湖、仙桃)逐月水質(zhì)數(shù)據(jù)、重點(diǎn)排污口和支流等監(jiān)測(cè)斷面的流量和水質(zhì)數(shù)據(jù)，以丹江口壩下斷面為上游邊界，漢川斷面為下游邊界，構(gòu)建MIKE21二維水動(dòng)力-水質(zhì)耦合模型，并對(duì)模型進(jìn)行率定驗(yàn)證。

3.2 計(jì)算區(qū)域網(wǎng)格劃分

計(jì)算區(qū)域網(wǎng)格劃分的形式采用三角形與四邊形網(wǎng)格的結(jié)合，由于漢江干流又長(zhǎng)又窄，用三角形網(wǎng)格需要大量網(wǎng)格才能反映地形，故沿河道主流采用四邊形網(wǎng)格，不僅可以減少網(wǎng)格量，減少計(jì)算時(shí)間，同時(shí)也能反映地形。本研究范圍共劃分為23765個(gè)單元網(wǎng)格，包含21588個(gè)節(jié)點(diǎn)(圖3)。流域DEM數(shù)字高程共包含215199組原始數(shù)據(jù)，將高程數(shù)據(jù)導(dǎo)入到坐標(biāo)系中，對(duì)網(wǎng)格進(jìn)行高程插值，得到漢江中下游地形圖(圖4)。

3.3 初始條件與邊界條件設(shè)定

由于河道太長(zhǎng)，地形差異較大，造成各河段之間水位相差較大，將模擬時(shí)期第一天的河道水位作為模型的初始水位。漢江中下游面積廣闊，單一污染物的模擬評(píng)價(jià)并不能準(zhǔn)確反應(yīng)流域內(nèi)水質(zhì)情況，因此選取CODMn、氨氮和TP作為水質(zhì)模型模擬的指標(biāo)，水質(zhì)初始條件采用模型計(jì)算第一天的各污染物濃度，其中CODMn為1.8 mg/L，氨氮為0.247 mg/L，TP為0.03 mg/L。

本研究分別以丹江口壩下斷面逐日平均流量和漢川斷面逐日平均水位作為水動(dòng)力模型的上、下邊界，以丹江口壩下逐月水質(zhì)作為水質(zhì)模型入流邊界。經(jīng)實(shí)地調(diào)研和查閱資料，發(fā)現(xiàn)漢江中下游存在多個(gè)污染源和支流匯入，由于缺乏地形數(shù)據(jù)，將污染源和支流的流量及污染物濃度概化為點(diǎn)源輸入到模型作為支流邊界，經(jīng)分析，該模型存在30個(gè)排污口和4條支流(唐白河、小清河、蠻河、竹皮河)匯入。

3.4 模型關(guān)鍵參數(shù)

3.4.1 水動(dòng)力模塊參數(shù)

底床糙率(底床摩擦力)表示水域底部對(duì)水流阻力的大小，是水力學(xué)模塊中重要的率定參數(shù)[10, 11]。漢江中下游地形情況復(fù)雜，底床糙率需根據(jù)不同地形設(shè)置不同數(shù)值。查閱相關(guān)資料，并經(jīng)過多次率定后，確定不同地形區(qū)域邊界處底床糙率采用分段常數(shù)插值方法生成(圖5)，其中丹江口壩下至襄陽(yáng)段糙率為0.05，襄陽(yáng)至余家湖段糙率為0.029，余家湖至宜城段糙率為0.05，宜城至皇莊段糙率為0.015，皇莊至沙洋段糙率為0.033，沙洋至漢川段糙率為0.015。

3.4.2 水質(zhì)模塊參數(shù)

降解系數(shù)是水質(zhì)模塊重要的率定參數(shù)[12]，表示污染物在水流的作用下衰減的快慢。查閱相關(guān)文獻(xiàn)，并經(jīng)過多次率定后，確定本研究水質(zhì)模塊CODMn的降解系數(shù)為0.045/d，氨氮的降解系數(shù)為0.062/d，TP的降解系數(shù)為0.042/d。

4 模型驗(yàn)證

4.1 水動(dòng)力模型驗(yàn)證

在模型的驗(yàn)證中，通常采用2個(gè)目標(biāo)函數(shù)來衡量模型結(jié)果，即絕對(duì)平均誤差(AME)和均方根誤差(RMSE)：

圖5 底床糙率分布

(1)

(2)

式(1)、(2)中，ai為實(shí)測(cè)值；bi為預(yù)測(cè)值；n為預(yù)測(cè)值個(gè)數(shù)。AME和RMSE值越接近0，說明模型模擬精度越高。

輸入2014年7～12月份上下游斷面水位流量和入?yún)R支流流量作為初始條件和邊界條件，得到各斷面水位(圖6)、流量模擬結(jié)果(圖7)，并計(jì)算實(shí)測(cè)水位與模擬值的平均誤差、絕對(duì)平均誤差及均方根誤差(表1)，實(shí)測(cè)流量與模擬值的平均相對(duì)誤差(表2)。

綜上可知，各斷面水位驗(yàn)證結(jié)果較好，實(shí)測(cè)水位與模擬值變化過程擬合程度較高，平均誤差均小于21 cm、絕對(duì)平均誤差均小于0.213 cm、均方根誤差均小于0.229，模擬精度較高；實(shí)測(cè)流量與模擬值變化過程擬合程度較高，流量平均相對(duì)誤差均小于13.8%。說明該水動(dòng)力模型能較好地模擬漢江中下游水動(dòng)力學(xué)變化過程。

4.2 水質(zhì)模型驗(yàn)證

輸入水質(zhì)模型的初始條件和邊界條件，即可模擬得到漢江中下游水質(zhì)監(jiān)測(cè)斷面的水質(zhì)動(dòng)態(tài)變化過程和污染物濃度模擬結(jié)果(圖8～10)，并計(jì)算各監(jiān)測(cè)斷面污染物實(shí)測(cè)濃度與模擬值的平均相對(duì)誤差(表3)。

由圖8～10、表3可知，各水質(zhì)監(jiān)測(cè)斷面CODMn實(shí)測(cè)值與模擬值平均相對(duì)誤差在13.3%～22.3%之間，模擬值與實(shí)測(cè)值比較接近，在可接受范圍內(nèi)。氨氮實(shí)測(cè)值與模擬值平均相對(duì)誤差在17.7%～29.7%之間，也在可接受范圍內(nèi)。然而TP實(shí)測(cè)值與模擬值平均相對(duì)誤差在12.7%～44.3%之間，模擬值與實(shí)測(cè)值變化趨勢(shì)相近、但相對(duì)誤差稍大，其主要原因是水體中水環(huán)境因子變化復(fù)雜、相互影響，且目前大壩調(diào)度數(shù)據(jù)、排污口和支流數(shù)據(jù)未經(jīng)完全統(tǒng)計(jì)，不能精確反映漢江中下游污染物濃度變化?？傮w而言，CODMn、氨氮和TP的誤差整體在可控范圍內(nèi)，表明構(gòu)建的水質(zhì)模型能夠較準(zhǔn)確的模擬流域內(nèi)CODMn、氨氮和TP的變化情況。

圖6 各斷面水位驗(yàn)證

圖7 各斷面流量驗(yàn)證

表1 各斷面水位驗(yàn)證結(jié)果統(tǒng)計(jì)

表2 各斷面流量驗(yàn)證結(jié)果統(tǒng)計(jì)

圖8 襄陽(yáng)站水質(zhì)驗(yàn)證

圖9 余家湖站水質(zhì)驗(yàn)證

圖10 仙桃站水質(zhì)驗(yàn)證

表3 各斷面水質(zhì)驗(yàn)證結(jié)果統(tǒng)計(jì)

5 突發(fā)水污染事故對(duì)干流水質(zhì)影響模擬研究

5.1 情景模擬

近年來，漢江中下游頻頻出現(xiàn)突發(fā)水污染事故，導(dǎo)致下游水質(zhì)在短期內(nèi)突然惡化，嚴(yán)重威脅沿江居民用水安全，制約生態(tài)平衡和社會(huì)經(jīng)濟(jì)的發(fā)展。假定余家湖段下游潛江市城北污水處理廠(112°52′23″E，30°28′26″N)入河排污口斷面于模擬期內(nèi)(2014年3月1日)突發(fā)水污染事故，污染物就近流入漢江，COD濃度為10000 mg/L，排放流量為0.2 m3/s，事故持續(xù)排放時(shí)間為1 d(3月1日0點(diǎn)至3月2日0點(diǎn))，其余時(shí)間采用城鎮(zhèn)污水處理廠一級(jí)A排放標(biāo)準(zhǔn)。模擬該河段突發(fā)水污染事故后需歷時(shí)多久才會(huì)基本恢復(fù)事故前背景值(Ⅱ類水)[13～15]。

5.2 邊界條件確定

5.2.1 流量設(shè)計(jì)

分別收集整理漢江中下游代表性斷面(襄陽(yáng)、余家湖、皇莊、仙桃)逐日流量和丹江口壩下逐月流量，并將其分別按10%、50%、90%水文頻率設(shè)置大、中、小3種流量，進(jìn)行水文頻率配線，得到P-Ⅲ型曲線下的統(tǒng)計(jì)參數(shù)(表4)。根據(jù)水文分析結(jié)果，襄陽(yáng)、余家湖、皇莊、仙桃斷面均選擇月最大流量P=10%時(shí)的計(jì)算值作為大流量設(shè)計(jì)值、均值流量P=50%時(shí)的計(jì)算值作為中流量設(shè)計(jì)值、最小流量P=90%時(shí)的計(jì)算值作為小流量設(shè)計(jì)值；丹江口壩下斷面選擇月均流量P=10%時(shí)的計(jì)算值作為大流量設(shè)計(jì)值、P=50%時(shí)的計(jì)算值作為中流量設(shè)計(jì)值、P=90%時(shí)的計(jì)算值作為小流量設(shè)計(jì)值。與前4個(gè)代表性斷面相比，丹江口壩下斷面大流量與中流量設(shè)計(jì)值偏小、而小流量設(shè)計(jì)值偏大的原因可能是丹江口水庫(kù)蓄水運(yùn)行調(diào)蓄，導(dǎo)致丹江口壩下流量過程均化，也可能是由于數(shù)據(jù)頻率較低，樣本數(shù)量不夠，導(dǎo)致水文頻率配線計(jì)算結(jié)果與實(shí)際存在偏差。

表4 各斷面流量水文配線參數(shù)

5.2.2 水動(dòng)力邊界條件

襄陽(yáng)、余家湖、皇莊、仙桃斷面計(jì)算結(jié)果中，大、中、小流量設(shè)計(jì)值結(jié)果相差不大，選取襄陽(yáng)斷面的流量計(jì)算結(jié)果作為模擬工況的入流特征流量，即選擇3926.63 m3/s為大流量設(shè)計(jì)值、816.05 m3/s為中流量設(shè)計(jì)值、387.23 m3/s為小流量設(shè)計(jì)值。另外由于丹江口壩下斷面數(shù)據(jù)頻率較低，水文頻率配線計(jì)算結(jié)果與其它斷面結(jié)果差異較大，為綜合考慮各種工況，同時(shí)選取丹江口壩下P=10%時(shí)的計(jì)算值1787.68 m3/s為流量工況設(shè)計(jì)值，共計(jì)4種入流流量設(shè)計(jì)工況。

下游邊界條件中由于缺少漢川斷面流量數(shù)據(jù)，故采用仙桃斷面流量與漢川斷面水位進(jìn)行擬合，首先將漢川日平均水位與仙桃日平均流量根據(jù)水位進(jìn)行排序，將相同水位對(duì)應(yīng)流量的平均值作為該水位對(duì)應(yīng)的流量值，將水位與對(duì)應(yīng)的流量進(jìn)行對(duì)數(shù)擬合(圖11)，水位(y)與流量(x)的關(guān)系為：y=4.5903In(x)-11.696，擬合的R2值為0.8967，計(jì)算工況下的漢川水位流量關(guān)系曲線見圖12。

5.2.3 水質(zhì)邊界條件

模型模擬的水質(zhì)指標(biāo)主要為CODMn、氨氮和TP，考慮到漢江中下游干流各監(jiān)測(cè)斷面以這3個(gè)指標(biāo)進(jìn)行評(píng)價(jià)時(shí)結(jié)果較好，故選擇多年監(jiān)測(cè)結(jié)果中最高濃度即最不利情景作為模型計(jì)算工況中的入流水質(zhì)濃度，故丹江口壩下斷面CODMn取4 mg/L，氨氮取0.546 mg/L，TP取0.06 mg/L。

圖11 漢川水位與對(duì)應(yīng)平均流量對(duì)數(shù)擬合結(jié)果

圖12 計(jì)算工況下漢川水位流量關(guān)系曲線

5.3 污染物擴(kuò)散影響研究

5.3.1 大流量設(shè)計(jì)值(Q=3926.63 m3/s)

此工況下，該突發(fā)事故并未對(duì)下游河道造成污染(圖13)，突發(fā)事故剛發(fā)生時(shí)形成的污染帶均為Ⅱ類水，其中有部分河段逐漸被稀釋降解為Ⅰ類水。

圖13 3月2日0點(diǎn)突發(fā)事故工況形成的污染帶

5.3.2 中流量設(shè)計(jì)值(Q=816.05 m3/s)

此工況下，突發(fā)事故剛發(fā)生時(shí)形成的污染帶下游部分河段被污染為Ⅲ類水(圖14)，造成的污染并非很嚴(yán)重；隨著水體流動(dòng)，歷經(jīng)24 h后污染帶范圍向下游擴(kuò)散(圖15)，上游逐漸被稀釋為Ⅱ類水；歷經(jīng)約55 h后污染帶基本恢復(fù)到Ⅱ類水(圖16)。

圖14 3月2日0點(diǎn)突發(fā)事故工況形成的污染帶

圖15 突發(fā)事故24 h后形成的污染帶

5.3.3 小流量設(shè)計(jì)值(Q=387.23 m3/s)

此工況下，突發(fā)事故剛發(fā)生時(shí)形成了顯著的污染帶，下游很長(zhǎng)一段河道被污染為Ⅳ類水(圖17)；隨著水體流動(dòng)，歷經(jīng)約71h后Ⅳ類水逐漸被稀釋降解，還存在較長(zhǎng)河段的Ⅲ類水(圖18)；歷經(jīng)97h后所受污染河段基本恢復(fù)為Ⅱ類水(圖19)。

圖16 突發(fā)事故約55 h后形成的污染帶

5.3.4 丹江口壩下斷面P=10%時(shí)流量設(shè)計(jì)值(Q=1787.68 m3/s)

此工況下，該突發(fā)事故并未對(duì)河段造成污染(圖20)，與大流量設(shè)計(jì)值時(shí)效果一致，突發(fā)事故剛發(fā)生時(shí)形成的污染帶水質(zhì)均滿足Ⅱ類水，并有部分河道逐漸被稀釋降解為Ⅰ類水。

圖17 3月2日0點(diǎn)突發(fā)事故工況形成的污染帶

圖18 突發(fā)事故約71 h后形成的污染帶

圖19 突發(fā)事故97 h后形成的污染帶

5.3.5 小結(jié)

當(dāng)流量設(shè)計(jì)值取3926.63 m3/s和1787.68 m3/s時(shí)，該突發(fā)事故下游河道仍為Ⅱ類水，可見并未對(duì)下游河道造成污染；當(dāng)流量設(shè)計(jì)值取816.05 m3/s時(shí)，該突發(fā)事故下游部分河段被污染為Ⅲ類水，并歷經(jīng)約55 h后該突發(fā)事故下游河道基本恢復(fù)到Ⅱ類水；當(dāng)流量設(shè)計(jì)值取387.23 m3/s時(shí)，該突發(fā)事故工況下游很長(zhǎng)一段河道被污染為Ⅳ類水，歷經(jīng)約71h后Ⅳ類水污染范圍被稀釋降解為Ⅲ類水，共計(jì)歷經(jīng)97 h后，污染帶基本恢復(fù)到Ⅱ類水。

圖20 3月2日0點(diǎn)突發(fā)事故工況形成的污染帶

6 結(jié)論

(1)根據(jù)多年實(shí)際監(jiān)測(cè)資料，漢江中下游干流各斷面水質(zhì)基本較好，以CODMn、氨氮和TP這3個(gè)指標(biāo)分析，各斷面基本滿足Ⅱ類水質(zhì)要求；而入?yún)R支流各斷面水質(zhì)較差，大部分監(jiān)測(cè)時(shí)段均無法保證滿足Ⅲ類水質(zhì)，尤其氮污染比較嚴(yán)重，需特別重視。

(2)構(gòu)建漢江中下游二維水動(dòng)力水質(zhì)耦合模型，上游邊界為丹江口壩下斷面，下游邊界為漢川斷面，沿程排污口及支流數(shù)據(jù)作為點(diǎn)源輸入模型。通過對(duì)該模型的率定驗(yàn)證，發(fā)現(xiàn)該水動(dòng)力水質(zhì)耦合模型能較好地模擬漢江中下游水動(dòng)力學(xué)和水質(zhì)動(dòng)態(tài)變化過程。

(3)確定工況模擬的邊界條件，模擬了4種不同流量設(shè)計(jì)值工況下突發(fā)水污染事故對(duì)干流水質(zhì)的影響。由結(jié)果可知，隨著丹江口壩下入流流量增大，污染物不持續(xù)排放且污染泄露量較小時(shí)，突發(fā)事故工況形成的污染帶范圍和污染程度都會(huì)明顯降低。因此，為減緩?fù)话l(fā)事故對(duì)干流水質(zhì)影響，可以適當(dāng)加大丹江口水庫(kù)的下泄流量。