，，,3

(1.河海大學(xué) 環(huán)境學(xué)院，南京 210000； 2.揚州市職業(yè)大學(xué) 資源與環(huán)境工程學(xué)院，江蘇 揚州 225000； 3.黃河水利科學(xué)研究院，鄭州 450003)

近年來在水庫周邊工業(yè)點源污染獲得有效治理后，如何控制周邊農(nóng)業(yè)面源污染對庫區(qū)水質(zhì)的影響已成為水環(huán)境治理研究重點。童曉霞等[1]分析了降雨量及降雨波動對農(nóng)田中氮磷流失的影響。陳會等[2]開展了主要面源污染在排水過程中的遷移轉(zhuǎn)化過程分析，確定了對面源污染形成起主要作用的排水形式。常蒲婷等[3]應(yīng)用確定性水質(zhì)模型及面源污染負(fù)荷模型對面源污染進(jìn)行風(fēng)險和污染負(fù)荷計算，模擬預(yù)測了突發(fā)面源污染對水質(zhì)的影響。張廣納等[4]在山峽庫區(qū)以區(qū)縣為研究單位，分析了面源污染對水質(zhì)影響的變化規(guī)律?，F(xiàn)有研究更偏重于農(nóng)業(yè)面源污染形成原因分析及對污染負(fù)荷的估算，而就降雨過程中農(nóng)業(yè)面源污染入庫與庫區(qū)水質(zhì)變化這一動態(tài)響應(yīng)過程的分析則明顯不足。開展這一過程的模擬分析可以為庫區(qū)面源污染的預(yù)報及治理提供技術(shù)依據(jù)。本次研究考慮入庫流量與風(fēng)應(yīng)力的綜合影響，構(gòu)建研究區(qū)三維水動力水質(zhì)模型，分析強降雨過程中特征污染物入庫后的濃度變化特征，所得結(jié)論可用于指導(dǎo)面源污染治理。

1 研究區(qū)概況

大伙房水庫是遼寧省重要水源地，為沈陽、盤錦等7市2 300萬人口提供水源。水庫東西長約35 km、水面最寬處達(dá)4 km、最窄處約0.3 km，水庫最大水深37 m，最大庫容量為21.87億m3，呈帶狀河谷型水庫特征。庫區(qū)降雨主要集中在夏季7—8月份，流域內(nèi)多年平均風(fēng)速為1.5～3.8 m/s。社河是重要入庫水源，全長53 km，控制流域面積516 km2，流域總?cè)丝?7 364人，流域現(xiàn)狀調(diào)查表明地表水污染源主要包括農(nóng)田地表徑流、畜禽養(yǎng)殖、農(nóng)村居民生活污染物等[5]。一旦發(fā)生強降雨，污染物將隨降雨徑流進(jìn)入河道最后匯入水庫，影響庫區(qū)水質(zhì)。如圖1所示，南彰黨水文站離社河入庫口距離為8 687 m，如不考慮污染物在這段距離上的降解稀釋，可將南彰黨水文站水文、水質(zhì)監(jiān)測值作為社河入庫口處的值。IP1-IP2為橫向垂線監(jiān)測斷面；IP3-IP4為縱向垂線監(jiān)測斷面。

圖1 大伙房水庫位置Fig.1 Location of Dahuofang Reservoir

2 典型污染因子及模擬時段選取

圖2 2012年8月4—8日 入庫流量變化曲線Fig.2 Curve of inflow of reservoir from August 4 to August 8, 2012

根據(jù)文獻(xiàn)的研究結(jié)論[6]，總磷是影響大伙房水庫水體富營養(yǎng)化的重要水化學(xué)因子，且大伙房水庫總磷濃度在近年內(nèi)變化較大，因此選取總磷作為典型污染因子。水庫總磷主要來源于上游入庫河流，尤其當(dāng)夏季豐水期，入庫流量特別大時，總磷污染會明顯加重，因此選擇豐水期作為研究時段。在2012年6月10日—8月10日這2個月的監(jiān)測期內(nèi)，8月4日有一次降水接近大暴雨過程，如圖2所示，社河流量于8月4日12時增加至最高的478 m3/s，TP通量也在8月4日16時達(dá)到峰值。本次研究選擇該次降雨過程開展分析。

