馬天海，孫 娟，顏劍波

（1.南京大學(xué) 金陵學(xué)院，南京 210089；2.中國電建集團(tuán)中南勘測(cè)設(shè)計(jì)研究院有限公司，長(zhǎng)沙 410014）

隨著經(jīng)濟(jì)快速發(fā)展和城市化進(jìn)程的加快，陽澄湖的有機(jī)污染處于緩慢上升的態(tài)勢(shì)，富營(yíng)養(yǎng)化水平不斷提高[1-2]，湖體和飲用水源富營(yíng)養(yǎng)化趨勢(shì)尚未得到有效控制，離城市供水水源地的要求更是存在著一定的差距。通過稀釋污染物和加快水循環(huán)過程是控制湖泊富營(yíng)養(yǎng)化、保障飲水水源安全的有效途徑之一[3-5]。近年來，蘇州市政府非常重視陽澄湖水污染問題，2008年蘇州市政府出臺(tái)了《陽澄湖水污染防治工作計(jì)劃》，正式將保護(hù)陽澄湖水源水質(zhì)，保障飲用水源和戰(zhàn)略備用飲用水源安全，改善陽澄湖水環(huán)境質(zhì)量作為工作的重點(diǎn)，并規(guī)劃從長(zhǎng)江經(jīng)楊林塘、七浦塘和永昌涇引長(zhǎng)江水入湖改善陽澄湖湖體水質(zhì)。

由于陽澄湖湖泊水動(dòng)力受周邊河道條件影響的復(fù)雜性，僅通過建立湖泊水力模型而不考慮湖泊邊界的互相影響，將帶來湖流模擬的不準(zhǔn)確因素。鑒此，本文利用MIKE FLOOD，建立了陽澄湖及周邊河網(wǎng)的一、二維水力耦合模型，模擬了陽澄湖的水位變化過程，為陽澄湖的水量調(diào)度提供了參考，亦為其水質(zhì)模擬提供了前提條件。

1 研究區(qū)域概況

研究區(qū)域是以陽澄湖為調(diào)蓄中心的蘇州陽澄淀泖水文片區(qū)。陽澄淀泖區(qū)位于蘇州市中部、太湖流域的東北部，該區(qū)北以望虞河為界；南以太浦河為界；東以長(zhǎng)江、上海為界；西以太湖東岸線為界。面積4780km2，占蘇州市陸地面積的75%左右。區(qū)內(nèi)有蘇州、吳縣、吳江、昆山、太倉和常熟等6個(gè)大中城市。陽澄湖位于蘇州古城東北10km，湖跨蘇州市相城區(qū)（原吳縣）和昆山市，是陽澄地區(qū)防洪、排澇、引水、灌溉的調(diào)蓄湖泊，同時(shí)也是蘇州市區(qū)和昆山城區(qū)主要飲用水水源地。

陽澄湖南北長(zhǎng)13~17km，東西寬2~4km，總面積119km2，湖中2條天然土埂貫穿南北，將湖面分為東、中、西湖，三湖之間有眾多港汊相通，湖體水深1.7~5.0m。

2 M IKE水動(dòng)力學(xué)模型

2.1 模型簡(jiǎn)介

MIKE模型體系主要包括MIKE11，MIKE21和MIKE3，是丹麥DHI公司開發(fā)的產(chǎn)品。MIKE11用于一維水動(dòng)力學(xué)、水質(zhì)、富營(yíng)養(yǎng)化和泥沙輸移計(jì)算及洪水預(yù)報(bào)等，適用于河流、湖庫及灌渠等。MIKE21用于二維水動(dòng)力學(xué)、水質(zhì)、富營(yíng)養(yǎng)化、石油泄漏等計(jì)算，適用于河流、湖庫、河口及海灣等。MIKE3用于三維水動(dòng)力學(xué)和水質(zhì)模擬，適用于河流、湖庫、河口、海洋等。MIKE21模型是MIKE11的姐妹模型，在全世界廣泛應(yīng)用，是一個(gè)極優(yōu)秀的模型，用來模擬在水質(zhì)預(yù)測(cè)中垂向變化常被忽略的湖泊、河口、海岸地區(qū)[6-7]。

根據(jù)河湖串聯(lián)河網(wǎng)區(qū)水動(dòng)力及其水環(huán)境演變特征，選取MIKE系統(tǒng)模型的MIKE11，MIKE21，MIKE FLOOD，構(gòu)建河湖串聯(lián)河網(wǎng)區(qū)水動(dòng)力系統(tǒng)模型。

2.2 水動(dòng)力控制方程

2.2.1 一維水動(dòng)力

