石佛寺水庫二維水動力及水質數值模擬與分析

胡 艷 海，周 林 飛

(沈陽農業大學 水利學院，遼寧 沈陽 110161)

長期以來，水庫水質狀況和庫區水生態環境受到廣大學者的關注，其中以水庫水動力條件為基礎進行水庫污染物質運移消散過程模擬是主要的研究方式。國外關于水質模型的研究起步于20世紀60～80年代，計算機技術的迅速發展促進了雙向耦合模型對水質模型變化所產生的影響，與此同時水質模型也伴隨著水動力系數、空間變量、主要的外界影響因素(例如溫度等)的引進而變得更加完善[1]。20世紀80年代后，多維數值模型求解伴隨二維水質模型的成熟得以實現。國內關于水質模型的研究比較晚，20世紀80年代，由于官廳水庫遭受嚴重的重金屬污染事件，由此引起了多個領域的研究學者對水庫水質問題的關注，并意識到利用水質數值模型解決水環境相關問題的重要性[2]。因此，眾多學者針對不同類型的湖泊、湖泊中各種水質問題、富營養化問題以及不同的水質因素構建水質模型，最終實現水質模型可行性驗證，同時為進一步實現模型的簡化、促進數值模型的應用進展創造條件[3-5]。龔春生等[6]通過模擬玄武湖混合流的水質動態變化過程為底泥污染和淺水湖泊混合流水質變化模擬提供新途徑。20世紀90年代，水質模型逐漸被廣泛應用到水質優化方案設計方面[7]。21世紀后，水質模型逐漸向多維、多層次方向發展[8-11]。例如，楊具瑞[12]構建的垂直分層的二維水質模型，計勇[8]構建的水動力-水質-底泥模型，廖臨毓[13]構建MIKE21水動力-水質耦合模型對COD、TN、TP三種污染物在3種調水方案下的流場和濃度場時空變化情況進行模擬，馬寧[14]構建水動力-水質耦合模型模擬水庫水動力和水質狀況并獲得水質改善方式。

石佛寺水庫位于遼河干流，承擔防洪、供水、生態等重任。因此，掌握庫區水動力條件以及水質狀況顯得尤為重要。本次研究以石佛寺水庫作為研究對象，通過總結前人經驗構建水庫的二維水動力-水質耦合模型，將石佛寺水庫來水量作為水庫水動力模型的主要影響因素，以水質模型濃度場變化和水動力模型流場變化作為運算基礎，模擬運算了庫區汛期和非汛期兩種工況下溶解氧(DO)、氨氮(NH3-N)、總氮(TN)、五日生化需氧量(BOD5)、化學需氧量(CODCr)、總磷(TN)6種水質指標的運移消散狀況。最后，根據分析模擬結果提出水環境綜合治理的具體措施，用以指導水庫運行，使遼河水流經水庫后能夠達到地表水Ⅲ類水質標準。

1 研究區域概況及水質監測

1.1 研究區域概況

石佛寺水庫位于遼寧省沈陽市沈北新區黃家鄉和法庫縣依牛堡鄉，地理坐標為東經123.427°～123.520°，北緯42.144°～42.188°，距沈陽市區約47 km，如圖1所示。石佛寺水庫具有防洪、供水、生態三大功能，其中最主要的是防洪，為典型的河道型平原水庫，是遼河干流上僅有的控制性水利工程。為改善區域生態環境，凈化遼河上游來水，實現其生態功能，2009年起在石佛寺水庫實施了生態建設工程，并在庫區內遼河左岸種植水生植物蘆葦、蒲草和荷花，同時進行生態蓄水，蓄水面積16.13 km2，蓄水位46.2 m。

1.2 水質監測

為了了解遼河水流經石佛寺水庫后水質凈化狀況，在水庫的入口、中間和出口布設水質監測點位，分別標記為采樣點1、采樣點2、采樣點3，如圖1所示。各監測點監測的水質指標為溶解氧(DO)、氨氮(NH3-N)、總氮(TN)、五日生化需氧量(BOD5)、化學需氧量(CODCr)、總磷(TP)等11個指標。DO利用多參數光譜水質分析儀測定，NH3-N采用水楊酸鈉法測定，TN采用全自動間斷化學分析儀測定，BOD5采用稀釋與接種法測定，CODCr采用重鉻酸鹽法測定，TP采用硫酸鉀—鉬銻抗法測定。監測頻率為每月2次，監測標準以GB 3838-2002《地表水環境質量標準》為依據。

