基于MODFLOW的地下水動態變化特征分析
——以位山灌區為例

王艷茹，徐征和，王 通

(濟南大學資源與環境學院，山東 濟南 250022)

0 引 言

MODFLOW是模塊化三維有限差分地下水流動模型(The modular finite difference groundwater flow model,簡稱MODFLOW)，由美國地質調查局開發、用來模擬地下水流動和污染物遷移等特性的計算機程序。國內外學者利用其構建地下水系統模型，從多個方面和角度研究地下水系統動態演變規律。比如，Kapil K. Narula等人在對喜馬拉雅流域研究的基礎上，采用SWAT、MODFLOW和MT3DMS軟件構建該流域的地下水系統模型，并研究分析在變化條件下地下水的演變規律[1]。在國內，婁華君等學者通過利用MODFLOW軟件模擬分析并提出解決地下水超采造成的地質災害的措施，為合理利用地下水資源提供科學依據[2]。綜合看來，當前針對地下水系統動態演變規律的區域研究仍較為匱乏，位山灌區也不例外。

位山灌區地處山東省魯西北平原，是黃河下游區域最大的引黃灌區。目前，僅吳傳貴等人依據位山灌區地質勘查和實測地下水位等數據，分析農業灌溉和地下水開采對當地地下水補給、徑流和排泄的影響[3]。綜合來看，已有研究中尚未應用MODFLIW軟件進行模型模擬，所涉及的數據空間過于單一，針對性較弱等問題，可見該區域相關研究尚顯薄弱。本次研究以位山灌區作為研究對象，根據灌區地質、水文地質條件調查等相關資料分析，構建灌區淺層地下水系統數值模型，通過采用翔實的地下水動態監測資料等進行模型調參校準識別，基于校準驗證后的數值模型對位山灌區現狀地下水變化進行水均衡狀態分析，為灌區地下水資源合理開發利用提供理論依據。

1 研究區概況

1.1 地理位置

位山灌區位于山東省聊城市境內，南臨黃河，地處東經115°16′～116°30′、北緯35°47′～37°03′之間，總面積位0.54 萬km2，其中設計灌溉面積為0.36 萬km2。詳見圖1。

圖1 研究區域地理位置圖Fig.1 Location of study area

1.2 水文地質條件

位山灌區地屬半干旱半濕潤大陸性季風氣候，四季變化明顯，據調查可知，灌區年均降水量是548.7 mm；位山灌區位于海河流域，水資源主要是由大氣降水、河流、地下水以及客水資源構成，根據統計資料灌區地下水資源總量為91 622.6 萬m3，其現狀年地下水開發利用量為59 985.1 萬m3；位山灌區地處黃河沖淤積平原，地勢變化幅度較小，灌區西南區域地勢較高，東北區域地勢較低，地表高程范圍為28.99～34.88 m，地面自然坡降為1/7 500～1/10 000；位山灌區屬于黃河下游魯西北平原，在地質構造上地處中朝準地臺，根據灌區地質資料可知區內相間分布有多種類型地質構造單元，詳見圖2。

圖2 位山灌區地質構造分布圖Fig.2 Geologic structure distribution of study area

位山灌區地下水動態分析中的主要水文地質參數包括給水度、滲透系數等。

(1)給水度μ。本文通過地下水動態監測資料，并結合位山灌區有關實驗研究成果確定了不同巖性的給水度，見表1。

表1 位山灌區給水度μ值表Tab.1 Specific yield μ of Weishan irrigation area

(2)滲透系數K。本文根據有關抽水試驗和同位素示蹤實驗成果進行綜合分析，確定了不同巖性的滲透系數K，見表2。

表2 位山灌區不同巖性K值表Tab.2 K of different lithology in Weishan irrigation area

根據以上水文地質參數范圍，對位山灌區進行水文地質參數分區，見表3，圖3。

2 研究方法

2.1 水文地質概念模型