MIKE 11 HD模塊通過求解明渠流完全非線性St．Venant方程，可以對(duì)包含多種水工建筑物，如堰、箱涵、橋梁和人工控制閘站等河道、渠道水流進(jìn)行模擬[8]。平原河網(wǎng)地區(qū)河道水流運(yùn)動(dòng)可以用垂向積分的連續(xù)性方程和動(dòng)量方程（圣維南方程）來描述：

式中Q為流量（m3/s）；A為斷面面積（m2）；q為側(cè)向入流流量（m3/s）；h為基準(zhǔn)面以上水深（m）；C為柯西阻力系數(shù)（m1/2/s）；R為水力半徑（m）；α為動(dòng)量分布系數(shù)；x計(jì)算空間坐標(biāo)（m）；t計(jì)算時(shí)間坐標(biāo)（s）；g重力加速度（m/s2）。

2.2.2 二維水動(dòng)力

MIKE 21水流模擬基于的控制方程是不可壓流三維雷諾Navier-Stokes平均方程沿水深積分的連續(xù)方程和動(dòng)量方程[9]，在笛卡爾坐標(biāo)系中可用如下方程表示：

連續(xù)性方程：

動(dòng)量方程：

式中t為時(shí)間（s）；x，y分別為笛卡爾坐標(biāo)（m）；h為總水深（m）；和分別為水深平均的值（m/s）；；S為點(diǎn)源的排放量（kg）；ρ為水的密度（kg/m3）；f=2Ωsinφ，表示Coriolis因子（Ω為地球自轉(zhuǎn)的角速度，φ為地理緯度）；g為重力加速度（m/s2）；η為水位（m）；sxx，sxv，和svv為radiation應(yīng)力張量（Pa）；τsx，τsy分別為水面風(fēng)應(yīng)力張量 （Pa）；τbx，τby分別為河床床面應(yīng)力張量（Pa）；pa為大氣壓（Pa）；ρ0為水的相對(duì)密度（kg/m3）；（us，vs）為外界排放到環(huán)境水體的速率（m/s）；橫向應(yīng)力Txx，Txy，Tyy包括黏滯阻力、紊流摩擦阻力和差動(dòng)平流摩擦阻力（Pa），用垂向流速平均的渦粘方程來計(jì)算：

3 陽澄湖一二維水力模型建立

3.1 一維河網(wǎng)概化

平原河網(wǎng)地區(qū)大小河道及湖泊數(shù)量眾多，需要先對(duì)計(jì)算區(qū)域的河道湖泊進(jìn)行概化，概化的范圍大體在陽澄淀泖區(qū)，西以太湖為界，北至望虞河，東至長(zhǎng)江，南至吳淞江—青陽港—瀏河邊界。經(jīng)概化的河網(wǎng)應(yīng)能夠反應(yīng)本地區(qū)天然河網(wǎng)的水動(dòng)力情況，根據(jù)過水能力等效原則及調(diào)蓄容積不變?cè)瓌t，對(duì)周邊河網(wǎng)進(jìn)行概化后，得到概化河道105條，河流節(jié)點(diǎn)464個(gè)，邊界節(jié)點(diǎn)22個(gè)，閘門35個(gè)。計(jì)算區(qū)域內(nèi)湖泊除陽澄湖做二維概化外，其他湖泊均作為零維調(diào)蓄節(jié)點(diǎn)來處理。計(jì)算區(qū)域一維河網(wǎng)概化如圖1。

圖1 蘇州市河網(wǎng)概化示意圖

3.2 二維水動(dòng)力模型地形的建立

根據(jù)蘇州市水文水資源勘測(cè)局測(cè)量的1∶25000 CAD湖形圖，陽澄湖水下地形導(dǎo)出x，y坐標(biāo)數(shù)據(jù)（北京54坐標(biāo)，吳淞高程系），生成閉合水陸邊界線（Land.xyz）和水深散點(diǎn)數(shù)據(jù)（Water.xyz），然后倒入網(wǎng)格生成器（mesh generator）生成mesh格式的湖泊網(wǎng)格。水下地形和網(wǎng)格分別如圖2和圖3。

圖2 陽澄湖湖區(qū)水下地形圖

圖3 陽澄湖湖區(qū)網(wǎng)格剖分示意圖

3.3 一、二維水力耦合條件處理

MIKE11對(duì)陽澄湖周邊河網(wǎng)的模擬采用基于Q，h交替網(wǎng)格的6點(diǎn)Abbott&Ionescu隱式有限差分法[9]，MIKE21對(duì)陽澄湖的水力模擬采用三角網(wǎng)格的有限體積法。根據(jù)連接處動(dòng)量守恒原則，利用MIKE FLOOD進(jìn)行一、二維模型的耦合，其連接如圖4。

圖4 一、二耦合連接處處理方式

