沈媛媛,辛寶東,郭高軒,紀軼群
(北京市水文地質工程地質大隊,北京 100195)
北京市人均水資源量不足300m3,水資源供需矛盾突出.北京市為應對缺水的緊張局面,2002年以后相繼建設了5處應急水源工程.其中房山張坊水源地,以巖溶水為應急水源.本文在詳細分析研究區地質及水文地質條件的基礎上,建立研究區水文地質概念模型和地下水流數值模型,采用MODFLOW求解,完成模型識別和驗證.并預測給定條件下水源地區地下水系統變化情況,為其合理開發利用提供科學依據.
研究區在構造單元上位于霞云嶺-龍門臺復向斜南翼,其東部為北嶺向斜,北部為大石河背斜(圖1).核部基本位于拒馬河與大石河分水嶺處,霧迷山組地層構成向斜兩翼,鐵嶺組、青白口系、寒武系在向斜核部形成一些列寬緩波狀褶皺.斷裂構造主要有北東向的黃莊-高麗營斷裂,牛口峪-長溝斷裂和霞云嶺-大峪溝斷裂.
研究區基巖含水層以霧迷山組白云巖分布最廣,在山前地帶有第四系砂礫石層,為較好的巖溶水含水層,單井涌水量可達3000m3/d.研究區北部還有鐵嶺組和洪水莊組巖溶含水層,其富水性較差.
分布較廣的的碳酸鹽巖裂隙發育,有利于降水直接入滲.地下水自西北向東南方向徑流,地下水位埋由山區到山前深逐漸減小.山前的地下水主要以泉水、側向徑流和人工開采.
本地區巖溶含水層雖然溶蝕裂隙有統一的地下水面,地下水系統可用多孔介質滲流理論模型來概化.
(1)含水層概化
根據研究區的水文地質條件確定了模型邊界,模擬面積為401km2.模擬對象為薊縣系霧迷山和鐵嶺組之碳酸鹽巖含水層.其巖溶含水層埋深400m左右,可將其概化為非均質各項同性含水層.
(2)含水層概化邊界條件
模擬含水層北部邊界為拒馬河與大石河的地表分水嶺,概化為流量邊界;東部為青白口系下馬嶺組板巖、千枚巖,為相對阻水層,概化為隔水邊界;南部與河北省相鄰,大峪溝斷裂成為地下水側向徑流的通道;另外,地下水可通過深層徑流向下游排泄,概化為流出邊界;西部以煌斑巖脈為界,概化為隔水邊界(圖2).


模型的上邊界為潛水水面,巖溶地下水系統通過該邊界與外界降雨入滲、地表水體滲漏等相關系.以含水層埋深400m處為模型底部概化為隔水邊界.
(3)地下水的水力特征
研究區含水層厚度大、分布廣.地下水以水平流動為主、垂向流動為輔,可將地下水流動概化為空間三維非穩定流.
(1)數學模型
概化后的巖溶地下水系統用如下數學公式表示[11]:

式中:Ω為滲流區域;H為巖溶含水層的水位標高(m);B為模型底界標高(m);Kx、Ky、Kz、Kn分別為x、y、z、邊界法向方向的滲透系數(m/d);μ為潛水含水層的給水度;w為含水層的源匯項(m/d);H0為含水層的初始水位分布(m);Γ0為滲流區域的上邊界,即地下水的自由表面;Γ2為滲流區域的側向邊界;Γ4為滲流區域的下邊界,即含水層底部的隔水邊界;n為邊界面的法線方向; q(x,y,z,t)為二類邊界的流量(m/d),設計流入為正,流出為負,隔水邊界值為0.
采用美國地質調查局發布的三維地下水有限差分計算軟件MODFLOW進行地下水數值模擬[12],將模擬區剖分為250mX250m的矩形網格,共剖分為120行、140列,有效單元格6416個.
模擬區分為6個滲透系數參數分區(圖2),其參數初值根據含水層巖性及以往研究資料設定.給水度根據以往試驗結果取為0.01.根據不同的巖性及以往研究成果,確定各巖組巖層的降雨入滲系數:霧迷山組為0.35;鐵嶺組為0.3,洪水莊組及下馬嶺組巖層為0.1.用達西公式計算其邊界側向徑流量.
在張坊以北拒馬河水接受地下水補給,在張坊以南變為河水帶補給地下水.采用MODFLOW河流計算子程序包,根據其河水位、底板高程、河流附近地下水位以及河道滲透系數,計算河流與含水層的交換水量.
泉水流量大小與附近地下水位呈正比,并與含水層滲透性能有關.采用與河流類似的模擬方式,以泉附近含水層滲透系數、地下水位與地表高程差值確定泉流量.
(2)模型計算結果
選取2004年6月統測地下水位作為模擬的初始地下水位,12月為模型識別期;2005年1月至6月為模型驗證期.采用有限差分法求解,以一個月為一個時間段,每個時間段分為10個時間步長,嚴格控制每次迭代的誤差.
選用"試錯法"調整模型參數,觀測井識別期和驗證期的模擬水位與實測水位殘差較小,給出張坊、長溝、廣潤莊和廣祿莊四處的觀測井地下水位擬合曲線對比圖(圖3).從圖中可以看出,除張坊觀測井的水位在個別應力期外,其他均擬合較好,滿足了精度要求,建立的模型較真實地刻畫了其研究區地下水系統特征,可以用于地下水預測評價.

