敦煌月牙泉地區(qū)人工回灌下的地下水動(dòng)態(tài)模擬

劉 暢，成建梅，蘇春利，楊俊倉(cāng)

(1.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北武漢 430074;2.甘肅省地質(zhì)環(huán)境監(jiān)測(cè)院，甘肅蘭州 730050)

甘肅省敦煌市月牙泉是具有世界影響力的旅游資源，然而近幾十年來(lái)泉湖水位持續(xù)下降，并呈現(xiàn)出逐漸加劇的趨勢(shì)［1-2］。為防止月牙泉干涸，人們嘗試了湖底清淤挖深、直接向泉湖注水等措施，但效果并不顯著，與此還伴生了湖岸失穩(wěn)坍塌、湖水渾濁等問(wèn)題，因此人們轉(zhuǎn)換思路，于2007 年開(kāi)展“月牙泉湖水位下降應(yīng)急治理工程”，在月牙泉附近建滲水池A 場(chǎng)和C 場(chǎng)進(jìn)行人工回灌，以提高局部地下水水位，改變?cè)卵廊苓呅》秶貐^(qū)地下水滲流場(chǎng)特征，從而達(dá)到抑制泉湖水位下降的目的［3］。該工程于2008 年建成并投入使用至今，泉湖水位已有明顯提升。

模擬技術(shù)是地下水流場(chǎng)變化研究與水資源管理的工具之一。近十多年來(lái)，已有不少學(xué)者運(yùn)用數(shù)值模擬技術(shù)進(jìn)行西北干旱半干旱區(qū)水資源研究。但是多數(shù)西北地下水模型研究?jī)A向于整體盆地或流域的模擬，對(duì)典型小范圍地區(qū)的精細(xì)模型研究相對(duì)較少［4-8］。特別是針對(duì)人工回灌、研究局部地下水流場(chǎng)變化的模型并不多見(jiàn)［9］，同時(shí)不同因素對(duì)地下水動(dòng)態(tài)變化的敏感性分析不多。

因此，筆者應(yīng)用FEFLOW 軟件，在研究敦煌盆地區(qū)域地下水流場(chǎng)的基礎(chǔ)上，建立月牙泉地區(qū)小區(qū)域地下水模型，重點(diǎn)分析人工回灌對(duì)月牙泉周邊小范圍地區(qū)地下水滲流場(chǎng)的影響，分析其提升月牙泉湖水位的機(jī)理，探討影響泉湖水位動(dòng)態(tài)的主要因素，并通過(guò)模型預(yù)測(cè)求解未來(lái)幾年人工回灌的最佳注入方式，為月牙泉保護(hù)工程合理實(shí)施提供理論依據(jù)。

1 水文地質(zhì)概況

敦煌盆地屬于典型大陸性干旱氣候，地處疏勒河流域下游西部。盆地南部卡拉塔什塔格山、北截山與北部北山在東部雙塔水庫(kù)相接，構(gòu)成楔形邊界，西部邊界位于甘肅、新疆交界的庫(kù)姆塔格沙漠邊緣，面積約12961 km2。天然條件下，敦煌盆地地下水接受河水滲漏和東部雙塔水庫(kù)處地下水的側(cè)向徑流補(bǔ)給，徑流自東向西，在西湖鄉(xiāng)、玉門(mén)關(guān)一帶與黨河沖洪積扇自南向北的地下徑流相匯，最終向西徑流，在下游水位淺埋深區(qū)以蒸發(fā)形式排泄。月牙泉位于黨河沖洪積扇南北向地下水徑流帶(圖1(據(jù)耿昕等［10］))。

第四紀(jì)以來(lái)敦煌盆地不斷沉降，盆地內(nèi)沉積了一套厚達(dá)數(shù)百米的第四系松散堆積物，在地貌上構(gòu)成了黨河與西水溝兩個(gè)規(guī)模較大的沖洪積扇，月牙泉位于扇間洼地處。黨河沖洪積扇第四系厚度一般大于300 m，且沉積物多以礫石、砂礫石為主，地下水屬單一結(jié)構(gòu)潛水含水層［11］。黨河是區(qū)內(nèi)主要常年性地表水系，距月牙泉約4.50 km。1975 年10 月黨河水庫(kù)修建后，河水全部被截流入庫(kù)，下游河床除暴雨期有暫時(shí)性水流外，一般不產(chǎn)流。黨河沖洪積扇地下水自南向北徑流，受微地形控制，局部地下水流向偏離主流向。主要表現(xiàn)在:三危山隆起阻隔黨河上游及南部地下徑流，山區(qū)接受降水匯集并補(bǔ)給其北部山前含水層，使得山前地下水向北東月牙泉方向徑流，如圖2(據(jù)黎志恒等［12］，略有修改)所示。因此，月牙泉地區(qū)地下水徑流方向自西南向北東，如圖3 所示。

