劉中峰,李 然,陳明千,黃 翔,曹玲玲

(四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,四川成都 610065)

近年來,二灘、溪洛渡、白鶴灘等一大批高壩工程已建、在建或待建。這些工程在帶來巨大經(jīng)濟(jì)效益和社會(huì)效益的同時(shí),不可避免地對(duì)河流水質(zhì)產(chǎn)生了一定的影響。因此,對(duì)大型水庫的水質(zhì)預(yù)測(cè)模擬具有重大的理論和工程意義。

大型水庫的水質(zhì)要素分布有著與天然河道不同的特點(diǎn)[1]。與天然河道相比,水庫水深顯著增大,流速明顯減小,水體水環(huán)境條件發(fā)生改變,同時(shí)由于地形、流場(chǎng)等的三維效應(yīng),溫度、水質(zhì)要素可能出現(xiàn)垂向分層現(xiàn)象。目前,一、二維水質(zhì)預(yù)測(cè)模型已經(jīng)比較成熟,三維模型的研究也正逐漸開展。聶晶[2]根據(jù)質(zhì)量守恒定律和Fick定律,推導(dǎo)出水庫水體總磷三維遷移轉(zhuǎn)化生態(tài)數(shù)學(xué)模型,并將其應(yīng)用到密云水庫,取得了較好的模擬結(jié)果;MAO等[3]建立了一種物理過程和生化過程耦合的三維模型,用來模擬太湖中營(yíng)養(yǎng)物質(zhì)的年變化情況;George等[4-5]針對(duì) Lake Washington提出了一個(gè)復(fù)雜的三維富營(yíng)養(yǎng)化模型,成功分析了浮游植物的動(dòng)力學(xué)特征。

以上模型能夠較好地模擬水體中水質(zhì)要素的基本分布規(guī)律,但它們大多針對(duì)特定的研究對(duì)象,且都沒有考慮溫度分層的影響,而在大型水庫中溫度分層對(duì)水質(zhì)要素的分布有較大的影響[6]。另外,上述模型都是非穩(wěn)態(tài)模型,這在提高模型預(yù)測(cè)精度的同時(shí),卻很大程度上提高了計(jì)算代價(jià),尤其是大型水庫的計(jì)算區(qū)域非常大,采用非穩(wěn)態(tài)模型計(jì)算所花費(fèi)的時(shí)間非常長(zhǎng)。在實(shí)際工程中,大型水庫在一個(gè)水文周期(豐、平、枯)內(nèi),水文條件和污染負(fù)荷往往變化不大,采用穩(wěn)態(tài)模型就能滿足工程要求。因此,筆者建立了一種考慮溫度分層影響的大型水庫穩(wěn)態(tài)三維水質(zhì)預(yù)測(cè)模型,并用二灘魚感魚河支庫的實(shí)測(cè)資料對(duì)模型進(jìn)行了驗(yàn)證。將模型應(yīng)用到某大型水庫,成功模擬出水質(zhì)要素的三維分布特性。

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

1.1 數(shù)學(xué)模型方程

數(shù)學(xué)模型方程包括連續(xù)方程、動(dòng)量方程、k方程、ε方程及污染物輸移方程。

連續(xù)方程

動(dòng)量方程

k方程

ε方程)

污染物輸移方程

其中

式中:i為坐標(biāo)方向,i=1,2,3;ρ為水的密度,kg/m3;Ui為i方向的水流速度值,m/s;P為壓強(qiáng),Pa;μ為動(dòng)力黏性系數(shù),m2/s;μt為紊動(dòng)渦黏系數(shù),m2/s;k為紊動(dòng)動(dòng)能,m2/s2;ε為紊動(dòng)動(dòng)能的耗散率,m2/s3;σk,σε和 σC分別為紊動(dòng)動(dòng)能、紊動(dòng)動(dòng)能的耗散率和水質(zhì)要素的Prandtl數(shù),對(duì)應(yīng)值分別為1.0,1.3和1.0;φ為DO,CODCr,TP,TN等水質(zhì)要素的質(zhì)量濃度,mg/L;Cμ,C1ε和 C2ε為模型常數(shù) ,取值由基本實(shí)驗(yàn)確定[7],分別為0.09,1.44和1.92;Sφ為污染物輸移方程的源項(xiàng),具體形式見第1.2節(jié)。

