初道忠,趙忠琦,2
(1.山東理工大學 資源與環境工程學院,山東 淄博 255049;2.山東黃金集團玲瓏金礦,山東 煙臺 265400)
隨著信息技術的飛速發展,利用數值法計算礦井涌水量的方法得到了廣泛的應用。這類方法能在考慮含水層的各向異性的同時,處理水文地質邊界條件較復雜的礦區,模擬礦山開采過程中礦井疏排水情況及涌水量變化特征。
地下水的數值模擬計算方法主要有有限單元法、有限差分法、邊界元法以及有限分析法等。本文利用有限差分法進行數值模擬計算,其基本思想為把連續定解區域用有限個離散點構成網格代替,把連續定解區域上的連續變量函數用離散變量函數近似表示。本文利用GMS軟件中的BOREHOLE模塊及MODFLOW模塊對大莊子礦區的礦井涌水量進行預測,并利用模擬結果對礦山開采進行指導。
礦區分布于華北板塊(Ⅰ級)、膠遼隆起區(Ⅱ級)、膠北隆起(Ⅲ級)、膠北斷隆(Ⅳ級)、明村-但山凸起(Ⅴ級)處。該區域屬河流沖洪積堆積地貌類型,大多為平原,地勢起伏情況較緩,總體呈現東南高西北低的地勢特點,最低侵蝕基準面標高位于礦體以上,為+10 m左右。
礦區屬于北溫帶半濕潤季風氣候,根據該區域氣象局1992~2014年氣象資料,該區域多年平均氣溫12~12.2 ℃,1998年8月13日達到高溫最值,為38.6 ℃,1992年1月19日達到低溫最值,為-18.3 ℃。區域降水情況如圖1所示,年平均降水量680 mm,年平均蒸發量1 711.1 mm,蒸發量大于降水量,年平均相對濕度69%。礦區地表排水設施完善,可輕松排走由于降水所存積的積水。該區域距離大規模水體較遠,區域西部有一條流入萊州灣的膠萊河,區域南部有一條匯入膠萊河的淄陽河,區域北部有一條流入萊州灣的沙河,三者都為間歇性河流,旱季大都干涸,僅雨季具有較大流量,不會影響到目標區域礦床的開采。

圖1 1981—2010年區域降水量柱形圖Fig.1 Column chart of regional precipitation from 1981 to 2010
依據地下水的賦存條件、水理性質、水力特征等,將礦區地下水大致劃分為三大類型:松散巖類孔隙水、碎屑巖類孔隙裂隙水和基巖裂隙水。各層的補給條件主要為大氣降水、基巖風化裂隙水以及相鄰的裂隙含水層的地下水靜儲量等。通過對鉆孔終孔地下水位埋深的研究發現,在枯水期時各鉆孔水位埋深大多低于第四系底板,而豐水期時水位埋深則超過第四系底板。故目標區域地下水不承壓,為潛水。
根據礦區含水層的主要分布特點,對復雜含水層分布情況進行概化[1-2],在不失真的前提下,將礦區含水層概化為四層,由上到下分別為第四系松散孔隙含水層、林家莊組砂巖含水層及荊山群基巖裂隙含水層聯合含水層、云山單元二長花崗巖隔水層以及荊山群基巖裂隙弱含水層。具體概化分層方法如圖2所示。

圖2 含水層概念模型Fig 2 Aquifer conceptual model
在四層概化地層之下為另一隔水層,由于礦體賦存于荊山群基巖裂隙弱含水層內,故分析四層概化地層之下的隔水層無意義。其中第四系松散孔隙含水層、林家莊組砂巖含水層受季節影響較大,降水補給量占比較多,地下水為潛水;由于有云山單元二長花崗巖隔水層遮擋,荊山群基巖裂隙含水層地下水具有微承壓性,但由于其基巖裂隙不發育,含水量較少,故承壓性質不明顯[3-4]。
礦區屬魯東低山丘陵水文地質區膠萊拗陷水文地質亞區膠萊盆地孔隙裂隙弱—強富水地段,地處濱海平原區,地勢比較平緩,最大比高10 m左右。礦區內最低侵蝕基準面標高為+10 m。南2 km有淄陽河自東向西流入澤河,為間歇性河流,大部分河段常年干涸,故地表水與礦區地下水水力聯系不密切。
根據礦區的水位等值線圖(圖3)及礦區所在區域方位,礦區西部邊緣主要為裸露的云山單元二長花崗巖,因二長花崗巖裂隙不發育,隔水性良好,故以其為礦區西部第二類隔水邊界;礦區南北部邊緣以其等水位線為邊界,由于水頭已知,故為第一類邊界;由于礦區東部為另一礦區的勘探區域,地質資料不詳,故以礦區東部邊界為界,根據原始流場模擬情況確定為第一類邊界[5]。