圖1 敦煌盆地平面位置

圖2 研究區(qū)水文地質(zhì)略圖

2 模擬模型

2.1 概念模型

圖3 模擬區(qū)范圍及邊界條件

模擬區(qū)范圍及邊界條件如圖3 所示，模擬區(qū)淺層含水系統(tǒng)概化為非均質(zhì)各向異性三維不穩(wěn)定流動(dòng)系統(tǒng)。為提高滲流場(chǎng)在人工回灌下變化的模擬精度，模擬區(qū)僅包括月牙泉地區(qū)約9 km2范圍。由于模擬時(shí)難以選取相對(duì)獨(dú)立的天然水文地質(zhì)單元，只能進(jìn)行人為邊界的處理。筆者在劃定人為邊界時(shí)充分考慮區(qū)域地下水流場(chǎng)在模擬區(qū)范圍的典型特征，參考研究區(qū)實(shí)測(cè)等水位線圖，劃定的人為隔水邊界總體與等水位線垂直(與流線一致)，并且在模型識(shí)別階段反復(fù)擬合長(zhǎng)觀孔水位變化，作為邊界控制與識(shí)別的依據(jù)，以保證計(jì)算的模擬區(qū)邊界特征符合區(qū)域流場(chǎng)特征。因此，模型邊界條件概化為:南、北邊界為人為隔水邊界，且南部邊界與模擬區(qū)南部地表分水嶺基本重合;西部和東部分別為側(cè)向流入(q1)、流出(q2)邊界;頂部為潛水面邊界;由于淺層含水系統(tǒng)底部存在多層弱透水粉砂及亞砂土層，模型忽略深部地下水系統(tǒng)與淺層地下水系統(tǒng)之間的水力聯(lián)系，將其底部300 m 處近似處理為隔水底板。需指出的是，黨河水庫(kù)規(guī)律性放水時(shí)河道滲漏對(duì)模擬區(qū)含水層局部位置的側(cè)向徑流補(bǔ)給由q'1刻畫(huà)。南部q″1刻畫(huà)受微地形控制局部地下水流向偏離區(qū)域地下水主流向而引起的徑流改變。

2.2 數(shù)學(xué)模型

上述水文地質(zhì)概念模型可用下列數(shù)學(xué)模型［13］描述:

式中:x，y，z 為笛卡爾坐標(biāo)系，L;H 為水頭，L;Kij為滲透系數(shù)張量(i，j =x，y，z)，L/T;t 為時(shí)間，T;H0為含水層的初始水頭，L;μs為單位儲(chǔ)水系數(shù)，L-1;μd為重力給水度，無(wú)量綱;q 為單位體積多孔介質(zhì)源(或匯)的流量，抽水時(shí)為負(fù)，L/T;W 為潛水面上的降水入滲和蒸發(fā)強(qiáng)度，L/T;n 為邊界外法線方向;Ω 為模擬滲流區(qū)域;Γ0為滲流區(qū)域的上邊界，即地下水的自由水面;Γ1為滲流區(qū)的第二類邊界;v(x，y，z，t)為t 時(shí)刻點(diǎn)(x，y，z)處的地下水滲透流速，L/T。

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

采用空間6 結(jié)點(diǎn)三棱柱為基本剖分單元，平面用三角形單元離散，結(jié)點(diǎn)數(shù)為24 544，三角形單元數(shù)為44 475。區(qū)內(nèi)含水層根據(jù)巖性水平和垂直方向均按由粗變細(xì)的遞變規(guī)律細(xì)劃為13 個(gè)子層。采用2008 年1 月1 日長(zhǎng)觀測(cè)孔水位數(shù)據(jù)，插值得到初始水位(圖4)。識(shí)別時(shí)段為2008 年1 月—2009 年12月，時(shí)間步長(zhǎng)為30 天。

