基于Visual MODFLOW的地下水?dāng)?shù)值模擬研究及水位預(yù)測

韓琳，顏翠翠，趙振偉，吳婧，馬鵬飛

(1.山東省地質(zhì)礦產(chǎn)勘查開發(fā)局八〇一水文地質(zhì)工程地質(zhì)大隊(duì)，山東 濟(jì)南 250014；2.山東省地勘局第二水文地質(zhì)工程地質(zhì)大隊(duì)(山東省魯北地質(zhì)工程勘察院)，山東 德州 253072)

0 引言

水是人類社會(huì)賴以生存和發(fā)展的不可代替的但可調(diào)控的重要自然資源，是構(gòu)成區(qū)域生態(tài)環(huán)境系統(tǒng)的基本要素[1-2]。目前，濟(jì)南市城市生活用水由地下水和地表水組成，占生活用水絕大部分的地表水的使用不僅要經(jīng)過多道工序，而且水質(zhì)、口感也與地下水有一定差距[3]。近年來，市政府下決心解決居民飲用優(yōu)質(zhì)地下水問題，已開發(fā)長孝水源地，設(shè)計(jì)供水能力10萬m3/d，目前已實(shí)現(xiàn)供水4萬m3/d。因此對(duì)長孝地區(qū)的水文地質(zhì)條件開展研究及時(shí)且必要，但在地下水的開采使用過程中，往往存在對(duì)地下水分布與動(dòng)態(tài)變化了解不清的盲目開采，因此了解地下水的分布與動(dòng)態(tài)變化就顯得格外重要。但地下水的儲(chǔ)藏條件特殊，常規(guī)的探測方法不能及時(shí)全面地了解地下水的賦存狀態(tài)[4-5]，而數(shù)值模擬方法因其直觀全面的優(yōu)點(diǎn)，非常適合用于地下水的研究中。地下水?dāng)?shù)值模擬的研究方面，不少研究者展開了深入研究。于翠翠等[6]應(yīng)用GMS軟件三維地下水流數(shù)值模擬模型，對(duì)濟(jì)南泉域內(nèi)巖溶地下水進(jìn)行數(shù)值模擬和水平衡分析，并應(yīng)用時(shí)間序列分析法對(duì)泉水水位動(dòng)態(tài)進(jìn)行了預(yù)測。朱君等[7]利用FEFLOW軟件模擬計(jì)算了丘陵山區(qū)地下水流動(dòng)特征下氚的遷移規(guī)律。不少研究者就地下水?dāng)?shù)值模擬方法進(jìn)行了系統(tǒng)性總結(jié)，分析了各種軟件的適用條件[8-15]。Visual MODFLOW軟件在水位預(yù)測方面應(yīng)用廣泛，不少學(xué)者應(yīng)用該軟件對(duì)不同地點(diǎn)的水位進(jìn)行了預(yù)測，預(yù)測結(jié)果達(dá)到了工程要求[16-20]。本文基于以上研究利用Visual MODFLOW軟件對(duì)模型進(jìn)行數(shù)值模擬計(jì)算，展開地下水水位預(yù)測相關(guān)研究。

1 研究區(qū)范圍及水文地質(zhì)背景

1.1 研究區(qū)邊界條件

濟(jì)南地區(qū)在大地構(gòu)造上橫跨新華夏第二隆起帶的魯西隆起及新華夏第二沉降帶的魯西北坳陷區(qū)。南部是北傾單斜構(gòu)造，北部是松散堆積物的凹陷區(qū)。南部單斜構(gòu)造中發(fā)育有多條規(guī)模較大的NNW向斷裂，如牛角店斷裂、馬山斷裂、千佛山斷裂、文化橋斷裂、東塢斷裂、文祖斷裂等,這些NNW向斷裂自東向西大致等距分布，將單斜構(gòu)造分割為若干個(gè)斷塊，這就形成了濟(jì)南地區(qū)地質(zhì)構(gòu)造的基本特征。根據(jù)區(qū)內(nèi)斷裂構(gòu)造的水理性質(zhì)將濟(jì)南地區(qū)巖溶地下水劃分為5個(gè)水文地質(zhì)單元(圖1)。

