新疆北部某煤礦持續(xù)采掘條件下礦井涌水量動態(tài)模擬預(yù)測

楊 晟，賈 超，袁 涵，張建配，陳 陽

(1.山東大學(xué)土建與水利學(xué)院，山東 濟(jì)南 250061；2.山東大學(xué)海洋研究院，山東 青島 266237)

礦井涌水量是衡量礦井開采方式合理性和安全性的重要指標(biāo)，直接關(guān)系到礦井的安全生產(chǎn)。因此礦井涌水量的準(zhǔn)確預(yù)測對于礦井開采和生產(chǎn)具有重要意義[1]。礦井開采是一個動態(tài)的過程，隨著開采面的變化，對應(yīng)的開采區(qū)疏水面積和疏水層位也發(fā)生變化，因此礦井的涌水量也是動態(tài)變化的。而以往預(yù)測礦井涌水量時往往較少考慮礦井的開采進(jìn)度和開采面的變化，因此導(dǎo)致預(yù)測礦井涌水量與實(shí)測值差距較大，且無法反映礦井涌水量的動態(tài)規(guī)律，因此考慮礦井采掘進(jìn)度的動態(tài)涌水量預(yù)測會更加真實(shí)準(zhǔn)確地反映礦井開采過程中涌水量實(shí)際情況。

目前礦井涌水量的預(yù)測方法主要有水文地質(zhì)比擬法[2]、時間序列法[3-4]、解析法[5]、大井法[6]、神經(jīng)網(wǎng)絡(luò)法[7-8]、數(shù)值法[11-15]等。其中數(shù)值法可以模擬礦區(qū)復(fù)雜的三維地質(zhì)結(jié)構(gòu)，并且可以概化仿真采掘進(jìn)度，使涌水量的模擬預(yù)測更為精確。因此以新疆北部某礦區(qū)為例，在利用數(shù)值法的基礎(chǔ)上考慮礦區(qū)的采掘進(jìn)度，運(yùn)用GMS數(shù)值模擬軟件對礦區(qū)涌水量進(jìn)行動態(tài)模擬預(yù)測[14]，預(yù)測結(jié)果與實(shí)際情況較為符合。

1 研究區(qū)概況

研究區(qū)位于新疆北部吐魯番盆地一帶，地形屬于西北邊緣低山丘陵地帶，礦區(qū)內(nèi)絕大部分為斜平地，坡度50‰，北高南低。區(qū)內(nèi)干旱少雨，降雨對地下水補(bǔ)給微弱，侏羅系基巖裂隙含水層主要接受來自北部山前雪山融水的深層補(bǔ)給，并向南部以地下徑流的方式排泄。地下水流動方向主要為東北向西南。區(qū)域廣泛覆蓋第四系松散堆積物，主要斷層導(dǎo)水性較弱，地層構(gòu)造較為簡單。區(qū)內(nèi)主要地層自下而上分別為石炭系中統(tǒng)、三疊系上統(tǒng)、侏羅系下統(tǒng)八道灣組、侏羅系下統(tǒng)三工河組、侏羅系中統(tǒng)西山窯組、侏羅系中統(tǒng)頭屯河組、第三系及第四系。礦區(qū)的含水層組主要為第四系松散孔隙含水層、下第三系陸相半膠結(jié)孔隙裂隙含水層、侏羅系基巖裂隙含水層和石炭系含水層。主采煤層賦存于侏羅系西山窯組含水層，含水層巖性為灰白色粗砂巖-細(xì)砂巖與灰色粉砂巖、泥巖互層，偶夾砂礫巖，含5號煤組，平均總厚度12.1 m，含煤系數(shù)5.38%。礦區(qū)主要分北一采區(qū)、北二采區(qū)和北三采區(qū)三個采區(qū)，采區(qū)內(nèi)開采順序?yàn)楹笸耸?，即由井田邊界向上山方向回采，區(qū)段及水平開采順序?yàn)橄滦惺?，即先采上水平區(qū)段，后采下水平區(qū)段。各采區(qū)分布如圖1所示。

圖1 采區(qū)分布圖Fig.1 Mining area distribution map

2 水文地質(zhì)模型的建立

2.1 模擬范圍及邊界條件的概化