模擬區(qū)內(nèi)有巖性剖面資料的鉆孔10 個(gè)，以淺孔為主，深孔2 個(gè)。模型總厚度約300 m，因此每一子層巖性分布主要通過(guò)深孔巖性柱狀圖并結(jié)合沖洪積扇扇間洼地巖性遞變規(guī)律外推得到，以保證模擬區(qū)各分層巖性變化符合區(qū)域巖性變化規(guī)律。相應(yīng)巖性對(duì)應(yīng)的水文地質(zhì)參數(shù)取值范圍根據(jù)前期抽水試驗(yàn)現(xiàn)場(chǎng)資料、室內(nèi)滲流實(shí)驗(yàn)以及水文地質(zhì)手冊(cè)經(jīng)驗(yàn)數(shù)據(jù)綜合分析確定(表1)。需要指出的是:由于模擬區(qū)主要含水層巖性為松散孔隙介質(zhì)，其滲透性能表現(xiàn)為非均質(zhì)二度各向異性，水平方向滲透系數(shù)Kxx=Kyy，而垂向上滲透系數(shù)約為水平向滲透系數(shù)的1/5［14］，即Kzz=0.2Kxx=0.2Kyy。

表1 模擬區(qū)各分層水文地質(zhì)參數(shù)取值范圍

模擬區(qū)地下水側(cè)向流入和流出量根據(jù)邊界過(guò)水?dāng)嗝婧穸?、水力梯度以及區(qū)域總體水量均衡初步確定。敦煌盆地地下水每年7—8 月因開(kāi)采而大幅下降，12 月至翌年3 月水位逐漸回升，存在年內(nèi)波動(dòng)。因此，以臨近流量邊界處水位長(zhǎng)觀測(cè)孔水位動(dòng)態(tài)序列為參照，賦予流量邊界相應(yīng)的波動(dòng)規(guī)律和強(qiáng)度。

圖4 2008 年1 月初始等水位線

模擬區(qū)西北側(cè)有沙山支渠通過(guò)，渠系水來(lái)源于黨河水庫(kù)泄水，其水量動(dòng)態(tài)基本與黨河水庫(kù)每月泄水動(dòng)態(tài)保持一致。將水庫(kù)泄水量乘以適當(dāng)?shù)恼蹨p系數(shù)，以此模擬渠系水滲漏對(duì)地下水的補(bǔ)給。

月牙泉水位下降應(yīng)急治理工程總體設(shè)計(jì)方案按其供水水源可分為地表水回灌和地下水回灌兩個(gè)系統(tǒng)。地表水源通過(guò)明渠引黨河水庫(kù)水至滲水池，地下水源將采自月牙泉西上豐口子北部黨河河漫灘的地下水通過(guò)專用管道引至滲水池。模擬區(qū)人工回灌A 場(chǎng)距月牙泉西1 000 m，C 場(chǎng)距月牙泉東500 m。A場(chǎng)在2008—2009 年2 年中滲水量很小，識(shí)別階段可忽略不計(jì)。C 場(chǎng)以面源方式補(bǔ)給注入，將每月注入總量平均分配到滲水區(qū)相應(yīng)結(jié)點(diǎn)。C 場(chǎng)總注入量歷時(shí)變化見(jiàn)圖5，模型識(shí)別過(guò)程水位擬合曲線見(jiàn)圖6，2009 年12 月計(jì)算等水位線見(jiàn)圖7。

圖5 C 場(chǎng)總注入量歷時(shí)變化

模型識(shí)別的部分參數(shù):①種植、非種植區(qū)降雨入滲系數(shù)分別為0.13、0.09;②種植、非種植區(qū)蒸發(fā)系數(shù)分別為0.03、0.01;③綠地回灌系數(shù)為0.2;④折減系數(shù)為9.09 ×10-6;⑤側(cè)向流入、流出邊界水力梯度分別為0.075%、1.85%。

圖6 模型識(shí)別過(guò)程水位擬合曲線

雖然建模中邊界選取的人為因素較多，但從擬合結(jié)果(圖6、表2)可知，其擬合效果較好，模擬流場(chǎng)可近似代替地下水實(shí)際流場(chǎng)，模型能夠用于預(yù)測(cè)。

圖7 2009 年12 月計(jì)算等水位線

表2 水位擬合誤差統(tǒng)計(jì)

通過(guò)模擬識(shí)別得到以下結(jié)論。

a. C 場(chǎng)人工回灌沒(méi)有直接入滲補(bǔ)給到月牙泉湖，而是使回灌點(diǎn)周圍局部潛水面隆起，水力梯度減小(等水位線變稀疏，如圖7)，造成月牙泉向下游排泄相對(duì)受阻，在上游來(lái)水量不變的情況下，令泉湖水位抬升。C 場(chǎng)現(xiàn)有回灌滲入量并未宏觀改變模擬區(qū)原有地下水滲流場(chǎng)補(bǔ)徑排條件。

