沈媛媛,辛寶東,郭高軒,紀(jì)軼群

(北京市水文地質(zhì)工程地質(zhì)大隊,北京 100195)

北京市人均水資源量不足300m3,水資源供需矛盾突出.北京市為應(yīng)對缺水的緊張局面,2002年以后相繼建設(shè)了5處應(yīng)急水源工程.其中房山張坊水源地,以巖溶水為應(yīng)急水源.本文在詳細(xì)分析研究區(qū)地質(zhì)及水文地質(zhì)條件的基礎(chǔ)上,建立研究區(qū)水文地質(zhì)概念模型和地下水流數(shù)值模型,采用MODFLOW求解,完成模型識別和驗證.并預(yù)測給定條件下水源地區(qū)地下水系統(tǒng)變化情況,為其合理開發(fā)利用提供科學(xué)依據(jù).

1 研究區(qū)地質(zhì)、水文地質(zhì)概況

1.1 地質(zhì)構(gòu)造

研究區(qū)在構(gòu)造單元上位于霞云嶺-龍門臺復(fù)向斜南翼,其東部為北嶺向斜,北部為大石河背斜(圖1).核部基本位于拒馬河與大石河分水嶺處,霧迷山組地層構(gòu)成向斜兩翼,鐵嶺組、青白口系、寒武系在向斜核部形成一些列寬緩波狀褶皺.斷裂構(gòu)造主要有北東向的黃莊-高麗營斷裂,牛口峪-長溝斷裂和霞云嶺-大峪溝斷裂.

1.2 水文地質(zhì)條件

研究區(qū)基巖含水層以霧迷山組白云巖分布最廣,在山前地帶有第四系砂礫石層,為較好的巖溶水含水層,單井涌水量可達(dá)3000m3/d.研究區(qū)北部還有鐵嶺組和洪水莊組巖溶含水層,其富水性較差.

分布較廣的的碳酸鹽巖裂隙發(fā)育,有利于降水直接入滲.地下水自西北向東南方向徑流,地下水位埋由山區(qū)到山前深逐漸減小.山前的地下水主要以泉水、側(cè)向徑流和人工開采.

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

本地區(qū)巖溶含水層雖然溶蝕裂隙有統(tǒng)一的地下水面,地下水系統(tǒng)可用多孔介質(zhì)滲流理論模型來概化.

(1)含水層概化

根據(jù)研究區(qū)的水文地質(zhì)條件確定了模型邊界,模擬面積為401km2.模擬對象為薊縣系霧迷山和鐵嶺組之碳酸鹽巖含水層.其巖溶含水層埋深400m左右,可將其概化為非均質(zhì)各項同性含水層.

(2)含水層概化邊界條件

模擬含水層北部邊界為拒馬河與大石河的地表分水嶺,概化為流量邊界;東部為青白口系下馬嶺組板巖、千枚巖,為相對阻水層,概化為隔水邊界;南部與河北省相鄰,大峪溝斷裂成為地下水側(cè)向徑流的通道;另外,地下水可通過深層徑流向下游排泄,概化為流出邊界;西部以煌斑巖脈為界,概化為隔水邊界(圖2).

模型的上邊界為潛水水面,巖溶地下水系統(tǒng)通過該邊界與外界降雨入滲、地表水體滲漏等相關(guān)系.以含水層埋深400m處為模型底部概化為隔水邊界.

(3)地下水的水力特征

研究區(qū)含水層厚度大、分布廣.地下水以水平流動為主、垂向流動為輔,可將地下水流動概化為空間三維非穩(wěn)定流.

3 巖溶地下水系數(shù)學(xué)模型及其計算結(jié)果

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

概化后的巖溶地下水系統(tǒng)用如下數(shù)學(xué)公式表示[11]:

式中:Ω為滲流區(qū)域;H為巖溶含水層的水位標(biāo)高(m);B為模型底界標(biāo)高(m);Kx、Ky、Kz、Kn分別為x、y、z、邊界法向方向的滲透系數(shù)(m/d);μ為潛水含水層的給水度;w為含水層的源匯項(m/d);H0為含水層的初始水位分布(m);Γ0為滲流區(qū)域的上邊界,即地下水的自由表面;Γ2為滲流區(qū)域的側(cè)向邊界;Γ4為滲流區(qū)域的下邊界,即含水層底部的隔水邊界;n為邊界面的法線方向; q(x,y,z,t)為二類邊界的流量(m/d),設(shè)計流入為正,流出為負(fù),隔水邊界值為0.

