李樹文，于虎廣

(1.河北工程大學(xué)，河北邯鄲 056038;2.河北省資源勘測(cè)研究重點(diǎn)實(shí)驗(yàn)室，河北邯鄲 056038)

MODFLOW在磁縣地下水?dāng)?shù)值模擬中的應(yīng)用

李樹文1，2，于虎廣1

(1.河北工程大學(xué)，河北邯鄲 056038;2.河北省資源勘測(cè)研究重點(diǎn)實(shí)驗(yàn)室，河北邯鄲 056038)

磁縣東部平原區(qū)水環(huán)境由于水資源利用不合理、地下水超采等原因發(fā)生惡化.通過分析研究區(qū)水文地質(zhì)條件，運(yùn)用MODFLOW進(jìn)行淺層地下水?dāng)?shù)值模擬.結(jié)果表明，在保持開采條件不變的情況下，至2020年，漏斗面積由18.59 km2增大為23.31 km2.提出合理調(diào)整和優(yōu)化地下水開采方案可以有效地避免降落漏斗影響范圍的增大.

地下水;數(shù)值模擬;MODFLOW

研究區(qū)位于磁縣東部平原，屬暖溫帶半濕潤(rùn)半干旱大陸性季風(fēng)氣候區(qū)，多年平均降水量531.4 mm，每年汛期(6—9月份)的降水量約占全年降水量的70% ～80%.

研究區(qū)是由漳河、滏陽(yáng)河和牤牛河的山前沖積物組成的太行山前沖積扇.其中南部為漳河沖積扇，北部是沙洺河和漳河沖積、洪積扇之間的交接帶，多是滏陽(yáng)河和牤牛河出山口的沖積和洪積物.

地下水賦存于第四系松散巖層中，主要受大氣降水和地表水體滲入補(bǔ)給及山前潛水的側(cè)向補(bǔ)給.含水砂層分布顯示了典型的山前沖積扇規(guī)律，主流部位，含水層砂粒粗、厚度大，向主流兩翼漸變?yōu)橹屑?xì)砂，地下水流向與地形基本一致.

近年來該地區(qū)地下水供水量約占區(qū)域總供水量的70%左右［1］，長(zhǎng)期超量開采致使地下水資源衰竭，造成地面沉降，當(dāng)?shù)氐墓┧踩艿絿?yán)重影響.

1 三維地下水流模擬

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

研究區(qū)西部為山區(qū)，接受山前側(cè)向徑流補(bǔ)給，將研究區(qū)西部邊界概化為二類流量邊界，流量以注水井的形式賦入模型，各段流量值利用達(dá)西定律分別計(jì)算.研究區(qū)東部有多個(gè)地下水位檢測(cè)井，故將東部邊界概化為一類水頭邊界.模擬區(qū)通過上部邊界接受補(bǔ)給和排泄，補(bǔ)給方式主要為地表水體入滲，排泄方式為垂直蒸發(fā)，故將上部平面邊界概化為二類流量邊界.下部邊界主要為幾乎不透水的泥巖，故概化為隔水邊界.

綜上所述，研究區(qū)可以概化為非均質(zhì)各向同性三維非穩(wěn)定地下水流系統(tǒng).

1.2 數(shù)學(xué)模擬模型

數(shù)學(xué)模型是在研究區(qū)水文地質(zhì)條件的基礎(chǔ)上建立的，如下式所示:

式中:K為含水層滲透系數(shù)，m/d;W為單位時(shí)間、單位體積上從垂直方向上流入或流出含水層的水量(流入為正，流出為負(fù))，m3/d;SS為貯水率;Ω為滲流區(qū)域;H為水位，m;t為時(shí)間，d;H0為初始時(shí)刻水位值，m.