1.2 污染物輸移方程源項(xiàng)表達(dá)式

污染物輸移方程源項(xiàng)表達(dá)式分為2種形式。

a.DO方程源項(xiàng)。DO方程源項(xiàng)分表層水體和下層水體2部分區(qū)別考慮。對(duì)于表層水體,DO方程源項(xiàng)包括表面復(fù)氧和降解耗氧2部分;對(duì)于下層水體,DO方程源項(xiàng)只有降解耗氧,而溶解氧的補(bǔ)充依靠輸移擴(kuò)散過程。表層水體和下層水體的DO方程源項(xiàng) Sφ表及Sφ下的表達(dá)式分別為

其中

式中:k1(T)為溫度T時(shí)CODCr的降解系數(shù),s-1;k2(T)為溫度T時(shí)的復(fù)氧系數(shù),s-1;CS為水體中飽和溶解氧質(zhì)量濃度(計(jì)算公式來源于文獻(xiàn)[8]),mg/L;CDO為水體DO質(zhì)量濃度,mg/L;LC為水體CODCr質(zhì)量濃度,mg/L;θ為溫度修正系數(shù);k2(20)為溫度在20℃時(shí)的復(fù)氧系數(shù)(計(jì)算公式來源于文獻(xiàn)[8]),s-1;U10為水面上10m處的風(fēng)速,m/s;h為表層水體厚度,計(jì)算中取0.5m。

b.一般污染物輸移方程源項(xiàng)。假定污染物隨時(shí)間變化為一階動(dòng)力學(xué)反應(yīng),源項(xiàng)表達(dá)式為

其中

式中:k(T)為溫度T時(shí)水質(zhì)要素的降解(衰減)系數(shù),s-1;k(20)為溫度在20℃時(shí)相應(yīng)水質(zhì)要素的降解(衰減)系數(shù),s-1。

1.3 模型求解

控制方程組的離散采用有限體積法。離散后的代數(shù)方程組采用SIMPLEC算法求解。

2 參數(shù)率定與模型驗(yàn)證

采用二灘魚感魚河支庫原型觀測(cè)資料對(duì)模型進(jìn)行參數(shù)率定與驗(yàn)證。

2.1 二灘水庫概況

二灘水庫位于金沙江支流雅礱江上。最大壩高240m,回水長(zhǎng)度約140km。魚感魚河支庫是二灘庫區(qū)右岸的主要支庫,回水長(zhǎng)度約30km。2008年7月對(duì)魚感魚河支庫水溫和DO分布進(jìn)行了原型觀測(cè),同時(shí)收集整理了有關(guān)部門的魚感魚河支庫水文資料、TP與TN定點(diǎn)取樣監(jiān)測(cè)資料及相關(guān)污染源資料,用于對(duì)三維水質(zhì)模型進(jìn)行驗(yàn)證。

2.2 計(jì)算區(qū)域及網(wǎng)格劃分

二灘魚感魚河支庫上游由永興河、惠民河和新坪河交匯而成。計(jì)算區(qū)域上游至永興河匯口,下游至魚感魚河入二灘主庫的匯口,全長(zhǎng)約30km。計(jì)算區(qū)域上表面水位與水庫正常蓄水位一致,為1200m。計(jì)算區(qū)域垂向最大水深為150m,橫向最大水面寬度為1165m。計(jì)算區(qū)域平面如圖1所示。

圖1 計(jì)算區(qū)域平面示意圖

計(jì)算區(qū)域共劃分163740個(gè)網(wǎng)格,其中橫向最大網(wǎng)格間距為55m,縱向最大網(wǎng)格間距為55m,垂向最大網(wǎng)格間距為10m。網(wǎng)格劃分情況如圖2所示。