(1)含水層的概化。根據灌區地質勘探資料和水文地質資料分析可知，灌區內含水介質的透水性具有空間分布特點，但是每處各方向行的滲透性能相同，因此將研究區域含水層系統概化為一層非均質、各向同性潛水含水層。依據灌區試驗站提供的地下水水位、流場分布等相關資料可以將灌區地下水流場概化為服從滲流連續性方程和達西定律的平面流動。通過分析最終確定把研究區域淺層地下水系統概化為二維非穩定平面流運動模型。

表3 位山灌區水文地質參數分區Tab.3 The Hydrogeological parameter zone of Weishan irrigation area

圖3 位山灌區水文地質參數分區圖Fig.3 The Hydrogeologic parameter zone of Weishan irrigation area

(2)邊界概化。位山灌區的北部局部區域以衛河為邊界，即從王莊電灌站到張官屯水庫區域；而南部局部區域以黃河為邊界，即從牛屯鄉到魚山鄉區域，因此將該區域含水層沿河邊界概化為第一類已知水頭邊界。依據位山灌區地下水流場分布特征將牛屯鄉到王莊電灌站和魚山鄉到張官屯水庫區域概化為第二類水頭邊界。灌區淺層地下水主要從灌區上部邊界接受水量補給，依據模型概化原則將灌區上部區域概化為水量交換邊界。通過水文地質資料可知灌區下邊界為滲透性較差的黏土組成的隔水層，潛水下滲補給量較少，據此將該層概化為研究區域的隔水邊界。

(3)水文地質參數和源匯項的處理。位山灌區的水文地質參數主要是由水文地質條件和抽水試驗結果確定。根據灌區水文地質勘察報告和抽水試驗資料，位山灌區的滲透系數為0.8～5.8 m/d，給水度為0.09～0.25。在建模過程中，根據水文地質條件給出參數初始值，通過水位擬合調整參數，確定各分區的水文參數值。

本次研究以灌區淺層地下水作為研究對象，源匯項是由補給和排泄兩項構成。灌區淺層地下水主要接受降水和農業灌溉水下滲補給、河流渠道滲漏以及含水層側向流入；而通過蒸散發、地下水開發利用以及含水層側向流出等方式排泄。

2.2 數學模型的建立

按照地下水運動定律如質量守恒原理和達西定律等原理以及對位山灌區水文地質條件、滲透系數等相關參數和地下水流場進行科學系統分析，采用下述的偏微分方程組和邊界條件表示前述已概化的位山灌區地下水流動系統。

(1)

式中：A為研究區域；x，y為空間坐標，m；k是滲透系數，m/d；S是飽和含水層的貯水系數，1/m；t是時間變量，d；h(x,y)是地下水待求水位，m；h(x,y)是初始水位，m；h1(x,y)是第一類邊界水位，m；q(x,y,t)是第二類邊界單寬流量，m3/d；n是第二類邊界條件內法線方向的單位向量；S1，S2是第一類和第二類邊界。

構建的地下水數學模型屬于二階非線性偏微分方程，在求解過程中需要先將待解方程轉化為線性方程再求解。本文選用以有限差分法為基礎求解方法的Visual MDFLOW軟件進行數學模型求解。該軟件可以根據研究區域不同區域功能進行模塊化設置，有利于進行模型參數和邊界條件等劃分和設定。除此之外，軟件自帶多種解算器如PCG2、SIP、SOR等，可選擇不同解算器進行地下水流動模型數值方程迭代計算以減少計算時間和過程，提高數值模擬效率。根據本文研究要求選擇采用WHS解算器進行模型模擬計算，需要對其中的一些相關參數進行初步設置。

圖4 WHS法主要參數設置圖Fig.4 Main parameter setting in WHS method

2.3 數值模型的建立

