殷曉曦， 陳陸望， 林曼利， 劉 鑫， 桂和榮

（1.合肥工業(yè)大學(xué) 資源與環(huán)境工程學(xué)院，安徽 合肥 230009；2.宿州學(xué)院 地球科學(xué)與工程學(xué)院，安徽 宿州 234000）

地下水系統(tǒng)是指受控于一定水文地質(zhì)條件，在某一空間域形成的與相鄰地質(zhì)體相對隔絕而內(nèi)部水力聯(lián)系密切的含水系統(tǒng)或單元，既包括地下水含水系統(tǒng)又包括地下水流系統(tǒng)［1－2］。水文地質(zhì)條件的復(fù)雜性（地形起伏多變、含水層厚度不一、非均質(zhì)、各向異性、多種形式源匯項）和水文地質(zhì)參數(shù)的時空變異性，在很大程度上限制了人們對地下水系統(tǒng)的認識和量化［3－4］。利用數(shù)值模擬軟件對地下水流等問題進行模擬的方法以其有效性、靈活性和相對廉價性逐漸成為地下水研究領(lǐng)域的一種不可缺少的重要方法，得到了廣泛應(yīng)用［5－6］。國外對地下水系統(tǒng)模擬研究較早，理論及軟件技術(shù)發(fā)展比較成熟，眾多學(xué)者對地下水系統(tǒng)數(shù)值模擬作了大量研究［7－9］。我國地下水數(shù)值模擬起步于20世紀70年代。隨著新技術(shù)、新方法的廣泛應(yīng)用，國內(nèi)多學(xué)者在理論和方法上不斷創(chuàng)新，并結(jié)合相關(guān)研究方向理論，提高了模擬結(jié)果的可靠性［10－15］。但從上述研究看，其所建模型基本上以一維流、平面二維流、擬三維流為主，未能完全真實地對含水層進行刻畫，尤其是在高水壓、構(gòu)造復(fù)雜且受采動影響的礦井，更難以真實地刻畫其含水層。

本文在系統(tǒng)分析煤礦地質(zhì)與水文地質(zhì)條件的基礎(chǔ)上，經(jīng)過對地質(zhì)體的科學(xué)概化，建立研究區(qū)水文地質(zhì)概念模型、數(shù)學(xué)模型與數(shù)值模型，并利用三維有限元地下水數(shù)值模擬軟件FEFLOW對模型進行求解；著重考慮采動影響以及煤礦特殊水文地質(zhì)條件，利用煤礦觀測孔的長期觀測資料對建立的模型進行驗證與標識；通過驗證與標識后的模型重新展現(xiàn)近期采動影響下煤礦主要充水含水層地下水滲流場的平面與空間分布規(guī)律，為煤礦深部地下水的防治提供理論依據(jù)。

1 研究區(qū)地質(zhì)及水文地質(zhì)概況

淮北煤田中南部平原地區(qū)均為新第三系和第四系松散層覆蓋；基底為太古代和早元古代深、中深變質(zhì)巖系及中元古代淺變質(zhì)巖系組成；蓋層是穩(wěn)定地臺型沉積，有上元古界青白口系至古生界二疊系；中生代侏羅系、白堊系及老第三系主要分布在斷陷盆地內(nèi)；區(qū)內(nèi)缺失晚奧陶系至下石炭系和三疊系地層。研究區(qū)任樓煤礦位于淮北煤田童亭背斜東南翼和南翼，如圖1所示，北起界溝斷層，南至F8斷層與許疃煤礦比鄰，東南部以F31斷層為界，淺部以11號煤層露頭為界，深部至31煤層－800等高線的平面投影，是一獨立的水文地質(zhì)單元，與鄰近礦井水力聯(lián)系微弱。任樓煤礦走向長9.8～14km，傾斜寬1.2～3.5km，總面積約43km2。任樓煤礦含水層從上至下為：①松散層第4含水層，簡稱“四含”，由含水層沉積環(huán)境與基巖面控制，埋深300m左右；② 二疊系煤系砂巖裂隙含水層，簡稱“煤系”，與開采煤層埋深有關(guān)，其中3煤至4煤、4煤至7煤、7煤至8煤、鋁質(zhì)泥巖至11煤之間砂巖裂隙較發(fā)育，可分別劃分為不同的含水組；③ 石炭系太原組巖溶含水層，簡稱“太灰”，埋深600m左右，一般劃分為14～15含水組；④ 及奧陶系巖溶含水層，簡稱“奧灰”，埋深900m左右，勘探不詳，埋藏較深，水文地質(zhì)與區(qū)域構(gòu)造復(fù)雜，為煤礦主要充水含水層。

圖1 任樓煤礦地質(zhì)概況

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

