呂情緒 ，狄軍貞 ，李 果 ，林 鑫

（1.神華神東煤炭集團(tuán)有限責(zé)任公司, 陜西 神木 719315；2.遼寧工程技術(shù)大學(xué) 土木工程學(xué)院, 遼寧 阜新 123000）

0 引 言

我國正處于工業(yè)化、城鎮(zhèn)化加快發(fā)展的歷史階段，能源需求總量仍有增長空間[1]。在相當(dāng)長時期內(nèi)，煤炭仍是保障我國能源安全穩(wěn)定的基礎(chǔ)能源，將持續(xù)為我國宏觀經(jīng)濟(jì)的高增長作出貢獻(xiàn)[2]。煤炭資源與地下水資源聯(lián)系緊密，在開采過程中會引發(fā)含水層疏干、河溝斷流等一系列環(huán)境問題，還可能發(fā)生礦井水害事故，嚴(yán)重威脅生產(chǎn)[3-6]。我國內(nèi)蒙古地區(qū)煤炭資源儲量豐富、賦存穩(wěn)定，易于進(jìn)行大規(guī)模、高效率、高強(qiáng)度的開采，但內(nèi)蒙古地區(qū)植被稀少、荒漠化嚴(yán)重，本就干旱的氣候遇上高強(qiáng)度的開采，水資源短缺形勢日益嚴(yán)峻，因此更應(yīng)重視該地區(qū)煤礦的地下水資源保護(hù)。采礦活動對地下水?dāng)_動研究歸根結(jié)底是對地下水流運(yùn)動的探究。目前，針對煤礦地下水流場變化的數(shù)值模擬研究已有大量報道[7-10]，ZHANG 等[11]認(rèn)為平頂山十礦-270 m 水平與-300 m水平同屬一個水文地質(zhì)單元，遂采用比擬法，以掘進(jìn)面積作為變量，來預(yù)測采掘-300 m 時礦井涌水量，預(yù)測結(jié)果符合實(shí)際。張鳳娥等[12]采用數(shù)值模擬法建立了神府礦區(qū)大柳塔礦區(qū)地下水平面滲流模型，模擬預(yù)測了采礦活動下地下水滲流場的動態(tài)變化過程，分析了地下水水位隨采掘范圍變化而波動的動態(tài)規(guī)律。白曉等[13]運(yùn)用地下水流數(shù)值模型和降水量預(yù)測模型，預(yù)測分析了峰峰礦區(qū)巖溶地下水資源量和地下水流場的變化。由此可見，數(shù)值模擬法應(yīng)用于礦區(qū)地下水演變具有重要意義。GMS(Groundwater Modeling System)是美國楊百翰大學(xué)環(huán)境模型實(shí)驗(yàn)室聯(lián)合美國軍隊(duì)排水工作站研發(fā)的三維可視化地下水模擬軟件。該軟件集成了多種所需的地下水模擬系統(tǒng)，將它們以模塊化的形式分類，主要包含了解決地下水三維流動問題的MODFLOW、解決地下水溶質(zhì)傳輸問題的MODPATH、MT3D、解決飽和非飽和地下水流動及溶質(zhì)傳輸問題的FEMWATER 模塊，還有MAP、TINs、GIS 等多種輔助模塊，對解決地下水問題的功能十分齊全。模擬過程中能根據(jù)研究區(qū)實(shí)際水文地質(zhì)條件插入并調(diào)節(jié)點(diǎn)、線、面等對象，解決研究區(qū)不規(guī)則邊界概化的難題；不僅如此，GMS 軟件還可以利用鉆孔數(shù)據(jù)生成Solid 模型，該建模方法能刻畫逼真的地層空間結(jié)構(gòu)，含、隔水層分布狀態(tài)一目了然，可視化效果明顯；數(shù)值軟件不僅自身功能強(qiáng)大，還可與其他系統(tǒng)提供數(shù)據(jù)兼容，達(dá)到聯(lián)合使用的效果，充分發(fā)揮其信息處理和分析功能[14-18]。