1—第四系；2—石炭系；3—奧陶系；4—寒武系；5—泰山巖群；6—侵入巖；7—巖脈；8—奧陶系頂板埋深等值線(m)；9—水文地質(zhì)單元分區(qū)界線；10—研究區(qū)范圍界線圖1 濟(jì)南地區(qū)水文地質(zhì)單元分區(qū)略圖

長清-孝里鋪水文地質(zhì)單元是濟(jì)南地區(qū)五大水文地質(zhì)單元之一，面積935.40km2，各地質(zhì)單元的面積見表1。市政府為解決居民飲用優(yōu)質(zhì)地下水問題開發(fā)的長孝水源地就位于長清-孝里鋪水文地質(zhì)單元內(nèi)。長清-孝里鋪水文地質(zhì)單元邊界條件及與周邊水文地質(zhì)單元水力聯(lián)系如下：

表1 濟(jì)南水文地質(zhì)單元分區(qū)及分區(qū)面積

東部邊界：馬山斷裂是研究區(qū)內(nèi)規(guī)模較大的NW向斷層，兩盤地層相對(duì)錯(cuò)動(dòng)較大，形成了2個(gè)斷塊地質(zhì)體，兩斷塊在新周莊以南具隔水性質(zhì)，分割為相互無水力聯(lián)系的不同水文地質(zhì)單元；斷層在新周莊至老屯段兩盤皆為灰?guī)r，斷層為弱透水性質(zhì)，認(rèn)為具備作為水文地質(zhì)邊界的條件；斷層在老屯以北至前隆段為透水性質(zhì)，但其影響是局部且有限的。因此，把馬山斷層作為長清-孝里鋪水文地質(zhì)單元與東部濟(jì)南泉域水文地質(zhì)單元的邊界。

南部邊界：本區(qū)的南部和東南部是地表分水嶺。在南部，地表分水嶺西起黃山巖脈，經(jīng)賢子峪，冷飯店、興隆鎮(zhèn)南，玉皇頂、木閣寨北等制高點(diǎn)至雙山東南出研究區(qū)。巖溶水南部補(bǔ)給區(qū)以崗辛起向東沿黃石崖、桃花峪分水嶺為界。根據(jù)水文地質(zhì)調(diào)查證實(shí)用地表分水嶺近似代替地下分水嶺是合理的，構(gòu)成長清-孝里鋪水文地質(zhì)單元南部的地下水分水嶺邊界。

西部邊界：牛角店斷裂為一NW向推測斷裂，東南端起自平陰縣安城鄉(xiāng)大官莊，向東北經(jīng)欒灣后被第四系覆蓋，繼續(xù)向北延伸。已有地質(zhì)勘探表明牛角店斷裂(黃河以北部分)為透水?dāng)嗔?，故將牛角店斷?黃河以北部分)作為本區(qū)西部的透水邊界；黃山巖脈(主支)及其分支阻水性質(zhì)極為明顯，為一阻水邊界?；谝陨戏治龉蕦⑴＝堑陻嗔?黃河以北段)及黃山巖脈作為長清-孝里鋪水文地質(zhì)單元的西邊界。

北部邊界：東阿斷裂為長清-孝里鋪水文地質(zhì)單元北部邊界，西自東阿縣高集鎮(zhèn)，向東北經(jīng)齊河縣趙官鎮(zhèn)、胡官屯鎮(zhèn)東部出境。長清-孝里鋪水文地質(zhì)單元邊界條件統(tǒng)計(jì)情況見表2。

表2 長清-孝里鋪水文地質(zhì)單元邊界條件統(tǒng)計(jì)

1.2 水文地質(zhì)背景

1.2.1 水位動(dòng)態(tài)特征

基于2007—2013年曹樓村(CK1)和西關(guān)村(CK15)的鉆孔地下水水位多年動(dòng)態(tài)變化曲線，可以發(fā)現(xiàn)，地下水水位對(duì)降水量很敏感，并隨降水量的“多—少”呈現(xiàn)“高—低”變化(圖2)。根據(jù)圖2可將區(qū)內(nèi)巖溶地下水水位多年動(dòng)態(tài)變化可分為3個(gè)階段：①2007年1月—2009年12月，水位先下降后上升穩(wěn)定階段；②2010年1月—2011年7月，水位先下降后上升再逐漸下降階段；③2011年8月—2013年12月，水位相對(duì)穩(wěn)定階段。同時(shí)發(fā)現(xiàn)地下水水位在枯水期水位總體下降，在豐水期水位總體有所回升。另外CK1、CX15兩鉆孔點(diǎn)自2007年1月至2013年12月，水位分別上升了0.19m和0.64m，水位呈上升趨勢。