2.1 含水層結(jié)構(gòu)概化

任樓煤礦四含地下水的動、靜儲量十分豐富，與煤系、太灰與奧灰含水層均呈角度不整合接觸，不僅能順層正常補給煤系含水層，而且也與太灰、奧灰含水層之間有著極為廣闊的互補關(guān)系。為了提高水文地質(zhì)數(shù)值模擬的精度，把整個煤礦地下水看作一個完整的系統(tǒng)。煤系、太灰與奧灰含水層巖層裂隙發(fā)育，斷裂構(gòu)造復(fù)雜，含水組間的連通性較好，可分別視為統(tǒng)一的含水子系統(tǒng)，將其如四含子系統(tǒng)一樣概化為單層結(jié)構(gòu)。地下水流整體上以水平運動為主、垂向運動為輔，地下水系統(tǒng)的垂向運動主要是層間的越流。地下水系統(tǒng)的輸入、輸出隨時空變化，故地下水為非穩(wěn)定流。參數(shù)隨空間變化，體現(xiàn)了系統(tǒng)的非均質(zhì)性。根據(jù)采動影響下煤礦水文地質(zhì)以及地下水補給和動態(tài)變化特征，將煤礦概化成非均質(zhì)、各向異性、空間三維結(jié)構(gòu)的地下水流系統(tǒng)。

2.2 邊界條件概化

任樓煤礦四含在側(cè)向上主要受其本身厚度及形成時各種沉積特征影響，表現(xiàn)為北部、南部和西部側(cè)向補給，屬于流量邊界；而在東部，四含厚度減小，富水性差，多數(shù)地段趨于尖滅，屬于隔水邊界。煤系地下水以靜儲量為主，邊界條件設(shè)為隔水邊界。太灰含水層在平面上北以界溝斷層為界，該斷層落差大，走向延伸長度大于3650m，斷層帶內(nèi)巖石具有擠壓、破碎現(xiàn)象，泥質(zhì)充填，可視為隔水邊界；南至F7斷層，該斷層在平面走向上呈“S”彎曲，鉆孔及斷層處巖心擠壓破碎，裂隙發(fā)育，具有滑面和擦痕，根據(jù)以往抽水試驗及勘探資料分析，將其作為流量邊界；西自11煤層露頭線，東至31煤層－800m走向的地表投影，東南角以F23斷層為界，均設(shè)為流量邊界。由于奧灰含水層距離開采煤層較遠，設(shè)置奧灰東部、南部與北部為隔水邊界，西部為透水邊界。

3 數(shù)學(xué)模型建立與地下水數(shù)值模擬模型

3.1 數(shù)學(xué)模型的建立

采用三維數(shù)學(xué)模型模擬研究區(qū)含水層系統(tǒng)地下水流動［6］，數(shù)學(xué)方程如下：

其中，kxx、kyy、kzz分別為沿x、y、z方向的滲透系數(shù)；H為地下水系統(tǒng)水位標高，是時間和空間上的變化函數(shù)，H＝H（x，y，z，t）；Ss為點（x，y，z）處的儲水率；W 為源匯項，是時間與空間的函數(shù)，W＝W（x，y，z，t）；t為時間；Ω 為滲流計算域；Γ2為流量邊界；nx、ny、nz分別為邊界Γ2沿x、y、z軸方向單位矢量；q（x，y，z，t）為流量在時間和空間上的變化函數(shù)。

3.2 地下水的數(shù)值模擬模型

3.2.1 模擬軟件選取及網(wǎng)格剖分

采用德國WASY公司開發(fā)的基于有限元法的FEFLOW軟件，建立研究區(qū)地下水流數(shù)值模擬模型。首先要對模擬區(qū)進行三角形剖分，剖分單元46060個，節(jié)點26992個，相應(yīng)空間上則是由6節(jié)點組成的三棱柱體構(gòu)成的空間地質(zhì)體。

3.2.2 模擬期的處理

本次研究以2008年1月1日流場作為地下水模擬的初始流場，2009年12月31日流場作為地下水模擬的擬合流場。擬合時段為2008年1月1日—2009年12月31日2個水文年，以1個月作為一個時間段，每個時間段內(nèi)包括若干時間步長，時間步長為模型自動控制。

3.2.3 定解條件的處理

以2008年1月1日觀測的地下水水位為基礎(chǔ)，采用內(nèi)插和外推法獲得各含水層的初始水位。各個流量邊界的參數(shù)主要考慮模擬初期與模擬末期的流場，有詳細資料的邊界，擬合邊界流入與流出量。時間步長由程序控制，每一步運算都嚴格控制誤差。在與采動影響區(qū)對應(yīng)的煤系含水層中分配若干單井，按照年平均礦井涌水量分配單井流量。井下太灰含水層施工的風(fēng)1、風(fēng)2、風(fēng)3與風(fēng)4按照年平均放水量平均分配井流量。