礦區(qū)北部邊界為山前地下水徑流補(bǔ)給邊界，作為定流量邊界；南部邊界為地下水排泄邊界作為定流量邊界；東西邊界為與流線垂直的邊界，作為隔水邊界處理；模型頂部覆蓋大面積第四系松散堆積物，接受大氣降雨入滲，概化為補(bǔ)給邊界；煤層所在西山窯組含水層以下緊鄰著三工河組黑色粉砂巖致密隔水層，巖層致密，隔水性能極好，將其作為模型的底部邊界，概化為隔水邊界。

2.2 模擬區(qū)源匯項(xiàng)的概化

1) 降雨蒸發(fā)。礦區(qū)降雨稀少，根據(jù)礦區(qū)附近氣象站的監(jiān)測數(shù)據(jù)，礦區(qū)年平均降水量為41.4 mm，地表第四系巖性主要為中粗砂，根據(jù)地下水流數(shù)值模擬技術(shù)要求，降雨入滲系數(shù)取經(jīng)驗(yàn)值0.3，通過Recharge模塊進(jìn)行概化(表1)。

礦區(qū)蒸發(fā)強(qiáng)度較大，年均蒸發(fā)量為1 500 mm，地下水水位埋深普遍大于5 m，根據(jù)地下水流數(shù)值模擬技術(shù)要求，蒸發(fā)系數(shù)取經(jīng)驗(yàn)值0.01，通過Evaporation模塊進(jìn)行概化(表2)。

表1 區(qū)域多年平均降雨量表Table 1 Regional average annual rainfall

表2 區(qū)域多年平均蒸發(fā)量表Table 2 Regional average annual evaporation

2) 邊界補(bǔ)給和排泄。研究區(qū)地下水補(bǔ)給主要來自于北部山前雪山融水的地下側(cè)滲補(bǔ)給，并由南部地下水流出邊界排泄。根據(jù)斷面法結(jié)合達(dá)西定律初步估計北部補(bǔ)給邊界流量約為q補(bǔ)給=1 489.64 m3/d；南部排泄邊界的流量約為q排泄=-1 242.15 m3/d，同時流量值在后期流場校核檢驗(yàn)過程中再不斷調(diào)整優(yōu)化以使其符合實(shí)際情況。

2.3 地層結(jié)構(gòu)概化

依據(jù)礦區(qū)水文地質(zhì)資料及鉆孔勘探數(shù)據(jù)，將礦區(qū)水文地質(zhì)模型地層結(jié)構(gòu)劃分為6層含(隔)水層(圖2)，自上而下主要為第四系松散孔隙含水層(H-I-1)、下第三系陸相半膠結(jié)孔隙裂隙含水層(H-I-2)、侏羅系頭屯河組承壓裂隙含水層(H-Ⅱ-1)、西山窯組上段5-3煤層上部至5-1煤層底板間承壓裂隙含水層(H-Ⅱ-2)、西山窯組下段5-1煤層底板以下粗粒相承壓裂隙含水層(H-Ⅱ-3)、三工河組粉砂巖致密隔水層(G-Ⅱ-1)。各地層在空間上由于巖性的差異，水文地質(zhì)參數(shù)表現(xiàn)出差異性，因此將各地層概化為非均質(zhì)各向異性地層，模型概化為6層非均質(zhì)各向異性三維非穩(wěn)定流模型。

圖2 區(qū)域地層綜合柱狀圖Fig.2 Regional stratigraphic comprehensive borehole histogram

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

根據(jù)礦區(qū)水文地質(zhì)概念模型，礦區(qū)地下水滲流符合達(dá)西定律，結(jié)合連續(xù)性方程可建立礦區(qū)地下水非穩(wěn)定流數(shù)學(xué)模型見式(1)[16]。

式中：H為水頭，m；Kxx、Kyy、Kzz為x、y、z方向上的滲透系數(shù)，m/d；n為法向向量；Ss為彈性釋水率，1/d；W為降雨入滲補(bǔ)給強(qiáng)度和蒸散發(fā)強(qiáng)度，m3/d；H0為模擬區(qū)初始流場，m；q為第二類定流量邊界流量，m3/d，流入為正，流出為負(fù)；Ω為模擬區(qū)。

