湖泊-河網(wǎng)耦合水動(dòng)力水質(zhì)模型研究

趙琰鑫，張萬順，湯 怡，吳 靜

（武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430079）

1 研究背景

河網(wǎng)水動(dòng)力水質(zhì)模型是描述河道水體中污染物遷移轉(zhuǎn)化規(guī)律的數(shù)學(xué)模型，是進(jìn)行河流水質(zhì)模擬和水污染控制的有力工具。近幾年，河網(wǎng)水質(zhì)模型得到長足的發(fā)展，如國外應(yīng)用較多的美國的QUAL-2E，WASP模型等［1］。國內(nèi)的彭虹等［2］建立了河流綜合水質(zhì)模型，李錦秀等［3］建立了三峽水庫整體一維水質(zhì)模型，儲君達(dá)［4］、韓龍喜等［5］建立了河網(wǎng)水質(zhì)模型并對模型求解方法進(jìn)行了改進(jìn)。但是以上研究仍存在一些不足，如其模型僅僅用于樹狀河網(wǎng)和單一河道，對于環(huán)狀的河網(wǎng)有限制。其次，其模型考慮的水質(zhì)變量及各個(gè)變量間的遷移轉(zhuǎn)化不夠全面。另外，針對復(fù)雜湖泊-河網(wǎng)區(qū)河道縱橫、水系呈網(wǎng)狀的特點(diǎn)建立的數(shù)值模型尚未見報(bào)道。

本文采用有限控制體積法對求解一、二維水動(dòng)力學(xué)和水質(zhì)模型控制方程，在一維、二維模型連接斷面處，采用一維模型模擬河道水位、流量變化，并作為隱式變量帶入二維湖泊模型進(jìn)行求解，實(shí)現(xiàn)了一、二維模型的耦合，并根據(jù)太湖典型流域河網(wǎng)區(qū)實(shí)測水文和水質(zhì)資料對耦合模型進(jìn)行率定和驗(yàn)證。

2 湖泊-河網(wǎng)區(qū)的一/二維水動(dòng)力水質(zhì)耦合模型

2.1 一維河網(wǎng)水動(dòng)力水質(zhì)模型

連續(xù)性方程：

動(dòng)量方程：

式中：uj為河網(wǎng)河道j的斷面平均流速；Qj為河網(wǎng)河道j的流量；Aj為河道j的過水面積；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ì)遷移控制方程充分考慮水體污染物的對流、擴(kuò)散降解作用。水質(zhì)遷移控制方程：

式中：Ci為水中污染物i的濃度；Kx、Ky為x、y方向上的擴(kuò)散系數(shù)。

2.3 一維、二維模型耦合方法 在一、二維模型連接斷面處，根據(jù)兩種模型模擬的水位、流量、濃度相等的條件，實(shí)現(xiàn)一、二維模型的耦合［6］。研究中通過設(shè)置過渡單元實(shí)現(xiàn)這一耦合，過渡單元為一維模型單元與二維模型單元的連接單元。圖1為一維模型單元和二維模型單元的過渡單元網(wǎng)格布置。通過在連接斷面處補(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ì)模型方程改寫成統(tǒng)一形式

采用非正交非交錯(cuò)網(wǎng)格，在控制體內(nèi)，對流項(xiàng)采用迎風(fēng)格式處理，對上式進(jìn)行積分和離散，得到對流擴(kuò)散方程的離散方程：

其中：ap、anb分別是系數(shù)。

采用SIMPLE正交算法，獲得自由表面η校正方程和速度修正方程，即

η校正方程：

速度修正方程：

表面和速度方程組屬于同一類對角型的代數(shù)方程組，可以應(yīng)用SIMPLE方法進(jìn)行快速求解。

3 耦合模型的應(yīng)用實(shí)例

3.1 區(qū)域概況 太湖流域位于長江三角洲，地跨江蘇、浙江和上海二省一市，流域面積36500km2。流域內(nèi)是我國經(jīng)濟(jì)最發(fā)達(dá)、目前經(jīng)濟(jì)發(fā)展最快的地區(qū)之一，水質(zhì)污染問題十分突出。流域內(nèi)河網(wǎng)縱橫交錯(cuò)，閘站眾多，水流運(yùn)動(dòng)十分復(fù)雜。

針對太湖流域復(fù)雜的河流水系網(wǎng)絡(luò)結(jié)構(gòu)特征，選取太湖湖西滆湖典型區(qū)域?yàn)檠芯糠秶瑓^(qū)域內(nèi)骨干河道及次級河道總長度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)河來水，輸向滆湖或太湖，區(qū)內(nèi)河流主流向自西往東，自北往南。受長江、太湖相對水位的影響和通江河口閘門控制，河道流向不穩(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)流量過程，下游太湖百瀆口站、大浦口站和宜興站的水位過程，作為水動(dòng)力計(jì)算的上下游邊界條件，采用研究區(qū)域2007年排污負(fù)荷作為水質(zhì)模型驗(yàn)證的計(jì)算條件，對模型進(jìn)行率定驗(yàn)證，采用水文站實(shí)測水溫作為水質(zhì)模型計(jì)算的溫度條件。

