高瑋岐，胡 煜，任華堂，Tomasz TYMISKI

(1.中央民族大學(xué)生命與環(huán)境科學(xué)學(xué)院，北京 100081；2.弗羅茨瓦夫環(huán)境與生命科學(xué)大學(xué)環(huán)境工程學(xué)院，弗羅茨瓦夫 50-363)

0 引 言

萬(wàn)家寨水庫(kù)位于黃河北干流上段托克托至龍口河段峽谷內(nèi)，是太原、大同、朔州3市的工業(yè)和生活用水主要水源地。近年來(lái)上游取水量增加，且上游工業(yè)排污持續(xù)增加，水環(huán)境污染呈現(xiàn)出日益嚴(yán)重之勢(shì)[1]，嚴(yán)重威脅山西人民的生產(chǎn)生活。開(kāi)展萬(wàn)家寨水庫(kù)水環(huán)境納污特性的研究，對(duì)于引水區(qū)域的供水安全具有重要的現(xiàn)實(shí)意義。

對(duì)于萬(wàn)家寨水庫(kù)的研究，現(xiàn)有成果主要集中于凌汛期水動(dòng)力特性的研究。冀鴻蘭[2]等根據(jù)歷史實(shí)測(cè)冰情資料，統(tǒng)計(jì)分析得出上游河段冰情特征規(guī)律，給出上游河段防凌調(diào)度原則；張寶森[3]等以水泥廠為觀測(cè)點(diǎn)，對(duì)凌汛期庫(kù)尾水位變幅增大受庫(kù)區(qū)水位因子及環(huán)境因子影響的比重進(jìn)行研究發(fā)現(xiàn)，岔河口形成冰塞、冰壩是上游水位突升的主要影響因子；梁貴生[4]等研究發(fā)現(xiàn)，庫(kù)區(qū)冰塞頭部形成位置取決于河道條件、水庫(kù)水位和河道流量。近年，對(duì)于萬(wàn)家寨水庫(kù)水環(huán)境污染的研究日漸重視，漆強(qiáng)強(qiáng)[5]選取具代表性斷面，進(jìn)行pH值、高錳酸鹽、化學(xué)需氧量、總氮等污染因子的監(jiān)測(cè)發(fā)現(xiàn)，大部分?jǐn)嗝嫠|(zhì)指標(biāo)均有不同程度的超標(biāo)現(xiàn)象；吳喜軍[6]等應(yīng)用改進(jìn)的內(nèi)梅羅污染指數(shù)法、灰色聚類法和單因子評(píng)價(jià)法對(duì)庫(kù)區(qū)水質(zhì)進(jìn)行評(píng)價(jià)，得出黃河干流整體水質(zhì)狀況較好的結(jié)論；王若儀[7]應(yīng)用營(yíng)養(yǎng)分級(jí)評(píng)分法和綜合營(yíng)養(yǎng)狀態(tài)指數(shù)法進(jìn)行總磷、總氮等5項(xiàng)污染因子的評(píng)價(jià)發(fā)現(xiàn)，萬(wàn)家寨水庫(kù)的葉綠素a、透明度已達(dá)到輕度富營(yíng)養(yǎng)水平，總氮達(dá)到中度富營(yíng)養(yǎng)水平。

目前，對(duì)于研究區(qū)域水環(huán)境污染的研究成果多為水質(zhì)總體評(píng)價(jià)，對(duì)于庫(kù)區(qū)水環(huán)境特性精細(xì)的時(shí)空分布研究相對(duì)匱乏。為此，本文應(yīng)用水環(huán)境數(shù)學(xué)模型，對(duì)萬(wàn)家寨水庫(kù)庫(kù)區(qū)水環(huán)境進(jìn)行模擬，分析水動(dòng)力特性和水環(huán)境特性，并計(jì)算典型工況下區(qū)域內(nèi)COD最大排放濃度分布及水環(huán)境容量。

1 數(shù)學(xué)模型

1.1 控制方程

采用二維水流水質(zhì)模型進(jìn)行求解[8]，即

(1)

(2)

(3)

