王艷茹，徐征和，王 通

(濟(jì)南大學(xué)資源與環(huán)境學(xué)院，山東 濟(jì)南 250022)

0 引 言

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

位山灌區(qū)地處山東省魯西北平原，是黃河下游區(qū)域最大的引黃灌區(qū)。目前，僅吳傳貴等人依據(jù)位山灌區(qū)地質(zhì)勘查和實(shí)測(cè)地下水位等數(shù)據(jù)，分析農(nóng)業(yè)灌溉和地下水開采對(duì)當(dāng)?shù)氐叵滤a(bǔ)給、徑流和排泄的影響[3]。綜合來看，已有研究中尚未應(yīng)用MODFLIW軟件進(jìn)行模型模擬，所涉及的數(shù)據(jù)空間過于單一，針對(duì)性較弱等問題，可見該區(qū)域相關(guān)研究尚顯薄弱。本次研究以位山灌區(qū)作為研究對(duì)象，根據(jù)灌區(qū)地質(zhì)、水文地質(zhì)條件調(diào)查等相關(guān)資料分析，構(gòu)建灌區(qū)淺層地下水系統(tǒng)數(shù)值模型，通過采用翔實(shí)的地下水動(dòng)態(tài)監(jiān)測(cè)資料等進(jìn)行模型調(diào)參校準(zhǔn)識(shí)別，基于校準(zhǔn)驗(yàn)證后的數(shù)值模型對(duì)位山灌區(qū)現(xiàn)狀地下水變化進(jìn)行水均衡狀態(tài)分析，為灌區(qū)地下水資源合理開發(fā)利用提供理論依據(jù)。

1 研究區(qū)概況

1.1 地理位置

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

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

1.2 水文地質(zhì)條件

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

圖2 位山灌區(qū)地質(zhì)構(gòu)造分布圖Fig.2 Geologic structure distribution of study area

位山灌區(qū)地下水動(dòng)態(tài)分析中的主要水文地質(zhì)參數(shù)包括給水度、滲透系數(shù)等。

(1)給水度μ。本文通過地下水動(dòng)態(tài)監(jiān)測(cè)資料，并結(jié)合位山灌區(qū)有關(guān)實(shí)驗(yàn)研究成果確定了不同巖性的給水度，見表1。

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

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

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

根據(jù)以上水文地質(zhì)參數(shù)范圍，對(duì)位山灌區(qū)進(jìn)行水文地質(zhì)參數(shù)分區(qū)，見表3，圖3。

2 研究方法

2.1 水文地質(zhì)概念模型

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

表3 位山灌區(qū)水文地質(zhì)參數(shù)分區(qū)Tab.3 The Hydrogeological parameter zone of Weishan irrigation area

圖3 位山灌區(qū)水文地質(zhì)參數(shù)分區(qū)圖Fig.3 The Hydrogeologic parameter zone of Weishan irrigation area

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

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

本次研究以灌區(qū)淺層地下水作為研究對(duì)象，源匯項(xiàng)是由補(bǔ)給和排泄兩項(xiàng)構(gòu)成。灌區(qū)淺層地下水主要接受降水和農(nóng)業(yè)灌溉水下滲補(bǔ)給、河流渠道滲漏以及含水層側(cè)向流入；而通過蒸散發(fā)、地下水開發(fā)利用以及含水層側(cè)向流出等方式排泄。

2.2 數(shù)學(xué)模型的建立

按照地下水運(yùn)動(dòng)定律如質(zhì)量守恒原理和達(dá)西定律等原理以及對(duì)位山灌區(qū)水文地質(zhì)條件、滲透系數(shù)等相關(guān)參數(shù)和地下水流場(chǎng)進(jìn)行科學(xué)系統(tǒng)分析，采用下述的偏微分方程組和邊界條件表示前述已概化的位山灌區(qū)地下水流動(dòng)系統(tǒng)。

(1)

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

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

