馬天海，孫 娟，顏劍波

（1.南京大學 金陵學院，南京 210089；2.中國電建集團中南勘測設(shè)計研究院有限公司，長沙 410014）

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

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

1 研究區(qū)域概況

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

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

2 M IKE水動力學模型

2.1 模型簡介

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

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

2.2 水動力控制方程

2.2.1 一維水動力

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

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

2.2.2 二維水動力

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

連續(xù)性方程：

動量方程：

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

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

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

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

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

3.2 二維水動力模型地形的建立

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

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

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

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

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

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

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

圖5 一二維耦合計算示意圖

4 陽澄湖一二維水力模型的率定與驗證

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

4.1.1 計算時段及時間空間步長的選取

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

4.1.2 計算邊界條件

4.1.2.1 外邊界條件

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

4.1.2.2 內(nèi)邊界條件

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

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

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

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

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

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

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

4.2 河道水位驗證

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

表1 實測水位與計算水位誤差統(tǒng)計

圖6 各站計算水位和實測水位過程擬合情況

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

4.3 湖區(qū)水位驗證

利用2008年實測湖泊水位資料，對建立的陽澄湖二維水動力模型進行率定，選擇湖區(qū)11個測點進行水位計算值與實測值對比分析，率定與驗證站點如圖7。

圖7 陽澄湖湖體水文測點分布

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

表2 測點水位率定結(jié)果

續(xù)表2

5 結(jié)語

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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