3 三維水動力、水質(zhì)數(shù)學(xué)模型

大伙房水庫中流速、污染物濃度沿垂向分布差異較大，因此，采用三維湖泊水動力水質(zhì)模型分析水動力參數(shù)、污染物濃度的空間分布特征。模型采用σ坐標(biāo)變換的相對分層法進(jìn)行求解。

3.1 σ坐標(biāo)下三維水動力數(shù)學(xué)模型

(1)

式中：H為總水深(m)；u,v分別為坐標(biāo)x,y方向上的流速分量(m/s)；w*指通過某一σ坐標(biāo)層的流速(m/s)；?x*,?y*均表示σ坐標(biāo)系下的坐標(biāo)方向;S為源項排放量(kg/(m3·s))；f為柯氏參數(shù)；ρ0為水流參考密度(kg/m3)；ρ為液體密度(kg/m3)；patm為壓強(Pa)；g為重力加速度(m/s2)；vz為垂向渦黏系數(shù)；us，vs為源項排放速度在x,y方向上的流速分量(m/s);HFX,HFY為擴散項。

3.2 σ坐標(biāo)下三維水質(zhì)數(shù)學(xué)模型

(2)

其中,

(3)

式中：C為污染物濃度(g/L)；kp為污染物線性衰減系數(shù)；Dv為垂向擴散系數(shù)；Cs為排放源項污染物濃度(g/L)；Dh為水平擴散系數(shù)。

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

4.1 強降雨過程中庫區(qū)水動力特征模擬分析

4.1.1 研究區(qū)網(wǎng)格劃分、定解條件及模擬參數(shù)取值

4.1.1.1 研究區(qū)網(wǎng)格布置

對研究區(qū)域開展水動力模擬，水平方向劃分為無結(jié)構(gòu)三角網(wǎng)格，網(wǎng)格邊長250 m,對社河入庫口位置進(jìn)行加密處理，網(wǎng)格邊長100 m，整個平面共分為3 652個網(wǎng)格單元。垂向上采用相對分層法剖分為10層。

4.1.1.2 定解條件及模擬參數(shù)取值

采用Kolmogorov-Prandl方程求解垂向渦黏系數(shù)，采用Smagorinsky亞網(wǎng)格尺度紊動模型計算水平渦黏系數(shù)。入流邊界根據(jù)豐水期實測平均入庫流量確定，蘇子河邊界值為138.7 m3/s，渾河邊界值為100.6 m3/s。社河入庫流量取流量的時間序列值。定解條件：初始流速均取為0；初始水位和出流邊界水位為129.85 m；依據(jù)不可入邊界條件，岸邊界法向流速為0 m/s。模擬參數(shù)的取值：水平與垂向擴散系數(shù)為1 m2/s；總磷(TP)降解系數(shù)為0.001 8 d-1；粗糙高度為0.01 m。

桌上的電話驟然響了。楊秉奎抓起聽筒:“對,是我,‘養(yǎng)病虧’站長……放心,我知道……哎,你說話客氣點嘛……我不管你是誰,給老子記著!”

4.1.2 研究區(qū)表層流場隨時間演變分析

以主導(dǎo)風(fēng)向(東北風(fēng))和平均風(fēng)速2.3 m/s作為自由表面條件，模擬分析豐水期水庫三維水動力特征。當(dāng)水動力模型計算穩(wěn)定時，取研究區(qū)表層流場進(jìn)行分析，如圖3所示。

圖3 水流入庫后研究區(qū)表層流場Fig.3 Surface flow field in the study area after water flows into the reservoir

