MODFLOW在新疆灌區地下水資源開發評價中的應用

劉 磊，董新光，吳 斌

(1.新疆農業大學 水利與土木工程學院，新疆 烏魯木齊830052;2.新疆水利廳，新疆 烏魯木齊830000)

新和縣位于新疆維吾爾自治區天山南麓，塔里木盆地北緣，東隔渭干河與庫車相望，北隔確勒塔格山與拜城縣相鄰，南連沙雅縣，東西長136 km，南北寬91 km，全縣總面積為8 823 km2，耕地面積49.35×m 104萬畝。地處歐亞大陸腹地，塔里木盆地北部，遠離海洋，北部和西部受天山屏障的阻隔，大西洋和印度洋暖濕氣流部分可翻越帕米爾高原或天山進入本縣，同時東南面又受塔克拉瑪干大沙漠的干熱氣流影響，從而使新和縣呈現暖溫帶大陸性干旱氣候的特征。

新疆的農業屬于灌溉型農業，尤其是在南疆地表水缺乏的情況下，有限的且不斷減少的地表水資源已無法保證農業經濟可持續發展。為改善新和縣灌區耕地農業用水季節性和高峰期水資源供需矛盾，提高灌溉保證程度，保證農業經濟可持續發展，科學合理的開發利用該區域地下水資源是一條有效途徑。

1 研究區水文地質概況

評價區處于渭干河沖洪積傾斜平原中上部，渭干河是該沖洪積平原區第四紀松散層孔隙地下水形成的唯一來源，渭干河在出山口后形成了巨大的沖洪積扇。渭干河沖洪積扇的中上部主要為地下水的徑流區，地下水徑流條件由北向南，由扇頂至扇前緣逐漸變差，隨著地層巖性的逐漸變細，地形坡度變緩，潛水水位埋深也逐漸變淺。在沖洪積扇的下部，由于地形坡度變緩，地層顆粒變細，從而減緩了地下水向下游排泄的速度，垂向排泄逐漸成為主要的排泄方式。

評價區地下水的補給方式主要有上游地下水側向流入補給和區內地表水的垂向入滲補給，大氣降水入滲補給鑒于區內降雨量小且次降雨量很少達到10 mm以上，在此可予以忽略。評價區地層結構簡單，地層顆粒較粗大，滲透性好，地下水在接受了渭干河及上游渠系水入滲補給后，以側向徑流流入的方式補給本區。

評價區地處渭干河西岸，渠系發育，尤魯都斯干渠、塔什艾日克干渠、依其艾日克干渠和沙雅總干渠等多條干渠均在工作區分布或穿過，同時，區內支、斗渠配套齊全，形成縱橫交錯的水系網，除干渠做防滲外，支、斗渠多防滲效果差或未防滲，從而使渠道滲漏成為測區地下水的另一主要補給源。此外，做為渭干河流域的古老灌區，評價區內農田分布廣泛，耕作層主要為粉土和粉砂，田間灌溉水大量入滲，從而形成了本區地下水的另一重要補給源。

項目區地層巖性在314國道以北主要以砂礫石夾中粗砂、粉土、粉質粘土為主，以南逐漸變為中粗砂夾粉土、粉質粘土層，第四系松散堆積物在工作區厚度大于150 m。巨厚的松散堆積層為地下水提供了良好的賦存空間，地下水類型主要為潛水，局部存在微承壓水。

在項目區東北部，含水層結構為單一的砂礫石潛水賦存區，水位埋深多為3～5 m，在縣城北部和東北部大于5 m，在西南部多在1～3 m，局部地段由于受地形影響在3～5 m之間。地下水流向，總體上與地形坡降一致，即向南、西南方向徑流，但局部因受引水渠及排水渠的影響而稍有改變。