1—降水量；2—CK1水位標(biāo)高；3—CX15水位標(biāo)高圖2 長清-孝里鋪水文地質(zhì)單元巖溶水水位及降水量變化關(guān)系曲線圖

1.2.2 巖溶發(fā)育特征及分布規(guī)律

(1)巖溶發(fā)育具有層控性

碳酸鹽巖是研究區(qū)內(nèi)主要的巖石，具有分布廣、厚度大的特點(diǎn)，但各層組碳酸鹽巖的生成時(shí)代，沉積環(huán)境不同，它們的化學(xué)成分，礦物成分及結(jié)構(gòu)構(gòu)造均有一定的差別，故其溶蝕特點(diǎn)及巖溶分布發(fā)育程度和特點(diǎn)不盡相同。巖溶最易發(fā)育的層位是寒武-奧陶紀(jì)三山子組和奧陶紀(jì)馬家溝群二段、四段和六段和寒武紀(jì)中統(tǒng)張夏組，巖溶是易發(fā)育的巖性是泥晶灰?guī)r和豹斑灰?guī)r，其次是鮞狀灰?guī)r類，第三為白云巖類。由于巖性與層位密切相關(guān)，巖溶發(fā)育具有一定的層狀特點(diǎn)。例如張夏組基本上是巨厚層鮞狀灰?guī)r組成，成為單獨(dú)的巖溶發(fā)育的層狀溶蝕孔洞-溶隙網(wǎng)絡(luò)系統(tǒng)，而寒武-奧陶紀(jì)三山子組和奧陶紀(jì)馬家溝群，巖性組成比較復(fù)雜多樣，雖在巖溶發(fā)育上相互聯(lián)系，但各層均有自己的特點(diǎn)，如其中泥晶和豹斑灰?guī)r以相對(duì)均勻性略差的寬溶隙系統(tǒng)為主，白云巖類及角礫灰?guī)r類則以較均勻的層狀溶孔及小型孔洞為主，其間又間隔一些巖溶不甚發(fā)育的層位。

(2)巖溶發(fā)育的分布特征

①水平分布特征。水平分布最主要的特征是由補(bǔ)給區(qū)到徑流排泄區(qū)，巖溶發(fā)育呈現(xiàn)由弱到強(qiáng)的變化規(guī)律。補(bǔ)給區(qū)巖溶發(fā)育差，為弱巖溶發(fā)育帶；補(bǔ)給徑流區(qū)巖溶發(fā)育，為中等巖溶發(fā)育帶；在徑流排泄區(qū)，由于地下徑流長期作用、構(gòu)造發(fā)育，巖溶十分發(fā)育，往往能形成具有開采價(jià)值的水源地。

②垂直分布特征。前面已經(jīng)提到不同巖層的生成時(shí)代，沉積環(huán)境不同，它們的化學(xué)成分，礦物成分及結(jié)構(gòu)構(gòu)造均有一定的差別，故將巖溶的垂直分布特征分述如下：

奧陶紀(jì)馬家溝群裂隙、巖溶垂直分布特征：長孝水文地質(zhì)單元發(fā)育深度一般在30～340m之間，其中長孝水文地質(zhì)單元140～230m這一深度區(qū)間發(fā)育最普遍，巖溶形態(tài)多為溶孔和溶蝕裂隙(圖3)。

圖3 長清-孝里鋪水文地質(zhì)單元奧陶紀(jì)馬家溝群裂隙-巖溶發(fā)育深度圖

寒武-奧陶紀(jì)三山子組裂隙-巖溶垂直分布特征：長孝水文地質(zhì)單元裂隙、巖溶的集中發(fā)育深度一般在30～100m和130～200m兩個(gè)區(qū)段，尤其在50～150m這一深度區(qū)間，發(fā)育最普遍。

