劉志鵬

（山西省水利水電科學(xué)研究院 山西太原 030002）

1 研究區(qū)概況

1.1 研究區(qū)斜井涌水概況

本次研究的對(duì)象為潘家窯煤礦整合前的原崔家?guī)X礦副斜井涌水問(wèn)題，研究區(qū)范圍選取潘家窯煤礦的東部原崔家?guī)X礦副斜井周?chē)欢ǚ秶?。?jù)現(xiàn)場(chǎng)調(diào)查，潘家窯煤礦原崔家?guī)X礦副斜井內(nèi)距井口約60 m處有三股地下水涌出，進(jìn)入斜井內(nèi)，其中兩股涌水點(diǎn)位于斜井壁腳處，標(biāo)高分別為1 233.34 m和1 230.29 m，分別距井口61.0 m和68.8 m，另一股位于斜井壁頂部，距井口68 m，三處涌水點(diǎn)涌水量分別為13.17 m3/h、40.71 m3/h和4.20 m3/h，總涌水量為58.08 m3/h。

1.2 研究區(qū)范圍及深度的確定

潘家窯煤礦原副斜井位于上窯村七里河谷，副斜井涌水補(bǔ)給屬淺層地下水補(bǔ)給。據(jù)查，淺層地下水的補(bǔ)給來(lái)源主要為附近的地表水和大氣降水，本次研究工作的具體范圍是：北部、東部以七里河為界，西部、南部考慮地下水影響到涌水點(diǎn)的距離，分別距涌水點(diǎn)126 m和200 m，研究區(qū)范圍總面積0.152 km2。

副斜井涌水點(diǎn)出露處的位置埋深約29.77～32.82 m，巖性為二疊系下石盒子組砂巖，其下為泥巖，即涌水出露點(diǎn)位置為砂巖與泥巖接觸帶附近，下石盒子組泥巖為淺層地下水的隔水層，故本次研究深度范圍為地表至二疊系下石盒子組泥巖，上部為孔隙潛水，下部為裂隙水，兩者之間沒(méi)有隔水層，具有統(tǒng)一的水力聯(lián)系，泥巖為隔水層。

1.3 水庫(kù)概況

七里河水庫(kù)位于七里河上游安家?guī)X礦附近干流上，研究區(qū)東南部，距潘家窯煤礦原副斜井約200 m。為小II型水庫(kù)，工程等級(jí)為V等。庫(kù)壩為土壩，壩頂高程為1 252.0 m。中隔墻將庫(kù)區(qū)分為兩區(qū)，右側(cè)行洪河槽總長(zhǎng)為852.42 m，高程由1 247.0 m降至泄洪閘進(jìn)口底板高程1239.0m，左側(cè)為蓄清區(qū)，面積約5.7萬(wàn)m2，高程由1 243.8 m降至1 239.0 m。

2 研究區(qū)水文地質(zhì)結(jié)構(gòu)及涌水來(lái)源、補(bǔ)給途徑分析

2.1 研究區(qū)水文地質(zhì)結(jié)構(gòu)

研究區(qū)淺層地下水含水層從巖性上來(lái)分有兩層：即第四系砂卵石層和二疊系砂巖，前者賦存孔隙水，后者賦存裂隙水。兩者之間沒(méi)有或只有局部隔水層，因而孔隙水和裂隙水具有統(tǒng)一的水力聯(lián)系。砂巖下伏泥巖為隔水層，是以上含水層的隔水底板。含水層、隔水層共同構(gòu)成研究區(qū)淺層地下水儲(chǔ)水結(jié)構(gòu)，雖然地下水含水介質(zhì)不同，但它們具有統(tǒng)一的水力聯(lián)系，從埋藏條件來(lái)分它們屬于潛水類(lèi)型。研究區(qū)淺層地下水具有統(tǒng)一的補(bǔ)給來(lái)源，主要來(lái)源于七里河水的下滲補(bǔ)給、大氣降水的入滲補(bǔ)給以及調(diào)蓄水庫(kù)高水位時(shí)水庫(kù)的滲漏補(bǔ)給等。

研究區(qū)北部、東南部為該淺層地下水的補(bǔ)給邊界，由于七里河上游截流后無(wú)大的洪水流經(jīng)，而平時(shí)各種排水基本常年不變，該補(bǔ)給邊界可按定水頭對(duì)待。地下水補(bǔ)給后向西南方向徑流，為流量邊界。

2.2 副斜井涌水來(lái)源、補(bǔ)給途徑分析