地下水的徑流條件主要受地形條件和地層巖性控制，測區地形平坦開闊，地勢北高南低，地形坡降1.0‰ ～2.5‰。含水層巖性主要為礫石、中粗砂、細砂和粉砂，自北向南，顆粒逐漸變細，徑流條件也相應逐漸變差。評價區地下水的排泄方式有潛水蒸發蒸騰排泄、地下水側向排泄、排渠排泄和人工開采等。

2 模型的建立

2.1 模型研究范圍

地下水模型的研究范圍與地下水調查范圍一致，在平面上呈半扇形沿東北－西南方向展布，總面積為295 km2，東北端角點坐標為(640000，4620000)，東南端角點坐標為(635500，4590000)，西北角點的坐標為(621700，4600000)(見圖1)。

圖1 地下水模型研究范圍

2.2 含水層的概化

根據模型范圍內含水層介質的成因、物質組成、空間分布規律，將模擬對象概化為非均質各項同性含水層。平面上根據含水層分布規律、結合勘探孔抽水試驗資料進行了參數分區;在垂直方向大部分地段含水層變化不明顯，僅作為一層，模擬深度為150 m。

2.3 數學模型

根據含水層埋藏分布、巖性及水流特征，模型區數學模型選用一層二維非穩定流數值模型，其具體數學模型表達式如下。

式中，Kxx、Kyy分別為 X、Y方向的滲透系數，L/T;H 為水頭值，L;M為含水層厚度、ε為源匯項，L/T;S為給水度或彈性釋水系數，潛水區取重力給水度，承壓區取彈性釋水系數;Ω為模擬范圍;n為邊界面的外法線方向;Γ為側邊界;B為底邊界。

對上述地下水流數學模型，采用正交網格三維有限差分數值方法求解。為保證計算結果的收斂性，采用交替方向隱式差分格式，利用強隱式(SIP)求解差分方程。采用美國地調局(USGS)三維地下水流模擬軟件包Processing MODFLOW進行數值計算。

研究區以2000年8月實測地下水位作為模擬預測模型的初始水位。利用MODFLOW自帶的工具Digitizer對實測地下水等水位線圖數字化，再用插值模塊Field Interpolator進行插值，對每一個單元格賦初始水頭值。

2.4 邊界條件的處理

依據模擬區水文地質條件，平面上邊界條件概化如下(見圖2)，將底邊界處理為隔水邊界，上邊界作為開放邊界，考慮渠系入滲和人工開采，分別用Processing Modflow中的Recharge(面狀補給)、Well(水井)程序包來處理;側向補給邊界作為定流量邊界，用Well程序包來處理;側向排泄邊界用General Head Boundary(一般水頭邊界)程序包來處理;排水渠道用 Drain(排水)程序包來處理;蒸發用Evapotranspiration(蒸發)程序包來處理。其中AB邊界為第二類補給流量邊界;BC和FA邊界為零流量邊界;CD、DE、EF邊界為第二類排泄流量邊界。

2.4.1 空間與時間離散

空間上采用等間距有限差分的離散方法，在地下水模型中采用自動剖分，模擬范圍內將含水層離散為1層、81行、69列，差分網格的大小為400 m×400 m，網格平均單元面積0.16 km2。模型計算區單元數(有效單元數)為1 842個，面積為 294.72 km2。

時間上以2000年1月作為初始時刻，根據地下水位的觀測時間，模擬時段以月為單位，將模型調試與參數率定過程的7年分成84個時段，將2000年的初始水位運行到2007年，再以模擬好的07年的流場作為規劃年的初始流場，規劃年模型調試與參數率定過程的10年分為120個時段，模擬計算以天為最小單位。

圖2 模型邊界條件概化及規劃開采井平面分布圖

2.4.2 模型中數據選用

地面高程數據從1:50000比例尺的地形圖上提取，然后利用MODFLOW自帶的插值模塊Field Interpolator插值，再賦值到各個單元格，作為地面高程。底部高程值用地面高程減去150 m。結合含水層成因、空間分布規律及勘探孔抽水試驗資料進行了參數分區。依據抽水試驗資料確定含水層的滲透系數和給水度。

