基于MIKE21的南四湖上級湖水量水質(zhì)響應模擬研究

宮雪亮 ，孫 蓉，蘆昌興，孫秀玲

(山東大學 土建與水利學院，山東 濟南 250061)

南水北調(diào)工程是一項促進水資源優(yōu)化配置的重大戰(zhàn)略性基礎工程。其中，南水北調(diào)東線項目成功的關鍵是污染控制。東線工程第一階段的主要供水要求滿足三類地表水質(zhì)量要求。南四湖獨特的地理位置和東線工程的調(diào)水方式?jīng)Q定了當?shù)匦×饔虻乃|(zhì)對保證調(diào)水質(zhì)量起著決定性的作用。

關于水域水質(zhì)水量模擬的研究起源于20世紀60年代。1965年，Hansne等建立了首個關于淺水海域水位變化的二維函數(shù)模型；20世紀70年代，國外學者提出了基于水質(zhì)水量的湖泊水動力模擬方案，日本、美國等科研機構(gòu)分別對琵琶湖(日本)、伊利湖(美國)等湖泊水流流態(tài)進行了數(shù)值模擬，并獲得了相應研究成果[1,2]。1975年，Gallagher與Liggerttatal等在前人基礎上考慮了氣象要素(以風為例)，構(gòu)建了二維風生流數(shù)學模型[3]；1985年，Simon等提出了建立三維模型的開發(fā)評價方案；1993年，Kersmiita Zic等改進了該模型，并以日本琵琶湖為例模擬了該地的內(nèi)波變化以及湖環(huán)流結(jié)構(gòu)[4]。2004年，Chao等計算了密西西比流域三維淺水牛扼湖的水動力數(shù)值模型，其結(jié)果與實際情況吻合度較高[5]。

由于南四湖湖泊屬于淺水湖泊，所以本文模擬選取二維模擬，基于丹麥水科所開發(fā)的MIKE21軟件，該軟件主要應用于港口、河流、湖泊、河口及海岸水動力、泥沙及水質(zhì)的模擬研究[6]。MIKE21中的水動力學模型在計算過程中可以綜合考慮各種因素的影響，如地形、氣象條件、下墊面情況等，在一定程度上保證了模擬的精確度。

1 構(gòu)建模型

1.1 研究區(qū)概況

南四湖是位于山東省西南部，淮河流域北部的一處狹長水域，處于東經(jīng) 116°34′～117°21′，北緯 34°27′～35°20′之間，由南陽湖、獨山湖、昭陽湖、微山湖這 4 個湖泊相連組成。

南四湖屬于淺水大型淡水湖泊，平均水深1.46 m，最大水深2.76 m。1960年在湖腰處建成攔湖大壩，稱為二級壩，壩上興建了節(jié)制閘和船閘，二級壩將南四湖分為上級湖和下級湖2個部分。二級壩以北為上級湖，入上級湖河流共29條，集流面積26 934 km2，占總集流面積的88.4%;二級壩以南為下級湖，入下級湖河流共24條，集流面積3 519 km2，僅占總集流面積的11.6%。南四湖流域位于半濕潤季風氣候區(qū)，平均年降水量在750 mm左右，年平均風速為2.7～3.8 m/s。

1.2 確定計算區(qū)域

上級湖湖區(qū)邊界根據(jù)經(jīng)緯度轉(zhuǎn)化為WGS_1984_UTM_ZONE_50N通用投影坐標。根據(jù)投影坐標在MIKE軟件中得到上級湖的湖區(qū)邊界，并將此區(qū)域劃定為計算區(qū)域進行網(wǎng)格劃分。

本次湖泊模擬所選取的單元網(wǎng)格劃分方式為三角形網(wǎng)格，相較于四邊形網(wǎng)格劃分方式，劃分區(qū)域更加細致穩(wěn)定。本次研究上級湖共劃分18 188個網(wǎng)格，其中湖中南陽島為高出水面湖心島，不進行劃分。圖1為上級湖計算區(qū)域網(wǎng)格詳細劃分以及局部網(wǎng)格展示。

圖1 計算區(qū)域網(wǎng)格劃分、網(wǎng)格細致圖Fig.1 Calculate area meshing, grid detailing

1.3 確定模型邊界