模型采用Visual MODFLOW 4.2進(jìn)行求解.利用MapGIS對(duì)研究區(qū)1∶50 000第四系水文地質(zhì)圖做矢量化處理，作為模擬區(qū)的底圖.對(duì)研究區(qū)進(jìn)行矩形網(wǎng)格剖分，網(wǎng)格大小為500 m×500 m，共剖分70行50列，去除無效單元格，有效單元體共計(jì)996個(gè).基于前期研究工作中各個(gè)階段的調(diào)查，并結(jié)合抽水試驗(yàn)結(jié)果，將研究區(qū)含水層參數(shù)分為3個(gè)分區(qū)，如圖1所示.給定滲透系數(shù)和給水度的初值，見表1.

圖1 研究區(qū)含水層參數(shù)分區(qū)

表1 水文地質(zhì)參數(shù)初值

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

根據(jù)水文地質(zhì)概念模型建立的數(shù)學(xué)模型必須反

表2 水文地質(zhì)參數(shù)識(shí)別結(jié)果

為了確保數(shù)學(xué)模型的可靠性，還需要進(jìn)行模型驗(yàn)證.模型驗(yàn)證數(shù)據(jù)為從已有的地下水位觀測(cè)資料中篩選出的其他時(shí)段數(shù)據(jù);選擇的模型驗(yàn)證時(shí)段為2009年8月30日至10月30日.降水入滲為主要補(bǔ)給來源，人工開采為主要排泄方式.潛水含水層絕對(duì)誤差的平均值為0.017 m，其中90%的觀測(cè)井?dāng)M合后絕對(duì)誤差小于0.032 m.通過模型識(shí)別和驗(yàn)證，所映實(shí)際流場(chǎng)的特點(diǎn)［2］.為此，首先對(duì)數(shù)學(xué)模型中的數(shù)學(xué)方程及邊界條件進(jìn)行校正，以確定是否能確切地反映實(shí)際水文地質(zhì)條件.模型識(shí)別時(shí)段選取2009年3月30日至6月30日.隨著春灌時(shí)節(jié)開采量持續(xù)增加，地下水水位會(huì)持續(xù)下降.通過持續(xù)觀測(cè)抽水井，可以獲得實(shí)際水位的變化數(shù)據(jù)，為下一步的參數(shù)識(shí)別提供依據(jù).

對(duì)模型采用PCG2(預(yù)調(diào)共軛梯度法)進(jìn)行求解，得到在給定水文地質(zhì)參數(shù)和各均衡項(xiàng)條件下地下水位的時(shí)空分布［3］.通過實(shí)際水位與計(jì)算水位的擬合分析［4］，反復(fù)地修改參數(shù)，直到達(dá)到文獻(xiàn)［5］中規(guī)定的“水位變化值較小(＜5 m)的情況下，水位擬合誤差一般應(yīng)小于0.5 m”時(shí)，認(rèn)為此時(shí)的參數(shù)值為含水層的實(shí)際參數(shù).由于是以實(shí)際的水文地質(zhì)條件作為依據(jù)劃分參數(shù)分區(qū)、選取參數(shù)初值，加上多次調(diào)參擬合，模型識(shí)別結(jié)果較為滿意.在模擬區(qū)分散選取4 個(gè)觀測(cè)井(S42，S73，S53，S61)，做出地下水位觀測(cè)值與計(jì)算值擬合圖，如圖2所示.絕對(duì)誤差的平均值為0.023 m，其中所有觀測(cè)井?dāng)M合絕對(duì)誤差均小于0.043 m.可見實(shí)測(cè)水位與計(jì)算水位相關(guān)程度較高，擬合效果較好［6］，說明建立的水文地質(zhì)概念模型以及數(shù)學(xué)模型是正確的.表2列出了具體的參數(shù)識(shí)別結(jié)果.建立的數(shù)學(xué)模型、選取的邊界條件及最終確定的水文地質(zhì)參數(shù)均正確、可靠，地下水流系統(tǒng)的預(yù)報(bào)可以運(yùn)用該模型進(jìn)行.

2 模型的應(yīng)用

在模型的應(yīng)用過程中，認(rèn)定開采條件不變，采用Monte-Carxo(蒙特卡羅)方法預(yù)報(bào)研究區(qū)內(nèi)的降水量，進(jìn)行隨機(jī)模擬和預(yù)報(bào).