寒武紀(jì)中統(tǒng)張夏組裂隙-巖溶垂直分布特征：長孝水文地質(zhì)單元裂隙、巖溶的集中發(fā)育深度一般在0～50m和270～320m兩個(gè)區(qū)段。

2 水文地質(zhì)模型概化

根據(jù)地質(zhì)資料，在考慮最大程度擬合實(shí)際水文地質(zhì)情況下，建立了長清-孝里鋪水文地質(zhì)單元的概念模型及結(jié)構(gòu)模型，利用Visual MODFLOW軟件對(duì)模型進(jìn)行數(shù)值計(jì)算，并根據(jù)觀測資料對(duì)模型進(jìn)行了為期3年的識(shí)別與驗(yàn)證，保證了模型的有效性，隨即展開了不同開采條件下的水位預(yù)測研究。

2.1 模擬范圍

基于前面已詳細(xì)描述的長清-孝里鋪水文地質(zhì)單元的邊界條件，結(jié)合前人的研究成果，表明長清-孝里鋪水文地質(zhì)單元基本為一具有系統(tǒng)補(bǔ)給、徑流、排泄條件的水文地質(zhì)單元。故將模擬區(qū)總面積確定為780km2，占整個(gè)地質(zhì)單元的83%(圖4)。

1—模擬范圍；2—透水?dāng)鄬樱?—弱透水?dāng)鄬樱?—阻水?dāng)鄬樱?—巖脈；6—一、二層與三層分界線；7—斷層及推測斷層界線圖4 長清-孝里鋪水文地質(zhì)單元模擬范圍示意圖

2.2 概念模型

水文地質(zhì)資料表明研究區(qū)主要有4個(gè)含水巖組：松散巖類孔隙含水巖組、碳酸鹽巖類裂隙巖溶含水巖組、碎屑巖夾碳酸鹽巖類巖溶裂隙含水巖組、塊狀巖類裂隙含水巖組。不同含水巖組的物理性質(zhì)存在差異，根據(jù)研究的重點(diǎn)及目的，將模擬區(qū)概化為三層：第一層為孔隙水含水層，為潛水，含水層總厚度10～40m；第二層弱透水層，最大厚度約200m；第三層為巖溶水含水層，具承壓性，最大深度為500m左右。

2.3 邊界條件概化

2.3.1 孔隙水含水層邊界概化

根據(jù)前面的地質(zhì)邊界條件可將東部邊界南段設(shè)為隔水邊界，南部、西部、北段邊界設(shè)為流量邊界。

2.3.2 巖溶水含水層邊界條件概化

東部邊界以馬山斷裂為界，南端為阻水邊界，北段為流量邊界。西部邊界以西部北端的牛角店斷裂為透水邊界，黃山巖脈為隔水邊界。南部以地下水(地表)分水嶺為邊界，即零通量邊界，地下分水嶺與地表分水嶺基本一致。北部邊界以奧陶紀(jì)灰?guī)r頂板埋深600m等值線為界，推測巖溶不發(fā)育，可視為阻水邊界。

2.3.3 垂向邊界

孔隙水含水層自由水面為系統(tǒng)的上邊界。通過該邊界，潛水與系統(tǒng)外界發(fā)生垂向水量交換?？紫端畬雍蛶r溶水含水層通過越流進(jìn)行水量交換，越流量由淺、深層的水頭差、弱透水層在垂向上的滲透系數(shù)和厚度決定。巖溶水底板作為系統(tǒng)的下邊界，底板以下巖層為不透水地層，定義為隔水邊界。

3 模型的建立及求解

根據(jù)對(duì)勘探區(qū)水文地質(zhì)條件的概化，建立數(shù)學(xué)模型，采用有限差分法進(jìn)行不同開采條件下的巖溶水開采方案計(jì)算，并進(jìn)行地下水位預(yù)報(bào)。

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

地下水流系統(tǒng)具有非均質(zhì)性、各向同性、空間三維結(jié)構(gòu)、非穩(wěn)定性，可用如下方程的定解問題來描述：