b. 影響模擬區(qū)地下水水位大幅度抬升和波動(dòng)的主要因素是沙山支渠滲漏和C 場(chǎng)人工回灌，而降水、蒸發(fā)、局部綠地灌溉對(duì)其影響不大。

c. 月牙泉湖水位持續(xù)下降由區(qū)域地下水負(fù)均衡引起，在模型中主要表現(xiàn)在降水量、局部綠地回灌補(bǔ)給、側(cè)向流入、渠系滲漏等的總量小于蒸發(fā)量與側(cè)向流出量的總量。

4 模型預(yù)測(cè)

4.1 預(yù)測(cè)方案設(shè)計(jì)

人工回灌在短期內(nèi)能夠顯著提升月牙泉湖水位，但在區(qū)域地下水負(fù)均衡造成區(qū)域地下水水位持續(xù)下降的背景下，它能否長(zhǎng)久有效是本文最為關(guān)注的問(wèn)題。因此，本節(jié)在模型識(shí)別的基礎(chǔ)上，預(yù)測(cè)A場(chǎng)、C 場(chǎng)在不同方案設(shè)計(jì)(如表3)下月牙泉湖水位5年內(nèi)的動(dòng)態(tài)，并求解最優(yōu)回灌方案。

表3 模型預(yù)測(cè)方案設(shè)計(jì)

預(yù)測(cè)中降雨蒸發(fā)量、渠系滲漏量、側(cè)向流入流出量等均采用多年平均值。預(yù)測(cè)時(shí)段為2010 年1 月至2014 年12 月，初始條件采用2009 年12 月模型識(shí)別結(jié)束時(shí)刻相應(yīng)值。

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

以2009 年12 月泉湖水位計(jì)算值1 134.34 m 為基準(zhǔn)，預(yù)測(cè)泉湖水位大于1 134.34 m 時(shí)取正值，小于1 134.34 m 時(shí)取負(fù)值，結(jié)果如圖8 所示。

圖8 不同回灌方案下泉湖水位變化

由于回灌量、回灌點(diǎn)位置差異，不同方案下的泉湖水位動(dòng)態(tài)各異。模擬預(yù)測(cè)結(jié)論如下:

圖9 2014 年12 月(方案6)預(yù)測(cè)等水位線

a. 截至2014 年12 月，使泉湖水位上升的方案包括方案2、3、6 和7，其中5 年后方案2 泉湖水位上升最多，達(dá)2.25 m，方案7 上升最小，約0.65 m;考慮預(yù)測(cè)時(shí)段區(qū)域地下水年降幅可能更為劇烈，以及綜合投入產(chǎn)出經(jīng)濟(jì)效益，最終確定方案6 為所求最優(yōu)解，即A 場(chǎng)回灌量2 500 m3/d。

b. A 場(chǎng)對(duì)泉湖水位的影響較C 場(chǎng)顯著，方案2和方案6 A 場(chǎng)日回灌量均為2500 m3/d，且方案2 比方案6 增開(kāi)C 場(chǎng)(5 000 m3/d)，但預(yù)測(cè)結(jié)束最終提升泉湖水位僅較方案6 多0.05 m。結(jié)果表明，為提高月牙泉湖水位，在其上游增加補(bǔ)給量?jī)?yōu)于在下游減少排泄量。

5 結(jié) 論

通過(guò)建立人工回灌下月牙泉周邊小范圍地區(qū)淺層地下水三維不穩(wěn)定流數(shù)值模型，探清C 場(chǎng)回灌補(bǔ)給和渠系滲漏是2008—2009 年影響模擬區(qū)水位動(dòng)態(tài)的主要因素;A、C 兩場(chǎng)回灌并未宏觀改變模擬區(qū)原有地下水滲流場(chǎng)補(bǔ)徑排條件;A 場(chǎng)增大上游補(bǔ)給量可提高泉湖水位，C 場(chǎng)減小下游水力梯度進(jìn)而抑制泉湖排泄;在提升月牙泉湖水位方面，A 場(chǎng)效果優(yōu)于C 場(chǎng)，A 場(chǎng)以回灌量2 500 m3/d 為最優(yōu)方案，能滿足泉湖水位在2014 年前穩(wěn)定的目標(biāo)。研究成果為月牙泉應(yīng)急治理工程合理實(shí)施提供理論依據(jù)。

