馮忠倫，曹 濱，王 維，刁維杰，林洪孝，王 剛

(山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東 泰安 271018)

0 引 言

南水北調(diào)工程的實(shí)施，可對(duì)沿線(xiàn)地下水進(jìn)行有效補(bǔ)給，這將從根本上改變黃淮海平原干旱缺水及生態(tài)環(huán)境惡化的局面[1,2]，沿線(xiàn)輸水河道是河流系統(tǒng)與地下水系統(tǒng)之間發(fā)生水力聯(lián)系的橋梁，也是區(qū)域水文循環(huán)的重要環(huán)節(jié)[3-4]。輸水可通過(guò)河道滲漏補(bǔ)給地下水[5]，引起地下水位的變化。目前，對(duì)南水北調(diào)通水對(duì)該區(qū)域地下水位的研究主要針對(duì)靠近南四湖的區(qū)域[1]，尚未對(duì)整個(gè)輸水沿線(xiàn)的地下水位進(jìn)行研究，且建立的模型未經(jīng)過(guò)輸水期實(shí)測(cè)資料驗(yàn)證。本文選取南水北調(diào)東線(xiàn)工程梁濟(jì)運(yùn)河區(qū)進(jìn)行研究，該運(yùn)河工程位于山東省西南部，作為輸水干線(xiàn)，連接江、淮、黃、海四大流域，賦有長(zhǎng)期輸水任務(wù)[6]。梁濟(jì)運(yùn)河北起梁山縣路那里村，南至濟(jì)寧市郊區(qū)李集村入南陽(yáng)湖，全長(zhǎng)87.8 km，流域面積達(dá)3 306 km2，在河道上建有三座泵站，自南向北分別為長(zhǎng)溝、鄧樓、八里灣泵站，見(jiàn)圖1。該輸水河道兩岸邊坡均采用混凝土襯砌，河底保留原有地貌，未做防滲處理。本文應(yīng)用Visual MODFLOW建立了梁濟(jì)運(yùn)河區(qū)地下水流場(chǎng)模型，對(duì)南水北調(diào)輸水后未來(lái)一段時(shí)間內(nèi)地下水位進(jìn)行預(yù)測(cè)，為有效管理輸水和新情景下的水資源管理提供了有力的技術(shù)依據(jù)。

1 研究區(qū)概況

根據(jù)梁濟(jì)運(yùn)河區(qū)的水文地質(zhì)狀況和南水北調(diào)輸水前后地下水流場(chǎng)的變化情況，研判其輸水影響范圍，以梁濟(jì)運(yùn)河兩側(cè)10 km處作為研究區(qū)的東、西部邊界，北部以梁山縣邊境為界，南部以梁濟(jì)運(yùn)河流域區(qū)邊界為界，擬定地下水流場(chǎng)模擬范圍，面積1 704.32 km2，區(qū)內(nèi)共有25眼地下水觀測(cè)井，見(jiàn)圖1。

2 研究區(qū)水文地質(zhì)概念模型

2.1 地下水系統(tǒng)分析

梁濟(jì)運(yùn)河沿岸地下水均為第四紀(jì)孔隙潛水，主要儲(chǔ)藏于壤土、黏土和亞砂土層中，地表透水性差，地下水運(yùn)動(dòng)滯緩，地下水主要補(bǔ)給來(lái)源為大氣降水，其次是引黃輸水側(cè)滲水及溝渠滲水。根據(jù)地質(zhì)鉆孔資料及含水層埋藏條件，將研究區(qū)內(nèi)地層概化為一層，即潛水含水層。研究區(qū)水文地質(zhì)分區(qū)主要包括湖西黃泛平原區(qū)、濱湖平原區(qū)和湖東山前井灌平原區(qū)。湖西黃泛平原區(qū)潛水含水層厚度約為20～40 m，以粉砂為主，個(gè)別有細(xì)砂，西部較粗，近湖較細(xì)。濱湖平原區(qū)含水層巖性以粉細(xì)砂為主，局部黏土裂隙發(fā)育及黏土夾礓石。湖東山前井灌平原區(qū)第四系含水砂層較厚，顆粒較粗。由東北向西南第四系逐漸加厚，含水砂層層次增多，厚度增大，顆粒變細(xì)，含水砂層主要由細(xì)砂、中砂組成。綜上所述，研究區(qū)內(nèi)含水層多孔介質(zhì)顆粒粗細(xì)不一，表現(xiàn)出非均質(zhì)性。故將研究區(qū)地下水含水系統(tǒng)概化為非均質(zhì)、各向異性介質(zhì)、三維非穩(wěn)定地下水流系統(tǒng)。

2.2 邊界條件處理