式中：Kx，Ky和Kz為x，y，z方向滲透系數(shù)(m/d)；h為水位標(biāo)高(m)；W為含水層的源匯項(xiàng)(1/d)；S為含水介質(zhì)的儲(chǔ)水率(1/m)；h0為初始水位(m)；Γ2為二類邊界；Kn為邊界面法向方向的滲透系數(shù)(m/d)；n為Γ2邊界的外法線方向；q(x,y,z,t)為二類邊界上已知流量函數(shù)，流入為正，流出為負(fù)，隔水邊界為0(m/d)；Ω為滲流區(qū)域。

3.2 數(shù)值模擬模型

根據(jù)模擬區(qū)的含水層結(jié)構(gòu)特征、邊界條件和地下水流場等，在垂向上分為三層，其中第一層為潛水含水層，二、三層為承壓含水層。在平面上剖分成大小為250m×250m的正方形單元，各層均剖分為172行，152列，其中一層有效單元格為6878個(gè)，二層有效單元格為6932個(gè)，三層有效單元格為12480個(gè)。網(wǎng)格剖分示意圖見圖5。

圖5 含水層模型剖分示意圖(左：孔隙水；右：巖溶水)

3.3 源匯項(xiàng)的處理

模擬區(qū)源匯項(xiàng)有點(diǎn)、線、面3種要素。點(diǎn)狀要素主要代表水廠、自備井、工農(nóng)業(yè)開采等源匯項(xiàng)；線狀要素主要代表河流(通過River模塊處理)、引水渠、灌渠等補(bǔ)給項(xiàng)；面狀要素由Recharge模塊給出，主要代表降雨入滲、灌溉入滲等補(bǔ)給項(xiàng)。其余源匯項(xiàng)資料通過Well模塊帶入模型計(jì)算，均處理成開采或者補(bǔ)給井。

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

為保證已建立模型的有效性，模型的識(shí)別與驗(yàn)證是必不可少的。本研究進(jìn)行了充分的模型識(shí)別與驗(yàn)證。識(shí)別期(2010.5.31—2011.6.11)；驗(yàn)證期(2011.6.2—2013.5.26)。模型的識(shí)別與檢驗(yàn)過程是整個(gè)模擬中極為重要的一步工作，本文采用試估—校正法進(jìn)行了反復(fù)修改參數(shù)和調(diào)整某些源匯項(xiàng)達(dá)到了較為理想的擬合結(jié)果。經(jīng)過多次調(diào)整參數(shù)擬合，得到了3層含水層的水文地質(zhì)參數(shù)。最終的模擬流場與實(shí)際流場總體流動(dòng)方向相同，等值線密間距接近，初步的擬合效果較好。

考慮模擬區(qū)水位觀測點(diǎn)的分布情況及觀測資料的連續(xù)性，在第一含水層選擇了5個(gè)觀測點(diǎn)、第三含水層選擇了10個(gè)觀測點(diǎn)，點(diǎn)位具體分布情況見圖6。進(jìn)行水位動(dòng)態(tài)過程曲線擬合，計(jì)算各擬合點(diǎn)的計(jì)算水位與觀測水位的絕對(duì)誤差和相對(duì)誤差，并對(duì)其進(jìn)行綜合分析，部分觀測孔過程線擬合情況見圖7，觀測值和模擬值的變化趨勢基本一致，2條曲線的偏離程度較低，曲線真實(shí)的反映了豐枯水期地下水水位的動(dòng)態(tài)變化，基本反映了地下水系統(tǒng)的水力特征，模型的有效性得到了驗(yàn)證。

1—巖脈；2—透水?dāng)鄬樱?—弱透水?dāng)鄬樱?—阻水?dāng)鄬樱?—斷層及推測斷層界線；6—孔隙水觀測井；7—巖溶水觀測井圖6 地下水位監(jiān)測孔分布示意圖

5 開采方案模擬與預(yù)測分析

5.1 開采方案模擬

在已建立的模型上，進(jìn)行不同開采方案的開采動(dòng)態(tài)響應(yīng)研究，開采的目的層位為巖溶含水層，在研究區(qū)內(nèi)選擇了一處富水地段，根據(jù)原有地質(zhì)資料顯示的巖溶水儲(chǔ)量情況確定了2種開采方案：

表3 開采方案具體實(shí)施方案