面狀補給量用Recharge程序包進行模擬，紅旗農場的面狀補給量為3 373×104m3/a，試驗林場的面狀補給量為1 685×104m3/a。側向補給量及側向排泄量用well程序包進行模擬，灌區邊界BC、CD和DE的側向補給量分別為316×104m3/a、773 ×104m3/a和621 ×104m3/a，灌區邊界 FI、IJ和JK的側向排泄量分別為、3 405×104m3/a、2 790×104m3/a和945×104m3/a。

地下水開采量用well程序包進行模擬。規劃機井開采期每月的地下水開采量依據灌區機井開采方案中的數據進行賦值。研究區排水渠道的排水量用Drain程序包進行模擬，排水量為 2419×104m3/a。研究區蒸發量用Evapotranspiration程序包進行模擬，蒸發量為689×104m3/a。

3 模型參數率定

3.1 率定依據和方法

2007年8月實測等水位線圖作為模型檢驗的流場圖，2000年的等水位線經模型運行7年后，與2007年的等水位線擬合較好的等水位線作為規劃年的初始流場。2000年到2007年評價區地下水動態觀測資料作為率定含水層參數的依據。以地下水均衡計算的地下水補給量、排泄量作為模型檢驗的輸入，以控制性勘探孔地下水位作為模型率定的檢驗標準。

采用自動與手動相結合的方法，通過計算水位和實際水位的擬合分析，不斷地修改反復調整參數，當兩者之間誤差“最小時”，即認為此時的參數值代表含水層的參數。

3.2 平面流場的擬合

經調試計算，2000年的等水位線經模型運行7年后，得到2007年8月流場與實測流場比較，趨勢基本一致，將2007年8月的流場作為規劃年的初始流場，模擬潛水初始流場與實測流場比較見圖3。計算值與實際值偏差最大為0.5 m，大部分小于0.1 m。地下水位擬合誤差為0.0～0.2 m的單元數占總單元數的82%，0.2～0.5 m單元數占總單元數的18%。故認為擬合效果較好。

圖3 07年8月模擬流場與實際流場擬合圖

3.3 數值模擬結果與地下水均衡結果對比分析

通過數值法計算的地下水補給量和排泄量，與地下水均衡法計算的補給量和排泄量結果比較分析(見表1)，可進一步提高地下水資源的評價精度。

從全區計算結果分析:地下水均衡法計算的地下水總補給量為12 671×104m3/a，模型計算為12 511×104m3/a，相差160×104m3/a，相對誤差為1.26%。地下水均衡法計算的總排泄量為13 224×104m3/a，模型計算為13 113×104m3/a，相差111×104m3/a，相對誤差為0.84%。地下水均衡法計算的均衡差為－552×104m3/a，模型計算為－602×104m3/a，相差 49 ×104m3/a。

地下水模型進行參數和各項補排量調試和調整的主要依據是將2000年流場運行7年后要與2007年實測的流場基本一致，以2007年的水均衡計算的各項地下水補排量作為模型調試的初始值。在率定參數同時，調整各項補排量，使地下水模擬的初始流場與2007年的實測流場擬合最佳為止，得出評價區模型模擬各項補給和排泄量與地下水均衡法算的量擬合最佳。本次運用兩種評價方法，從流場上講，研究區地下水模擬的初始流場與2007年的實測流場基本一致;從補排關系上，水均衡法與地下水模型計算的總補給量、總排泄量和均衡結果差異不大。故認為此計算成果是可信的。

表1 地下水均衡法與數值法計算結果比較 104m3

4 應用

4.1 開采方案布置