式中，η為水位高度；D為單元的全水深；u、v分別為流速矢量U在x、y方向的分量；g為重力加速度；Cb為拖曳力系數(shù)，與水面的粗糙度有關(guān)；ρ為水體密度；AM為渦粘系數(shù)；f=2Ωsinφ為科氏力參數(shù)，Ω為地球自轉(zhuǎn)角速度，φ為地理緯度。

采用對(duì)流擴(kuò)散方程模擬區(qū)域污染物輸移擴(kuò)散規(guī)律，即

(4)

式中，C為污染物垂向平均濃度；AH為水平擴(kuò)散系數(shù)；S為源項(xiàng)。

應(yīng)用有限體積法(FVM)對(duì)于上述物質(zhì)連續(xù)方程、動(dòng)量方程以及濃度方程進(jìn)行離散，主要包括對(duì)對(duì)流及擴(kuò)散通量做體積積分和時(shí)間積分處理，根據(jù)對(duì)流通量、擴(kuò)散通量、源(匯)項(xiàng)的表達(dá)式即可建立各單元的物質(zhì)守恒方程，將所有單元的濃度守恒方程進(jìn)行聯(lián)立求解，進(jìn)而得知全場(chǎng)污染物的濃度分布特性。

1.2 網(wǎng)格劃分及計(jì)算條件

計(jì)算區(qū)域?yàn)?9°34′4″N～40°13′05″N，111°10′23″E～111°27′06″E，上游邊界取在拐上斷面，下游邊界取在壩上斷面。研究區(qū)域網(wǎng)格劃分見(jiàn)圖1。采用二維水動(dòng)力模型進(jìn)行計(jì)算，庫(kù)區(qū)平面應(yīng)用非結(jié)構(gòu)三角形網(wǎng)格進(jìn)行貼體劃分，其劃分網(wǎng)格10 476個(gè)，三角形網(wǎng)格節(jié)點(diǎn)數(shù)為6 109個(gè)，三角形界面數(shù)16 584個(gè)。沿河寬方向網(wǎng)格間距為33.0～66.5 m，平均間距為51.1 m；沿河道方向網(wǎng)格間距為38.7～62.2 m，平均間距為53.6 m。

圖1 萬(wàn)家寨庫(kù)區(qū)網(wǎng)格劃分

1.3 計(jì)算條件

根據(jù)已掌握的萬(wàn)家寨水庫(kù)長(zhǎng)期監(jiān)測(cè)數(shù)據(jù)，采用上游入庫(kù)流量、下游壩前水位邊界進(jìn)行控制，將拐上斷面多年月均監(jiān)測(cè)資料作為模型入庫(kù)的水質(zhì)邊界條件。模型計(jì)算時(shí)間步長(zhǎng)為0.1 s，底部拖曳力系數(shù)為0.025。參考研究庫(kù)區(qū)相似河段及相關(guān)水庫(kù)水質(zhì)模型中有關(guān)各污染物衰減系數(shù)的研究成果[9-10]，根據(jù)季節(jié)的差異，非冰封期化學(xué)需氧量(COD)的降解系數(shù)設(shè)為0.240 d-1，氨氮設(shè)置為0.105 d-1。

2 模型驗(yàn)證

從水動(dòng)力和水質(zhì)2個(gè)方面進(jìn)行驗(yàn)證。采用水庫(kù)庫(kù)容日變化過(guò)程對(duì)水動(dòng)力結(jié)果進(jìn)行驗(yàn)證，見(jiàn)圖2。從圖2可知，實(shí)測(cè)值和計(jì)算值平均相對(duì)誤差為1.08%，具有較好的一致性。

圖2 模擬庫(kù)容與實(shí)測(cè)庫(kù)容驗(yàn)證

選取COD與氨氮濃度進(jìn)行不同水文年條件下水質(zhì)模型壩前污染物濃度驗(yàn)證，見(jiàn)圖3、4。利用平均對(duì)數(shù)比偏差Dr計(jì)算計(jì)算值與實(shí)測(cè)值的差異性，公式如下

(5)

圖3 不同水文年條件下水質(zhì)模型COD濃度驗(yàn)證