根據(jù)南四湖上級湖入湖河流，選取匯水面積較大的河流，分為湖東湖西兩個區(qū)域。湖東區(qū)選取6條入湖河流，包括洸府河、泗河、白馬河、界河、北沙河、城郭河；湖西區(qū)選取7條入湖河流，包括洙水河、洙趙新河、萬福河、老萬福河、蔡河、惠河、東魚河。下級湖通過二級壩入上級湖。圖2為入南四湖流域劃分圖。

圖2 入南四湖流域劃分Fig.2 Watershed division into Nansi Lake

根據(jù)入湖河流流域劃分情況以及流域匯流面積情況將入湖河流編號整理，形成如表1，最終形成輸入模型的邊界條件。

1.4 模型初始數(shù)據(jù)及輸入數(shù)據(jù)

模型初始數(shù)據(jù)包括處理后的上級湖DEM數(shù)字高程數(shù)據(jù)、上級湖湖區(qū)初始水質(zhì)和湖區(qū)初始水位。模型輸入項包括各入湖河流的入流流量及污染物數(shù)據(jù)，出湖水位數(shù)據(jù)，湖內(nèi)風速風向以及降水蒸發(fā)量。

(1)南四湖上級湖DEM數(shù)字高程數(shù)據(jù)：高程點數(shù)據(jù)一共包含41 152個原始數(shù)據(jù)，將高程點數(shù)據(jù)經(jīng)緯度轉(zhuǎn)化為WGS_1984_UTM_ZONE_50N通用投影坐標，并對數(shù)據(jù)中不符合實際情況的個別高程點數(shù)據(jù)進行相對處理，得到真實可靠的上級湖基本地形數(shù)據(jù)。

表1 南四湖上級湖入湖河流及編號 km2Tab.1 Incoming rivers and numbers of the upstream section of Nansi Lake

(2)入湖流量數(shù)據(jù)：假設流域各河流入湖流量總量為100,150和 200 m3/s，由于很多河流未設立水文站，所以使用水文比擬法分配至各入湖河流，進行計算。

(3)風向風速：根據(jù)文獻資料設置模型上級湖內(nèi)風速風向，冬季盛行東北風和西北風，夏季盛行東南風和西南風，多年平均風速2.9 m/s，最大風速12.5 m/s。本次研究選取多年平均風速。

(4)湖面蒸發(fā)：南四湖湖區(qū)面積較大，蒸發(fā)量不作為主要變量進行考慮，根據(jù)文獻資料[7]設置模型上級湖內(nèi)降水蒸發(fā)資料，如表2所示。根據(jù)多年平均資料設置不同月份的蒸發(fā)數(shù)據(jù)的時間序列文件。

表2 南四湖多年平均月水面蒸發(fā)量表Tab.2 Annual average monthly water surface evaporation of Nansi Lake

(5)本次研究分析湖內(nèi)水質(zhì)污染物包括化學需氧量和總磷濃度含量。上級湖湖區(qū)初始水質(zhì)由于監(jiān)測數(shù)據(jù)時間尺度較大，加之近幾年湖內(nèi)水質(zhì)有明顯改善，大部分湖區(qū)可達三類水水質(zhì)，所以湖內(nèi)初始水質(zhì)定為三類水標準，即COD=20 mg/L，TP=0.05 mg/L。

1.5 模型參數(shù)選取

本次湖泊模擬分為水動力模擬與水質(zhì)模擬兩個部分，其中水動力模擬需要控制參數(shù)為湖區(qū)介質(zhì)糙率，水質(zhì)模擬需要控制參數(shù)為污染物水平擴散系數(shù)以及污染物降解系數(shù)。通過對以上參數(shù)進行率定，確定最終模型參數(shù)。

南四湖湖內(nèi)情況復雜，航道、深槽、湖草、蘆葦、魚池等遍布湖區(qū)，相互交錯，且疏密程度有所差異，所以湖區(qū)內(nèi)不同地區(qū)的糙率也需要進行概化分區(qū)處理。經(jīng)查閱資料[8]以及分析，最終確定采用4種情況的糙率為：航道(深槽)n=0.025，蘆葦n=0.619，湖草n=0.176，明湖n=0.084，,圖3為南四湖上級湖湖內(nèi)粗糙率mesh文件。

圖3 湖內(nèi)介質(zhì)糙率分布Fig.3 Distribution of medium roughness in the lake

參考文獻[9]得出降解系數(shù)在0.001～0.1 d-1之間，基于此數(shù)值范圍參考近年來對于南四湖污染物降解的研究報告，最終確定數(shù)值：COD降解系數(shù)為0.010 6 d-1，TP降解系數(shù)為0.003 1 d-1。