因研究區(qū)西北部邊界與地下水等水位線(xiàn)近乎垂直，與外界幾乎無(wú)水力聯(lián)系，因此該區(qū)可視為零通量邊界。而在其余邊界處，由于無(wú)明顯自然邊界，且與研究區(qū)內(nèi)存在較強(qiáng)水力聯(lián)系，故將其確定為第二類(lèi)流量邊界。在垂向上，把潛水含水層自由水面作為模型上邊界，第四系松散堆積物下邊界為模型的底部邊界。通過(guò)上邊界，地下水流系統(tǒng)與外界發(fā)生垂向水量交換，如接受大氣降水及灌溉回歸補(bǔ)給、潛水蒸發(fā)及農(nóng)業(yè)灌溉開(kāi)采等[7]。

3 地下水流數(shù)值模擬模型

根據(jù)上述水文地質(zhì)概念模型，研究區(qū)地下水流的數(shù)學(xué)模型可表達(dá)如下[8-10]。

(1)

式中：Ω1為潛水滲流域；Γ1、Γ2為滲流域的第一類(lèi)、第二類(lèi)邊界；h為潛水位高程，m；h0為潛水含水層初始水位高程，m；b為含水層地板高程，m；h1為Γ2邊界上水頭分布，m；q為含水層二類(lèi)邊界單位面積流量，m3/(d·m2)；kn為邊界面法線(xiàn)方向的滲透系數(shù)，m/s；Kx、Ky、Kz分別為x,y,z方向上的滲透系數(shù)，m/s；μ為潛水含水層給水度；qs為源匯項(xiàng)，mm/d。

3.1 地下水流模型的求解

將研究區(qū)平面剖分成75列，20行，共計(jì)1500個(gè)單元格。求解過(guò)程中，首先根據(jù)研究區(qū)水文地質(zhì)條件及已有的基于抽水試驗(yàn)得出的參數(shù)值，給出模型中各區(qū)的參數(shù)初值。

3.2 源匯項(xiàng)的確定

研究區(qū)地下水補(bǔ)給量主要包括降水入滲、河道入滲、側(cè)向徑流、灌溉回歸補(bǔ)給；排泄項(xiàng)主要有潛水蒸發(fā)、農(nóng)業(yè)灌溉開(kāi)采、側(cè)向徑流排泄等。模型中，將降水入滲、灌溉回歸以及農(nóng)業(yè)灌溉開(kāi)采三者按各自分區(qū)進(jìn)行疊加處理，流入為正，流出為負(fù)。側(cè)向流入和流出通過(guò)邊界條件控制。梁濟(jì)運(yùn)河則按河流模塊分段進(jìn)行賦值河水位、河底高程、河寬、河床厚度、滲透系數(shù)等參數(shù)。

3.3 模型的識(shí)別與驗(yàn)證

選取2007年1月1日至2010年12月31日作為模型識(shí)別期，將邊界條件、初始條件以及源匯項(xiàng)輸入模型中運(yùn)行，用該時(shí)段25眼觀測(cè)井的實(shí)際觀測(cè)水位對(duì)模型進(jìn)行識(shí)別、調(diào)參。通過(guò)擬合期的地下水位，識(shí)別水文地質(zhì)參數(shù)、邊界條件及均衡項(xiàng)，使建立的模型更加符合實(shí)際水文地質(zhì)條件。經(jīng)多次運(yùn)行求解，所得地下水流場(chǎng)變化基本一致(見(jiàn)圖2)，說(shuō)明所建立的模型基本可行。識(shí)別后的水文地質(zhì)參數(shù)值及分區(qū)分別見(jiàn)表1、圖3。

為進(jìn)一步驗(yàn)證模型建立的可靠性，用4次輸水期地下水位觀測(cè)資料對(duì)數(shù)值模型進(jìn)行檢驗(yàn)，其中初始條件和邊界條件不變，其余均采用實(shí)測(cè)資料進(jìn)行輸入。通過(guò)對(duì)模型運(yùn)算結(jié)果誤差參數(shù)統(tǒng)計(jì)(見(jiàn)表2)可見(jiàn)，各輸水期內(nèi)模型運(yùn)算結(jié)果均滿(mǎn)足模型精度要求，且歷次輸水中地下水位觀測(cè)值和計(jì)算值的絕對(duì)誤差<0.5 m的觀測(cè)井均占總數(shù)的70%以上[11]。通過(guò)繪制歷次輸水末期流場(chǎng)擬合圖(見(jiàn)圖4)，可見(jiàn)歷次輸水末期流場(chǎng)擬合程度均較好。表明模型擬合效果較好，說(shuō)明含水層結(jié)構(gòu)、邊界條件的概化及水文地質(zhì)參數(shù)的選取是合理的，該模型能夠真實(shí)反映當(dāng)?shù)氐乃牡刭|(zhì)條件，可進(jìn)行下一步分析工作。