4 數(shù)值模型建立

4.1 模擬區(qū)域網(wǎng)格剖分及時間步長的確定

在滿足模型精度要求前提下，將模型在橫向剖分300個網(wǎng)格，縱向剖分280個網(wǎng)格，垂向有6層含(隔)水層，因此按6層剖分，最終得到300×280×6個單元格(圖3)。

模擬過程中選取2007年1月～2007年12月作為模型的識別期，識別調(diào)整模型水文地質(zhì)參數(shù)；以2008年1月～2008年12月作為模型的檢驗(yàn)期，檢驗(yàn)計算流場和水位的準(zhǔn)確性；礦區(qū)三個采區(qū)采掘期為2009年1月～2019年1月，將2009年1月～2020年1月作為涌水量動態(tài)模擬期，時間劃分以月為步長。

4.2 水文地質(zhì)參數(shù)分區(qū)

根據(jù)礦區(qū)水文地質(zhì)資料，將各地層參數(shù)分區(qū)概化如圖4所示，同時各分區(qū)的水文地質(zhì)參數(shù)初始值由鉆孔抽水試驗(yàn)等資料初步給定，后期根據(jù)各井點(diǎn)地下水位實(shí)測值進(jìn)行校核調(diào)參確定最終校正值。

圖3 模型網(wǎng)格剖分圖Fig.3 Mesh of the numerical simulation model

圖4 各含(隔)水層水文地質(zhì)參數(shù)分區(qū)圖Fig.4 Aquifers hydrogeological parameter zoning map

4.3 采掘空間的概化

為了較為準(zhǔn)確地刻畫出礦井采掘過程中的涌水，運(yùn)用drain模塊刻畫礦井的排水[17-18]，三個采區(qū)分為10個工作面，每個工作面按計劃在一年內(nèi)完成開采。開采前先對工作面單元格進(jìn)行短期疏水，將水位迅速下降至工作面的安全開采平面，然后進(jìn)行開采。對于已開采和正在開采的工作面均設(shè)置排水(drain)邊界，保持地下水位低于該工作面的安全開采平面，未開采的工作面則不進(jìn)行排水，通過這樣的方式可以較為合理準(zhǔn)確地刻畫礦區(qū)開采過程中地下水流場的變化。

模型排水(drain)邊界的計算見式(2)。

(2)

式中：QD為排水量，m3/d，即由含水層流入drain的水量；CD為采取滲透性系數(shù)，m2/d；HD為排水標(biāo)高，m。

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

運(yùn)用礦區(qū)水位觀測孔2007年的水位監(jiān)測數(shù)據(jù)對模型計算的流場進(jìn)行識別，最終調(diào)整得到礦區(qū)水文地質(zhì)模型各分區(qū)水文地質(zhì)參數(shù)見表3。

將識別后的水文地質(zhì)參數(shù)和概化的水文地質(zhì)條件代入模型，最終計算得到天然狀態(tài)下地下水流場計算結(jié)果如圖5所示，計算流場和實(shí)測流場在水位、流向、流態(tài)上基本一致，計算流場能基本反映實(shí)際流場的宏觀形態(tài)。

表3 模型水文地質(zhì)參數(shù)一覽表Table 3 Hydrogeological parameters of the hydrogeology model

圖5 計算流場與實(shí)測流場對比圖Fig.5 Comparison diagram of calculated and observed flow field

選取礦區(qū)6號觀測孔為例，用2008年檢驗(yàn)期地下水位動態(tài)監(jiān)測數(shù)據(jù)對模型計算水位進(jìn)行動態(tài)檢驗(yàn)，水位擬合曲線如圖6所示，結(jié)果顯示模型計算水位與觀測孔實(shí)測動態(tài)水位基本擬合，計算水位可以基本反映實(shí)際水位的動態(tài)變化規(guī)律，模型整體可靠度較高，可作下一步預(yù)測分析。

圖6 地下水水位計算值與實(shí)測值擬合圖Fig.6 Fitting curves of observed and calculated groundwater level

5 礦井涌水量的動態(tài)預(yù)測

5.1 考慮采掘進(jìn)度的礦井涌水量動態(tài)模擬預(yù)測