采用美國地質(zhì)調(diào)查局發(fā)布的三維地下水有限差分計算軟件MODFLOW進(jìn)行地下水?dāng)?shù)值模擬[12],將模擬區(qū)剖分為250mX250m的矩形網(wǎng)格,共剖分為120行、140列,有效單元格6416個.

模擬區(qū)分為6個滲透系數(shù)參數(shù)分區(qū)(圖2),其參數(shù)初值根據(jù)含水層巖性及以往研究資料設(shè)定.給水度根據(jù)以往試驗結(jié)果取為0.01.根據(jù)不同的巖性及以往研究成果,確定各巖組巖層的降雨入滲系數(shù):霧迷山組為0.35;鐵嶺組為0.3,洪水莊組及下馬嶺組巖層為0.1.用達(dá)西公式計算其邊界側(cè)向徑流量.

在張坊以北拒馬河水接受地下水補(bǔ)給,在張坊以南變?yōu)楹铀畮аa(bǔ)給地下水.采用MODFLOW河流計算子程序包,根據(jù)其河水位、底板高程、河流附近地下水位以及河道滲透系數(shù),計算河流與含水層的交換水量.

泉水流量大小與附近地下水位呈正比,并與含水層滲透性能有關(guān).采用與河流類似的模擬方式,以泉附近含水層滲透系數(shù)、地下水位與地表高程差值確定泉流量.

(2)模型計算結(jié)果

選取2004年6月統(tǒng)測地下水位作為模擬的初始地下水位,12月為模型識別期;2005年1月至6月為模型驗證期.采用有限差分法求解,以一個月為一個時間段,每個時間段分為10個時間步長,嚴(yán)格控制每次迭代的誤差.

選用"試錯法"調(diào)整模型參數(shù),觀測井識別期和驗證期的模擬水位與實測水位殘差較小,給出張坊、長溝、廣潤莊和廣祿莊四處的觀測井地下水位擬合曲線對比圖(圖3).從圖中可以看出,除張坊觀測井的水位在個別應(yīng)力期外,其他均擬合較好,滿足了精度要求,建立的模型較真實地刻畫了其研究區(qū)地下水系統(tǒng)特征,可以用于地下水預(yù)測評價.

4 應(yīng)急水源地開采對河流滲漏的影響

(1)應(yīng)急開采前地下水系統(tǒng)均衡情況

模擬期內(nèi)含水層系統(tǒng)水量均衡結(jié)果見表1,表中補(bǔ)給量包括降雨入滲補(bǔ)給量和邊界流入量以及灌溉滲漏量.其中降雨入滲量占總補(bǔ)給量的89.1%,是其最重要的補(bǔ)給源,邊界流入量占總補(bǔ)給量的9.0%,灌溉滲漏補(bǔ)給量為1.5%.區(qū)域地下水主要以工農(nóng)業(yè)開采、河流和泉水排泄以及邊界側(cè)向流出.其中工農(nóng)業(yè)開采占總排泄量的26.0%,泉排泄38.4%,邊界流出量33.3%.

模擬結(jié)果表明,在房山應(yīng)急水源地開采前,地下水向河流排泄量為139.4X104m3/a.

表1 模擬期內(nèi)地下水系統(tǒng)均衡結(jié)果表(X104m3/a)

(2)應(yīng)急開采后地下水系統(tǒng)均衡情況

為評價應(yīng)急水源地開采后的地下水系統(tǒng)與河流之補(bǔ)排關(guān)系,選擇2009年8月實測的地下水流場作為初始流場,預(yù)測枯水年(年均降雨量513mm)和應(yīng)急水源地現(xiàn)狀開采條件下、5年后的地下水系統(tǒng)均衡結(jié)果見表2.

表2 應(yīng)急水源地開采條件下地下水系統(tǒng)均衡結(jié)果表(X104m3/a)

預(yù)測結(jié)果,應(yīng)急水源地開采后,河流與地下水系統(tǒng)總的補(bǔ)排關(guān)系由應(yīng)急水源地開采前地下水補(bǔ)給河流轉(zhuǎn)變?yōu)楹恿鳚B漏補(bǔ)給地下水,河流滲漏補(bǔ)給量為398.0X104m3/a,占總補(bǔ)給量的5.1%.河流與地下水系統(tǒng)的交換量的變化量為537.4X104m3/a,占應(yīng)急水源地總開采量23%.