參考文獻[10,11]得出擴散系數(shù)在0.01～20 m2/s之間，基于此數(shù)值范圍參考近年來對于南四湖污染物擴散的研究報告，最終確定擴散系數(shù)為0.15 m2/s。

1.6 模型驗證

1.6.1 水位驗證

選取2010年數(shù)據(jù)進行模型水量水質(zhì)參數(shù)的驗證。模型給定計算邊界條件及參數(shù)選擇：梁濟運河作為自然入湖河流，上級湖通過二級壩出水渠向下級湖泄水，作為上級湖的出湖口，出湖水位以二級壩實測水位數(shù)據(jù)為準。

根據(jù)上級湖內(nèi)實測水位資料，選取上級湖內(nèi)測站南陽站作為模型驗證站點，將實測數(shù)據(jù)與模擬數(shù)據(jù)進行比較分析。模擬尺度為1 h，輸出數(shù)據(jù)尺度為1 d。

南陽站位于東經(jīng)116.54°北緯35.01°，處于湖內(nèi)中游位置，經(jīng)過對照發(fā)現(xiàn)模擬數(shù)值與實際水位相差不大，如圖4(a)所示，曲線趨勢基本相似，波峰與波峰相對，波谷與波谷相對，7-9月份呈現(xiàn)總體上升趨勢，10月份之后，非汛期時段內(nèi)呈現(xiàn)總體下降趨勢，與實際情況相符合。南陽站水位相對誤差圖見圖4(b)，整體誤差均在可以接受的誤差范圍之內(nèi)。實測全年水位平均值為34.305 m，模擬全年水位平均值為34.368 m，相差0.063 m，全年平均水位較為符合。說明建立的模型能夠較好地模擬湖區(qū)水動力的變化，可以用于各種情景分析。

圖4 南陽站水位Fig.4 Water level of Nanyang Station

1.6.2 水質(zhì)模型驗證

水質(zhì)模型驗證是在基于水動力模型模擬基礎上，對模型的水質(zhì)參數(shù)進行驗證。驗證站選取南陽站作為水質(zhì)模型驗證的參證站。選取2010年6月-8月進行模擬，前期為模型預熱，最終根據(jù)入湖流量以及水質(zhì)測量時間(7月10日左右)選取7月7日-13日的水質(zhì)數(shù)據(jù)與實測數(shù)據(jù)進行驗證。選取COD和TP兩種污染物作為水質(zhì)模擬污染物指標。下面就這兩種污染物指標，對模型進行驗證分析。

圖5為南陽站2010年7月7日-13日的COD和TP濃度模擬值與實測值對比圖。污染物濃度總體變化幅度不大，由于汛期有持續(xù)入湖徑流，所以湖內(nèi)污染物濃度有上升的趨勢，實測值在模擬值變化范圍之內(nèi)。由于總磷易聚集，較難降解的特性，導致在短時間內(nèi)上升速度較快。綜合兩種污染物的情況，水質(zhì)數(shù)據(jù)均能由模型很好地體現(xiàn)出來，所以認為模型有代表性并且可應用于后續(xù)的湖內(nèi)水質(zhì)模擬。

圖5 南陽站污染物濃度模擬值與實測值對比Fig.5 Comparison between simulated and measured pollutant concentration of Nanyang Station

2 水量水質(zhì)模擬

2.1 模擬方案設定

二維水質(zhì)模型是建立在水動力模型的基礎上進行模擬分析，本文旨在研究南四湖上級湖的水量水質(zhì)響應結(jié)果，所以方案設置以此為重點。南水北調(diào)調(diào)水由下級湖向上級湖經(jīng)二級壩入湖至梁濟運河出湖口調(diào)水出湖，本文僅對特殊干旱時期調(diào)水量在80 m3/s情況進行模擬分析。

湖內(nèi)初始水位設為興利水位34.5 m，風向為東北風(315°)，初始水質(zhì)為三類水。由于二級壩向上級湖調(diào)水已經(jīng)經(jīng)過下級湖的調(diào)蓄過程，所以假設調(diào)水水質(zhì)也已達到三類水，即COD=20 mg/L，TP=0.05 mg/L。

假設流域各河流入湖流量總量為100,150和200 m3/s，使用水文比擬法分配至各入湖河流。入湖徑流水質(zhì)選擇較差的Ⅴ類水，即COD=40 mg/L，TP=0.2 mg/L。為了模擬在及時關閘后，湖內(nèi)總體的污染物變化情況，模型運行時長設定為180 d，其中前20 d入湖河流有入流。