圖4 不同水文年條件下水質(zhì)模型氨氮濃度驗(yàn)證

式中，DLm為實(shí)測(cè)值；DLp為計(jì)算值；n為樣本個(gè)數(shù)。

經(jīng)計(jì)算，各種條件下2種水質(zhì)指標(biāo)計(jì)算值與實(shí)測(cè)值之間的平均對(duì)數(shù)比偏差均在0.1～0.3范圍內(nèi)，最大為0.297，滿足工程應(yīng)用精度要求。

圖5 不同水文年下夏季壩前30 km范圍內(nèi)流速分布狀態(tài)

3 結(jié)果及分析

3.1 水動(dòng)力場(chǎng)分析

圖5給出了典型豐水年、平水年及枯水年萬(wàn)家寨水庫(kù)壩前30 km區(qū)域內(nèi)夏季流速分布狀態(tài)。從圖5可知，不同水文年流速差異較大，豐水年流速>平水年流速>枯水年流速，這是由于不同水文條件下的流量差異所致。此外，由于紊流作用，在曲率較大的岸邊界處會(huì)形成漩渦，由于離心慣性力和重力導(dǎo)致局部水體出現(xiàn)摻混現(xiàn)象，有利于水體交換。

3.2 水力停留時(shí)間

水力停留時(shí)間(HRT)是反映水體交換能力的重要參數(shù)，對(duì)水體自凈有著重要影響[11]，計(jì)算公式如下

(6)

式中，HRT為水力停留時(shí)間；V為水庫(kù)庫(kù)容；Q為入流流量。

庫(kù)區(qū)水力停留時(shí)間計(jì)算統(tǒng)計(jì)見(jiàn)表1。從表1可知，在典型豐水年、平水年、枯水年3種水文年條件下，水庫(kù)段水力停留時(shí)間在1.76～11.78 d范圍內(nèi)，萬(wàn)家寨水庫(kù)具備典型的河道水庫(kù)型水動(dòng)力特征。具體而言，夏季水力停留時(shí)間小于冬季水力停留時(shí)間，最短水力停留時(shí)間出現(xiàn)在豐水年夏季，為1.76 d；最長(zhǎng)水力停留時(shí)間出現(xiàn)在枯水年冬季，為11.78 d。產(chǎn)生此結(jié)果的原因是水庫(kù)在豐水年夏季具有較常年各季節(jié)偏高的徑流量，水力停留時(shí)間較其他時(shí)期短?？菟甓揪哂休^常年各季節(jié)偏低的徑流量，水力停留時(shí)間較其他時(shí)期長(zhǎng)。

表1 庫(kù)區(qū)水力停留時(shí)間計(jì)算統(tǒng)計(jì)

3.3 COD的最大排放濃度及水環(huán)境容量計(jì)算

3.3.1計(jì)算工況

為了探究水庫(kù)水環(huán)境容量，假設(shè)萬(wàn)家寨水庫(kù)上游入流斷面處以固定流量連續(xù)排放一定濃度的有機(jī)污水，設(shè)計(jì)3種不同的計(jì)算工況，見(jiàn)表2。

根據(jù)庫(kù)區(qū)HJ 338—2018《飲用水水源保護(hù)區(qū)劃分技術(shù)規(guī)范》，飲用水源一級(jí)保護(hù)區(qū)的水質(zhì)應(yīng)滿足二類水質(zhì)標(biāo)準(zhǔn)，其中水質(zhì)指標(biāo)COD應(yīng)小于15 mg/L。在此要求下，本文采用控制斷面法，將下游飲水水源一級(jí)保護(hù)區(qū)的3個(gè)監(jiān)測(cè)斷面作為控制斷面(見(jiàn)圖6)，根據(jù)控制斷面的水質(zhì)要求確定水環(huán)境容量。

表2 連續(xù)源污染排放計(jì)算設(shè)計(jì)工況

圖6 飲用水源保護(hù)區(qū)控制斷面

3.3.2污染物濃度空間分布