依據《新和縣布喀塔木水源地(擴建)工程報告》，水源地擬布機井50眼，計算區原水源地現狀開采量為6 786.14×104m3/a，原水源地如按設計開采量開采時計算區總開采量為(6 786.14～1 880) ×104m3/a=4 906.14 ×104m3/a，計算區可開采量為6 351.81×104m3/a，則尚可開(6 351.81～4 906.14) ×104m3/a=1 445.67×104m3/a。即計算區在今后的地下水開發利用過程中，尚可開采1 445.67×104m3/a的地下水量。

設計開采量是在可開采量滿足的條件下綜合灌區缺水量、缺水時間及考慮本水源地置換地表水方案進行確定的。根據這50眼機井的不同的規劃開采量，擬定了3種方案對評價區在項目運行5年后和10年后的水量與水位進行預測。

方案一:按照原水源地的開采井數和開采量繼續運行。

方案二:根據《新和縣布喀塔木水源地(擴建)工程報告》，布喀塔木水源地預計新建機井50眼，總開采量為1 000×104m3/a。原有機井的開采量保持不變，仍然為6 786.14×104m3/a，總的開采量為 7 786.14 ×104m3/a。

方案三:根據《新和縣布喀塔木水源地(擴建)工程報告》，布喀塔木水源地預計新建機井50眼，總開采量為1 000×104m3/a。原有機井的開采量減小到設計開采量4 906.14×104m3/a，總的開采量為 5 906.14 ×104m3/a。

4.2 各方案模型預測

通過不同開采方案預測，方案一與方案二的水井開采量大于水源地的可開采量，運行5年后(2012年)潛水區大部分地段地下水位在現狀年基礎上下降0.3～2 m，集中開采區下降1～2.5 m。從方案一到方案二，隨開采量增大，地下水位降幅增大，同等降幅面積在擴大。方案三開采量小于水源地的可開采量，原水源地水位在現狀年基礎上有所上升0.1～0.5 m，新建開采井在現狀年基礎上水位下降0.1～1.1 m。

各方案運行10年后(2017年)方案一與方案二的水井潛水區大部分地段地下水位在現狀年基礎上下降0.4～3 m，集中開采區下降2～3 m。從方案2到方案1，隨開采量增大，地下水位降幅增大。方案三開采量小于水源地的可開采量，原水源地水位在現狀年基礎上有所上升0.2～1.0 m，新建開采井在現狀年基礎上水位下降0.2～1.2 m。

5 結語

對新和縣老灌區通過MODFLOW軟件建立的地下水數值模擬通過驗證符合研究區地下水含水層的邊界條件和實際真實情況，模型運行穩定可靠。利用模型對新灌區的地下水的三個開采方案進行了預測分析評價，得出按照原方案或加大開采量都超過了改研究區的可開采量，運行5年和10年后該區地下水的水量和水位都有下降趨勢，故在增加開采井的情況下減少開采量是比較符合地下水資源合理利用的方案。模擬說明這一優化方案可以較好的緩解該區域地下水水資源與農業發展的矛盾，為該區域的農業經濟發展和地下水開采提供了參考。

［1］Mc Donald G，Michael and Arlen W，Harbaugh.A modular threedimensional finite－difference ground－water flow model［M］.United State Government Printing Office，Washington，1998.

［2］盧文喜.地下水運動數值模擬過程中邊界條件問題探討［J］.水利學報，2003(3):33－36.

［3］吳劍鋒，朱學愚.由MODFLIOW淺談地下水流數值模擬軟件的發展趨勢［J］.工程勘察，2000，(2):12－15.

［4］王智，高瑾，董新光.MODFLOW在三工河流域地下水資源評價中的應用［J］.新疆農業大學報，2002，25(4):29 －34.

［5］卞玉梅，盧文喜，馬洪云.Visual MODFLOW在水源地地下水數值模擬中的應用［J］.東北水利水電，2006，(3).

［6］Sun Minzheng，Chunmiao.Long－term groundwater management by a MODFLOW based dynamic optimization tool［J］.Journal of the American Water Resources Association，1999，35(1):99－111.