圖2 魚感魚河支庫網(wǎng)格劃分示意圖

2.3 污染源

魚感魚河支庫的點(diǎn)源污染主要來自漁門鎮(zhèn)的生活污水排放。該點(diǎn)源的污水排放量為24.97萬t/a,污水中各污染物的質(zhì)量濃度分別為:ρ(CODCr)=363mg/L,ρ(TP)=4.05mg/L,ρ(TN)=61.6mg/L 。

經(jīng)調(diào)查,庫區(qū)所在位置為山區(qū),耕地稀少,工農(nóng)業(yè)很不發(fā)達(dá),面源污染主要來自庫區(qū)內(nèi)的網(wǎng)箱養(yǎng)魚。收集了庫區(qū)漁業(yè)產(chǎn)量、飼料用量等漁業(yè)養(yǎng)殖資料,并根據(jù)文獻(xiàn)調(diào)研成果[9-10],計(jì)算得到魚感魚河支庫面源污染物排放量。在計(jì)算中,面源污染物隨沿程匯流(末斷面流量與起始斷面流量之差,認(rèn)為從兩岸沿程均勻匯入)從兩岸進(jìn)入水體。面源污染排放特征值見表1,其中質(zhì)量濃度值為污染物排放量與沿程匯流量之比。

2.4 邊界條件

入流邊界包括永興河、惠民河和新坪河3條支流入口、漁門鎮(zhèn)生活污水排放口及兩岸沿程匯流,設(shè)置為速度入口。在入流邊界給定速度分布和水質(zhì)要素濃度分布,紊動(dòng)特征量(紊動(dòng)動(dòng)能 k和紊動(dòng)能耗散率ε)根據(jù)入口流速 u0和入口水深 H0按下述經(jīng)驗(yàn)公式確定[11]:

表1 魚感魚河支庫面源污染排放特征值

出流邊界為魚感魚河支庫進(jìn)入二灘庫區(qū)的匯口斷面,設(shè)置為充分發(fā)展的紊流邊界,即所有標(biāo)量(k,ε,φ)法向梯度為零,平面的法向速度分量梯度為零。在水面上給定對(duì)稱面邊界條件,即沒有垂向的速度分量。水底為固壁邊界,水動(dòng)力學(xué)條件按標(biāo)準(zhǔn)壁函數(shù)法處理,并認(rèn)為通過壁面的污染物通量為零。計(jì)算中,水溫采用魚感魚河支庫2008年7月的實(shí)測(cè)數(shù)據(jù)。

2.5 驗(yàn)證結(jié)果與分析

模型中各水質(zhì)要素的 k(20)及θ初值參照已有的研究成果[12-13]取值,見表2。

表2 各水質(zhì)要素k(20)及 θ取值

采用表2中模型參數(shù)初值,模擬得到二灘水庫各水質(zhì)要素的空間分布。表3給出了Ⅰ斷面(新坪河匯口)、Ⅱ斷面(新坪河匯口下游3000m)、Ⅲ斷面(新坪河匯口下游6000m)自由水面中間點(diǎn)上TP和TN模擬結(jié)果與實(shí)測(cè)結(jié)果的比較。由表中數(shù)據(jù)可以看出,數(shù)值模擬結(jié)果與實(shí)測(cè)結(jié)果吻合較好。

表3 自由水面中間點(diǎn)TP和TN模擬結(jié)果與實(shí)測(cè)結(jié)果對(duì)比

圖3給出了惠民河匯口斷面和新坪河匯口斷面中垂線水面下15m范圍內(nèi)DO模擬結(jié)果與實(shí)測(cè)結(jié)果的比較。從圖3可以看出,模擬結(jié)果與實(shí)測(cè)結(jié)果總體趨勢(shì)接近,呈現(xiàn)上大下小的趨勢(shì)。計(jì)算可得2個(gè)斷面中垂線DO質(zhì)量濃度的計(jì)算值和測(cè)量值的均方根誤差分別為0.78mg/L和0.31mg/L,表明兩者吻合較好。

