朱 君,安瑞瑞
(1.中國輻射防護(hù)研究院,山西 太原 030006; 2.山西晉環(huán)科源環(huán)境資源科技有限公司,山西 太原 030024)
隨著地下水?dāng)?shù)值模擬技術(shù)的發(fā)展,應(yīng)用該技術(shù)反映礦區(qū)地下水動態(tài)變化已是目前最常用的方法,國際上通用的商業(yè)軟件有Visual Modflow、GMS、FEMWATER、FEFLOW等。如何提高數(shù)值模擬的精準(zhǔn)度使之更加符合實(shí)際情況是亟待解決的問題。礦區(qū)斷層構(gòu)造復(fù)雜,有的區(qū)域甚至縱橫交錯,斷層的概化和處理成為影響模擬結(jié)果的一個重要因素。目前,國內(nèi)在地下水?dāng)?shù)值模擬中能概化三維斷層并應(yīng)用到模型中的實(shí)例很少,普遍的做法是將斷層兩側(cè)地層的滲透系數(shù)分區(qū)取值,或者與相鄰含水層視為一個系統(tǒng)。而考慮斷層因素的情況也是應(yīng)用Visual Modflow、GMS針對單斷層的簡單構(gòu)造區(qū),Visual Modflow和GMS采用有限差分法,缺點(diǎn)是網(wǎng)格剖分方向必須與斷層走向一致,見圖1中的Fa、Fb;否則概化斷層的網(wǎng)格呈“鋸齒”狀,見圖1中的Fc,影響計算結(jié)果的收斂性和精度,顯然在復(fù)雜斷層區(qū)域采用有限差分法已無法滿足精度要求[1-5]。FEMWATER和FEFLOW采用有限節(jié)點(diǎn)法,能夠消除“鋸齒”狀網(wǎng)格的不利影響,更好地描述復(fù)雜的斷層構(gòu)造,模擬計算結(jié)果更加真實(shí)[6-8]。國外應(yīng)用FEFLOW在地下水模型中對復(fù)雜斷層的概化已經(jīng)能夠做得非常精準(zhǔn),因此國內(nèi)亟須彌補(bǔ)差距[9-13]。本文以內(nèi)蒙古某煤礦內(nèi)的10條斷層為例,建立復(fù)雜斷層區(qū)的地下水三維數(shù)值模型。
圖1 有限差分法斷層概化示意圖
研究礦區(qū)位于內(nèi)蒙古自治區(qū)與寧夏回族自治區(qū)交界處,屬于華北地臺、鄂爾多斯盆地西緣褶皺沖斷帶的北段,即賀蘭山—橫山段。褶皺、斷裂較發(fā)育,褶皺軸向及主要斷裂多呈近南北向展布,其次為北東—北東東向斷裂及北西向斷裂。總體構(gòu)造形態(tài)為受斷層切割的單斜,地層總體東傾,走向近南北,在單斜基礎(chǔ)上發(fā)育次級褶曲和少量斷層。
礦區(qū)內(nèi)共有10條斷層,其中逆斷層6條(F1、F2、F3、F5、F6、F8),正斷層4條(F4、F7、F9、F10),均屬于傾斜斷層,錯斷了新近系以下地層,斷層特征見表1。
表1 斷層特征
礦區(qū)開采煤層位于二疊系山西組,上覆含水層從上至下為第四系松散孔隙含水層,厚2.05~59.50 m,主要以大氣降水補(bǔ)給,受下伏新近系黏土隔水層控制,厚度82.60~397.32 m,由東北向西南徑流。新近系下部砂礫石層含水層,厚度5.45~60.48 m,平均厚19.78 m,孔隙發(fā)育,透水性好,為該地區(qū)富水性最佳的含水層,主要接受區(qū)域側(cè)向補(bǔ)給,由西北向東南徑流。二疊系石盒子組砂巖裂隙含水層,厚128.06~638.76 m,主要接受上覆新近系砂礫石層含水層補(bǔ)給,由西北向東南徑流。二疊系山西組砂巖裂隙含水層,厚24.42~121.69 m,主要接受上覆新近系砂礫石層含水層越流補(bǔ)給和二疊系石盒子組砂巖裂隙含水層直接補(bǔ)給,由西北向東南徑流。因此二疊系砂巖裂隙含水層與新近系砂礫石層含水層地下水有著相同的變化規(guī)律。
模型的東邊界(AB段),以蔥溝斷層為界,概化為流量邊界。模型的北邊界(AD段)、南邊界(BC段),大致以地下水流線為界,概化為流量邊界。模型的西邊界(CD段),以已知水頭為界,概化為定水頭邊界。模擬區(qū)面積約130 km2,基本概況見圖2。
圖2 模擬區(qū)基本概況示意圖
開采煤層上覆4含水層,分別為第四系松散孔隙含水層、新近系砂礫石含水層、二疊系石盒子組砂巖裂隙含水層和二疊系山西組砂巖裂隙含水層。另外第四系松散孔隙含水層與新近系砂礫石含水層之間有82.60~397.32 m的黏土層,可視為弱透水層,因此整個模擬區(qū)自上而下劃分為5層(表2)。
表2 目標(biāo)含水層概化情況
礦區(qū)內(nèi)的10條斷層均屬于傾斜斷層,傾角在70°~75°之間,縱橫交錯,較為復(fù)雜,故將傾斜斷層處理為垂直斷層??臻g離散時,垂向上分為5層,斷層處網(wǎng)格加密,共生成三角網(wǎng)格370 060個,節(jié)點(diǎn)217 537個。網(wǎng)絡(luò)剖分示意圖見圖3,斷層概化三維示意圖見圖4。
圖5 各含水層抽水歷時曲線
圖3 網(wǎng)格剖分示意圖
圖4 斷層概化三維示意圖
水文地質(zhì)參數(shù)的取值,包括含水層、弱透水層、斷層x、y、z方向的主滲透系數(shù)Kxx、Kyy、Kzz、彈性釋水系數(shù)Sy、單位貯存量Ss。滲透系數(shù)的取值主要是參考礦區(qū)及相鄰礦區(qū)的抽水試驗(yàn)資料,抽水歷時曲線見圖5,圖中S為水位降深,Q為流量,得到第四系松散孔隙含水層、新近系砂礫石含水層、二疊系石盒子組砂巖裂隙含水層、二疊系山西組砂巖裂隙含水層水平方向的滲透系數(shù)Kxx、Kyy,垂直方向的滲透系數(shù)Kzz一般為水平方向滲透系數(shù)的1/10[9-12]。另外,根據(jù)礦區(qū)《水文地質(zhì)勘查成果報告》,斷層透水,滲透系數(shù)約為相鄰含水層的10倍,經(jīng)模型校核后確定該參數(shù)取值。彈性釋水系數(shù)Sy、單位貯存量Ss主要取經(jīng)驗(yàn)值。另外,第四系松散孔隙含水層根據(jù)水文地質(zhì)參數(shù)不同分為A、B 2個區(qū)。具體取值情況見圖6和表3。
圖6 第四系松散孔隙含水層水文地質(zhì)參數(shù)分區(qū)
地層概化層分區(qū)Kxx/(m·d-1)Kyy/(m·d-1)Kzz/(m·d-1)SySs/m-1Q松散孔隙含水層A6.56.50.650.180.00013B10.8910.891.0890.250.00049N2弱透水層8.64×10-68.64×10-68.64×10-70.00082.6×10-5砂礫石含水層2.32722.32720.23720.350.001Psh砂巖裂隙含水層0.00340.00340.000340.026.9×10-5斷層0.0340.0340.0340.026.9×10-5P1s砂巖裂隙含水層0.003730.003730.0003730.036.9×10-5斷層0.03730.03730.03730.036.9×10-5
模擬區(qū)內(nèi)共有12個第四系松散孔隙含水層地下水水位監(jiān)測點(diǎn),其位置分布見圖2。在2014年8月—2015年3月的一個連續(xù)水文年對該含水層的水位進(jìn)行了監(jiān)測,具體監(jiān)測值見表4,從監(jiān)測結(jié)果分析,在一個連續(xù)水文年中,所有點(diǎn)位水位變化幅度在0.56~0.61 m,可視為穩(wěn)定流,將3期水位的平均值作為模型校核資料。
表4 第四系松散孔隙含水層水位監(jiān)測值 m
圖7 第四系松散孔隙含水層水位擬合線
圖8 二疊系石盒子組砂巖裂隙含水層水位擬合線
模型建立后,對比計算水位與監(jiān)測水位的吻合程度(圖7~8),驗(yàn)證模型的合理性和準(zhǔn)確性。在加入疏干排水條件前,以穩(wěn)定流模式運(yùn)行模型,計算得到第四系松散孔隙含水層的初始水位分布值。該地區(qū)年平均降雨量為270.4 mm,降雨入滲系數(shù)經(jīng)過模型驗(yàn)證和校核后取0.12,輸入模型的降雨入滲量為32.85 mm/d。通過穩(wěn)定流模型計算得到第四系松散孔隙含水層初始水位與實(shí)際監(jiān)測水位的偏差范圍在0.39~1.83 m,其中偏差較大的有3處,分別是SMJ-11為1.69 m,SMJ-35為1.83 m,SMJ-44為1.74 m。其余監(jiān)測點(diǎn)的偏差均小于0.72 m(表5)。
表5 擬合結(jié)果 m
表6 二疊系石盒子組砂巖裂隙含水層監(jiān)測水位與計算水位 m
表7 二疊系山西組砂巖裂隙含水層和新近系砂礫石含水層之間的水位關(guān)系
另外根據(jù)《水文地質(zhì)勘查成果報告》,礦區(qū)有9個水文地質(zhì)鉆孔監(jiān)測過二疊系石盒子組砂巖裂隙含水層水位(表6);2個水文地質(zhì)鉆孔監(jiān)測過二疊系山西組砂巖裂隙含水層和新近系砂礫石含水層之間的水位關(guān)系(表7)。由鉆孔ZK2803、ZK2805可知,新近系砂礫石含水層(E)、二疊系山西組砂巖裂隙含水層(P1s)水位相差不大,間接說明新近系砂礫石層含水層以下地層沒有穩(wěn)定隔水層,二疊系石盒子組砂巖裂隙含水層和二疊系山西組砂巖裂隙含水層主要靠新近系砂礫石含水層補(bǔ)給,地下水有著相同的變化規(guī)律。
二疊系石盒子組砂巖裂隙含水層初始水位與實(shí)際監(jiān)測水位的偏差范圍在0.33~2.51 m,根據(jù)以上2含水層水位擬合程度可知,建立的模型和水文地質(zhì)參數(shù)的取值基本上能夠反映實(shí)際情況。第四系松散孔隙含水層計算水位等值線圖和二疊系石盒子組砂巖裂隙含水層計算水位等值線見圖9~10。
圖9 第四系松散孔隙含水層計算水位等值線
圖10 二疊系石盒子組砂巖裂隙含水層計算水位等值線
根據(jù)冒落裂隙帶和導(dǎo)水裂隙帶計算,在煤層隱伏露頭處附近最大高度能導(dǎo)通新近系底部砂礫石含水層,礦區(qū)開采后正常涌水量為13 800 m3/d,垂向上地下水流速度加劇,將涌水量分采區(qū)、分時段輸入模型,并以非穩(wěn)定流模型計算礦區(qū)開采30年后對地下水資源的影響,得到礦區(qū)開采30年水位等值線見圖11。預(yù)測結(jié)果顯示,開采后水位變化部分集中在新近系砂礫石含水層及以下區(qū)域,第四系松散孔隙含水層的水位等值線基本上沒有變化,分析原因是第四系松散孔隙含水層底部有厚度82.60~397.32 m的完整且連續(xù)的棕紅色黏土層,有效阻隔了地下水的滲漏和越流,對該含水層影響較小。
圖11 礦區(qū)開采30年后水位等值線
礦區(qū)開采30年后,二疊系山西組砂巖裂隙含水層、二疊系石盒子組砂巖裂隙含水層、新近系砂礫石含水層形成以開采區(qū)為中心的降落漏斗,最大降深及影響面積詳見表8。分析原因是二疊系山西組煤層開采后,“兩帶”直接破壞二疊系山西組砂巖裂隙含水層、二疊系石盒子組砂巖裂隙含水層下部以及煤層露頭處的新近系底部砂礫石含水層,加大了地下水的滲漏和越流。另外礦區(qū)內(nèi)的10條透水?dāng)鄬?形成滲水通道,加劇了上述含水層地下水的下滲。
表8 含水層水位降深、影響面積
圖12 開采30年后新近系砂礫石含水層水位及降深等值線
從預(yù)測結(jié)果來看(圖12~14),二疊系石盒子組砂巖裂隙含水層與新近系砂礫石含水層的水位及水位降深等值線基本相似。原因是該地區(qū)二疊系石盒子組砂巖裂隙含水層主要接受上覆新近系砂礫石層含水層補(bǔ)給,且兩含水層之間沒有隔水層,因此地下水有著相同的變化規(guī)律。
圖13 開采30年后二疊系石盒子組砂巖裂隙含水層水位及降深等值線
圖14 開采30年后二疊系山西組砂巖裂隙含水層水位及降深等值線
復(fù)雜斷層的概化和處理方式是影響數(shù)值模擬結(jié)果的重要因素,本文以內(nèi)蒙古某礦區(qū)中的10條錯綜復(fù)雜的斷層為研究對象,應(yīng)用FEFLOW軟件將傾斜斷層簡化為垂直透水?dāng)鄬?建立三維地下水?dāng)?shù)值模型,克服了以往Visual Modflow、GMS建模過程中斷層處“鋸齒”狀網(wǎng)格對計算結(jié)果收斂性和精度的影響,并預(yù)測了礦區(qū)開采30年后,對上覆各含水層地下水資源的破壞程度,二疊系山西組砂巖裂隙含水層、二疊系石盒子組砂巖裂隙含水層、新近系砂礫石含水層受到了不同程度的影響,而第四系松散孔隙含水層由于底部較厚且完整、連續(xù)的黏土層阻隔幾乎沒有受到影響。
[1] 王博,劉耀煒,孫小龍,等.斷層對地下水滲流場特征影響的數(shù)值模擬[J].地震,2008,28(3):115-123.(WANG Bo,LIU Yaowei,SUN Xiaolong,et al.Numerical simulation of the influence of fault on the groundwater seepage field[J].Journal of Earthquake,2008,28(3):115-123.(in Chinese))
[2] 梁世川,徐明,王磊,等.GMS在地下水?dāng)?shù)值模擬及斷層處理中的應(yīng)用:以蓋孜河水源地為例[J].地下水,2013,35(6):53-54,125.(LIANG Shichuan,XU Ming,WANG Lei,et al.Application of GMS in numerical simulation of groundwater and fault generalizability: as an example for the water source location of Gez river[J].Ground Water,2013,35(6):53-54,125.(in Chinese))
[3] 武強(qiáng),朱斌,徐華,等.MODFLOW在淮北地下水?dāng)?shù)值模擬中的應(yīng)用[J].遼寧工程技術(shù)大學(xué)學(xué)報(自然科學(xué)版),2005,24(4):500-503.(WU Qiang,ZHU Bin,XU Hua,et al.Application of MODFLOW in groundwater numerical simulation in Huaibei City[J].Journal of Liaoning Technical University (Natural Science),2005,24(4):500-503.(in Chinese))
[4] 牛建立.地下水?dāng)?shù)值計算中斷層處理的“切割—導(dǎo)通法”[J].煤田地質(zhì)與勘探,2002,30(5):39-40.(NIU Jianli.Cutting-conduction method for fault generalizability in groundwater numerical calculation[J].Coal Geology & Exploration,2002,30(5):39-40.(in Chinese))
[5] 朱斌,武強(qiáng).斷層影響下的地下水流數(shù)值模擬[J].桂林理工大學(xué)學(xué)報,2005,25(1):31-35.(ZHU Bin,WU Qiang.Numerical simulation of groundwater flow under the influence of faults[J].Journal of Guilin University of Technology,2005,25(1):31-35.(in Chinese))
[6] 殷曉曦,陳陸望,林曼利,等.采動影響下任樓煤礦地下水流三維數(shù)值模擬[J].合肥工業(yè)大學(xué)學(xué)報(自然科學(xué)版),2013,36(1):93-98.(YIN Xiaoxi,CHEN Luwang,LIN Manli,et al.Three-dimensional numerical simulation of groundwater flow in Renlou Mine under mining disturbance[J].Journal of Hefei University of Technology(Natural Science),2013,36(1):93-98.(in Chinese))
[7] 韓巍,李國敏,黎明,等.大武水源地巖溶地下水開采動態(tài)數(shù)值模擬分析[J].中國巖溶,2008,27(2):182-188.(HAN Wei,LI Guomin,LI Ming,et al.Dynamic numerical simulation of groundwater in Dawu karst water source location[J].Carsologica Sinica,2008,27(2):182-188.(in Chinese))
[8] 劉志峰,柳碩林,王琳琳,等.數(shù)值模擬法在西龍河嶧山斷層帶巖溶水允許開采量評價中的應(yīng)用[J].地質(zhì)學(xué)刊,2010,34(3):271-274.(LIU Zhifeng,LIU Shuolin,WANG Linlin,et al.The application of numerical simulation method to evaluation the allowable exploitation amount of cavern water in Xilong River Yishan fault zone[J].Journal of Geology,2010,34(3):271-274.(in Chinese))
[11]LUO J.Hydrogeological model Jura Os[R].Bern:Arbeitsbericht NAB 13-26,2014.