同時,在地下水排泄量中,泉水流量和邊界流出量均有所減小,泉流量從應(yīng)急開采前的2333.0X104m3/a減小為1953.7X104m3/a,減小約20%;邊界流出量從2020.0X 104m3/a減小為1391.0X104m3/a,減小了30%.

應(yīng)急水源地開采后,地下水系統(tǒng)仍為正均衡狀態(tài),但地下水系統(tǒng)與地表河流的補(bǔ)排關(guān)系總體發(fā)生了轉(zhuǎn)變,開采進(jìn)一步激發(fā)了河流滲漏量,并且襲奪了部分泉流量和側(cè)向排出量.

5 結(jié)論

本研究利用地下水流數(shù)值模擬法建立了北京房山巖溶水應(yīng)急水源地地下水流數(shù)值模型,確定了水文地質(zhì)參數(shù),并進(jìn)行了地下水系統(tǒng)均衡分析,模擬了地下水系統(tǒng)變化過程.模擬和預(yù)測結(jié)果表明,地下水系統(tǒng)與地表河流和泉水關(guān)系密切.在枯水年份,應(yīng)急水源地持續(xù)開采使得地下水系統(tǒng)與地表河流總體補(bǔ)排關(guān)系發(fā)生改變,由地下水向河流排泄轉(zhuǎn)變?yōu)楹恿鳚B漏補(bǔ)給地下水.其研究結(jié)果能夠為應(yīng)急水源地建設(shè)和運(yùn)行提供科學(xué)依據(jù).

[1]謝振華,張兆吉,邢國章等.華北平原典型城市地下水供水安全保障分析[J].資源科學(xué),2009,31(3):400～405.

[2]郭高軒,劉文臣,辛寶東等.北京巖溶水勘查開發(fā)的現(xiàn)狀與思考[J].南水北調(diào)與水利科技,2011(2).

[3]郭高軒,辛寶東,劉文臣.電測深在北京-張坊應(yīng)急巖溶水源地勘查中的應(yīng)用[J].物探與化探,2010,34(2):225～228.

[4]北京市水文地質(zhì)工程地質(zhì)大隊,北京京燕水利管理有限公司.北京應(yīng)急水源-房山巖溶地下水供水水文地質(zhì)勘察與評價[R].2007.

[5]辛寶東.北京市房山區(qū)巖溶地下水水文地球化學(xué)特征[J].水文地質(zhì)工程地質(zhì),2005,(3):74～75.(XIN Bao-dong.The Geochemistry Characteristics of Karst Groundwater in Fangshan Region, SW Beijing[J].Hydrogeology and Engineering Geology, 2005,(3):74～75.(in Chinese))

[6]劉記來.北京應(yīng)急地下水源地開采極限研究[D].中國科學(xué)院研究生院博士學(xué)位論文,2010.

[7]錢家忠,吳劍鋒,董洪信等.徐州市張集水源地裂隙巖溶水三維等參有限元數(shù)值模擬[J].水利學(xué)報,2003,(3):37～41.

[8] Dufresne D P, C W Drake.Regional groundwaterflow model construction and wellfield site selection in a karst area, Lake City, Florida [J].Engineering Geology, 1999,52(1):129～139.

[9] 韓 巍,李國敏,黎 明等.大武水源地巖溶地下水開采動態(tài)數(shù)值模擬分析[J].中國巖溶,2008,27(2):182～188.

[10] 周 訓(xùn),陳明佑,方 斌等.埋藏型巖溶地下水源地的三維數(shù)值模擬-以天津市寧河北巖溶地下水源地為例[J].中國巖溶,2006,25(1):6～11.

[11] 李俊亭.地下水流數(shù)值模擬[M].北京:地質(zhì)出版社,1989.

[12]McDonaldM G, Harbaugh A W.A modular threedimensional finite-difference ground-water flow model: U.S.Geological Survey Techniques of Water-Resources Investigations[R].Book6, Chapter A1, 1988:586.

[13]盧文喜.地下水系統(tǒng)的模擬預(yù)測和優(yōu)化管理[M].北京:科學(xué)出版社,1999.

[14]陳 喜,劉傳杰,胡忠明等.泉域地下水?dāng)?shù)值模擬及泉流量動態(tài)變化預(yù)測[J].水文地質(zhì)工程地質(zhì),2006,(2): 36～40.