現(xiàn)以內(nèi)蒙古神東礦區(qū)某井田為例，利用數(shù)值模擬法，預(yù)測礦區(qū)高強(qiáng)度開采下地下水流場的動態(tài)變化，探究高強(qiáng)度采動條件下地下水均衡狀況和多含水層系統(tǒng)水流場的演變規(guī)律。

1 研究區(qū)概況

1.1 自然地理概況

研究區(qū)如圖1 所示，位于內(nèi)蒙古與陜西交界處，屬溫帶大陸性半干旱氣候，年降雨量小、蒸發(fā)量大，研究區(qū)位于烏蘭木倫河一級階地的西緣，總體地勢呈西北高、東南低的楔形特征。井田內(nèi)徑流主要有烏蘭木倫河、補(bǔ)連溝和呼和烏素溝。

圖1 研究區(qū)位置與等高線地形Fig.1 Location and contour map of study area

1.2 采掘概況

目前該煤礦開采設(shè)計能力為28 Mt/a，工作面采高達(dá)5.5 m，深度達(dá)500 m，工作面長度300 m，推進(jìn)長度4 500 m。根據(jù)范立民[19]對高強(qiáng)度開采的定義（表1），該煤礦屬于空間高強(qiáng)度開采。采用多盤區(qū)開采方式和全部垮落、一次采全高的長壁綜合機(jī)械化采煤方法。

表1 煤炭資源開采強(qiáng)度劃分指標(biāo)Table 1 Coal resources exploitation intensity classification index

2 地下水?dāng)?shù)值模型的建立

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

2.1.1 含、隔水層概化

通過對研究區(qū)地層結(jié)構(gòu)和巖性進(jìn)行剖析，將含水層自上而下劃分為第四系孔隙潛水含水層、第三系紅土隔水層、白堊系志丹群含水層、侏羅系安定組隔水層、侏羅系直羅組含水層和侏羅系延安組隔水層，共3 個含水層和3 個隔水層。研究選取礦井的主要充水來源志丹群和直羅組含水層為模擬的目標(biāo)含水層，地層模型與剖面情況如圖2 所示。

圖2 地層結(jié)構(gòu)與剖面Fig.2 Stratigraphic structure and section

2.1.2 邊界條件概化

水平方向上：根據(jù)水文地質(zhì)資料分析，研究區(qū)北面的呼和烏素溝、南面的活雞兔溝、東面的烏蘭木倫河切割出一個較為獨(dú)立、完整的水文地質(zhì)單元，活雞兔溝和呼和烏素溝為季節(jié)性河流，故設(shè)定為定流量邊界；東側(cè)烏蘭木倫河設(shè)置為定水頭邊界；西側(cè)邊界為人為設(shè)定的邊界，采用軟件中的通用水頭模塊進(jìn)行模擬，流量為軟件系統(tǒng)自動運(yùn)算調(diào)整。

垂直方向上：上方松散含水層可接受大氣降水入滲和灌溉回歸入滲補(bǔ)給，主要排泄方式為人工開采和側(cè)向徑流。底部邊界為12 煤至22 煤底板，滲透能力較差，可作為垂向上的隔水邊界。

2.1.3 地下水流動特征概化

區(qū)內(nèi)含水層巖性各不相同，地層介質(zhì)的滲透性能隨空間位置和水流方向的不同而變化，可將含水層概化為非均質(zhì)各向異性。地下水系統(tǒng)各運(yùn)動要素隨時間發(fā)生改變，且含水層在垂向上存在相互的水力聯(lián)系，因此將研究區(qū)含水層概化為三維非穩(wěn)定流。

2.2 水流模型模擬

2.2.1 模型離散

模擬區(qū)面積約276 km2，根據(jù)水文地質(zhì)條件的概化和研究精度要求，在水平方向和垂直方向上離散為有限數(shù)量的長方體，圈定的平面網(wǎng)格尺寸為327.8 m ×286 m，剖分75 行75 列，共5 625 個單元格；垂向上6 層共33 750 個單元格，其中研究區(qū)活動單元格為18 028 個。

2.2.2 模型識別與驗(yàn)證