大氣降水和河水滲入地下后首先補(bǔ)給第四系松散層孔隙水，然后隨地下水徑流。其下伏砂巖由于存在裂隙，特別是構(gòu)造裂隙，孔隙水可以下滲進(jìn)入砂巖裂隙中，使裂隙水獲得補(bǔ)充，隨統(tǒng)一的徑流場(chǎng)流動(dòng)。由于研究區(qū)副斜井的掘進(jìn)，正好通過(guò)砂巖中裂隙發(fā)育帶，而副斜井通過(guò)淺層地下水底板掘進(jìn)至煤層，故砂巖裂隙中的地下水就進(jìn)入到副斜井中，形成涌水。由于砂巖裂隙的導(dǎo)水性較強(qiáng)，于是在裂隙發(fā)育帶和副斜井涌水點(diǎn)一帶形成強(qiáng)徑流帶，地下水除順裂隙向西南方向徑流外，很大部分水量涌入到副斜井中。

七里河水庫(kù)滲漏可能進(jìn)入副斜井的途徑與以上類(lèi)似，為水庫(kù)水滲漏進(jìn)入第四系松散層，然后隨地下水徑流通過(guò)導(dǎo)水裂隙進(jìn)入副斜井中。

3 礦井涌水?dāng)?shù)值模擬研究

3.1 水文地質(zhì)條件的概化

邊界條件表示的是研究區(qū)域的地下水系統(tǒng)與其周?chē)吔绛h(huán)境之間的關(guān)系，周?chē)h(huán)境的變化對(duì)研究系統(tǒng)有影響[1]。模擬區(qū)范圍，東部和北部以七里河及水庫(kù)庫(kù)岸為界，西北部、西部和南部為人為邊界，研究區(qū)東西長(zhǎng)約330 m，南北寬約445 m，面積約為118 145 m2。

根據(jù)區(qū)內(nèi)水文地質(zhì)鉆孔資料及地下水存儲(chǔ)條件，將模擬區(qū)內(nèi)含水層概化為一層。在對(duì)基本條件進(jìn)行系統(tǒng)分析的基礎(chǔ)之上，借鑒前人的工作成果進(jìn)行參數(shù)分區(qū)，統(tǒng)一參數(shù)分區(qū)內(nèi)含水層可視為均質(zhì)，水流服從達(dá)西定律。根據(jù)區(qū)內(nèi)的地質(zhì)、水文地質(zhì)條件，結(jié)合水動(dòng)力場(chǎng)特征，將七里河及水庫(kù)邊界概化為第一類(lèi)邊界，其他邊界概化為第二類(lèi)邊界，其中七里河及水庫(kù)庫(kù)岸為補(bǔ)給邊界，其余邊界為排泄邊界。整個(gè)地下水流態(tài)概化為三維非穩(wěn)定流，四周均按通用水頭邊界處理，系統(tǒng)的底部為隔水邊界。根據(jù)實(shí)際長(zhǎng)期觀測(cè)的系列資料，選擇模擬時(shí)段為2010年9月5日至2010年9月30日，以2010年9月5日的水頭分布作為研究區(qū)的初始水位。

3.2 地下水系統(tǒng)數(shù)學(xué)模型的求解原理與方法

3.2.1 求解原理

建立地下水水流模型之后，運(yùn)用Visual Modflow軟件模擬地下水流場(chǎng)時(shí)，首先要對(duì)研究區(qū)進(jìn)行離散化處理，下面將通過(guò)對(duì)含水層的離散化處理來(lái)說(shuō)明水流模型的求解原理。

首先對(duì)含水層進(jìn)行空間離散，即網(wǎng)格剖分。將一個(gè)三維含水層系統(tǒng)劃分為一個(gè)三維的網(wǎng)格系統(tǒng)，整個(gè)實(shí)際含水層被劃分為若干層，每一層又剖分為若干行和若干列。這樣，含水層被剖分成若干個(gè)長(zhǎng)方體，稱(chēng)為計(jì)算單元格，每個(gè)計(jì)算單元格的位置用該計(jì)算單元格所在的行號(hào)（i），列號(hào)（j）和層號(hào)（k）來(lái)表示。每個(gè)剖分出來(lái)的小長(zhǎng)方體的中心稱(chēng)為節(jié)點(diǎn)，一個(gè)計(jì)算單元格的水頭由水頭在該節(jié)點(diǎn)的值所表示[2]。

