萬家寨水庫水環境特性分析

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

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

0 引 言

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

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

目前，對于研究區域水環境污染的研究成果多為水質總體評價，對于庫區水環境特性精細的時空分布研究相對匱乏。為此，本文應用水環境數學模型，對萬家寨水庫庫區水環境進行模擬，分析水動力特性和水環境特性，并計算典型工況下區域內COD最大排放濃度分布及水環境容量。

1 數學模型

1.1 控制方程

采用二維水流水質模型進行求解[8]，即

(1)

(2)

(3)

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

采用對流擴散方程模擬區域污染物輸移擴散規律，即

(4)

式中，C為污染物垂向平均濃度；AH為水平擴散系數；S為源項。

應用有限體積法(FVM)對于上述物質連續方程、動量方程以及濃度方程進行離散，主要包括對對流及擴散通量做體積積分和時間積分處理，根據對流通量、擴散通量、源(匯)項的表達式即可建立各單元的物質守恒方程，將所有單元的濃度守恒方程進行聯立求解，進而得知全場污染物的濃度分布特性。

1.2 網格劃分及計算條件

計算區域為39°34′4″N～40°13′05″N，111°10′23″E～111°27′06″E，上游邊界取在拐上斷面，下游邊界取在壩上斷面。研究區域網格劃分見圖1。采用二維水動力模型進行計算，庫區平面應用非結構三角形網格進行貼體劃分，其劃分網格10 476個，三角形網格節點數為6 109個，三角形界面數16 584個。沿河寬方向網格間距為33.0～66.5 m，平均間距為51.1 m；沿河道方向網格間距為38.7～62.2 m，平均間距為53.6 m。

圖1 萬家寨庫區網格劃分

1.3 計算條件

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

2 模型驗證

從水動力和水質2個方面進行驗證。采用水庫庫容日變化過程對水動力結果進行驗證，見圖2。從圖2可知，實測值和計算值平均相對誤差為1.08%，具有較好的一致性。

圖2 模擬庫容與實測庫容驗證

選取COD與氨氮濃度進行不同水文年條件下水質模型壩前污染物濃度驗證，見圖3、4。利用平均對數比偏差Dr計算計算值與實測值的差異性，公式如下

(5)

圖3 不同水文年條件下水質模型COD濃度驗證

圖4 不同水文年條件下水質模型氨氮濃度驗證

式中，DLm為實測值；DLp為計算值；n為樣本個數。

經計算，各種條件下2種水質指標計算值與實測值之間的平均對數比偏差均在0.1～0.3范圍內，最大為0.297，滿足工程應用精度要求。

圖5 不同水文年下夏季壩前30 km范圍內流速分布狀態

3 結果及分析

3.1 水動力場分析

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

3.2 水力停留時間

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

(6)

式中，HRT為水力停留時間；V為水庫庫容；Q為入流流量。

庫區水力停留時間計算統計見表1。從表1可知，在典型豐水年、平水年、枯水年3種水文年條件下，水庫段水力停留時間在1.76～11.78 d范圍內，萬家寨水庫具備典型的河道水庫型水動力特征。具體而言，夏季水力停留時間小于冬季水力停留時間，最短水力停留時間出現在豐水年夏季，為1.76 d；最長水力停留時間出現在枯水年冬季，為11.78 d。產生此結果的原因是水庫在豐水年夏季具有較常年各季節偏高的徑流量，水力停留時間較其他時期短。枯水年冬季具有較常年各季節偏低的徑流量，水力停留時間較其他時期長。

表1 庫區水力停留時間計算統計

3.3 COD的最大排放濃度及水環境容量計算

3.3.1計算工況

為了探究水庫水環境容量，假設萬家寨水庫上游入流斷面處以固定流量連續排放一定濃度的有機污水，設計3種不同的計算工況，見表2。

根據庫區HJ 338—2018《飲用水水源保護區劃分技術規范》，飲用水源一級保護區的水質應滿足二類水質標準，其中水質指標COD應小于15 mg/L。在此要求下，本文采用控制斷面法，將下游飲水水源一級保護區的3個監測斷面作為控制斷面(見圖6)，根據控制斷面的水質要求確定水環境容量。

表2 連續源污染排放計算設計工況

圖6 飲用水源保護區控制斷面

3.3.2污染物濃度空間分布

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

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

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

3.3.3水環境容量

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

(7)

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

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

表3 庫區水環境容量

4 結 語

本文采用水環境數學模型，對萬家寨水庫庫區水動力特性和水環境特性進行了分析，計算出典型工況下區域內污染物濃度空間分布和水環境容量。在固定源連續排放下，本文設定典型的水文條件，利用數學模型計算，得到連續污染源的COD最大允許排放濃度與相應的COD負荷水環境容量，為庫區污染物總量消減提供參考依據。