運(yùn)用識別檢驗(yàn)后的數(shù)值模型對礦區(qū)三個采區(qū)的涌水量進(jìn)行下一步的模擬預(yù)測。三個采區(qū)按照北一采區(qū)、北二采區(qū)、北三采區(qū)順序依次開采，北一采區(qū)安全開采面底板高度+450 m，北二采區(qū)安全開采面底板高度+430 m，北三采區(qū)安全開采面底板高度+400 m，按照由高到低、自上而下的原則開采。

北一采區(qū)分三個工作面，按計劃三年開采完畢；北二采區(qū)分為四個工作面，按計劃四年開采完畢；北三采區(qū)分為三個工作面，按計劃三年開采完畢。對于已開采及正在開采的工作面，地下水位疏干至安全開采平面以下，還未開采的工作面則不設(shè)置排水。按照采掘先后的順序模擬預(yù)測，煤層所在西山窯組含水層的水位變化如圖7～9所示。

圖7 北一采區(qū)疏干后水位等值線圖Fig.7 Contour map of groundwater after dredging in north 1 mining area

圖8 北一采區(qū)、北二采區(qū)疏干后水位等值線圖Fig.8 Contour map of groundwater after dredging in north 1 and north 2 mining area

圖9 三采區(qū)全疏干后水位等值線圖Fig.9 Contour map of groundwater after dredging in three mining areas

通過Flow budget模塊計算出了三個采區(qū)的礦井涌水量隨時間的動態(tài)變化曲線如圖10所示。從計算結(jié)果來看，在一個采區(qū)的開采過程中，當(dāng)一個工作面開采完成，開始開采下一個工作面時，由于開采新的工作面需要迅速將水位疏降至安全開采面以下，因此短時間大量排水，涌水量突增，隨后將水位維持在安全開采面以下，涌水量會慢慢降至平穩(wěn)。最終隨著開采面的不斷擴(kuò)大，總體上涌水量在不斷增大。此外由于北二采區(qū)、北三采區(qū)鄰近北一采區(qū)，因此其疏水過程中將分擔(dān)北一采區(qū)部分涌水量，從而使得北一采區(qū)涌水量進(jìn)一步減小。而由于北三采區(qū)與北二采區(qū)間隔一定距離，因此北三采區(qū)疏水對北二采區(qū)的涌水量影響較小，北二采區(qū)涌水量下降幅度也較小。由此說明相鄰采區(qū)在疏水過程中由于流場的相互干擾和影響，因此涌水量也會互相影響，鄰近采區(qū)的疏水會分擔(dān)一部分先前開采區(qū)的靜儲水量，使得先采區(qū)的涌水量減小。而隨著流場的穩(wěn)定，各采區(qū)的涌水量也逐漸趨于穩(wěn)定。最終計算得到模擬期末北一采區(qū)涌水量為211 m3/d、北二采區(qū)涌水量為613 m3/d、北三采區(qū)涌水量為563 m3/d。

圖10 采掘期內(nèi)三采區(qū)計算涌水量動態(tài)變化曲線Fig.10 Curves of calculated dynamic mine inflow in three mining areas during the mining period

表4 采掘期逐年平均涌水量統(tǒng)計表
Table 4 Statistical table of average annual mine inflow during the mining period

年份20092010201120122013平均涌水量/(m3/d)309457556763899年份20142015201620172018平均涌水量/(m3/d)9761 0531 2791 4011 456

同時為了探究礦區(qū)開采過程中總體涌水量的變化規(guī)律，列出開采過程中逐年的平均涌水量見表4，并作變化曲線如圖11所示，結(jié)果表明礦區(qū)涌水量隨著開采的深入，總體呈現(xiàn)上升趨勢，但中間各個采區(qū)過渡時涌水量有較大的增幅，呈現(xiàn)階梯狀上升，原因是由于北一采區(qū)至北三采區(qū)的安全開采平面高度呈梯度下降，因此由已采區(qū)向未采區(qū)過渡時，疏水降深也隨之增大，涌水量會有短期較大幅度的增長，但隨著新采區(qū)的逐步開采，涌水量的上升幅度逐漸變緩，涌水量慢慢趨于穩(wěn)定。