其次對(duì)時(shí)間進(jìn)行離散。時(shí)間可以根據(jù)實(shí)際情況分為若干個(gè)應(yīng)力期。在每個(gè)應(yīng)力期內(nèi)，所有外部源匯項(xiàng)（如抽水量、入滲量等）的強(qiáng)度保持不變。每個(gè)應(yīng)力期又可以再分為若干個(gè)時(shí)段，每個(gè)時(shí)段內(nèi)的時(shí)間步長(zhǎng)可根據(jù)實(shí)際變化情況選擇。

對(duì)于計(jì)算單元格，從地下水流動(dòng)的連續(xù)性方程出發(fā)，根據(jù)流入和流出計(jì)算單元格的水流量之差等于單元格中貯水量的變化的水均衡原理，導(dǎo)出其六個(gè)界面上的有限差分公式，同時(shí)考慮外部源匯項(xiàng)對(duì)其影響，獲得計(jì)算單元格的地下水滲流計(jì)算有限差分公式。依次列出所有有效單元格的向后差分格式，這些方程構(gòu)成了一個(gè)線性代數(shù)方程組。利用初始條件對(duì)其進(jìn)行迭代求解，就可以獲得每一個(gè)時(shí)段末地下水位的平面分布。

3.2.2 求解方法

由于所計(jì)算的水頭值，既是空間坐標(biāo)的函數(shù)，也是時(shí)間坐標(biāo)的函數(shù)，因此不僅要將數(shù)學(xué)模型進(jìn)行空間上的離散化，同時(shí)也要進(jìn)行時(shí)間的離散化。

地下水運(yùn)動(dòng)的有限差分公式實(shí)際上是根據(jù)地下水流動(dòng)的連續(xù)性方程進(jìn)行的。按照連續(xù)性方程，流入和流出某個(gè)計(jì)算單元的水量之差應(yīng)等于該單元格中貯水量的變化。當(dāng)不考慮地下水密度的變化時(shí)，連續(xù)性方程可簡(jiǎn)單的表示為：

式中：∑Qi——單位時(shí)間內(nèi)流進(jìn)或流出該計(jì)算單元的水量

μ——含水層的給水度(無(wú)量綱)；

△x，△y——計(jì)算單元的空間間距(L)；

△h——某一時(shí)間段內(nèi)水頭的變化(L)；

△t——時(shí)間步長(zhǎng)(T)。

若考慮垂向水量交換，對(duì)某個(gè)計(jì)算單元（i,j)而言，其垂向水量增量可以用下式表示：

式中：Wij——流入計(jì)算單元（i,j)的水量(L3／T)；

ε1i,j——垂向入滲補(bǔ)給量（L3／T)；

ε2i,j——垂向排泄量(L3／T)。

將模擬區(qū)每個(gè)計(jì)算網(wǎng)格的水頭hi,j進(jìn)行表示，可得到n個(gè)線性方程，將它們組成一個(gè)n維的線性方程組，用矩陣表示為：

式中：[A]——水頭的系數(shù)矩陣；

{h}——所求的水頭矩陣；

{q}——表示各個(gè)方程中所包含的常數(shù)項(xiàng)和已知項(xiàng)，也稱(chēng)右端項(xiàng)。

在Visual Modflow中，系數(shù)矩陣和右端項(xiàng)是通過(guò)各個(gè)軟件包來(lái)逐步計(jì)算的，最后通過(guò)迭代法即可求得水頭矩陣的系列解。

3.3 模型的設(shè)計(jì)與實(shí)現(xiàn)

3.3.1 計(jì)算區(qū)剖分

將研究區(qū)在平面上剖分為80行，80列，其網(wǎng)格為矩形，共計(jì)6 400個(gè)單元格,每個(gè)單元格的平面面積約為22.9 m2。

3.3.2 定解條件的確定

根據(jù)2010年9月5日研究區(qū)內(nèi)長(zhǎng)觀孔的監(jiān)測(cè)資料，經(jīng)過(guò)差值計(jì)算求得含水層的滲流場(chǎng)，將其作為研究區(qū)的初始流場(chǎng)，見(jiàn)圖1。

圖1 初始流場(chǎng)圖

3.3.3 地下水均衡要素的確定

在處理地下水系統(tǒng)的源匯項(xiàng)時(shí)，包括補(bǔ)給項(xiàng)和排泄項(xiàng)兩項(xiàng)內(nèi)容。在均衡區(qū)內(nèi)，地下水補(bǔ)給項(xiàng)主要包括大氣降水入滲補(bǔ)給、地下水側(cè)向徑流補(bǔ)給和河道滲漏補(bǔ)給；地下水排泄項(xiàng)主要包括潛水蒸發(fā)和徑流排泄。計(jì)算研究區(qū)內(nèi)各均衡要素的目的在于確定地下水的各個(gè)補(bǔ)給排泄隨時(shí)間和空間的變化規(guī)律，為建立地下水?dāng)?shù)值模擬模型準(zhǔn)備數(shù)據(jù)。