圖4 WHS法主要參數(shù)設(shè)置圖Fig.4 Main parameter setting in WHS method

2.3 數(shù)值模型的建立

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

圖5 研究區(qū)域網(wǎng)格剖分圖Fig.5 Mesh generation of study area

(2)模型初始參數(shù)設(shè)置。根據(jù)位山灌區(qū)試驗(yàn)站提供的灌區(qū)水文地質(zhì)勘察報(bào)告和抽水試驗(yàn)資料，結(jié)合前人所得經(jīng)驗(yàn)參數(shù)值作為模型相關(guān)參數(shù)初始值，詳見表3。

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

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

圖6 灌區(qū)模型邊界條件設(shè)置圖Fig.6 Boundary condition set of irrigation area model

(4)模擬期的確定。依據(jù)位山灌區(qū)試驗(yàn)站提供的地下水位監(jiān)測(cè)資料以及聊城市水文局提供的灌區(qū)降水、蒸發(fā)等水文資料，分別選擇2010.1.1-2011.12.31期間的地下水?dāng)?shù)據(jù)和水文資料作為地下水?dāng)?shù)值模型調(diào)參校準(zhǔn)和驗(yàn)證的基礎(chǔ)數(shù)據(jù)。其中選擇2010.1.1-2010.12.31期間的相關(guān)數(shù)據(jù)作為模型構(gòu)建數(shù)據(jù)，主要用于模型參數(shù)的識(shí)別和調(diào)整；而2011.1.1-2011.12.31期間的數(shù)據(jù)主要用于數(shù)值模型的可靠性驗(yàn)證分析。根據(jù)研究數(shù)據(jù)完整程度，模型識(shí)別和驗(yàn)證時(shí)段均為被劃分為12個(gè)時(shí)段，每個(gè)時(shí)段內(nèi)模型的降水灌溉等補(bǔ)給項(xiàng)和蒸發(fā)等排泄項(xiàng)是一致的。

3 結(jié)果與分析

3.1 灌區(qū)模型的識(shí)別

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

圖7 研究區(qū)監(jiān)測(cè)井實(shí)測(cè)水位與模擬地下水位擬合曲線Fig.7 Fitting curve of measured and simulated water level of monitoring well in study area

位山灌區(qū)內(nèi)共布設(shè)了25眼地下水位監(jiān)測(cè)井，而且具有長(zhǎng)期的地下水位實(shí)測(cè)數(shù)據(jù)。選擇2010年實(shí)測(cè)地下水實(shí)測(cè)資料作為數(shù)值模型識(shí)別數(shù)據(jù)，同時(shí)選擇2010年6月7日灌區(qū)地下水流場(chǎng)和模擬流場(chǎng)進(jìn)行對(duì)比分析。經(jīng)過不斷調(diào)整校準(zhǔn)模型中的滲透系數(shù)等參數(shù)，使地下水位模擬結(jié)果與實(shí)測(cè)水位擬合程度較高。

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

圖8 2010年6月7日模擬流場(chǎng)與實(shí)際流場(chǎng)對(duì)照?qǐng)DFig.8 Comparison of simulated flow field and actual flow field on 7, June,2010

綜上所述，地下水位模擬值和實(shí)測(cè)值存有一定誤差，但是均分布在95%的置信區(qū)間內(nèi)，模擬結(jié)果未出現(xiàn)誤差累計(jì)和誤差擴(kuò)大問題，表明該模型的相關(guān)參數(shù)設(shè)置合理(詳見表4)，可以真實(shí)反映研究區(qū)域水文地質(zhì)概念模型的水文地質(zhì)參數(shù)，識(shí)別后的數(shù)值模型可以用于研究灌區(qū)地下水位的動(dòng)態(tài)變化和水量轉(zhuǎn)化關(guān)系。

表4 模型相關(guān)參數(shù)值Tab.4 The Correlation parameter values of model