在社河入庫后12 h，流量為478 m3/s，達(dá)到最大值。此時吞吐流對水庫流場的影響遠(yuǎn)大于風(fēng)生流，水流呈現(xiàn)受吞吐流支配的水動力特征。同時受地形條件影響，水流呈現(xiàn)由入庫口向壩址方向顯著的順岸流特征，流速相對較大。在社河入庫后50,80 h，流量分別下降為86,31 m3/s，此時吞吐流對水流場的影響逐步減弱，水流結(jié)構(gòu)呈現(xiàn)由風(fēng)生流和吞吐流共同支配逐步轉(zhuǎn)變?yōu)橛娠L(fēng)生流支配的變化特征。與此同時，受上游渾河、蘇子河入庫水流影響，庫區(qū)中的部分水流向壩址流動，部分水流向社河入庫口流動，與社河的入庫河流形成順時針大環(huán)流。

4.1.3 研究區(qū)三維流場隨時間演變分析

當(dāng)水動力模型計算穩(wěn)定時，選擇10個垂向水層中的表層、中層(第6層)及底層開展三維水流場特征分析。選擇3個模擬時刻研究區(qū)分層流場特征，見圖4—圖6。

圖4 水流入庫后12 h研究區(qū)分層流場Fig.4 Stratified flow field in the study area 12 hours after water flows into the reservoir

圖5 水流入庫后50 h研究區(qū)分層流場Fig.5 Stratified flow field in the study area 50 hours after water flows into the reservoir

圖6 水流入庫后80 h研究區(qū)分層流場Fig.6 Stratified flow field in the study area 80 hours after water flows into the reservoir

由圖4—圖6可知：由于水流場受風(fēng)應(yīng)力、社河入庫水流和上游渾河、蘇子河入庫水流的多重影響，而主導(dǎo)風(fēng)向、上游水流方向與社河入庫水流方向不同，表層水流受此影響最為明顯，因此表層水流速較慢，底層水流速較快。

4.2 強降雨過程中庫區(qū)水質(zhì)變化模擬分析

4.2.1 水平濃度場分布特征分析

水平方向上，取表層的TP分布狀態(tài)進(jìn)行特征分析。圖7為3個時刻研究區(qū)表層TP濃度場分布圖。

圖7 水流入庫后研究區(qū)表層TP濃度分布Fig.7 Distribution of TP concentration in surface water after water flows into the reservoir

由圖7可知，TP入庫12 h時，社河入庫流量較大，吞吐流對水庫流場的影響較大，污染團(tuán)在水流輸移作用下，沿水流方向往下游擴散，最大濃度范圍(濃度>0.4 mg/L)逐漸增大，隨后由于TP的擴散和降解作用，最大濃度范圍又逐步減小。TP入庫50，80 h時，社河入庫流量逐步減少，污染團(tuán)向下游遷移擴散時受主導(dǎo)風(fēng)向東北風(fēng)影響，一部分TP污染團(tuán)流向西南岸，向壩址方向擴散。隨著時間的推移，污染團(tuán)超標(biāo)面積(濃度>0.025 mg/L)不斷增大，即污染團(tuán)對受納水體的影響范圍不斷擴大。

圖8 水流入庫后特征橫向立面TP濃度場Fig.8 TP concentration field in specific horizontal elevation after water flows into the reservoir

4.2.2 面源污染空間分布特征分析

4.2.2.1 橫向立面分布特征分析

橫向立面方向上，取圖1中IP1-IP2線所示橫向垂線斷面濃度分布特征進(jìn)行分析。由于TP在社河中垂向擴散已經(jīng)達(dá)到穩(wěn)定，因此在入庫后，TP在垂向上的擴散也達(dá)到穩(wěn)定。圖8為TP自社河入庫12,50,80 h時庫區(qū)橫向立面TP濃度場。

由圖8可知，TP入庫12 h時，污染團(tuán)剛到達(dá)垂線斷面，橫向垂線斷面最大濃度僅為1.510-8mg/L,遠(yuǎn)小于TP標(biāo)準(zhǔn)限值0.025 mg/L，對垂線斷面區(qū)域的水質(zhì)影響不大；TP入庫50 h時，橫向垂線斷面最大濃度為0.4 mg/L，大于TP標(biāo)準(zhǔn)限值，此時污染團(tuán)對垂線斷面區(qū)域的水質(zhì)影響較大，并且此時水流場受風(fēng)生流影響較大，污染團(tuán)向下游擴散的同時逐漸流向西南岸，因此垂線斷面西側(cè)濃度逐漸大于東側(cè)濃度；TP入庫80 h時，橫向垂線斷面最大濃度為0.2 mg/L，雖有所下降，但仍大于TP標(biāo)準(zhǔn)限值，說明污染團(tuán)在此較長時間內(nèi)對橫向垂線斷面區(qū)域的水質(zhì)有一定程度的影響，同樣受風(fēng)生流的影響，橫向垂線斷面濃度表現(xiàn)出西側(cè)水域濃度比東側(cè)水域濃度大的特征。