模型識別驗(yàn)證過程中需要遵循以下原則[20]：①以流場為校正依據(jù)，需保證模擬地下水流場與實(shí)際流場走向與形態(tài)基本接近；②以水位為校正依據(jù)，需保證模擬地下水水位與實(shí)測水位變化趨勢一致，兩者差值滿足精度要求；③區(qū)域地下水均衡量應(yīng)該相差不大，接近于公式計算補(bǔ)給排泄變化量；④校正的水文地質(zhì)參數(shù)要與實(shí)際情況對應(yīng)，不能出入太多。模擬校正過程以2016 年全年地下水位觀測資料為基礎(chǔ)，1 個月為1 個應(yīng)力期，每個應(yīng)力期內(nèi)劃分若干時間步長；再以2019 年全年地下水位觀測資料為依據(jù)，對模型參數(shù)進(jìn)行再校正。采用試算法，不斷地調(diào)整試算水文地質(zhì)參數(shù)，使計算值落在95%置信區(qū)間范圍內(nèi)。由圖3、圖4 可看出模型水位擬合程度較好，校正后的水文地質(zhì)參數(shù)見表2。

表2 校正后水文地質(zhì)參數(shù)Table 2 Corrected hydrogeological parameters

圖3 白堊系含水層觀測孔校核情況Fig.3 Check of observation holes in the Cretaceous aquifer

圖4 直羅組含水層觀測孔校核情況Fig.4 Check of observation holes in Zhiluo Formation

再以2016 年12 月模擬水位為初始水位，在模型中輸入2017-2019 年各項(xiàng)模擬要素，一個月為一個時間步長，共36 個時間步長；根據(jù) 2019 年全年觀測水位，對模型進(jìn)行驗(yàn)證，觀測孔水位擬合狀態(tài)如圖5 所示。模擬結(jié)果誤差均在允許區(qū)間內(nèi)；誤差值緊密地散落在y=x曲線兩側(cè)，說明擬合精度較高，模擬結(jié)果可信度高，水流模型校核情況較為理想。

圖5 2019 年1 月志丹群和直羅組含水層校核結(jié)果Fig.5 Check results of Zhidan Group and Zhiluo Formation aquifer in January 2019

3 煤礦開采對含水層水流場的影響

根據(jù)礦井開采計劃，通過軟件對研究區(qū)2020-2029 年這10 年的地下水非穩(wěn)定流模型進(jìn)行計算，開采盤區(qū)疏干排水量以2019 年各盤區(qū)涌水平均值為基礎(chǔ)，每2 年增加3%。預(yù)測得到直羅組含水層水位流場變化如圖6 所示?？梢钥闯觯孩倜旱V開采后，在盤區(qū)內(nèi)形成了多個沿抽排水孔為中心的水位降落漏斗，第5 年末地下水徑流條件發(fā)生改變，由原先的西北向東南改變?yōu)樗闹芟蚵┒分行膮R聚，在三、四盤區(qū)內(nèi)表現(xiàn)得最為顯著。越靠近漏斗中心水力梯度越大、地下水徑流速度越快。②隨著采空區(qū)面積的逐年增大，與初始流場相比，直羅組降落漏斗面積逐漸增加，影響半徑不斷擴(kuò)大，漏斗中心處水位持續(xù)下降，水位下降幅度從第1 年末的105 m 到第5 年末的337 m，導(dǎo)致周邊地下水向漏斗中心補(bǔ)給加劇，該層地下水水位整體下降，至10 年末，該層地下水水位才趨于穩(wěn)定，水位下降幅度在351 m 左右。因此，開采活動對直羅組含水層影響程度較大。

圖6 直羅組含水層水流場Fig.6 Flow field of Zhiluo Formation aquifer