一二維耦合根據(jù)“水位～流量”銜接關(guān)系，一維水流流向二維區(qū)域時(shí)，由一維計(jì)算出連接河道末端（第一個(gè)Q點(diǎn)）的流量，并作為源項(xiàng)提供給二維連接網(wǎng)格單元；二維水流流向一維河道時(shí)，由二維模型計(jì)算出連接網(wǎng)格單元的水位，提供給一維河道連接節(jié)點(diǎn)，并作為該節(jié)點(diǎn)的水位邊界。在連接處的河道端點(diǎn)需要給定一個(gè)虛擬的水位邊界，作為耦合模型的啟動(dòng)條件，該虛擬水位不影響后續(xù)計(jì)算。一、二維耦合計(jì)算如圖5。

圖5 一二維耦合計(jì)算示意圖

4 陽澄湖一二維水力模型的率定與驗(yàn)證

4.1 模型邊界條件設(shè)置

4.1.1 計(jì)算時(shí)段及時(shí)間空間步長(zhǎng)的選取

根據(jù)現(xiàn)有資料情況，耦合模型率定及驗(yàn)證的時(shí)間段選2002年1月1日至12月31日?？紤]到和MIKE21的耦合，水力耦合模型時(shí)間步長(zhǎng)取60s?？臻g步長(zhǎng)由模型根據(jù)河長(zhǎng)自行給定，步長(zhǎng)范圍50～1500m。

4.1.2 計(jì)算邊界條件

4.1.2.1 外邊界條件

長(zhǎng)江沿線潮位過程；望虞河、太湖、瀏河—吳淞江沿線水位過程；整個(gè)陽澄淀泖區(qū)的氣象條件 （降雨、蒸發(fā)）及風(fēng)場(chǎng)條件。其中風(fēng)場(chǎng)條件參考相近太湖地區(qū)風(fēng)場(chǎng)[10]。

4.1.2.2 內(nèi)邊界條件

（1）望虞河?xùn)|岸控制線：當(dāng)湘城水位低于3.5m時(shí)，望虞河?xùn)|岸控制線除琳橋港閘控制50m3/s外，其余所有口門全線控制；當(dāng)湘城水位超過3.5m時(shí)，琳橋港閘關(guān)閉。

（2）陽澄區(qū)沿長(zhǎng)江水閘：瀏河、楊林、七浦、白茆、滸浦閘以陽澄湖湘城站為代表，當(dāng)湘城水位低于3.0m時(shí)全力引水；當(dāng)湘城水位在3.0～3.2m時(shí)不引不排；當(dāng)湘城水位超過3.2m時(shí)排水。沿江其他小型水閘按最高水位3.5m、最低水位3.0m控制運(yùn)行。

（3）沿太湖水閘：在楓橋站水位高于4.2m時(shí)關(guān)閉，低于4.2m時(shí)敞開。

4.1.3 河道水力模型參數(shù)

計(jì)算區(qū)域河道糙率取值0.02～0.025。

4.1.4 湖泊水力模型參數(shù)

曼寧系數(shù)是影響水流計(jì)算的關(guān)鍵參數(shù)。研究淺水湖泊湖底糙率一般在0.02～0.025之間，本文通過率定確定曼寧數(shù)取46～48m1/3/s。水平渦粘系數(shù)計(jì)算參數(shù)Cs參考相關(guān)文獻(xiàn)[11-12]進(jìn)行計(jì)算，其他參數(shù)設(shè)為默認(rèn)值。

4.2 河道水位驗(yàn)證

利用2002年水文實(shí)測(cè)資料進(jìn)行耦合水力模型的驗(yàn)證，常熟、昆山（二）站、直塘、湘城和巴城的計(jì)算水位和實(shí)測(cè)水位過程擬合情況如表1及圖6。

表1 實(shí)測(cè)水位與計(jì)算水位誤差統(tǒng)計(jì)

圖6 各站計(jì)算水位和實(shí)測(cè)水位過程擬合情況

由河道水位率定結(jié)果可知，內(nèi)部水文站實(shí)測(cè)水位與計(jì)算水位相對(duì)誤差控制范圍在-7.0%~9.0%，水位平均誤差在0.04～0.08之間，評(píng)價(jià)相對(duì)誤差在3%以下，河道水位計(jì)算具有一定的模擬精度。

4.3 湖區(qū)水位驗(yàn)證

利用2008年實(shí)測(cè)湖泊水位資料，對(duì)建立的陽澄湖二維水動(dòng)力模型進(jìn)行率定，選擇湖區(qū)11個(gè)測(cè)點(diǎn)進(jìn)行水位計(jì)算值與實(shí)測(cè)值對(duì)比分析，率定與驗(yàn)證站點(diǎn)如圖7。