圖2 識(shí)別期末期地下水流場(chǎng)擬合圖(單位：m)Fig.2 Chart of groundwater flow field at the end of the recognition period

參數(shù)分區(qū)滲透系數(shù)K/(m·s-1)給水度μ降水入滲補(bǔ)給系數(shù)αⅠ2.4×10-50.330.24Ⅱ5.0×10-30.400.28Ⅲ8.0×10-30.400.22Ⅳ5.0×10-40.460.28Ⅴ7.0×10-30.330.28Ⅵ1.9×10-40.440.24

圖3 參數(shù)分區(qū)圖(單位：m)Fig.3 Map of parameter zones

表2 模型驗(yàn)證期誤差參數(shù)統(tǒng)計(jì)表Tab.2 Statistical table of error parameters in themodel validation phase

注：模型精度達(dá)到要求標(biāo)準(zhǔn)為RMS<1 m，標(biāo)準(zhǔn)RMS<10%，相關(guān)系數(shù)越趨近于1越好[12]。

圖4 歷次輸水末期流場(chǎng)擬合圖(單位：m)Fig.4 Fittings of flow field in the last stage of water delivery

4 研究區(qū)地下水位預(yù)測(cè)

4.1 建立預(yù)測(cè)模型

本研究選取2015年12月31日研究區(qū)的地下水流場(chǎng)作為水位預(yù)測(cè)的初始流場(chǎng)，將梁濟(jì)運(yùn)河河道水位按設(shè)計(jì)水位相應(yīng)提高，降水、蒸發(fā)以及農(nóng)業(yè)灌溉開(kāi)采均采用多年平均數(shù)據(jù)，運(yùn)用經(jīng)過(guò)經(jīng)識(shí)別和驗(yàn)證過(guò)的模型預(yù)測(cè)未來(lái)5年內(nèi)研究區(qū)內(nèi)的地下水位變化情況。

4.2 預(yù)測(cè)結(jié)果分析

(1)地下水流場(chǎng)和水位變化。對(duì)比圖5、圖6所繪制的2015年和2020年末的地下水位流場(chǎng)圖發(fā)現(xiàn)，研究區(qū)內(nèi)地下水流場(chǎng)變化趨勢(shì)較小，總體流向仍由東北、東部向西南方向匯集。而2015年末研究區(qū)地下水位等值線(xiàn)相對(duì)密集，2020年末區(qū)域內(nèi)地下水位等值線(xiàn)相對(duì)稀疏，說(shuō)明2020年末地下水水力坡度變緩，水位得到整體抬升。由圖7可知，水位抬升主要集中在鄧樓～長(zhǎng)溝泵站區(qū)間內(nèi)靠近長(zhǎng)溝泵站的區(qū)域，最大抬升高度達(dá)6.0 m，在八里灣～鄧樓泵站區(qū)間內(nèi)水位抬升約為0.5～1.5 m，長(zhǎng)溝泵站以南至湖口區(qū)間內(nèi)水位抬升約為1.5～2.0 m，而在研究區(qū)西北部區(qū)域，由于遠(yuǎn)離輸水主河道，水位基本無(wú)變化。

圖5 地下水位初始流場(chǎng)(單位：m)Fig.5 Groundwater initial flow field

圖6 地下水位預(yù)測(cè)流場(chǎng)(單位：m)Fig.6 Groundwater predicted flow field

(2)地下水位年內(nèi)變化分析。選取M11801710、M119091A、M126047A、M1260920共4眼代表井繪制其2015年和2020年內(nèi)地下水位變化動(dòng)態(tài)，見(jiàn)圖8。由圖可知，2020年4眼代表井所呈現(xiàn)出的地下水位波動(dòng)狀態(tài)仍然是先降后升的趨勢(shì)，且最低水位仍出現(xiàn)在7月份左右，與2015年較為一致。但在2020年地下水的波動(dòng)狀態(tài)明顯要比2015年趨于穩(wěn)定，變幅也明顯減小。在2015年地下水年變幅在1.5～2.0 m之間，而在2020年地下水年變幅在0.5～1.0 m之間。說(shuō)明南水北調(diào)輸水后，由于河道長(zhǎng)時(shí)段高水位輸水，對(duì)區(qū)域地下水系統(tǒng)有實(shí)時(shí)補(bǔ)給，從而影響了地下水的波動(dòng)狀態(tài)。

圖7 地下水抬升等值線(xiàn)圖(單位：m)Fig.7 Contour map of groundwater table rising