圖3 斷面中垂線DO質(zhì)量濃度模擬結(jié)果與實(shí)測(cè)結(jié)果對(duì)比

綜合分析認(rèn)為,基于上述模型參數(shù)的預(yù)測(cè)結(jié)果與監(jiān)測(cè)結(jié)果誤差較小,能夠滿足實(shí)際要求,可以用于對(duì)TP,TN,CODCr等污染物的數(shù)值模擬。

3 模型誤差與適用性分析

在二灘魚感魚河支庫水質(zhì)模擬中,基于對(duì)網(wǎng)箱養(yǎng)魚規(guī)模的調(diào)查以及國(guó)內(nèi)關(guān)于網(wǎng)箱養(yǎng)魚污染的研究成果,進(jìn)行了面源計(jì)算。由于目前關(guān)于網(wǎng)箱養(yǎng)魚污染物排放特征的研究資料和水平有限,因此面源污染負(fù)荷水平的計(jì)算誤差可能直接影響到水質(zhì)模型預(yù)測(cè)結(jié)果的精度。

分析DO模擬結(jié)果與觀測(cè)結(jié)果的差別,認(rèn)為除本模型中考慮的污染物耗氧和表面復(fù)氧過程外,實(shí)際水庫中還存在底泥耗氧、光合作用產(chǎn)氧等過程。這些過程的綜合作用可能造成DO沿水深的非線性變化,而筆者建立的模型中忽略了這些作用,因此存在一定誤差。另外,本模型沒有考慮風(fēng)力作用對(duì)水庫垂向混合的影響,因此造成表面DO濃度較觀測(cè)值偏低。模型僅考慮了溫度分層對(duì)模型中降解系數(shù)、復(fù)氧系數(shù)的影響,沒有考慮秋、冬季節(jié)由溫度差引起的垂向混合作用,因此模型對(duì)秋、冬季節(jié)DO模擬誤差可能較大。

對(duì)于水庫中TP和TN等污染物的模擬,由于它們?cè)谒畮熘械拇瓜蚍謱有Ч幻黠@,同時(shí)相對(duì)于DO來說,它們受表層復(fù)氧過程的影響較小,因此忽略風(fēng)力混合作用和溫度差產(chǎn)生的誤差相對(duì)于DO來說要小得多。另外,水質(zhì)取樣中的隨機(jī)性以及測(cè)點(diǎn)代表性也可能影響到對(duì)模型精度的判斷。

4 模型應(yīng)用

將上述模型應(yīng)用于國(guó)內(nèi)某在建大型水庫的水質(zhì)預(yù)測(cè)研究,成功模擬出該水庫BOD,DO,CODCr等水質(zhì)要素的三維分布情況。受篇幅所限,在此僅給出該水庫某庫灣段豐水期CODCr的模擬結(jié)果。該庫灣長(zhǎng)約11km,共劃分190×31×27(縱向×橫向×垂向)個(gè)網(wǎng)格。

圖4給出了自由水面CODCr質(zhì)量濃度分布的模擬結(jié)果。可以看出,受污水排放影響,石灰窯溝匯口出現(xiàn)污染帶,污染帶橫向?qū)挾茸畲筮_(dá)到70m,縱向延伸至庫灣出口。污染帶橫向影響范圍較小,這主要是由于溝口下游河道束窄而形成回流區(qū),抑制了CODCr在橫向的擴(kuò)散。

圖4 豐水期自由水面ρ(CODCr)分布示意圖(單位:mg/L)

圖5 給出了黑水河匯口上游斷面CODCr質(zhì)量濃度分布的模擬結(jié)果。從圖5中可以看出,CODCr質(zhì)量濃度分布在垂向沒有出現(xiàn)明顯的分層現(xiàn)象。但是由于石灰窯溝污水匯入的影響,右岸CODCr質(zhì)量濃度比左岸高。

