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

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

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

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

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

1 構(gòu)建模型

1.1 研究區(qū)概況

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

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

1.2 確定計(jì)算區(qū)域

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

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

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

1.3 確定模型邊界

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

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

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

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

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

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

表1 南四湖上級(jí)湖入湖河流及編號(hào) km2Tab.1 Incoming rivers and numbers of the upstream section of Nansi Lake

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

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

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

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

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

1.5 模型參數(shù)選取

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

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

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

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

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

1.6 模型驗(yàn)證

1.6.1 水位驗(yàn)證

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

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

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

圖4 南陽(yáng)站水位Fig.4 Water level of Nanyang Station

1.6.2 水質(zhì)模型驗(yàn)證

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

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

圖5 南陽(yáng)站污染物濃度模擬值與實(shí)測(cè)值對(duì)比Fig.5 Comparison between simulated and measured pollutant concentration of Nanyang Station

2 水量水質(zhì)模擬

2.1 模擬方案設(shè)定

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

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

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

2.2 模擬結(jié)果

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

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

3 結(jié) 論

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

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

圖6 不同入湖流量下梁濟(jì)運(yùn)河斷面污染物濃度變化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 各方案不同時(shí)間三類水面積占比Tab.3 Water area proportion of third type in different time of each scheme

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

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

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

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

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