建立涌水量預(yù)測(cè)模型時(shí)涉及的主要參數(shù)為單元的滲透系數(shù)K、給水度μ、儲(chǔ)水系數(shù)S和孔隙比。其中最為關(guān)鍵的是滲透系數(shù)K以及儲(chǔ)水系數(shù)S。由于受構(gòu)造及水文地質(zhì)條件的影響，這些參數(shù)在空間上表現(xiàn)出較強(qiáng)的差異性，因此有必要先進(jìn)行參數(shù)分區(qū)[3]。滲透系數(shù)分區(qū)主要是根據(jù)已知資料和某些已知點(diǎn)以及區(qū)域水文地質(zhì)條件確定的，分區(qū)的取值受人為因素影響比較大[4]。水文地質(zhì)參數(shù)分區(qū)原則上需綜合考慮以下要素：含水層的富水性規(guī)律和補(bǔ)徑排條件、斷裂構(gòu)造的特點(diǎn)以及觀測(cè)孔的數(shù)量和分布[5]。

最后，運(yùn)用Visual Modflow軟件進(jìn)行地下水的數(shù)值模擬，輸入或修改各類(lèi)水文地質(zhì)參數(shù)和幾何參數(shù)，運(yùn)行模型，反演校正參數(shù)，直到顯示輸出結(jié)果。

3.4 模型的識(shí)別與驗(yàn)證

本文運(yùn)用了試錯(cuò)法對(duì)模型參數(shù)進(jìn)行確定，經(jīng)過(guò)反復(fù)調(diào)參，得到較為合理的模型識(shí)別結(jié)果和合理的參數(shù)組合。模型識(shí)別的過(guò)程實(shí)際上是反求參數(shù)以使其與水文地質(zhì)原型相仿真的過(guò)程[6、7]。

3.4.1模型的識(shí)別

根據(jù)研究區(qū)地下水位觀測(cè)資料的實(shí)際情況，選取2010年9月5日至2010年9月30日作為模型識(shí)別階段，以五天作為一個(gè)步長(zhǎng)，分六個(gè)步長(zhǎng)來(lái)計(jì)算研究區(qū)內(nèi)的水位。依據(jù)含水層的巖性、導(dǎo)水性和埋藏條件、厚度以及水位動(dòng)態(tài)特征，按照不同的水文地質(zhì)條件將研究區(qū)分為六個(gè)參數(shù)區(qū)，見(jiàn)圖2，具體分區(qū)方法本文不再贅述。

圖2 參數(shù)分區(qū)圖

將計(jì)算值與已有資料的觀測(cè)值比較，進(jìn)行水位擬合。盡管觀測(cè)區(qū)內(nèi)的觀測(cè)點(diǎn)較多，但是因部分觀測(cè)點(diǎn)的觀測(cè)資料不全，本次試驗(yàn)只選用兩個(gè)代表性水位觀測(cè)孔SK10及SK7進(jìn)行水位擬合，擬合曲線見(jiàn)圖3、圖4。

圖3 SK10孔水位擬合曲線圖

圖4 SK7孔水位擬合曲線圖

由圖3、4可知，本次識(shí)別各時(shí)段觀測(cè)孔水位計(jì)算值與實(shí)測(cè)值的擬合誤差符合模型誤差要求，模擬的結(jié)果表明系統(tǒng)穩(wěn)定性較好，模擬的精度和效果是理想的,建立的研究區(qū)水文地質(zhì)模型是比較符合實(shí)際的。

3.4.2 模型的驗(yàn)證

為進(jìn)一步驗(yàn)證識(shí)別后的模型和水文地質(zhì)參數(shù)的可靠性，對(duì)模型進(jìn)行驗(yàn)證，用校正后的模型及參數(shù)組合計(jì)算出各觀測(cè)點(diǎn)的水位，將計(jì)算結(jié)果與實(shí)測(cè)水位比較。根據(jù)計(jì)算區(qū)地下水位觀測(cè)資料的實(shí)際情況，選擇2010年10月5日至2010年10月25日期間的水位觀測(cè)資料對(duì)模型進(jìn)行驗(yàn)證。模型的結(jié)果見(jiàn)圖5、6，可知本次驗(yàn)證各時(shí)段水位擬合誤差符合驗(yàn)證要求。