圖3 礦區水位等值線Fig.3 Water level contour of mining area
根據礦區的含水層情況,結合地下水流補給及排泄情況,評估區視為三維非穩定流系統[6],利用三維非穩定流微分方程建立數學模型如下所示:

(1)
H(x,y,z,t)E,S,N=φ(x,y,z,t)
(2)
(3)
H(x,y,z,t)|t=0=H0(x,y,z)
(4)
其中:K為含水層滲透系數;H為評估區水頭分布情況;Ss為含水層貯水率;H0為t=0時的水頭分布情況;W為單位時間單位體積含水層補給或排泄水量;(x,y,z)為空間坐標,E,W,S,N代表評估區東西南北邊界;t代表時間。
滲流數學模型是地下水模擬的理論指導,是模型的概念基礎。該數學模型中,公式(1)表示評估區域各含水層采用潛水三維均質非穩定流微分方程計算方法進行模擬;公式(2)、(3)表示評估區域各個方位的邊界條件情況;公式(4)表示評估區域各位置的初始水頭情況。
根據已經建立的滲流模型,對各含水層賦予初始滲透系數,按照抽水試驗中抽水井及觀測井的位置,在滲流模型中同樣布置觀測井,不斷修改各含水層的滲透系數,對滲流模型進行反演迭代計算,使模擬結果與抽水試驗實測結果相一致,進而確定各含水層滲透系數。
利用觀測孔JG07、JG11、JG24在2016年2月19日~2016年2月27日時間段內的抽水試驗水位數據,對各個測量水位的時間節點進行編號,結合滲流模型模擬結果,其擬合曲線如圖4所示。

(a)JG07觀測孔水位擬合曲線
求得含水層1、含水層2、含水層3、含水層4的滲透系數分別為:3.54 m·d-1、0.81 m·d-1、0.000 52 m·d-1、0.13 m·d-1。
以地下水模擬軟件GMS為平臺,利用已知的鉆孔數據結合borehole模塊[7]生成鉆孔模型以及概化地層模型如圖5、6所示。
利用GMS的MODFLOW模塊[8-9]對含水層模型賦以邊界條件及降水條件見表1。
根據原始流場情況可知,流場與邊界水位反映情況一致,但由于降水緣故,流場中潛水水位面高于邊界水位面,說明降水對于含水層1、2層影響較大,而含水層3、4受到降水影響較小,主要受邊界徑流以及泄流影響。
模擬非抽水條件下原始流場如圖7所示。
在原始流場的基礎上,設立假想抽水井,不斷調試假想抽水井的日抽水量,觀察流場的變化情況,進而得到最優日抽水量。假想抽水井的最優日抽水量見表2。根據評估區域的氣象資料可知,評估區域降水量最大的時間段為六到八月份,故以該時間段評估礦井的最大涌水量。

表2 日抽水量匯總表Tab.2 Daily pumping summary table
根據最優日抽水量模擬流場情況如圖8、9所示。

圖8 抽水數值模擬結果Fig.8 Pumping numerical simulation results
根據抽水模擬結果,紅色區域代表地下水疏干區域,為四個抽水井形成的水位降落漏斗范圍,周圍彩色單元表示由于邊界水頭條件及降水條件仍然存留的地下水含水單元。由于礦體賦存于荊山群野頭組基巖裂隙弱含水層的圍巖蝕變帶中,故以第四含水層的疏水情況作為排水標準。

(a) 含水層1模擬結果
由以上紅色區域分布情況可知,紅色區域分布于四個含水層中的每一層,且至少占每層面積的80%以上。由此可知地下水日抽水量為1 715~2 498m3/d,疏水時間為六月到九月。
對比抽水數值模擬結果與該礦區地質詳查報告中利用大井法[10]預測抽水結構可知,大井法預測最大抽水量為1 532 m3/d,與數值模擬結果基本吻合。兩種方法存在偏差主要是由于大井法計算過程較為簡單,故精確度較低造成的。
通過對大莊子區域礦井涌水量的數值模擬過程可得出以下結論:(1)該礦區的地下水日疏水量為1 715~2 498 m3/d,礦井最大涌水量為2 498 m3/d;(2)利用GMS對地下水進行數值模擬的方法需要考慮地層結構、地下水的徑流影響、降水情況、邊界條件以及滲流模型的選取等問題,適用于地質資料較為詳細時使用,故其計算結果更為可信;(3)數值模擬過程中,如果地層情況較為復雜,概化地層難以表達,可將復雜地層詳細劃分后分別進行概化。