根據(jù)數(shù)值軟件對研究區(qū)未來10 年地下水流變化的動態(tài)模擬，將志丹群含水層第1 年末、3 年末、5年末、10 年末的預(yù)測水位流場與2016 年初始水位流場相比較（圖7）。①志丹群含水層流場并未發(fā)生明顯變化，水位下降情況較為緩和，下降最大處為三盤區(qū)，水位下降了20 m。②隨著開采年份、開采強(qiáng)度的增加，志丹群含水層水流場發(fā)生改變，第3 年末，水位降落漏斗面積持續(xù)增大，漏斗中心水位下降了33 m。③到第5 年末，降落漏斗面積進(jìn)一步擴(kuò)大，四周地下水逐漸向漏斗中心處補(bǔ)給，由于上部第四系潛水含水層滲透性較好，對下層的補(bǔ)給能力較強(qiáng)，因此，遠(yuǎn)離采掘區(qū)上方的地下水受擾動程度較弱，靠近采掘區(qū)上方的地下水流場及水位受影響程度相對較大。④第10 年末，志丹群降落漏斗下降高度達(dá)116 m，面積擴(kuò)大速率逐漸下降，地下水流場變化程度趨于平緩。

圖7 志丹群含水層水流場Fig.7 Flow field diagram of Zhidan Group aquifer

利用GMS 軟件中的Flow Budget 模塊得到地下水區(qū)域均衡量，地下水總補(bǔ)給量為357 589.74 m3/d，總排泄量為357 563.62 m3/d，誤差為0.007 3%，說明研究區(qū)處于均衡狀態(tài)，模擬計算補(bǔ)給排泄量差距小，說明所建立的水流模型較為合理。

由表3 水均衡統(tǒng)計可知，降雨入滲是研究區(qū)的主要補(bǔ)給來源，占比51.90%，其次為河流補(bǔ)給，占比45.21%；河流排泄為主要排泄方式，占比達(dá)78.24%，其次為潛水蒸發(fā)，占比9.82%，礦井排水占比4.52%。

表3 地下水均衡情況Table 3 Groundwater balance

補(bǔ)給總量與排泄總量均為42 903 m3/d，占區(qū)域總補(bǔ)給項(xiàng)和總排泄項(xiàng)的12.0%，補(bǔ)給排泄處于均衡狀態(tài)。補(bǔ)給總量與排泄總量均為20 297 m3/d，占區(qū)域總補(bǔ)給項(xiàng)和總排泄項(xiàng)的5.7%，補(bǔ)給排泄處于均衡狀態(tài)（圖8、圖9）。

圖8 志丹群含水層的區(qū)域水均衡計算Fig.8 Regional water balance calculation chart of Zhidan Group aquifer

圖9 直羅組含水層的區(qū)域水均衡計算Fig.9 Regional water balance calculation chart of Zhiluo Formation aquifer

4 結(jié) 論

1）在采礦活動中，第1 年末，漏斗中心水位下降了105 m；隨著采掘強(qiáng)度的加大，水位降落漏斗的面積逐漸擴(kuò)大，第3 年末下降了176 m；第5 年末地下水徑流條件發(fā)生改變，降落漏斗中心水位下降幅度增加到337 m；至第10 年末，降落漏斗面積增大不明顯、水位下降不顯著，穩(wěn)定在351 m 左右。

2）志丹群含水層地下水流場并未發(fā)生較大改變，僅在礦井上方的局部區(qū)域水流場發(fā)生改變。第1 年末水位降落漏斗逐漸出現(xiàn)，漏斗中心水位下降20 m，至第10 年末，采掘中心處水位下降了116 m。

3）根據(jù)區(qū)域均衡與層間均衡分析，區(qū)域地下水總補(bǔ)給量為357 589.74 m3/d，總排泄量為357 563.62 m3/d，誤差為0.007 3%，說明研究區(qū)處于均衡狀態(tài)，模擬計算補(bǔ)給排泄量差距小，說明所建立的水流模型較為合理。降雨入滲是研究區(qū)的主要補(bǔ)給來源，占比51.90%，其次為河流補(bǔ)給，占比45.21%；河流排泄為主要排泄方式，占比達(dá)78.24%，其次為潛水蒸發(fā)，占比9.82%，礦井排水占比4.52%。志丹群和直羅組含水層補(bǔ)給排泄分別占區(qū)域總均衡量的12.0%和5.7%，其中礦井排水量為直羅組含水層的主要排泄方式，占比71.3%。