綜上，模擬結(jié)果表明通過(guò)水文地質(zhì)條件的概化、邊界條件的確定、源匯項(xiàng)的處理、水文地質(zhì)參數(shù)的選取以及參數(shù)調(diào)整后所建模型能夠較好地反映研究區(qū)的水文地質(zhì)特征，且誤差均在工程允許范圍之內(nèi)，所以該模型可以用于地下水位動(dòng)態(tài)預(yù)報(bào)。

圖5 SK10計(jì)算結(jié)果與實(shí)測(cè)值對(duì)比圖

圖6 SK7計(jì)算結(jié)果與實(shí)測(cè)值對(duì)比圖

4 模型結(jié)果

本次模擬計(jì)算將研究區(qū)作為一個(gè)水文地質(zhì)系統(tǒng)，計(jì)算的目的層為淺層地下水，包含二疊系下石盒子組砂巖裂隙水及第四系松散層孔隙水，概化為非均質(zhì)各向異性。整個(gè)地下水流態(tài)概化為平面非穩(wěn)定流，四周均按通用水頭邊界處理，系統(tǒng)的底部為隔水邊界。通過(guò)模型運(yùn)轉(zhuǎn)，按水庫(kù)高水位和低水位，能夠預(yù)測(cè)水庫(kù)滲漏對(duì)副斜井涌水的影響。

在水庫(kù)不同蓄水水位條件下，研究區(qū)的側(cè)向補(bǔ)給量隨著水庫(kù)蓄水水位的升高和降低而增大和減少。目前水庫(kù)水位為1 245.1 m，屬一般高水位，據(jù)調(diào)查，水庫(kù)平時(shí)蓄水高度也在此高程附近，前面模型建立和驗(yàn)證時(shí)即按此作為初始水位。而水庫(kù)庫(kù)底高程為1 240 m左右，水庫(kù)設(shè)計(jì)高水位為1 250 m，故本次工作模擬預(yù)測(cè)水庫(kù)蓄水在1 240 m和1 250 m時(shí)對(duì)研究區(qū)地下水水位的影響。通過(guò)模型運(yùn)行，得出研究區(qū)滲流場(chǎng)，其流場(chǎng)模擬結(jié)果如圖7、圖8所示，具體運(yùn)算過(guò)程本文不做詳細(xì)敘述。

通過(guò)數(shù)值模擬可以看出，水庫(kù)滲漏對(duì)副斜井涌水有一定的影響，但其影響程度有限。目前庫(kù)水位為1245m，副斜井涌水量為58.083 m3/h；當(dāng)庫(kù)水位為1 240 m時(shí)，基本上為空庫(kù)，這時(shí)不會(huì)產(chǎn)生滲漏，副斜井涌水量為49.125 m3/h；當(dāng)庫(kù)水位為1 250 m時(shí)，副斜井涌水量增加為66.25m3/h，水庫(kù)滲漏對(duì)副斜井產(chǎn)生了影響。

圖7 水庫(kù)蓄水在1240m時(shí)滲流場(chǎng)模擬結(jié)果

圖8 水庫(kù)蓄水在1250m時(shí)滲流場(chǎng)模擬結(jié)果

［1］薛禹群.地下水動(dòng)力學(xué)［M］.北京:地質(zhì)出版社,1997.

［2］武強(qiáng),董東林等.水資源評(píng)價(jià)的可視化專(zhuān)業(yè)軟件(Visual Modflow)與應(yīng)用潛力[J].水文地質(zhì)工程地質(zhì),1999,26(5):21-23.

［3］宋西平,寇廣潮.地下水?dāng)?shù)值模擬面參數(shù)和水文參數(shù)研究的重要性［J］.地下水,2001,23(2):63-64.

［4］魏加華,陳良程,張遠(yuǎn)東等.地下水?dāng)?shù)值模型三維可視化研究［J］.煤田地質(zhì)與勘探,2003,31(4):33-36.

［5］鄭紅梅,劉明柱.Visual modflow在天津市地下水?dāng)?shù)值模擬中的應(yīng)用［J］.華北水利水電學(xué)院學(xué)報(bào),2007,28(2):8-1l.

［6］陳崇希.“防止模擬失真，提高仿真性”是數(shù)值模擬的核心［J］.水文地質(zhì)工程地,2003,47(2):1-5.

［7］蔣玲.和縣龍?zhí)裂罔F礦礦坑涌水量預(yù)測(cè)研究［D］.合肥:合肥工業(yè)大學(xué),2009.