3.2.4 水文地質(zhì)參數(shù)概化

根據(jù)研究區(qū)地質(zhì)與水文地質(zhì)條件及以往抽（注）水資料，結(jié)合數(shù)值模擬模型中的含隔水層結(jié)構(gòu)及其與采動的關(guān)系，對四含與太灰含水層視各向異性取參數(shù)，劃分為多個區(qū)，如圖2所示；相對于四含與太灰，奧灰含水層富水性強、煤系含水層與隔水層富水性弱，各自參數(shù)的微小變化并不影響數(shù)值模擬結(jié)果，故把煤系與奧灰以及隔水層視各向同性取參數(shù)，分別看作1個區(qū)。對應(yīng)的主要參數(shù)有滲透系數(shù)與儲水率。煤系頂板隔水層如受采動影響，相應(yīng)水文參數(shù)擴大為原來的5倍或10倍。根據(jù)研究區(qū)斷層含水性極弱、導(dǎo)水性差的特點，選取落差大于30m的典型斷層（邊界斷層除外），其滲透系數(shù)為斷層所貫穿巖層滲透系數(shù)中的最小值。2個陷落柱用增加陷落范圍內(nèi)垂向上的滲透系數(shù)的方法模擬（受陷落柱影響的隔水層滲透系數(shù)一般增加到原來的10倍，含水層滲透系數(shù)一般增加到原來的2倍）。采動影響下研究區(qū)水文地質(zhì)參數(shù)見表1所列。

圖2 四含與太灰參數(shù)分區(qū)

表1 研究區(qū)水文地質(zhì)參數(shù)

3.2.5 模型的識別與驗證

模型的識別與驗證過程是把數(shù)值模擬得到的結(jié)果和地面觀測孔長期水位動態(tài)觀測資料進行比較，看兩者所反映的動態(tài)變化規(guī)律是否相近。若不相近，就要修正模型及相關(guān)參數(shù)，直到滿意為止［3］。典型地面觀測孔計算水位和實測水位擬合過程曲線如圖3所示，其中橫軸1～12分別表示日 期：2008－01、2008－03、2008－05、2008－07、2008－09、2008－11、2009－01、2009－03、2009－05、2009－07、2009－09、2009－11。從圖3可知，所建立的模擬模型基本達到模型精度要求，符合研究區(qū)水文地質(zhì)條件與采動特征，基本反映了地下水系統(tǒng)動態(tài)。

圖3 地面觀測孔水位曲線擬合圖

4 數(shù)值模擬結(jié)果分析

4.1 主要充水含水層平面滲流場特征

2009年12月底研究區(qū)主要充水含水層平面滲流場特征，如圖4所示（箭頭表示滲流速度矢量，云圖表示水位標高）。從圖4可以看出，四含水由于受淺部煤層開采影響，靠近東部邊界出現(xiàn)明顯的低水位分布區(qū)，采動影響區(qū)地下水滲流速度大，礦井涌水量主要來自研究區(qū)的西部邊界。煤系水滲流速度小，在研究區(qū)內(nèi)難形成低水位分布區(qū)，但在采動影響范圍內(nèi)典型斷層局部附近出現(xiàn)較高的滲流速度。太灰水滲流速度大，但受采動影響小，地下水滲流方向與井下放水密切相關(guān)，其中部地下水滲流方向指向井下放水孔的影響區(qū)域，南部與北部地下水滲流方向從西向東，斷層與陷落柱是影響其特征的主要因素。奧灰水滲流速度最大，幾乎不受采動影響，在研究區(qū)的南部與北部由于受典型斷層控制滲流方向近似由西向東，但在研究區(qū)中部，由于受2個陷落柱的控制，滲流方向由北向南。

圖4 主要充水含水層平面滲流特征圖

4.2 采動影響下主要充水含水層空間滲流場

主要充水含水層空間滲流特征如圖5所示。從圖5a空間地下水位等值線可以看出，奧灰與太灰之間的隔水層空間水位等值線與頂?shù)装迤叫星颐芏却?，層間越流顯著，典型斷層與陷落柱是其控制的主因。太灰與奧灰含水層空間水位等值線均以垂直分布為主，等值線密度奧灰大于太灰，2個含水層在研究區(qū)范圍內(nèi)主要以平面滲流為主，滲流速度奧灰大于太灰。太灰與煤系含水層之間的隔水層空間水位等值線與頂?shù)装褰谄叫星颐芏容^小，但在采動區(qū)（特別是采動區(qū)內(nèi)有典型斷層）等值線密度加大，水力聯(lián)系顯著增強。煤系含水層垂向等值線密度稀疏，含水層平面滲流速度緩慢，但在采動影響區(qū)內(nèi)垂向等值線發(fā)生偏轉(zhuǎn)，與含水層頂?shù)装褰谄叫?，煤系水直接涌入礦坑。四含與煤系之間隔水層空間水位等值線也近于水平，但在采動影響區(qū)內(nèi)水平等值線密度最大，采動影響顯著，四含水大量涌入礦坑。從圖5b滲流速度空間矢量可以看出，受采動影響含水層主要是四含與煤系，太灰與奧灰影響較小。研究區(qū)內(nèi)典型斷層或陷落柱雖然一定程度地控制著含水層滲透速度矢量的分布，但四含與煤系水流場主要取決于采動影響范圍。