監測時間為2009～2017年。石佛寺水庫作為水源地，水標準執行地表水Ⅲ類。

2 基本原理

2.1 水動力數值模擬基本原理

MIKE21水動力模塊是在二維數值求解淺水方程的基礎上建立的[15-18]。二維水動力數值模擬分析基于以下4項基本假定展開計算：不可壓縮性假定、靜水壓力假定、Boussinesq假定、Reynolds值均布假定。水動力模塊控制方程應用沿水深積分的不可壓縮的Navier-Stokes方程，其主要包括二維平面流的連續性方程、動量方程[19]，如式(1)～(3) 所示。本文選用單元中心的有限體積法作為數值解法，并采用MIKE21FM的非結構化三角形網格進行模型求解。

(1)

(2)

(3)

圖1 石佛寺水庫位置及水質監測點位置Fig.1 Location and water quality monitoring points of Shifosi Reservoir

(4)

(5)

(6)

2.2 水質模型基本原理

庫湖水質控制方程是全面考慮污染物質水動力學、物理運輸和對流擴散等過程并將質量守恒原理作為前提建立起來的。該控制方程的數學表達式如式(7)～(8) 所示[20]。

(7)

S=Qs(Cs-C)

(8)

二維水動力模型作為水質模型的基礎，只有建立水動力模型并滿足模型參數率定和模型驗證，方能加載水質模型，其主要原理在于水動力模型能夠提供水質模型水質濃度輸出的流速等水動力條件。

3 石佛寺水庫水動力及水質數值模擬

3.1 水動力數值模擬

水庫水動力影響因素包括：進出水條件、風場、地形、庫內水工構筑物以及水庫的岸線形態等。其中，水庫的地形和岸線形態作為相對固定的因素，來水量為主要因素。

3.1.1模擬范圍

本文選用1∶5 000石佛寺水庫電子版地形圖作為地形資料，利用Auto CAD和南方CASS提取高程點并生成.xyz格式陸地land和水深water數據文件用于網格生成。采用Beijing_1954_GK_Zone_21N投影坐標，控制最大三角形面積20 000 m2、最小三角形角度30°以及最大節點個數150 000個，利用MIKE ZERO網格生成器生成網格，計算網格如圖2所示。完成水深water文件插值并利用網格分析器進行網格分析，獲得合理準確網格文件并由MIKE Animator生成模擬區域的三維地形，如圖3所示。

圖2 石佛寺水庫計算網格Fig.2 Computational grid of Shifosi Reservoir

圖3 石佛寺水庫三維地形Fig.3 Three-dimensional topography of Shifosi Reservoir

3.1.2定解條件

(1) 初始條件。模擬時長為61 d，即2017年6月1日08：00：00至2017年7月31日08：00：00，時間步長為3 600 s。由于模型在短時間內計算趨于穩定，故在初始時刻t=0時，庫區初始流速u和v設置為0，初始水位為δ=46.200 m。其他數據為實測值。

(2) 邊界條件。根據MIKE軟件前處理方法，上游邊界為石佛寺水庫入口處實測流量數據，下游邊界為石佛寺水庫出口處實測水位數據；水質模型利用岸壁法確定閉邊界，即沿法線方向上的濃度值和流速值均為零；主要根據底床摩擦應力確定底床邊界條件；采用可避免模型發散的“干濕水深判別法”確定移動邊界，即水深大于濕水深0.1 m時按水域處理，水深小于干水深0.005 m時按陸域處理。

3.1.3模型的率定與驗證

(1) 模型率定。選擇2017年6，7月石佛寺水庫入口處實測水位率定模型參數，率定結果為：曼寧系數45 m1/3/s(MIKE軟件的曼寧系數與糙率值互為倒數)、渦黏系數0.28。

圖4 出口合成流速模擬值與實測值對比Fig.4 Comparison of the simulated value and measured value of the exit flow velocity

圖5 入口水位模擬值與實測值對比Fig.5 Comparison of the simulated value and measured value of the entrance water level

3.2 水質數值模擬

3.2.1水質模板選擇與水質指標選取

本文選用改進的Eco Lab內置水質模板WQ with nutrients對石佛寺水庫水質進行模擬。該水質模板包含狀態變量7個，常量39個，作用力6個，輔助變量9個，過程21個。

綜合考慮實際水質監測數據的超標情況和生態數值實驗室(ECO Lab)內置模板中狀態變量的可模擬情況確定水質指標為：DO，NH3-N，BOD5，TN，CODCr，TP，其中TN和TP作為衍生結果輸出。

3.2.2定解條件與模型參數