3.4.1 水動(dòng)力模型的驗(yàn)證 采用太滆運(yùn)河黃埝橋站和漕橋河漕橋站2007年逐日實(shí)測流量、水位過程對模型進(jìn)行水動(dòng)力學(xué)驗(yàn)證。水位驗(yàn)證結(jié)果如圖2和圖3所示，由圖可知，水位計(jì)算值與實(shí)測值擬合程度較高，兩驗(yàn)證斷面絕對誤差均小于0.55m，相對誤差均小于12.97%，模型對河網(wǎng)水位模擬具有較高精度。

流量驗(yàn)證結(jié)果如圖4和圖5所示，由圖可知，流量計(jì)算值整體變化規(guī)律與實(shí)測值變化規(guī)律具有較好的一致性。由于研究區(qū)域?yàn)槠皆泳W(wǎng)區(qū)域，河道具有低流量、小流速的特點(diǎn)，部分時(shí)段河道流量及小，導(dǎo)致驗(yàn)證的相對誤差較大。模擬值與實(shí)測值相對誤差在20%以內(nèi)的天數(shù)占全年的82.5%，模型精度較高。

二維湖泊水動(dòng)力模型采用滆湖坊前站2007年逐日實(shí)測水位數(shù)據(jù)進(jìn)行驗(yàn)證，驗(yàn)證結(jié)果見圖6，由圖可以看出，計(jì)算值與實(shí)測值規(guī)律呈現(xiàn)較好一致性，模擬的絕對誤差控制在0.33m以內(nèi)，相對誤差控制在7.36%以內(nèi)，模擬效果較好。

3.4.2 水質(zhì)模型驗(yàn)證 在研究區(qū)域常規(guī)水質(zhì)監(jiān)測斷面中選取與河網(wǎng)計(jì)算斷面重合的黃埝橋、雪埝橋和漕橋等3個(gè)斷面，采用2007年逐月實(shí)測水質(zhì)數(shù)據(jù)對一維河網(wǎng)水質(zhì)模型進(jìn)行驗(yàn)證，模型計(jì)算結(jié)果與實(shí)測值對比，由表1黃埝橋、雪埝橋及漕橋斷面氨氮、COD、TN和TP的模擬值和實(shí)測值的平均誤差對比可以看出其相對誤差多在30%以內(nèi)，平均誤差在25%左右，表明模型模擬值與實(shí)測值吻合較好，能夠滿足河網(wǎng)一維水質(zhì)模擬的要求。

表1 各斷面指標(biāo)平均誤差值（%）

表2 各指標(biāo)相對誤差值（%）

選取滆湖內(nèi)的滆湖北常規(guī)水質(zhì)監(jiān)測點(diǎn)，采用2007年1月、4月、7月和10月實(shí)測水質(zhì)數(shù)據(jù)對二維湖泊水質(zhì)模型進(jìn)行驗(yàn)證，由表2滆湖北氨氮、COD、TN和TP的模擬值和實(shí)測值的對比可以看出，其相對誤差多在30%以內(nèi)，平均誤差在24%左右，表明模型模擬值與實(shí)測值吻合較好，能夠滿足湖泊二維水質(zhì)模擬的要求。

4 結(jié)論

本文提出了適用于復(fù)雜湖泊-河網(wǎng)區(qū)域的一維、二維水動(dòng)力和水質(zhì)耦合數(shù)學(xué)模型。采用有限控制體積法對求解一、二維水動(dòng)力學(xué)和水質(zhì)模型控制方程，在一維、二維模型連接斷面處，利用兩種模型模擬的水位、流量相同的條件，將一維與二維模型有機(jī)連接為一個(gè)整體，實(shí)現(xiàn)一維模型與二維模型的耦合。應(yīng)用太湖典型流域河網(wǎng)區(qū)實(shí)測水文和水質(zhì)資料對耦合模型進(jìn)行率定和驗(yàn)證，表明模型計(jì)算值與實(shí)測資料吻合較好，因此一、二維耦合模型的設(shè)計(jì)是合理的，適用于復(fù)雜湖泊-河網(wǎng)區(qū)的水動(dòng)力和水質(zhì)變化的模擬和研究。模型可用于深入分析滆湖與周邊河網(wǎng)中水質(zhì)濃度整體分布狀況、響應(yīng)關(guān)系，對流域不同污染物消減和水利調(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］彭虹，張萬順，夏軍.河流綜合水質(zhì)生態(tài)數(shù)值模型［J］.武漢大學(xué)學(xué)報(bào)，2002，35（4）：56-59.

［3］李錦秀，廖文根，黃真理.三峽水庫整體一維水質(zhì)數(shù)學(xué)模擬研究［J］.水利學(xué)報(bào)，2002（12）：7-10.

［4］褚君達(dá).河網(wǎng)對流輸移問題的求解及應(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.