胡 艷 海，周 林 飛

(沈陽農(nóng)業(yè)大學 水利學院，遼寧 沈陽 110161)

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

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

1 研究區(qū)域概況及水質(zhì)監(jiān)測

1.1 研究區(qū)域概況

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

1.2 水質(zhì)監(jiān)測

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

監(jiān)測時間為2009～2017年。石佛寺水庫作為水源地，水標準執(zhí)行地表水Ⅲ類。

2 基本原理

2.1 水動力數(shù)值模擬基本原理

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

(1)

(2)

(3)

圖1 石佛寺水庫位置及水質(zhì)監(jiān)測點位置Fig.1 Location and water quality monitoring points of Shifosi Reservoir

(4)

(5)

(6)

2.2 水質(zhì)模型基本原理

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

(7)

S=Qs(Cs-C)

(8)

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

3 石佛寺水庫水動力及水質(zhì)數(shù)值模擬

3.1 水動力數(shù)值模擬

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

3.1.1模擬范圍

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

圖2 石佛寺水庫計算網(wǎng)格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。由于模型在短時間內(nèi)計算趨于穩(wěn)定，故在初始時刻t=0時，庫區(qū)初始流速u和v設置為0，初始水位為δ=46.200 m。其他數(shù)據(jù)為實測值。

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

3.1.3模型的率定與驗證

(1) 模型率定。選擇2017年6，7月石佛寺水庫入口處實測水位率定模型參數(shù)，率定結(jié)果為：曼寧系數(shù)45 m1/3/s(MIKE軟件的曼寧系數(shù)與糙率值互為倒數(shù))、渦黏系數(shù)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 水質(zhì)數(shù)值模擬

3.2.1水質(zhì)模板選擇與水質(zhì)指標選取

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

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

3.2.2定解條件與模型參數(shù)

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

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

(3) 模型參數(shù)。建立生態(tài)模型時污染物的運移消散過程必然會發(fā)生，因此水質(zhì)的擴散系數(shù)以及衰減系數(shù)的設置顯得至關(guān)重要[21]。其中，擴散系數(shù)在各個方向上的上下限通常根據(jù)其與各網(wǎng)格速度分量之間的線性比例關(guān)系而確定；降解系數(shù)則通常根據(jù)原型觀測法[22]獲得，然后利用環(huán)境部門檢測結(jié)果和相關(guān)研究成果進行修正和再次修正，最后確定各水質(zhì)指標的降解系數(shù)如下：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水質(zhì)模型率定與驗證

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

(2) 模型驗證?；谀Ｐ蛥?shù)的成功率定，將對模型進行驗證，模型驗證時選用2017年5，7月出口處實測生化需氧量、化學需氧量、氨氮、總氮、總磷數(shù)據(jù)驗證模型，驗證結(jié)果如圖6所示。采用統(tǒng)計學中的偏差統(tǒng)計法，得到五日生化需氧量、化學需氧量、氨氮、溶解氧、總磷、總氮實測數(shù)據(jù)與計算模擬數(shù)據(jù)之間平均誤差分別為4.04%，7.39%，5.10%，3.80%，3.05%，1.40%，均在可接受范圍內(nèi)，表明實測值與模擬計算值擬合程度較好，完全能夠滿足精度要求，所以本文所構(gòu)建的水質(zhì)模型是有效可行的。

4 模擬結(jié)果與分析

4.1 水動力模擬結(jié)果與分析

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

4.1.1模擬方案

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

4.1.2模擬結(jié)果分析

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

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

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

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

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

4.2 水質(zhì)模擬結(jié)果與分析

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

4.2.1模擬結(jié)果

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

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

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

4.2.2分 析

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

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

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

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

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

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

5 結(jié) 論

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

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