(1)應急開采前地下水系統均衡情況
模擬期內含水層系統水量均衡結果見表1,表中補給量包括降雨入滲補給量和邊界流入量以及灌溉滲漏量.其中降雨入滲量占總補給量的89.1%,是其最重要的補給源,邊界流入量占總補給量的9.0%,灌溉滲漏補給量為1.5%.區域地下水主要以工農業開采、河流和泉水排泄以及邊界側向流出.其中工農業開采占總排泄量的26.0%,泉排泄38.4%,邊界流出量33.3%.
模擬結果表明,在房山應急水源地開采前,地下水向河流排泄量為139.4X104m3/a.

表1 模擬期內地下水系統均衡結果表(X104m3/a)
(2)應急開采后地下水系統均衡情況
為評價應急水源地開采后的地下水系統與河流之補排關系,選擇2009年8月實測的地下水流場作為初始流場,預測枯水年(年均降雨量513mm)和應急水源地現狀開采條件下、5年后的地下水系統均衡結果見表2.

表2 應急水源地開采條件下地下水系統均衡結果表(X104m3/a)
預測結果,應急水源地開采后,河流與地下水系統總的補排關系由應急水源地開采前地下水補給河流轉變為河流滲漏補給地下水,河流滲漏補給量為398.0X104m3/a,占總補給量的5.1%.河流與地下水系統的交換量的變化量為537.4X104m3/a,占應急水源地總開采量23%.
同時,在地下水排泄量中,泉水流量和邊界流出量均有所減小,泉流量從應急開采前的2333.0X104m3/a減小為1953.7X104m3/a,減小約20%;邊界流出量從2020.0X 104m3/a減小為1391.0X104m3/a,減小了30%.
應急水源地開采后,地下水系統仍為正均衡狀態,但地下水系統與地表河流的補排關系總體發生了轉變,開采進一步激發了河流滲漏量,并且襲奪了部分泉流量和側向排出量.
本研究利用地下水流數值模擬法建立了北京房山巖溶水應急水源地地下水流數值模型,確定了水文地質參數,并進行了地下水系統均衡分析,模擬了地下水系統變化過程.模擬和預測結果表明,地下水系統與地表河流和泉水關系密切.在枯水年份,應急水源地持續開采使得地下水系統與地表河流總體補排關系發生改變,由地下水向河流排泄轉變為河流滲漏補給地下水.其研究結果能夠為應急水源地建設和運行提供科學依據.
[1]謝振華,張兆吉,邢國章等.華北平原典型城市地下水供水安全保障分析[J].資源科學,2009,31(3):400~405.
[2]郭高軒,劉文臣,辛寶東等.北京巖溶水勘查開發的現狀與思考[J].南水北調與水利科技,2011(2).
[3]郭高軒,辛寶東,劉文臣.電測深在北京-張坊應急巖溶水源地勘查中的應用[J].物探與化探,2010,34(2):225~228.
[4]北京市水文地質工程地質大隊,北京京燕水利管理有限公司.北京應急水源-房山巖溶地下水供水水文地質勘察與評價[R].2007.
[5]辛寶東.北京市房山區巖溶地下水水文地球化學特征[J].水文地質工程地質,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]劉記來.北京應急地下水源地開采極限研究[D].中國科學院研究生院博士學位論文,2010.
[7]錢家忠,吳劍鋒,董洪信等.徐州市張集水源地裂隙巖溶水三維等參有限元數值模擬[J].水利學報,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] 韓 巍,李國敏,黎 明等.大武水源地巖溶地下水開采動態數值模擬分析[J].中國巖溶,2008,27(2):182~188.
[10] 周 訓,陳明佑,方 斌等.埋藏型巖溶地下水源地的三維數值模擬-以天津市寧河北巖溶地下水源地為例[J].中國巖溶,2006,25(1):6~11.
[11] 李俊亭.地下水流數值模擬[M].北京:地質出版社,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]盧文喜.地下水系統的模擬預測和優化管理[M].北京:科學出版社,1999.
[14]陳 喜,劉傳杰,胡忠明等.泉域地下水數值模擬及泉流量動態變化預測[J].水文地質工程地質,2006,(2): 36~40.