圖5 主要充水含水層空間滲流特征圖

5 結(jié) 論

從平面角度分析，四含由于受淺部煤層開采影響，東部出現(xiàn)低水位分布區(qū)；煤系含水層采動影響范圍內(nèi)典型斷層局部附近滲流異常；太灰與奧灰含水層受采動影響小，斷層與陷落柱是影響太灰與奧灰水滲流的主要因素。從空間角度分析，太灰與奧灰含水層間越流顯著，典型斷層與陷落柱是其控制的主因；煤系含水層空間水位垂向等值線密度稀疏，但在采動影響區(qū)垂向等值線發(fā)生偏轉(zhuǎn)，與含水層頂?shù)装褰谄叫?；四含與煤系之間隔水層空間水位等值線也近于水平，但在采動區(qū)內(nèi)水平等值線密度最大。采動對研究區(qū)四含與煤系2個含水層影響最為顯著。

［1］王大純，張人權(quán)，史毅虹，等.水文地質(zhì)學(xué)基礎(chǔ)［M］.北京：地質(zhì)出版社，1995：81－86.

［2］方向清，傅耀軍，華解明，等.北方巖溶地下水系統(tǒng)模式及特征［J］.合肥工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2011，34（2）：286－291.

［3］張繼成，張旭昇，胡雅杰，等.基于GIS技術(shù)與FEFLOW的秦川王盆地南部地下水數(shù)值模擬［J］.蘭州大學(xué)學(xué)報：自然科學(xué)版，2010，46（5）：31－38.

［4］陳 舟，錢家忠，姜常讓，等.多孔介質(zhì)地下水非達西滲流研究進展［J］.合肥工業(yè)大學(xué)學(xué)報：自然科學(xué)版，2008，31（10）：1539－1543.

［5］賀國平，邵景力，崔亞莉，等.FEFLOW在地下水流模擬方面的應(yīng)用［J］.成都理工大學(xué)學(xué)報：自然科學(xué)版，2003，30（4）：356－361.

［6］董東林，王存社，陳書客，等.典型煤礦地下水運動及污染數(shù)值模擬：Feflow及 Modflow應(yīng)用［M］.北京：地質(zhì)出版社，2010：64－68.

［7］Scheibe T，Yabusaki S.Scaling of flow and transport behavior in heterogeneous groundwater systems［J］.Advances in Water Resources，1998，22（3）：223－238.

［8］Sarwar A，Helmut E.Development of a conjunctive use model to evaluate alternative management options for sur－face and groundwater resources［J］.Hydrogeology Journal，2006，14：1676－1678.

［9］Disersch H J，Kolditz O.Coupled groundwater flow and transport：2.thermohaline and 3Dconverction systems［J］.Advances in Water Resources，1998，21（5）：401－425.

［10］周念清，朱 蓉，朱學(xué)愚.MODFLOW在宿遷市地下水資源評價中的應(yīng)用［J］.水文地質(zhì)工程地質(zhì)，2000，27（6）：9－13.

［11］楊青春，盧文喜，馬洪云.Visual－Modflow在吉林省西部地下水數(shù)值模擬中的應(yīng)用［J］.水文地質(zhì)工程地質(zhì)，2005，32（3）：67－69.

［12］任印國，柳華武，李明良，等.石家莊市東部平原FEFLOW地下水數(shù)值模擬與研究［J］.水文，2009，29（5）：59－62.

［13］張國偉，孫亞軍，徐智敏，等.FEFLOW在區(qū)域地下水流場演化數(shù)值模擬中的應(yīng)用［J］.地下水，2010，32（1）：4－33.

［14］孫志浩，王開章，孔凡亮.基于FEFLOW的淺層地下水水環(huán)境演化模擬［J］.中國農(nóng)村水利水電，2009（10）：5－11.

［15］魏 新，楊國勇，胡林煒.基于FEFLOW的采區(qū)放水實驗數(shù)值模擬與模型評估［J］.地下水，2010，32（5）：20－22.