［1］丁宏偉，龔開(kāi)誠(chéng).敦煌月牙泉湖水位持續(xù)下降原因及對(duì)策分析［J］. 水文地質(zhì)工程地質(zhì)，2004(6):74-77.(DING Hongwei，GONG Kaicheng. Analyses of the reasons and countermeasures for the decline in the water level of the Yueya Spring Lake near Dunhuang［J］.Hydrogeology and Engineering Geology，2004(6):74-77.(in Chinese))

［2］尹念文，侯燕軍，李彬.月牙泉水位變化的原因分析［J］.地下水，2010，32(3):33-34. (YIN Nianwen，HOU Yanjun，LI Bin. Analyses of the reasons for changes in the water level of the Yueya Spring［J］. Ground water，2010，32(3):33-34.(in Chinese))

［3］黎志恒，邱蔓，楊俊倉(cāng)，等.甘肅省敦煌市月牙泉湖水位下降治理綜合研究報(bào)告［R］.蘭州:甘肅省地質(zhì)環(huán)境監(jiān)測(cè)院，2004:62.

［4］朱高峰，馬金珠，李自珍，等.民勤綠洲地下水系統(tǒng)數(shù)值模擬［J］.蘭州大學(xué)學(xué)報(bào):自然科學(xué)版，2005，41(5):5-9. (ZHU Gaofeng，MA Jinzhu，LI Zizhen，et al.Groundwater simulation of Minqin Oasis［J］. Journal of Lanzhou University:Natural Sciences，2005，41(5):5-9.(in Chinese))

［5］董少剛，唐仲華，馬騰，等. 太原盆地地下水?dāng)?shù)值模擬［J］. 水資源保護(hù)，2009，25 (3):25-27. (DONG Shaogang，TANG Zhonghua，MA Teng，et al. Numerical simulation of groundwater flow in Taiyuan Basin［J］.Water Resources Protection，2009，25(3):25-27. (in Chinese))

［6］喬曉英，王文科，姜桂華，等.西北干旱內(nèi)陸盆地地下水生態(tài)功能的探討［J］.水資源保護(hù)，2005，21(5):6-10.(QIAO Xiaoying，WANG Wenke，JIANG Guihua，et al.Study on ecological function of groundwater in northwest arid inland basin ［J］. Water Resources Protection，2005，21(5):6-10. (in Chinese))

［7］SAYIT A P，YAZICIGIL H. Assessment of artificial aquifer recharge potential in the Kucuk Menderes River Basin，Turkey［J］. Hydrogeology Journal，2012，20(4):755-766.

［8］AL-ASSA’D T A，ABDULLA F A. Artificial groundwater recharge to a semi-arid basin:case study of Mujib aquifer，Jordan［J］. Environmental Earth Sciences，2010，60(4):845-859.

［9］李恒太，石萍，武海霞.地下水人工回灌技術(shù)綜述［J］.中國(guó)國(guó)土資源經(jīng)濟(jì)，2008，21(3):41-48. (LI Hengtai，SHI Ping，WU Haixia. Artificial recharge technology of groundwater［J］. Natural Resource Economics of China，2008，21(3):41-48. (in Chinese))

［10］耿昕，黎濤，田春艷，等.河西走廊疏勒河流域地下水動(dòng)態(tài)數(shù)值模型研究［R］. 武漢:中國(guó)地質(zhì)大學(xué)(武漢)，2010:9.

［11］桑學(xué)鋒，張明泉，王浩，等.敦煌盆地地下水?dāng)?shù)值模擬及可視化與管理［J］.蘭州大學(xué)學(xué)報(bào):自然科學(xué)版，2007，43(3):18-22. (SANG Xuefeng，ZHANG Mingquan，WANG Hao，et al. Visual simulation and management of groundwater in Dunhuang Basin［J］. Journal of Lanzhou University:Natural Sciences，2007，43(3):18-22. (in Chinese))

［12］黎志恒，邱蔓.甘肅敦煌月牙泉地區(qū)水文地質(zhì)與環(huán)境地質(zhì)綜合勘查報(bào)告［R］. 蘭州:甘肅省地質(zhì)環(huán)境監(jiān)測(cè)院，2008.

［13］陳崇希，林敏.地下水動(dòng)力學(xué)［M］.武漢:中國(guó)地質(zhì)大學(xué)出版社，2006:13.

［14］鄧惠森.巖土層垂向滲透系數(shù)與徑向滲透系數(shù)［J］.工程勘察，1993(1):44-45. (DENG Huisen. The vertical permeability coefficient and radial permeability coefficient of strata［J］. Geotechnical Investigation and Surveying，1993(1):44-45. (in Chinese))