在滿足水庫(kù)下游3條控制斷面水質(zhì)指標(biāo)COD濃度均達(dá)到二類水質(zhì)標(biāo)準(zhǔn)(15 mg/L)的條件下，利用數(shù)學(xué)模型計(jì)算了枯水期、平水期、豐水期3種水文條件下的最大允許排放濃度，結(jié)果分別為426、310 mg/L與252 mg/L。在該污染負(fù)荷條件下，COD濃度的空間分布見(jiàn)圖7。從圖7可知，總體而言，污染物濃度隨著水流方向縱向濃度逐漸降低，庫(kù)尾處污染物COD濃度最大，入庫(kù)后污染物隨水流向下游遷移，但隨著污染物運(yùn)動(dòng)過(guò)程中遷移擴(kuò)散及降解作用，沿水流輸運(yùn)方向污染物濃度會(huì)逐漸減小。

圖7 不同水文條件下區(qū)域COD濃度空間分布

豐水期、平水期、枯水期污染帶上下游濃度空間差異分別為230.59、279.65 mg/L和289.48 mg/L。豐水期入庫(kù)流量較大，由于推移作用強(qiáng)，縱向遷移距離較長(zhǎng)；枯水期入庫(kù)流量較小，推動(dòng)作用弱，污染物拉伸范圍較小，濃度空間分布更不均勻，濃度空間差異較大。

3.3.3水環(huán)境容量

(1)計(jì)算區(qū)域。選擇入庫(kù)處(拐上斷面)到壩前共73 km為研究水功能區(qū)域。水環(huán)境容量W采用下式計(jì)算

(7)

式中，Q為設(shè)計(jì)流量；C0為段首污染物濃度；(t2-t1)為時(shí)間間隔。

(2)不同水文條件下水環(huán)境容量分析。設(shè)計(jì)工況下水環(huán)境容量在不同水文條件下的計(jì)算結(jié)果見(jiàn)表3。其中，入庫(kù)COD濃度分別為在枯水期、平水期、豐水期3種水文條件下，連續(xù)污染源COD的最大臨界排放濃度。從表3可知，在枯水期、平水期、豐水期設(shè)計(jì)流量下，下游3個(gè)控制斷面COD均能達(dá)到二類水質(zhì)標(biāo)準(zhǔn)(15 mg/L)時(shí)，入庫(kù)污染物COD日均最大負(fù)荷總量分別為14 060、15 959 t和18 137 t。在等長(zhǎng)河段內(nèi)，水質(zhì)目標(biāo)相同的前提下，不同水文條件下COD的水環(huán)境容量差異明顯，豐水期>平水期>枯水期。產(chǎn)生此差異的原因是污染物的遷移受流速的影響，流量越大，流速越大，污染物的對(duì)流擴(kuò)散及遷移能力越強(qiáng)；枯水期流量小，水流流速較小，污染物遷移速度較為緩慢，容易聚集在庫(kù)區(qū)內(nèi)；豐水期水流流速較大，污染物擴(kuò)散迅速，水力停留時(shí)間短，在庫(kù)區(qū)聚集的可能性小，所以豐水期污染物COD日均負(fù)荷總量最大。該水環(huán)境容量?jī)H根據(jù)下游斷面的水質(zhì)要求計(jì)算得到，在水庫(kù)管理中還需根據(jù)整個(gè)庫(kù)區(qū)的水質(zhì)要求確定可容納的污染負(fù)荷。

表3 庫(kù)區(qū)水環(huán)境容量

4 結(jié) 語(yǔ)

本文采用水環(huán)境數(shù)學(xué)模型，對(duì)萬(wàn)家寨水庫(kù)庫(kù)區(qū)水動(dòng)力特性和水環(huán)境特性進(jìn)行了分析，計(jì)算出典型工況下區(qū)域內(nèi)污染物濃度空間分布和水環(huán)境容量。在固定源連續(xù)排放下，本文設(shè)定典型的水文條件，利用數(shù)學(xué)模型計(jì)算，得到連續(xù)污染源的COD最大允許排放濃度與相應(yīng)的COD負(fù)荷水環(huán)境容量，為庫(kù)區(qū)污染物總量消減提供參考依據(jù)。