圖7 陽澄湖湖體水文測(cè)點(diǎn)分布

陽澄湖測(cè)定水位計(jì)算值與實(shí)測(cè)值對(duì)比情況如表2，湖區(qū)水位計(jì)算值與實(shí)測(cè)值誤差范圍在-0.28%~0.28%。因此，一、二維河網(wǎng)湖泊耦合水動(dòng)力模型的，耦合模型應(yīng)用于陽澄湖水動(dòng)力模擬是可行的。

表2 測(cè)點(diǎn)水位率定結(jié)果

續(xù)表2

5 結(jié)語

（1）MIKE模型是個(gè)模擬內(nèi)容豐富的模型，其MIKE FLOOD可方便地將一維河網(wǎng)和二維湖泊水力耦合起來計(jì)算，然其自帶的產(chǎn)匯流模塊不適合平原地區(qū)特殊的圩區(qū)匯流方式。本文建立了平原圩區(qū)的產(chǎn)匯流模型替換改進(jìn)了MIKE自身的產(chǎn)匯流模型，并利用MIKE模型自帶的分布源接口，將產(chǎn)匯流模型計(jì)算的結(jié)果作為MIKE11河道旁側(cè)入流條件，實(shí)現(xiàn)產(chǎn)匯流和河網(wǎng)水力的耦合計(jì)算。

（2）MIKE21水力模型的建立應(yīng)用了效率和精度較高的三角網(wǎng)格有限體積法，較好地?cái)M合了陽澄湖復(fù)雜的湖形邊界，建立陽澄湖及周邊河網(wǎng)一、二維水動(dòng)力模型。模型驗(yàn)證結(jié)果表明，河道水位計(jì)算值與實(shí)測(cè)值誤差范圍在-7.0%~9.0%，湖區(qū)水位計(jì)算值與實(shí)測(cè)值誤差范圍在-0.28%~0.28%。因此，一、二維河網(wǎng)湖泊耦合水動(dòng)力模型的具有一定的模擬精度，耦合模型應(yīng)用于陽澄湖水動(dòng)力模擬是可行的。

（3）本文中水量耦合模型的構(gòu)建概念清楚，但計(jì)算區(qū)域內(nèi)沒有實(shí)測(cè)流量站點(diǎn)，水量模型僅作水位率定，因此模型精度可能受到一定影響。

[1]桂智凡，薛濱，姚書，魏文佳.陽澄湖水質(zhì)現(xiàn)狀及原因探討[J].地理科學(xué)，2011，31（12），1487-1491.

[2]趙凌宇，翁建中.陽澄湖浮游藻類現(xiàn)狀調(diào)查及水質(zhì)評(píng)價(jià)[J].環(huán)境監(jiān)測(cè)管理與技術(shù)，2013，25（1）：27-29.

[3]田向榮，馬巍，廖文根，等.調(diào)水對(duì)梅梁湖、五里湖水環(huán)境影響研究[J].人民長(zhǎng)江，2007（2）：69-72.

[4]郜會(huì)彩，李義天，何用，等.改善漢陽湖群水環(huán)境的調(diào)水方案研究[J].水資源保護(hù)，2006（22）：41-44.

[5]馬巍，廖文根，李錦秀，等.引水調(diào)控改善太湖湖灣水環(huán)境及其效果預(yù)測(cè)[J].長(zhǎng)江流域資源與環(huán)境，2007（16）：52-56.

[6]程海云，黃艷.丹麥水力研究所河流數(shù)學(xué)模擬系統(tǒng)[J].水利水電快報(bào)，1996，17（19）：24-27.

[7]許婷.MIKE21HD 計(jì)算原理及應(yīng)用實(shí)例[J].港工技術(shù)，2010，47（5）：1-5.

[8]盧士強(qiáng).平原感潮河網(wǎng)水環(huán)境數(shù)學(xué)模型研究與應(yīng)用[D].上海：同濟(jì)大學(xué)，2003.

[9]Danish Hydraulic Institute（DHI）.MIKE 21&MIKE3 FLOW MODEL FM Hydrodynamic and Transport Module Scientific Documentation[R]．2005．

[10]洪曉瑜.太湖水環(huán)境數(shù)學(xué)模型建立及排污總量控制研究[D].南京：河海大學(xué)，2005.

[11]田向榮，馬巍，廖文根.調(diào)水對(duì)梅梁湖、五里湖水環(huán)境影響研究[J].人民長(zhǎng)江，2007（2）：69-72.

[12]龔春生，姚琪，趙棣華.玄武湖風(fēng)生流數(shù)值模擬研究[J].河海大學(xué)學(xué)報(bào)（自然科學(xué)版），2005（1）：72-75.