在進(jìn)行多方案預(yù)測之前，進(jìn)行了自然均衡狀態(tài)下曹樓抽水試驗(yàn)和相同地點(diǎn)的少量增采試驗(yàn)，對(duì)兩者的地下水流場分析，證明了長孝水文地質(zhì)單元有一定增采潛力。在將每年分為枯水期與豐水期2個(gè)應(yīng)力期的基礎(chǔ)上，分別預(yù)測分析了2個(gè)方案運(yùn)行10年后地下水流場變化情況。預(yù)測期間，2個(gè)方案降水量均采用1995—2004年枯水系列年資料，模型中以559m平均降雨量進(jìn)行預(yù)測。

5.2 開采方案對(duì)地下水位的影響

運(yùn)行方案1后，地下水位產(chǎn)生明顯變化，在布井地段出現(xiàn)了一個(gè)明顯的地下水降落漏斗，在曹樓附近出現(xiàn)最大降深約9.02m，年降幅為0.902m；運(yùn)行方案2后，地下水位降落漏斗進(jìn)一步擴(kuò)大，曹樓最大降深為12.25m，年降幅達(dá)1.225(圖8)。隨著開采量增加，巖溶水得不到有效補(bǔ)給，導(dǎo)致降落漏斗增大，水位降深加大。運(yùn)行兩方案后，由模擬區(qū)馬山斷裂長清北排泄區(qū)觀測孔水位變化情況可以看出，觀測孔水位有持續(xù)下降趨勢，雖然下降幅度不大，但改變了北邊界的補(bǔ)排方向。增加的開采方案不僅增加了地下水補(bǔ)給量，同時(shí)也減少了排泄量減少，但補(bǔ)給增量是有限度的，故導(dǎo)致了地下水水位的降低。

1—觀測值；2—模擬值圖7 觀測孔水位過程線擬合圖(左：識(shí)別期；右：驗(yàn)證期)

圖8 巖溶水預(yù)測流場示意圖(左：方案1；右：方案2)

通過數(shù)值法對(duì)巖溶水開采方案模擬分析表明，實(shí)施開采方案后，地下水水位變化明顯，在開采井附近產(chǎn)生地下水降落漏斗，對(duì)水文地質(zhì)單元的地下水位也有一定影響。枯水系列年的2種開采方案改變了地下水均衡狀態(tài)。自然情況下，馬山斷裂為巖溶水的排泄邊界，而實(shí)施開采方案后，激化增加補(bǔ)給量，襲奪部分邊界排泄量，但激化與襲奪總量是有限度的，結(jié)合該區(qū)的實(shí)際情況，各方案說明在連續(xù)枯水系列年條件下，長清地區(qū)開采量不宜超過14萬m3/d(含分散開采)。若考慮東邊界以東外圍水源地的影響，應(yīng)建立長孝-濟(jì)南數(shù)值模型，統(tǒng)一進(jìn)行地下水資源分析與預(yù)測，才能進(jìn)行精度更高的模型研究。

6 結(jié)論與思考

本文以濟(jì)南地區(qū)的長清-孝里鋪水文地質(zhì)單元為研究對(duì)象，根據(jù)該區(qū)的水文地質(zhì)條件確定了數(shù)值模型研究的范圍，然后對(duì)模型進(jìn)行了區(qū)分并進(jìn)行了邊界條件概化，建立了用于地下水位預(yù)測的數(shù)值模型。經(jīng)過為期3年的模型識(shí)別與驗(yàn)證，保證的模型的有效性。隨后開展了兩種不同巖溶水開采方案下的水位預(yù)測數(shù)值模擬得出了以下結(jié)論：

(1)本文建立的地質(zhì)結(jié)構(gòu)影響的地下水預(yù)測模型是可行的，可為相關(guān)領(lǐng)域的研究提供一種有效的分析方法。

(2)在連續(xù)枯水系列年條件下，結(jié)合該區(qū)的實(shí)際情況，2種開采方案預(yù)測說明，長清地區(qū)巖溶水開采量不宜超過14萬m3/d，為后期的水源地建設(shè)提供了可靠依據(jù)。

(3)在考慮相鄰?fù)鈬吹氐挠绊憰r(shí)，應(yīng)建立范圍更廣的數(shù)值模型，統(tǒng)一進(jìn)行地下水資源分析與預(yù)測，數(shù)值模擬的結(jié)果才能與實(shí)際情況更加接近。