1)預(yù)測(cè)至2020年，降落漏斗的變化情況.結(jié)果表明，在保持開采量不變，隨機(jī)預(yù)報(bào)降水量情況下，降落漏斗持續(xù)增大，漏斗面積由18.59 km2增大至23.31 km2，如圖3所示.漏斗中心一直在磁州鎮(zhèn)內(nèi).

2)為了控制降落漏斗，須規(guī)劃研究區(qū)工農(nóng)業(yè)及城鎮(zhèn)生活用水.隨著南水北調(diào)中線工程的逐步實(shí)施，每年可為研究區(qū)補(bǔ)給水量945×104m3，再加上合理的工農(nóng)業(yè)節(jié)水措施，可以做到減少地下水開采量，將開采量由現(xiàn)在的6 032×104m3減少至2020年的4 848×104m3.漏斗中心水位將回升到模擬初始狀態(tài)，影響范圍基本未擴(kuò)大，如圖4所示.

圖3 開采條件不變時(shí)的降落漏斗模擬圖

3 結(jié)語(yǔ)

1)數(shù)值模擬結(jié)果表明，通過合理選取水文地質(zhì)參數(shù)以及源匯項(xiàng)，建立了合理的水文地質(zhì)概念模型，并最終建立了正確的數(shù)學(xué)模型.研究區(qū)地下水的實(shí)際情況可以得到正確的反映，且溶質(zhì)運(yùn)移模擬也可以應(yīng)用該模型.

2)采取有效的節(jié)水措施，可以使地下水位有所回升，控制降落漏斗影響范圍的擴(kuò)大.

［1］張兆吉，費(fèi)宇紅，陳宗宇，等.華北平原地下水可持續(xù)利用調(diào)查評(píng)價(jià)［M］.北京:地質(zhì)出版社，2009.

圖4 采取有效措施后降落漏斗模擬圖

［2］李平，盧文喜，馬洪云，等.Visual MODFLOW在地下水?dāng)?shù)值模擬中的應(yīng)用——以公主嶺市黃龍工業(yè)園水源地為例［J］.工程勘察，2006(3):24－28.

［3］San Juan Canna，Kolm Kenneth E.Concept ualization，characterization and numerical modeling of the Jackson Hole alluvial aquifer using ARC/INFO and MODFLOW［J］.Engineering Geology，1996，42(2 －3):119 －137.

［4］Peeve A S，Wariochaa J.Regional groundwater flow modeling of the Glacial Lake Agassiz Peatlands，Minnesota［J］.Journal of Hydrology，2001，243(1 －2):94 －100.

［5］中國(guó)地質(zhì)調(diào)查局.GWID1地下水流數(shù)值模擬技術(shù)要求［S］.北京:地質(zhì)出版社，2004.

［6］康敏娟.基于Visual MODFLOW的邯鄲縣淺層地下水量評(píng)價(jià)研究［D］.邯鄲:河北工程大學(xué)，2010.

The Application of MODFLOW in Numerical Simulation of Groundwater in Cixian

LI Shu-wen1，2，YU Hu-guang1
(1.Hebei University of Engineering，Handan 056038，China;2.Hebei Key Laboratory of Resources Survey and Research，Handan 056038，China)

At present，the water environment in the east plain of Cixian is deteriorating due to the irrational use of water and the over exploitation of the groundwater.Through the analysis of the study area’s hydrogeological condition，the shallow groundwater is simulated by using MODFLOW.The simulation results show that the cone of depression will increase from 18.59 km2to 3.31 km2until 2020 under the invariable working conditions.The paper proposes that the expansion of depression cone can be avoided by rationally modifying and optimizing the groundwater exploitation plan.

groundwater;numerical simulation;MODFLOW

1002－5634(2012)02－0053－03

2011－12－22

李樹文(1955—)，男，河北邯鄲人，教授，博士，主要從事水文水資源方面的研究.

(責(zé)任編輯:喬翠平)