3.2 灌區(qū)模型的驗(yàn)證

模型識(shí)別確定后需要根據(jù)研究區(qū)域識(shí)別期以外的數(shù)據(jù)進(jìn)行模型驗(yàn)證以確定模型是否真的可以準(zhǔn)確反映研究區(qū)域地下水動(dòng)態(tài)變化?；谝炎R(shí)別模型，以模型運(yùn)行結(jié)果作為初始條件，其他參數(shù)、邊界等條件不變，模擬時(shí)間為1年，運(yùn)行模型。最后選擇2011年灌區(qū)水位實(shí)測(cè)數(shù)據(jù)和模擬結(jié)果對(duì)比分析以確定數(shù)值模型的準(zhǔn)確性和可靠性。

根據(jù)模型模擬結(jié)果選擇4個(gè)典型觀測(cè)井(異于識(shí)別時(shí)選擇的觀測(cè)井)的實(shí)測(cè)水位進(jìn)行對(duì)比分析，詳見圖9。通過分析四個(gè)觀測(cè)井實(shí)測(cè)水位和模擬水位擬合曲線可知，模擬水位和實(shí)測(cè)水位擬合度較高，一些時(shí)間段模擬水位與實(shí)際水位存在偏差，但是水位動(dòng)態(tài)變化趨勢(shì)基本一致。

圖9 研究區(qū)監(jiān)測(cè)井實(shí)測(cè)水位與模擬地下水位擬合曲線Fig.9 Fitting curve of measured and simulated water level of monitoring well in study area

通過建立模型對(duì)2011年1月-2011年12月的位山灌區(qū)進(jìn)行水均衡模擬，結(jié)果詳見表5。依據(jù)《位山灌區(qū)地下水資源評(píng)價(jià)報(bào)告》可知在采用水量均衡法計(jì)算得出淺層地下水補(bǔ)給量是103 512 萬m3，淺層地下水過度開采量是3 225.4 萬m3，與模型模擬結(jié)果分別相差4.2和10.2 萬m3，可見水均衡模擬結(jié)果符合灌區(qū)實(shí)際情況，地下水?dāng)?shù)值模型模擬結(jié)果可靠。

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

3.3 地下水量均衡與水位波動(dòng)分析

(1)地下水量均衡分析。區(qū)域地下水量均衡是指在一定時(shí)段內(nèi)研究區(qū)域內(nèi)地下水補(bǔ)給量和排泄量之間的相對(duì)關(guān)系，即補(bǔ)給量大于排泄量稱之為正均衡反之為負(fù)均衡。本文依據(jù)翔實(shí)的地下水位監(jiān)測(cè)資料以及識(shí)別驗(yàn)證后的數(shù)值模型，對(duì)研究區(qū)域進(jìn)行地下水量均衡計(jì)算，詳見表6。

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

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

(2)地下水位波動(dòng)分析。大氣降水是淺層地下水的主要補(bǔ)給來源，降水與地下水位回升具有明顯的相關(guān)性。文中取16B號(hào)監(jiān)測(cè)井進(jìn)行試驗(yàn)，觀測(cè)地下水位在不同保證率下的變化，見圖10。

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

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

4 結(jié) 論

(1)MODFLOW模型程序結(jié)構(gòu)合理，離散方法簡(jiǎn)單，求解方法多樣化，易于操作，是一種適合于模擬地下水流動(dòng)態(tài)的數(shù)值模擬軟件，且擬合度較高，較為準(zhǔn)確，具有廣泛的應(yīng)用價(jià)值。

(2)由模擬結(jié)果可知灌區(qū)在模擬期內(nèi)地下水補(bǔ)給總量是103 508.2 萬m3，排泄總量是106 743.8 萬m3，二者相差-3 235.6 萬m3，由此可知位山灌區(qū)目前地下水系統(tǒng)屬于負(fù)平衡狀態(tài)。

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

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

□