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

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

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

南水北調工程是一項促進水資源優化配置的重大戰略性基礎工程。其中，南水北調東線項目成功的關鍵是污染控制。東線工程第一階段的主要供水要求滿足三類地表水質量要求。南四湖獨特的地理位置和東線工程的調水方式決定了當地小流域的水質對保證調水質量起著決定性的作用。

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

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

1 構建模型

1.1 研究區概況

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

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

1.2 確定計算區域

上級湖湖區邊界根據經緯度轉化為WGS_1984_UTM_ZONE_50N通用投影坐標。根據投影坐標在MIKE軟件中得到上級湖的湖區邊界，并將此區域劃定為計算區域進行網格劃分。

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

圖1 計算區域網格劃分、網格細致圖Fig.1 Calculate area meshing, grid detailing

1.3 確定模型邊界

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

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

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

1.4 模型初始數據及輸入數據

模型初始數據包括處理后的上級湖DEM數字高程數據、上級湖湖區初始水質和湖區初始水位。模型輸入項包括各入湖河流的入流流量及污染物數據，出湖水位數據，湖內風速風向以及降水蒸發量。

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

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

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

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

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

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

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

1.5 模型參數選取

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

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

圖3 湖內介質糙率分布Fig.3 Distribution of medium roughness in the lake

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

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

1.6 模型驗證

1.6.1 水位驗證

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

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

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

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

1.6.2 水質模型驗證

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

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

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

2 水量水質模擬

2.1 模擬方案設定

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

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

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

2.2 模擬結果

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

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

3 結 論

根據模擬結果,可以得到以下結論:

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

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

圖7 不同入湖流量下湖內流速分布及流向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時不同入湖流量下湖內污染物濃度Fig.8 Pollutant concentration in the lake under different inflow on day 120

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

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

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

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