(3)補(bǔ)給量分析。在模型中針對(duì)梁濟(jì)運(yùn)河輸水渠道沿線(xiàn)設(shè)置水均衡區(qū)，統(tǒng)計(jì)2020年河道滲漏量計(jì)算結(jié)果，見(jiàn)圖9。由圖可知，2020年研究區(qū)地下水系統(tǒng)得到河道輸水滲漏補(bǔ)給量為3 592.3 萬(wàn)m3，平均9.8萬(wàn)m3/d，且滲漏補(bǔ)給主要集中在7-9月份，占總滲漏量的77%。

5 結(jié) 論

(1)南水北調(diào)輸水后，到2020年，預(yù)測(cè)梁濟(jì)運(yùn)河區(qū)地下水流場(chǎng)變化趨勢(shì)較小，總體流向仍由東北、東部向西南方向匯集。

(2)2020年末地下水水力坡度變緩，水位得到整體抬升。水位抬升主要集中在鄧樓～長(zhǎng)溝泵站區(qū)間內(nèi)靠近長(zhǎng)溝泵站的區(qū)域，最大抬升高度達(dá)6.0 m，在八里灣～鄧樓泵站區(qū)間內(nèi)水位抬升約為0.5～1.5 m，長(zhǎng)溝泵站以南至湖口區(qū)間內(nèi)水位抬升約為1.5～2.0 m，而在研究區(qū)西北部區(qū)域，由于遠(yuǎn)離輸水主河道，水位基本無(wú)變化。

圖8 代表井地下水位年內(nèi)變化動(dòng)態(tài)Fig.8 Groundwater table dynamic change of representative wells in year

圖9 2020年河道逐月滲漏量Fig.9 River leakage by month in 2020

(3)預(yù)測(cè)得出2020年區(qū)域年內(nèi)地下水位波動(dòng)狀態(tài)仍然是先降后升的趨勢(shì)，且最低水位仍出現(xiàn)在7月份左右，與2015年較為一致。

(4)2020年地下水的波動(dòng)狀態(tài)明顯要比2015年趨于穩(wěn)定，變幅也明顯減小。在2015年地下水年變幅在1.5～2.0 m之間，而在2020年地下水年變幅在0.5～1.0 m之間。

(5)2020年河道輸水滲漏補(bǔ)給量為3 592.3 萬(wàn)m3，平均9.8 萬(wàn)m3/d，且有77%的滲漏量都集中在7-9月份。

□

[1] 祁元昊, 張 鑫, 焦 瑩,等. 南水北調(diào)供水前后京杭運(yùn)河梁濟(jì)段沿岸地下水埋深變化預(yù)測(cè)研究[J]. 水利與建筑工程學(xué)報(bào), 2012,10(3):145-149.

[2] 方 生, 陳秀玲. 關(guān)于黃淮海東部平原南水北調(diào)受水區(qū)水資源優(yōu)化配置與生態(tài)環(huán)境綜合治理的建議[J]. 南水北調(diào)與水利科技, 2004,24(3):131-133.

[3] Rejani R, Jha M K, Panda S N, et al. Simulation modeling for efficient groundwater management in balasore coastal basin, india[J]. Water Resources Management, 2007,22(1):23-50.

[4] Thomas T, Jaiswal R K, Galkate R, et al. Development of a rainfall-recharge relationship for a fractured basaltic aquifer in central india[J]. Water Resources Management, 2009,23(15):3 101-3 119.

[5] 席海洋, 馮 起, 司建華,等. 黑河下游額濟(jì)納三角洲河道滲漏對(duì)地下水補(bǔ)給研究綜述[J]. 冰川凍土, 2012,34(5):1 241-1 247.

[6] 楊安邦, 朱順初, 周家貴,等. 南水北調(diào)東線(xiàn)一期工程對(duì)梁濟(jì)運(yùn)河段地下水影響研究[J]. 南水北調(diào)與水利科技, 2009,7(6):245-249.

[7] 王剛, 周啟友, 魏國(guó)孝,等. 酒泉盆地地下水系統(tǒng)數(shù)值模擬與預(yù)測(cè)[J]. 工程勘察, 2009,37(2):37-41.

[8] 薛禹群. 地下水動(dòng)力學(xué)[M]. 北京:地質(zhì)出版社, 1997.

[9] 鄭春苗. 地下水污染物遷移模擬[M]. 北京:高等教育出版社, 2009.

[10] 薛禹群, 謝春紅. 地下水?dāng)?shù)值模擬[M]. 北京:科學(xué)出版社, 2007.

[11] GB/T 14497-1993, 地下水資源管理模型工作要求[S].

[12] 廖梓龍. 包頭市地下水動(dòng)態(tài)模擬與調(diào)控研究[D]. 北京：中國(guó)水利水電科學(xué)研究院, 2013.