4.2.2.2 縱向立面分布特征分析

縱向立面方向上，取如圖1中IP3-IP4所示縱向垂線斷面濃度分布特征進(jìn)行分析。圖9為TP入庫12,50,80 h時庫區(qū)縱向立面TP濃度場。

由圖9可知，TP入庫12 h時，縱向垂線斷面最大濃度為0.4 mg/L，大于TP標(biāo)準(zhǔn)限值0.025 mg/L，污染團(tuán)入庫時間不久，僅對社河入庫口附近水域水質(zhì)影響較大，此時TP污染團(tuán)中心處到入庫口的距離為0.52 km；TP入庫50 h時，縱向垂線斷面最大濃度為0.4 mg/L，大于TP標(biāo)準(zhǔn)限值，污染團(tuán)受東北風(fēng)向的影響，逐漸流向西南岸，因此TP濃度出現(xiàn)中間一段濃度較低而此段兩邊水域濃度較高的情況，但是靠近入庫口水域的濃度更高，此時TP污染團(tuán)中心處到社河入庫口的距離為2.40 km；TP入庫80 h時，縱向垂線斷面最大濃度為0.3 mg/L，雖然仍有部分區(qū)域TP濃度大于標(biāo)準(zhǔn)限值，但超標(biāo)范圍只有很小一部分，說明污染團(tuán)沿著水流方向，不斷向下游遷移，TP污染物濃度沿流程降低，對區(qū)域水質(zhì)影響較小。

5 結(jié) 論

針對面源污染入庫與庫區(qū)水質(zhì)變化的動態(tài)響應(yīng)過程開展過程分析?；跇?gòu)建的研究區(qū)三維水動力及水質(zhì)模型，分析了強降雨過程中庫區(qū)水動力及污染團(tuán)濃度變化與吞吐流、風(fēng)生流影響之間的關(guān)系。呈現(xiàn)了污染團(tuán)濃度在水平向、橫向和縱向上隨時間的動態(tài)變化過程。

分析結(jié)論可以實時預(yù)報污染團(tuán)位置的遷移過程,這為制定更為切實可行的面源污染治理手段提供了有效的技術(shù)參考。在后續(xù)研究中，可以關(guān)注數(shù)值模擬技術(shù)與地理信息系統(tǒng)技術(shù)的結(jié)合實現(xiàn)數(shù)值模擬的可視化。

參考文獻(xiàn)：

[1] 童曉霞,崔遠(yuǎn)來,史 偉. 降雨對灌區(qū)農(nóng)業(yè)面源污染影響規(guī)律的分布式模擬[J].中國農(nóng)村水利水電,2010,(9)：33-35.

[2] 陳 會,王 康,周祖昊. 基于排水過程分析的水稻灌區(qū)農(nóng)田面源污染模擬[J].農(nóng)業(yè)工程學(xué)報,2012,28(6):112-119.

[3] 常蒲婷,楊 侃,侯學(xué)勇. 突發(fā)性水污染事件模擬分析研究[J].水電能源科學(xué),2009,27(3):38-41.

[4] 張廣納,邵景安,王金亮. 基于農(nóng)業(yè)面源污染的三峽庫區(qū)重慶段水質(zhì)時空格局演變特征[J].自然資源學(xué)報,2015,30(11):1872-1881.

[5] 遼寧省大伙房水庫管理局.大伙房水庫志[M].北京:中國水利水電出版社，2006.

[6] 焉鴻啟,趙 文,郭 凱,等.遼寧省6座水源水庫富營養(yǎng)化狀況的分析與評價[J].大連海洋大學(xué)學(xué)報,2016,31(2):180-184.