(1) 初始條件。模擬范圍、模擬網格系統及模擬時間同水動力模型?？紤]模型整體一致，故在初始時刻t=0時，水質指標的初始濃度值設置為0。

(2) 邊界條件。閉邊界法線方向上設定流速和濃度值均為零；開邊界，將水庫進出口處水質指標實測濃度生成的時間序列文件作為上下游邊界條件。

(3) 模型參數。建立生態模型時污染物的運移消散過程必然會發生，因此水質的擴散系數以及衰減系數的設置顯得至關重要[21]。其中，擴散系數在各個方向上的上下限通常根據其與各網格速度分量之間的線性比例關系而確定；降解系數則通常根據原型觀測法[22]獲得，然后利用環境部門檢測結果和相關研究成果進行修正和再次修正，最后確定各水質指標的降解系數如下：DO為0.075 14 mg/d；NH3-N為0.113 1 mg/d；BOD5為0.140 53 mg/d；TN為0.014 95 mg/d；CODCr為0.102 05 mg/d；TP為0.042 12 mg/d。

3.2.3水質模型率定與驗證

(1) 模型率定。選擇2017年6～7月石佛寺水庫主槽區實測五日生化需氧量、化學需氧量、氨氮、總氮、總磷率定模型參數，率定結果為：橫向擴散系數5 m2/s、縱向擴散系數5 m2/s，化學需氧量衰減系數0.078 5 mg/d，五日生化需氧量衰減系數0.108 1 mg/d，氨氮衰減系數0.087 mg/d，溶解氧衰減系數0.057 8 mg/d，總氮衰減系數0.011 5 mg/d，總磷衰減系數0.032 4 mg/d。

(2) 模型驗證?；谀Ｐ蛥档某晒β识ǎ瑢δＰ瓦M行驗證，模型驗證時選用2017年5，7月出口處實測生化需氧量、化學需氧量、氨氮、總氮、總磷數據驗證模型，驗證結果如圖6所示。采用統計學中的偏差統計法，得到五日生化需氧量、化學需氧量、氨氮、溶解氧、總磷、總氮實測數據與計算模擬數據之間平均誤差分別為4.04%，7.39%，5.10%，3.80%，3.05%，1.40%，均在可接受范圍內，表明實測值與模擬計算值擬合程度較好，完全能夠滿足精度要求，所以本文所構建的水質模型是有效可行的。

4 模擬結果與分析

4.1 水動力模擬結果與分析

利用驗證合理的模型進行模擬，將1號、2號和3號監測點作為特征點，用于反映石佛寺水庫流速和流場特點，監測點布置如圖1所示。

4.1.1模擬方案

本次水動力數值模擬共分為兩種工況：即汛期和非汛期。根據石佛寺水庫多年運行經驗，6～9月為汛期，非汛期不考慮冰凍期，故非汛期僅模擬為10～11月和4～5月。汛期和非汛期初始流速均設置為0，初始水位均設置為46.200 m。汛期和非汛期上游邊界為石佛寺水庫入口處實測流量數據，下游邊界為石佛寺水庫出口處實測水位數據。模擬時長非汛期為2016年10月1日～11月30日、2017年4月1日～5月31日，共122 d；汛期為2017年6月1日～9月30日，共122 d?，F選取汛期和非汛期水動力模型運行末期時典型流場分布情況進行分析。

4.1.2模擬結果分析

(1) 根據圖7可知，汛期水位大部分區域介于46.280～46.256 m之間，非汛期水位大部分區域介于46.190～46.185 m之間，即汛期庫區水位大部分僅略高于非汛期。原因是雖然汛期水量較大，但會通過下游閘門控制水庫水位，使其保持在46.200 m左右，保障生態需水。

圖6 各水質指標模擬值與實測值對比Fig.6 Comparison of the simulated value and measured value of water quality indicators

(2) 圖8為汛期和非汛期不同水量條件下庫區典型流場流速分布情況，結果表明兩種工況下水庫具有較好的流動狀態且未出現死水區(流速≤0.3 cm/s時為死水區)和回水區。對比可以看出汛期的流速整體大于非汛期，說明不同的來水量水庫內不同區域內的流速變化明顯。兩種工況下，流速較大區域均出現在上游入口區域和下游出口區域，這是因為過水斷面面積小從而加速了水體運動；過水斷面面積較大的區域，流速比較緩慢，占庫區面積也大；流速較小區域出現在左右岸岸邊，與岸邊地形相關?？傮w來看，庫區內水體流動性良好、水流不急，有利于水體的交換以及污染物的去除。