(1)模型剖分。依據位山灌區試驗站提供的灌區區域圖和DEM圖，采用ArcGIS軟件進行圖形校準并進行矢量化處理和灌區高程提取，其分別作為構建數值模型的底圖和模型地表高程數據。根據灌區面積和數據質量以及模擬結果精度要求，使用Visual Modflow軟件對灌區模型進行網格劃分，如圖5。灌區模型共被劃分為200×200(行×列)，其中有效單元格(白色部分)為24 057個，無效單元格(綠色部分)為15 943個；藍色區域為灌區的河流補給邊界；由于研究對象為淺層含水層，因此在垂直方向上設定一層，而灌區含水層的底板標高設置為60 m。圖5中顯示為白色的區域均為有效研究單元格，而綠色區域的單元格為無效，模型計算時不進行模擬計算，計算節點位于剖分單元中心。

圖5 研究區域網格剖分圖Fig.5 Mesh generation of study area

(2)模型初始參數設置。根據位山灌區試驗站提供的灌區水文地質勘察報告和抽水試驗資料，結合前人所得經驗參數值作為模型相關參數初始值，詳見表3。

(3)邊界條件的處理。模型邊界主要由研究區域的四周邊界和河流邊界組成。根據上述邊界條件，依據灌區地下水流場分布特征確定該邊界參數，采用Visual Modflow軟件中的GHB程序包對該邊界賦值通過調節邊界水頭控制邊界對模型的補給和排泄。

基于聊城市提供的相關等資料，結合相關文獻，確定不同河流在各河段的滲透系數，利用軟件自帶的river程序設定相關參數賦值計算。模型底層設定為不透水邊界。詳見圖6。

圖6 灌區模型邊界條件設置圖Fig.6 Boundary condition set of irrigation area model

(4)模擬期的確定。依據位山灌區試驗站提供的地下水位監測資料以及聊城市水文局提供的灌區降水、蒸發等水文資料，分別選擇2010.1.1-2011.12.31期間的地下水數據和水文資料作為地下水數值模型調參校準和驗證的基礎數據。其中選擇2010.1.1-2010.12.31期間的相關數據作為模型構建數據，主要用于模型參數的識別和調整；而2011.1.1-2011.12.31期間的數據主要用于數值模型的可靠性驗證分析。根據研究數據完整程度，模型識別和驗證時段均為被劃分為12個時段，每個時段內模型的降水灌溉等補給項和蒸發等排泄項是一致的。

3 結果與分析

3.1 灌區模型的識別

模型識別是指通過不斷地調整模型中的相關參數，每調整一次使用模型運行模擬地下水位，然后分析其與該時段實測水位的誤差，如果誤差較大則需要再一次調整，直至模擬水位和實測水位誤差滿足要求為止，即模擬曲線和實測水位能夠最佳擬合。

圖7 研究區監測井實測水位與模擬地下水位擬合曲線Fig.7 Fitting curve of measured and simulated water level of monitoring well in study area

位山灌區內共布設了25眼地下水位監測井，而且具有長期的地下水位實測數據。選擇2010年實測地下水實測資料作為數值模型識別數據，同時選擇2010年6月7日灌區地下水流場和模擬流場進行對比分析。經過不斷調整校準模型中的滲透系數等參數，使地下水位模擬結果與實測水位擬合程度較高。

在模型識別期內，在研究區域內選擇4個典型地點和兩個對比流場進行模擬擬合分析。由圖7可知，四個典型地點的地下水位模擬結果和地下水位實測結果擬合程度較好，二者動態變化趨勢基本一致。選擇2010年6月7日的模擬地下水流場分布圖和實際流場分布圖基本相符，詳見圖8。按照《地下水資源管理模型工作要求》(GB/T14497-93)對比分析灌區所有監測井的模擬結果和實測數據可知，擬合誤差小于1 m的監測井有21個，占有所有觀測孔的84%，同時計算各觀測井水位實測值和模擬值的殘差均值為0.37 m，二者的相關系數為0.965。

圖8 2010年6月7日模擬流場與實際流場對照圖Fig.8 Comparison of simulated flow field and actual flow field on 7, June,2010

綜上所述，地下水位模擬值和實測值存有一定誤差，但是均分布在95%的置信區間內，模擬結果未出現誤差累計和誤差擴大問題，表明該模型的相關參數設置合理(詳見表4)，可以真實反映研究區域水文地質概念模型的水文地質參數，識別后的數值模型可以用于研究灌區地下水位的動態變化和水量轉化關系。

