趙琰鑫,張萬(wàn)順,湯 怡,吳 靜
(武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院,湖北 武漢 430079)
河網(wǎng)水動(dòng)力水質(zhì)模型是描述河道水體中污染物遷移轉(zhuǎn)化規(guī)律的數(shù)學(xué)模型,是進(jìn)行河流水質(zhì)模擬和水污染控制的有力工具。近幾年,河網(wǎng)水質(zhì)模型得到長(zhǎng)足的發(fā)展,如國(guó)外應(yīng)用較多的美國(guó)的QUAL-2E,WASP模型等[1]。國(guó)內(nèi)的彭虹等[2]建立了河流綜合水質(zhì)模型,李錦秀等[3]建立了三峽水庫(kù)整體一維水質(zhì)模型,儲(chǔ)君達(dá)[4]、韓龍喜等[5]建立了河網(wǎng)水質(zhì)模型并對(duì)模型求解方法進(jìn)行了改進(jìn)。但是以上研究仍存在一些不足,如其模型僅僅用于樹(shù)狀河網(wǎng)和單一河道,對(duì)于環(huán)狀的河網(wǎng)有限制。其次,其模型考慮的水質(zhì)變量及各個(gè)變量間的遷移轉(zhuǎn)化不夠全面。另外,針對(duì)復(fù)雜湖泊-河網(wǎng)區(qū)河道縱橫、水系呈網(wǎng)狀的特點(diǎn)建立的數(shù)值模型尚未見(jiàn)報(bào)道。
本文采用有限控制體積法對(duì)求解一、二維水動(dòng)力學(xué)和水質(zhì)模型控制方程,在一維、二維模型連接斷面處,采用一維模型模擬河道水位、流量變化,并作為隱式變量帶入二維湖泊模型進(jìn)行求解,實(shí)現(xiàn)了一、二維模型的耦合,并根據(jù)太湖典型流域河網(wǎng)區(qū)實(shí)測(cè)水文和水質(zhì)資料對(duì)耦合模型進(jìn)行率定和驗(yàn)證。
2.1 一維河網(wǎng)水動(dòng)力水質(zhì)模型
連續(xù)性方程:
動(dòng)量方程:
式中:uj為河網(wǎng)河道j的斷面平均流速;Qj為河網(wǎng)河道j的流量;Aj為河道j的過(guò)水面積;t為時(shí)間;qj為河段j的側(cè)流匯流流量;g為重力加速度;zj為河網(wǎng)河道j的水位;n為河道糙率;R為水力半徑。
污染物遷移轉(zhuǎn)換方程:
式中:Ci為污染物濃度;hj為河網(wǎng)河道j的水位;uj為河網(wǎng)河道j的流速;E為擴(kuò)散系數(shù);Sm為污染物排放量;kd為污染物的降解系數(shù),kd=k0θT-20,T為水溫,θ為系數(shù),取1~1.08,k0為常溫下的降解系數(shù)。
河網(wǎng)節(jié)點(diǎn)方程:
式中:cs是節(jié)點(diǎn)j的水質(zhì)濃度;ciin為入流的水質(zhì)濃度;Qiin為入流流量。
2.2 二維湖泊水動(dòng)力水質(zhì)模型
連續(xù)方程:
X方向動(dòng)量方程:
Y方向動(dòng)量方程:
式中:u、v為x、y方向的垂向平均速度;z為水面高程;h為水深;f為科氏力系數(shù)f=2Ωsinθ,Ω為地球旋轉(zhuǎn)的角頻率,θ為當(dāng)?shù)氐木暥?;γi為紊動(dòng)黏性系數(shù);ρa(bǔ)和ρw分別是空氣和水密度;fw為風(fēng)應(yīng)力系數(shù);wx、wy分別為x、y方向的風(fēng)速。
水質(zhì)遷移控制方程充分考慮水體污染物的對(duì)流、擴(kuò)散降解作用。水質(zhì)遷移控制方程:
式中:Ci為水中污染物i的濃度;Kx、Ky為x、y方向上的擴(kuò)散系數(shù)。
2.3 一維、二維模型耦合方法 在一、二維模型連接斷面處,根據(jù)兩種模型模擬的水位、流量、濃度相等的條件,實(shí)現(xiàn)一、二維模型的耦合[6]。研究中通過(guò)設(shè)置過(guò)渡單元實(shí)現(xiàn)這一耦合,過(guò)渡單元為一維模型單元與二維模型單元的連接單元。圖1為一維模型單元和二維模型單元的過(guò)渡單元網(wǎng)格布置。通過(guò)在連接斷面處補(bǔ)充物理量之間的關(guān)系(水位、流量、濃度相等),實(shí)現(xiàn)了一維模型與二維模型的耦合,即:水位連接條件:Z1=Z2;流量連接條件:Q1=∫Uεhεdε;水質(zhì)連接條件:C1=C2。其中:Z1、Z2分別為一、二維模型在連接斷面處的水位;C1、C2分別為一、二維模型在連接斷面處的水質(zhì)濃度;Q1為一維模型在連接斷面上的流量;Uε為二維模型在連接斷面法向上的流速;h為水深;ε為一、二維模型連接斷面坐標(biāo)。
2.4 數(shù)值離散格式 將水動(dòng)力和水質(zhì)模型方程改寫(xiě)成統(tǒng)一形式
采用非正交非交錯(cuò)網(wǎng)格,在控制體內(nèi),對(duì)流項(xiàng)采用迎風(fēng)格式處理,對(duì)上式進(jìn)行積分和離散,得到對(duì)流擴(kuò)散方程的離散方程:
其中:ap、anb分別是系數(shù)。
采用SIMPLE正交算法,獲得自由表面η校正方程和速度修正方程,即
η校正方程:
速度修正方程:
表面和速度方程組屬于同一類(lèi)對(duì)角型的代數(shù)方程組,可以應(yīng)用SIMPLE方法進(jìn)行快速求解。
3.1 區(qū)域概況 太湖流域位于長(zhǎng)江三角洲,地跨江蘇、浙江和上海二省一市,流域面積36500km2。流域內(nèi)是我國(guó)經(jīng)濟(jì)最發(fā)達(dá)、目前經(jīng)濟(jì)發(fā)展最快的地區(qū)之一,水質(zhì)污染問(wèn)題十分突出。流域內(nèi)河網(wǎng)縱橫交錯(cuò),閘站眾多,水流運(yùn)動(dòng)十分復(fù)雜。
針對(duì)太湖流域復(fù)雜的河流水系網(wǎng)絡(luò)結(jié)構(gòu)特征,選取太湖湖西滆湖典型區(qū)域?yàn)檠芯糠秶?,區(qū)域內(nèi)骨干河道及次級(jí)河道總長(zhǎng)度1642km,承擔(dān)著防洪排澇、農(nóng)田灌溉、城鄉(xiāng)供水保障和交通航運(yùn)等重要功能。區(qū)域內(nèi)主要河道有:武宜運(yùn)河、武進(jìn)港、采菱港、太滆運(yùn)河、扁擔(dān)河、湟里河、蕪申運(yùn)河、孟津河、中干河等。其中,扁擔(dān)河、南運(yùn)河、采菱港、武進(jìn)港等運(yùn)河南部的水網(wǎng)承接運(yùn)河來(lái)水,輸向滆湖或太湖,區(qū)內(nèi)河流主流向自西往東,自北往南。受長(zhǎng)江、太湖相對(duì)水位的影響和通江河口閘門(mén)控制,河道流向不穩(wěn),常有滯流、倒流現(xiàn)象;整體河道呈現(xiàn)平原河網(wǎng)低流速、小流量的特性。
3.2 河網(wǎng)概化 太湖地區(qū)典型流域模擬計(jì)算中的河網(wǎng)、湖泊是在天然河網(wǎng)湖泊的基礎(chǔ)上根據(jù)河道輸水能力相等的原理進(jìn)行合并、概化,概化河道的斷面為梯形。依據(jù)河網(wǎng)結(jié)構(gòu)和河道匯流特點(diǎn),將河網(wǎng)劃分成110個(gè)河段,共91個(gè)計(jì)算節(jié)點(diǎn)和519個(gè)計(jì)算斷面。
3.3 模型參數(shù)
(1)糙率。參考相關(guān)研究報(bào)告,太湖河網(wǎng)區(qū)河道糙率取0.02-0.03,湖底糙率取0.002~0.025。
(2)降解系數(shù)。根據(jù)太湖流域河網(wǎng)區(qū)水質(zhì)計(jì)算的經(jīng)驗(yàn)值氨氮降解系數(shù)取0.05~0.20d-1,COD降解系數(shù)取 0.08~0.25d-1, TN降解系數(shù)取0.06~0.15 d-1,TP降解系數(shù)取0.05~0.08 d-1,。
3.4 模型的率定驗(yàn)證 采用2007年上游入流斷面夏溪河夏溪橋站、湟里河湟里站、北干河?xùn)|安橋站、江南運(yùn)河常州站等水文站點(diǎn)流量過(guò)程,下游太湖百瀆口站、大浦口站和宜興站的水位過(guò)程,作為水動(dòng)力計(jì)算的上下游邊界條件,采用研究區(qū)域2007年排污負(fù)荷作為水質(zhì)模型驗(yàn)證的計(jì)算條件,對(duì)模型進(jìn)行率定驗(yàn)證,采用水文站實(shí)測(cè)水溫作為水質(zhì)模型計(jì)算的溫度條件。
3.4.1 水動(dòng)力模型的驗(yàn)證 采用太滆運(yùn)河黃埝橋站和漕橋河漕橋站2007年逐日實(shí)測(cè)流量、水位過(guò)程對(duì)模型進(jìn)行水動(dòng)力學(xué)驗(yàn)證。水位驗(yàn)證結(jié)果如圖2和圖3所示,由圖可知,水位計(jì)算值與實(shí)測(cè)值擬合程度較高,兩驗(yàn)證斷面絕對(duì)誤差均小于0.55m,相對(duì)誤差均小于12.97%,模型對(duì)河網(wǎng)水位模擬具有較高精度。
流量驗(yàn)證結(jié)果如圖4和圖5所示,由圖可知,流量計(jì)算值整體變化規(guī)律與實(shí)測(cè)值變化規(guī)律具有較好的一致性。由于研究區(qū)域?yàn)槠皆泳W(wǎng)區(qū)域,河道具有低流量、小流速的特點(diǎn),部分時(shí)段河道流量及小,導(dǎo)致驗(yàn)證的相對(duì)誤差較大。模擬值與實(shí)測(cè)值相對(duì)誤差在20%以內(nèi)的天數(shù)占全年的82.5%,模型精度較高。
二維湖泊水動(dòng)力模型采用滆湖坊前站2007年逐日實(shí)測(cè)水位數(shù)據(jù)進(jìn)行驗(yàn)證,驗(yàn)證結(jié)果見(jiàn)圖6,由圖可以看出,計(jì)算值與實(shí)測(cè)值規(guī)律呈現(xiàn)較好一致性,模擬的絕對(duì)誤差控制在0.33m以內(nèi),相對(duì)誤差控制在7.36%以內(nèi),模擬效果較好。
3.4.2 水質(zhì)模型驗(yàn)證 在研究區(qū)域常規(guī)水質(zhì)監(jiān)測(cè)斷面中選取與河網(wǎng)計(jì)算斷面重合的黃埝橋、雪埝橋和漕橋等3個(gè)斷面,采用2007年逐月實(shí)測(cè)水質(zhì)數(shù)據(jù)對(duì)一維河網(wǎng)水質(zhì)模型進(jìn)行驗(yàn)證,模型計(jì)算結(jié)果與實(shí)測(cè)值對(duì)比,由表1黃埝橋、雪埝橋及漕橋斷面氨氮、COD、TN和TP的模擬值和實(shí)測(cè)值的平均誤差對(duì)比可以看出其相對(duì)誤差多在30%以內(nèi),平均誤差在25%左右,表明模型模擬值與實(shí)測(cè)值吻合較好,能夠滿足河網(wǎng)一維水質(zhì)模擬的要求。
表1 各斷面指標(biāo)平均誤差值(%)
表2 各指標(biāo)相對(duì)誤差值(%)
選取滆湖內(nèi)的滆湖北常規(guī)水質(zhì)監(jiān)測(cè)點(diǎn),采用2007年1月、4月、7月和10月實(shí)測(cè)水質(zhì)數(shù)據(jù)對(duì)二維湖泊水質(zhì)模型進(jìn)行驗(yàn)證,由表2滆湖北氨氮、COD、TN和TP的模擬值和實(shí)測(cè)值的對(duì)比可以看出,其相對(duì)誤差多在30%以內(nèi),平均誤差在24%左右,表明模型模擬值與實(shí)測(cè)值吻合較好,能夠滿足湖泊二維水質(zhì)模擬的要求。
本文提出了適用于復(fù)雜湖泊-河網(wǎng)區(qū)域的一維、二維水動(dòng)力和水質(zhì)耦合數(shù)學(xué)模型。采用有限控制體積法對(duì)求解一、二維水動(dòng)力學(xué)和水質(zhì)模型控制方程,在一維、二維模型連接斷面處,利用兩種模型模擬的水位、流量相同的條件,將一維與二維模型有機(jī)連接為一個(gè)整體,實(shí)現(xiàn)一維模型與二維模型的耦合。應(yīng)用太湖典型流域河網(wǎng)區(qū)實(shí)測(cè)水文和水質(zhì)資料對(duì)耦合模型進(jìn)行率定和驗(yàn)證,表明模型計(jì)算值與實(shí)測(cè)資料吻合較好,因此一、二維耦合模型的設(shè)計(jì)是合理的,適用于復(fù)雜湖泊-河網(wǎng)區(qū)的水動(dòng)力和水質(zhì)變化的模擬和研究。模型可用于深入分析滆湖與周邊河網(wǎng)中水質(zhì)濃度整體分布狀況、響應(yīng)關(guān)系,對(duì)流域不同污染物消減和水利調(diào)度調(diào)度處理等方案等進(jìn)行情景模擬,并為太湖河網(wǎng)區(qū)的水污染治理工程提供技術(shù)支撐。
[1]鄭孝宇,褚君達(dá),朱維斌.河網(wǎng)非穩(wěn)態(tài)水環(huán)境容量研究[J].水科學(xué)進(jìn)展,1997,8(1):25-31.
[2]彭虹,張萬(wàn)順,夏軍.河流綜合水質(zhì)生態(tài)數(shù)值模型[J].武漢大學(xué)學(xué)報(bào),2002,35(4):56-59.
[3]李錦秀,廖文根,黃真理.三峽水庫(kù)整體一維水質(zhì)數(shù)學(xué)模擬研究[J].水利學(xué)報(bào),2002(12):7-10.
[4]褚君達(dá).河網(wǎng)對(duì)流輸移問(wèn)題的求解及應(yīng)用[J].水利學(xué)報(bào),1992(7):30-34.
[5]韓龍喜,金忠青.三角聯(lián)解法水力水質(zhì)模型的糙率反演及面污染源計(jì)算[J].水利學(xué)報(bào),1998(7):30-34.
[6]Wanshun Zhang,Yan Wang,Hong Peng.A coupled water quantity-quality model for water allocation analysis[J].Water Resource Manage,2010(24):485-511.