圖7 水位分布(單位：m)Fig.7 Water level distribution

圖8 流場流速分布(單位：m/s)Fig.8 Flow field velocity distribution

4.2 水質模擬結果與分析

來水量是石佛寺水庫水動力的主要影響因素，因此選擇非汛期和汛期兩種不同情景對其進行水質模擬分析。具體模擬時間段為：非汛期為2016年10月1日至11月30日，2017年4月1日至5月31日，共122 d；汛期2017年6月1日至9月30日，共122 d。初始時刻水質指標濃度設置為0，采用水質指標濃度值生成的時間序列文件作為上下邊界條件。

4.2.1模擬結果

本次數值模擬獲得兩種情景下水庫DO，NH3-N，BOD5，TN，CODCr和TP在不同時間段的濃度分布情況。汛期情景下各污染物在模擬末期時分布圖如圖9所示，非汛期情景下各污染物在模擬末期時分布如圖10所示。

圖9 汛期各水質指標濃度分布(單位： mg/L)Fig.9 Distribution of water quality indicators in flood season

圖10 非汛期各水質指標濃度分布(單位： mg/L)Fig.10 Distribution of water quality indicators in non-flood season

4.2.2分 析

汛期和非汛期除DO以外各污染物濃度整體分布情況為：入口>中間>出口，原因是石佛寺水庫具有凈化污染物質的作用，其中在非汛期盡管水流流速減緩，但是水體中的某些污染物質會伴隨泥沙發生沉降，污染物質被存儲在底泥中，進而實現凈化水質的作用。DO為出口>中間>入口，流經水庫復氧能力增強，水質得到改善。

在汛期NH3-N，TN，CODCr和BOD5濃度小于非汛期的濃度，原因一是汛期庫區內水生植物長勢茂盛，微生物具有很強的活性，有很好的水質凈化作用；原因二是汛期有較大的水量，流速也較大，自然增強了水體的自凈能力。在汛期DO的濃度大于非汛期的濃度，原因是水體復氧途徑由原來的的大氣復氧，增加了水生植物向水體中輸送氧氣和藻類等沉水植物的光合作用產生氧氣，而且汛期流速快增加了水體的復氧能力。TP的濃度汛期和非汛期變化不明顯。

水質從入口到出口逐漸得到凈化，說明遼河水流經水庫后水質得到了改善。但分析出口水質數據，仍有指標在一定時間段內達不到地表水Ⅲ類水質要求，因此，仍需加強庫區內的污染綜合控制，可以從以下3個方面進行考慮：

(1) 控制外源污染。進水污染物濃度為導致水庫污染的主要因素，因此可通過控制面源污染和點源污染，盡量減少經降雨、下滲、徑流過程形成的污染源匯入水庫。在外源污染超過水庫的凈化能力時，可根據石佛寺水庫水動力水質模擬結果進行調水，進一步改善水質。

(2) 減少內源污染。水庫在非汛期水流流速會減緩，泥沙會攜帶一些污染物質沉降而儲存在水庫底泥中，尤其在水流相對緩慢的非汛期，為避免再次污染可以通過對庫底底泥疏浚來控制內源的匯入。

(3) 種植水生植物。汛期的流場和水質指標濃度分布圖疊加后得出，生長水生植物的區域水流受阻流速變緩，水生植物的凈化作用使污染物濃度降低；在無水生植物的區域，水流流速較大的區域水質指標的濃度值相對較小些。因此，應該在水流流速較小的區域種植凈化能力較強且阻水較弱的水生植物，在高流速區域種植阻水能力較強的水生植物。

5 結 論

(1) 石佛寺水庫水動力數值模擬結果表明：來水量是水庫水流狀態的主要影響因素。汛期和非汛期兩種工況下，汛期的流速整體大于非汛期，說明不同的來水量水庫內不同區域內的流速變化明顯；兩種工況下水庫均具有較好的流動狀態，未出現死水區和回水區，且水流不急，有利于水體的交換以及污染物的去除。

(2) 石佛寺水庫水質數值模擬結果表明：遼河水流經水庫后逐級得到了凈化，而且凈化效果汛期優于非汛期。但在一定時間段內仍有指標達不到地表水Ⅲ類水質要求，需加強庫區內的污染綜合控制，具體措施如下：控制外源污染和調水凈化；通過底泥疏浚減少內源污染；種植具有凈化功能的水生植物，在水流流速較小的區域種植阻水較弱的水生植物，在高流速區域種植阻水能力較強的水生植物。