圖5 豐水期黑水河匯口上游斷面ρ(CODCr)分布示意圖(單位:mg/L)

5 結(jié)論與展望

建立了大型水庫三維水質(zhì)預(yù)測(cè)模型,并在污染物輸移方程中采用了考慮溫度分層影響的源項(xiàng)表達(dá)式。采用二灘魚感魚河支庫原型觀測(cè)資料對(duì)模型進(jìn)行驗(yàn)證,模擬結(jié)果和監(jiān)測(cè)數(shù)據(jù)吻合較好,表明該模型是可靠的。將模型應(yīng)用于某在建大型水庫的水質(zhì)預(yù)測(cè),成功模擬出水質(zhì)三維分布特性,模擬結(jié)果對(duì)于開展水環(huán)境影響的定量預(yù)測(cè)評(píng)價(jià)、制定水庫水污染防治規(guī)劃、保護(hù)庫區(qū)水環(huán)境具有一定的理論意義和應(yīng)用價(jià)值。

需要指出的是,本文工作是對(duì)三維水質(zhì)預(yù)測(cè)模型研究的初步探索。模型在污染源特別是面源估算、水溫分布引起的密度流影響、水庫底泥及光合作用影響等方面尚存不足。此外,沒有模擬水溫,而是將水溫?cái)?shù)據(jù)作為邊界條件輸入模型。因此,在深入污染源研究基礎(chǔ)上,綜合考慮各種影響因素,建立水溫、水質(zhì)耦合求解的三維水質(zhì)預(yù)測(cè)模型值得進(jìn)一步研究和探索。

[1]LI Ran,FU Xiao-ying,LI Jia,et al.Three-dimensional numerical simulation of BOD-DO for large reservoir[J].Journal of Hydrodynamics:Ser B,2002,3:50-53.

[2]聶晶.水庫水體總磷三維數(shù)學(xué)模型及其應(yīng)用[D].長(zhǎng)春:吉林大學(xué),2004.

[3]MAO Jing-qiao,CHEN Qiu-wen,CHEN Yong-can.Threedimensional eutrophication model and application to taihu lake,China[J].Journal of Environmental Sciences,2008,20(3):278-284.

[4]ARHONDITSIS G B,BRETT M T.Eutrophication model for Lake Washington(USA)PartⅠ:model description and sensitivity analysis[J].Ecological Modelling,2005,187(2/3):140-178.

[5]AR HONDITSIS G B,BRETT M T.Eutrophication model for lake washington(USA)PartⅡ:model calibration and system dynamic analysis[J].Ecological Modelling,2005,187(2/3):179-200.

[6]王冠,韓龍喜,常文婷.基于立面二維水動(dòng)力-水溫耦合模型的水庫水溫分布[J].水資源保護(hù),2009,25(2):59-63.

[7]金忠青.N-S方程的數(shù)值解和紊流模型[M].南京:河海大學(xué)出版社,1989:54.

[8]JEAN-MARC T.Simulation of a mesotrophic reservoir(Lake Pareloup)over a long period(1983-1998)using ASTER2000 biological model[J].Water Research,2004,38(2):393-403.

[9]陳丁,鄭愛榕.網(wǎng)箱養(yǎng)殖的氮、磷和有機(jī)物的污染及估算[J].福建農(nóng)業(yè)學(xué)報(bào),2005,20(B12):57-62.

[10]陳德春.投餌式網(wǎng)箱養(yǎng)魚對(duì)水質(zhì)的影響及保護(hù)措施[J].水利漁業(yè),1993(2):23-26.

[11]鄧云.大型深水庫的水溫預(yù)測(cè)研究[D].成都:四川大學(xué),2003.

[12]饒群.大型水體富營(yíng)養(yǎng)化數(shù)學(xué)模擬的研究[D].南京:河海大學(xué),2002.

[13]黃真理,李玉梁,李錦秀,等.三峽水庫水質(zhì)預(yù)測(cè)和環(huán)境容量計(jì)算[M].北京:中國(guó)水利水電出版社,2006:211-221.