5.2 未考慮采掘進(jìn)度的礦井涌水量預(yù)測

不考慮采掘進(jìn)度的情況下，則不考慮開采的實(shí)際順序，將三個采區(qū)的地下水位同時疏降至各采區(qū)的安全開采平面以下，三個采區(qū)涌水量動態(tài)變化曲線如圖12所示。三個采區(qū)初期均大量排水，后期隨著疏水的穩(wěn)定，涌水量慢慢減小并趨于穩(wěn)定。最終計算得到模擬期末北一采區(qū)涌水量為236 m3/d、北二采區(qū)涌水量為671 m3/d、北三采區(qū)涌水量為611 m3/d。

圖11 采掘期內(nèi)逐年平均涌水量變化曲線Fig.11 Curve of annual average mine inflow during the mining period

圖12 采掘期內(nèi)三采區(qū)計算涌水量動態(tài)變化曲線Fig.12 Curves of calculated dynamic mine inflow in three mining areas during the mining period

5.3 大井法計算礦井涌水量

根據(jù)礦井的水文地質(zhì)資料，采用大井法計算三個采區(qū)涌水量。將各采區(qū)巷道系統(tǒng)的排水視為理想大井，周界為統(tǒng)一的降落漏斗，由于侏羅系西山窯組含水層本為承壓含水層，隨著開采疏水作用，采區(qū)水位慢慢下降，含水層水位下降至頂板以下，開采區(qū)域由承壓轉(zhuǎn)無壓，因此采用承壓轉(zhuǎn)無壓完整井公式(式(3))。最終得到大井法計算的礦井涌水量,計算成果見表5。

(3)

式中：K為滲透系數(shù)，m/d；H為自然水位至含水層底板距離，m；M為含水層厚度，m；h0為疏干水位至含水層底板距離，m；R0為引用影響半徑，m；r0為引用半徑，m；S為水位降深，即自然水位至疏干水位的距離，m；P為預(yù)測區(qū)面積，km2。

5.4 分析結(jié)果

最終將考慮采掘進(jìn)度預(yù)測礦井涌水量、未考慮采掘進(jìn)度預(yù)測礦井涌水量、大井法計算礦井涌水量和實(shí)測涌水量進(jìn)行對比，將各種方法計算涌水量及其與實(shí)測涌水量的相對誤差列于表6。由對比結(jié)果分析，考慮采掘進(jìn)度模擬礦井涌水量最為接近礦井的實(shí)測涌水量，相對誤差最小，因此考慮采掘進(jìn)度模擬礦井涌水量相對其他方法精度更高，并且可以反映礦區(qū)涌水量的動態(tài)變化規(guī)律，相比其他方法具有一定的優(yōu)越性。

表5 大井法礦井涌水量計算成果表Table 5 Table of mine inflow calculated by large well method

表6 各方法預(yù)測礦井涌水量統(tǒng)計對比表Table 6 Statistical comparison table of the mine inflow predicted by various methods

6 結(jié) 論

1) 開采新工作面時由于需短期將水位疏排并維持在安全開采面以下，涌水量先期會突增，而后保持工作面水位不變進(jìn)行排水，涌水量慢慢降至平緩。新采區(qū)的開采疏水會分擔(dān)相鄰采區(qū)的靜儲水量，導(dǎo)致鄰近采區(qū)的涌水量減小，采區(qū)距離越近，互相影響越為明顯。

2) 礦區(qū)涌水量隨著開采的進(jìn)行總體上是逐年上升的。在舊采區(qū)向新采區(qū)過渡時由于安全開采面的下降，涌水量會有短期較大幅度的增長，后期隨著該新采區(qū)的開采，涌水量增長變緩，慢慢趨于穩(wěn)定。

3) 最終將各方法計算的涌水量和實(shí)測涌水量進(jìn)行對比分析，發(fā)現(xiàn)考慮采掘進(jìn)度計算得到的礦井涌水量與實(shí)測值相對誤差最小，相比于其他方法準(zhǔn)確性要高，能更真實(shí)反映礦井涌水量情況，為礦井后期的水災(zāi)防治和安全生產(chǎn)提供一定的參考依據(jù)。