表4 模型相關參數值Tab.4 The Correlation parameter values of model

3.2 灌區模型的驗證

模型識別確定后需要根據研究區域識別期以外的數據進行模型驗證以確定模型是否真的可以準確反映研究區域地下水動態變化。基于已識別模型，以模型運行結果作為初始條件，其他參數、邊界等條件不變，模擬時間為1年，運行模型。最后選擇2011年灌區水位實測數據和模擬結果對比分析以確定數值模型的準確性和可靠性。

根據模型模擬結果選擇4個典型觀測井(異于識別時選擇的觀測井)的實測水位進行對比分析，詳見圖9。通過分析四個觀測井實測水位和模擬水位擬合曲線可知，模擬水位和實測水位擬合度較高，一些時間段模擬水位與實際水位存在偏差，但是水位動態變化趨勢基本一致。

圖9 研究區監測井實測水位與模擬地下水位擬合曲線Fig.9 Fitting curve of measured and simulated water level of monitoring well in study area

通過建立模型對2011年1月-2011年12月的位山灌區進行水均衡模擬，結果詳見表5。依據《位山灌區地下水資源評價報告》可知在采用水量均衡法計算得出淺層地下水補給量是103 512 萬m3，淺層地下水過度開采量是3 225.4 萬m3，與模型模擬結果分別相差4.2和10.2 萬m3，可見水均衡模擬結果符合灌區實際情況，地下水數值模型模擬結果可靠。

表5 研究區域淺層地下水水量均衡分析表Tab.5 Water equilibrium analysis of shallow groundwater in irrigation area

3.3 地下水量均衡與水位波動分析

(1)地下水量均衡分析。區域地下水量均衡是指在一定時段內研究區域內地下水補給量和排泄量之間的相對關系，即補給量大于排泄量稱之為正均衡反之為負均衡。本文依據翔實的地下水位監測資料以及識別驗證后的數值模型，對研究區域進行地下水量均衡計算，詳見表6。

由表6可知，在不同保證率下，降雨量不同，補給總量隨著降雨量的減少而減少；補給項主要包括降水入滲補給量、農業灌溉水下滲補給量等，補給項受降雨影響較大，而排泄項主要是蒸散發量、地下水開采量等，受降雨影響較小，所以表6中，隨著降雨量的較少，補給總量減少明顯，排泄項量變化不大，均衡差增大。

表6 不同保證率下淺層地下水量均衡分析表Tab.6 Water equilibrium analysis of shallow groundwater under different guaranteed rates

(2)地下水位波動分析。大氣降水是淺層地下水的主要補給來源，降水與地下水位回升具有明顯的相關性。文中取16B號監測井進行試驗，觀測地下水位在不同保證率下的變化，見圖10。

圖10 不同保證率下監測井水位分析圖Fig.10 Analysis of underground water level of monitoring well under different guaranteed rates

由圖10可知，在不同保證率下，水位變化趨勢大致相同。但是，隨著降雨量的減少，地下水位降低，水位埋深越來越低，可見降雨量與地下水位呈正相關。

4 結 論

(1)MODFLOW模型程序結構合理，離散方法簡單，求解方法多樣化，易于操作，是一種適合于模擬地下水流動態的數值模擬軟件，且擬合度較高，較為準確，具有廣泛的應用價值。

(2)由模擬結果可知灌區在模擬期內地下水補給總量是103 508.2 萬m3，排泄總量是106 743.8 萬m3，二者相差-3 235.6 萬m3，由此可知位山灌區目前地下水系統屬于負平衡狀態。

(3)位山灌區淺層地下水過度開采量為3 225.4 萬m3，與模型模擬結果分別相差10.2 萬m3，表明灌區地下水均衡模擬結果可靠，符合灌區實際情況。

(4)大氣降水是影響水量均衡、地下水位波動的重要因素，且大氣降水與水量均衡、地下水位呈正相關。

□