2.2 模擬結(jié)果

圖6為入湖流量100,150和200 m3/s時梁濟運河調(diào)水出口斷面污染物濃度變化圖?？梢钥闯觯簼\河出口斷面污染物濃度在開始是持續(xù)下降的，這是由于湖內(nèi)自凈使得湖內(nèi)本身污染物濃度下降；十幾天之后會有一個上升，是因為臨近梁濟運河出口斷面的4條河流的污染物逐步擴散至出口斷面，使得斷面污染物濃度有明顯上升趨勢；在30天前后達到峰值，之后便持續(xù)下降，在50～70 d又開始有小幅度抬升，是由于較遠入湖口攜帶的污染物經(jīng)過與湖水的充分混合，由湖內(nèi)水流運動流至出口斷面使出口斷面有一個較小的峰值，但由于這部分污染物已經(jīng)經(jīng)過湖水的充分混合，所以這個峰值會比初始峰值小。

圖7為不同情況下上級湖內(nèi)平均流速以及流向圖，在調(diào)水情況下，水流流向由南向北，湖內(nèi)流向為由二級壩調(diào)水入湖，以及各入湖口入湖并向湖內(nèi)擴散，最終經(jīng)由梁濟運河抽水繼續(xù)北調(diào)。流域入湖流量越大，水流運動越快。各方案不同時間水質(zhì)達到三類水面積占比情況見表3。圖8為120 d時，不同入湖流量下湖內(nèi)污染物濃度情況。

3 結(jié) 論

根據(jù)模擬結(jié)果,可以得到以下結(jié)論:

(1)圖6、7可以看出,盡管調(diào)水出口斷面污染物濃度變化過程相似,但流域入湖流量越大，出口斷面污染物濃度峰值出現(xiàn)越早,說明在入湖河流來水水質(zhì)情況一定(Ⅴ類水)的情況下，湖內(nèi)流速的快慢決定了湖內(nèi)污染物擴散運移速度。流域入湖流量越大，流速越快，污染物運移速度也就越快，污染物越容易擴散。

圖6 不同入湖流量下梁濟運河斷面污染物濃度變化Fig.6 Pollutant concentration change of Liangji Canal section under different inflow

圖7 不同入湖流量下湖內(nèi)流速分布及流向Fig.7 The velocity distribution and flow direction in the lake under different inflow

表3 各方案不同時間三類水面積占比Tab.3 Water area proportion of third type in different time of each scheme

圖8 120 d時不同入湖流量下湖內(nèi)污染物濃度Fig.8 Pollutant concentration in the lake under different inflow on day 120

(2)雖然入湖流量越大污染物越容易擴散,但來水量越大，也就代表污染物攜帶量越多,如表3所示,同一時間同一污染物水質(zhì)達到三類水面積占比,隨著入湖流量增大而減小。

(3)當流域入湖流量與入湖水質(zhì)一定時，不同污染物所表現(xiàn)出的情況也不盡相同。COD降解速率快，最終上級湖內(nèi)COD達到三類水水質(zhì)的面積也較多；TP降解速率較慢，水質(zhì)達到三類水的面積與COD相差較大。

(4)考慮地形因素，由圖8可以看出，最終湖內(nèi)污染物較為聚集的地方均為湖岸邊。由于湖岸邊水深較淺，水流流速較緩，污染物容易聚集，所以相較于湖內(nèi)或者是輸水航道線上的水來說，最終得到的水質(zhì)較差。在后續(xù)的湖內(nèi)水資源治理中，可以重點考慮湖岸帶的污染治理，會使湖內(nèi)整體水質(zhì)上升。

(5)入湖污染物濃度必須嚴格控制，否則一旦進入到湖內(nèi)，如果沒有大量流域入湖流量或調(diào)水流量的沖洗，污染物很難擴散降解，最終導致湖內(nèi)整體水質(zhì)下降。尤其是非調(diào)水期，由于調(diào)水期過后湖內(nèi)整體初始水位不高，所以需要入上級湖河流開閘放水以補充湖內(nèi)水量，此時沒有調(diào)水流量沖洗，就更需要對入湖徑流水質(zhì)進行嚴格控制，否則入湖徑流雖然可以加速湖內(nèi)污染物擴散運移，但自身所攜帶的大量污染物質(zhì)反而對湖內(nèi)水質(zhì)改善起到反作用。