煤礦礦井水水質(zhì)形成及演化的水動力場驅(qū)動作用及數(shù)學模型構(gòu)建

孫亞軍 ， 趙先鳴 ， 徐智敏 ， 張 莉 ， 陳 歌 ， 馮 琳 ， 李 鑫 ， 陳天賜 ， 袁慧卿 ，劉 琪 ， 郭 娟 ， 張 健 ， 劉嘉欣 ， 熊小鋒

(1.中國礦業(yè)大學 資源與地球科學學院, 江蘇 徐州 221116；2.礦山水害防治技術(shù)基礎(chǔ)研究國家級專業(yè)中心實驗室, 江蘇 徐州 221116)

我國煤炭資源的開采方式以井工開采為主[1]，由于采掘活動對頂?shù)装搴畬拥挠绊懀ǔa(chǎn)生大量的礦井水[2-3]。同時，煤礦開采過程還必然引起區(qū)域地下水動力場的改變，這也是礦井水水質(zhì)演化最主要的驅(qū)動力[4]。煤礦開采過程中，一方面，受頂?shù)装宀蓜悠茐膸У男纬蒣5-6]、控水構(gòu)造活化失穩(wěn)[7]和封閉不良鉆孔[8]、采空區(qū)積水[9]等影響，使礦區(qū)地下水動力場發(fā)生重大變化并進一步影響礦井水的水量及流速，這直接決定了參與水–巖、水–水混合作用等反應(yīng)的水量及反應(yīng)的時間，從而影響水質(zhì)的形成與演化；另一方面，原始地層在人為擾動以后，其溫度[10]、壓力[11]、介質(zhì)條件[12-14]、氧化還原[15]等理化條件也會發(fā)生明顯的改變，由此引起礦井水中某些特征組分及微生物群落結(jié)構(gòu)發(fā)生變化，打破了原有地下水中的物理–化學–生物反應(yīng)平衡，進而進一步影響礦井水水質(zhì)[16]。因此，礦井水水質(zhì)的形成與演化是由煤礦區(qū)水動力場驅(qū)動，水化學場、微生物場[17]等多場協(xié)同作用的結(jié)果。目前，既能體現(xiàn)煤礦采空區(qū)特殊水動力場又能夠表現(xiàn)多場作用的研究鮮有報道，并且煤礦采空區(qū)水動力場的特殊性及多場作用的復(fù)雜性也進一步限制了相關(guān)理論研究的推進。

筆者團隊在國家重點研發(fā)項目的支持下，在大量分析我國煤礦水質(zhì)演化的水動力場、水化學場、微生物場多場資料的基礎(chǔ)上，提出了礦井水水質(zhì)演化的垂向“三帶”概念[4]：第1 帶的水質(zhì)演化受開采活動影響最小，主要與地下水自然流動與蒸發(fā)濃縮作用相關(guān)；第2 帶多為煤礦常規(guī)生產(chǎn)階段產(chǎn)生的礦井水，其水質(zhì)形成與演化與不同含水層間的水–水混合作用以及水巖作用相關(guān)。根據(jù)《煤礦防治水細則》有關(guān)要求，生產(chǎn)礦井排水能力應(yīng)滿足能在20 h 內(nèi)排出礦井24 h 的正常涌水量[18]，故正常生產(chǎn)礦井產(chǎn)生的礦井涌水一般會在24 h 內(nèi)由井下排出，因此影響第2 帶礦井水水質(zhì)的水–水混合作用以及水–巖作用進行的時間相對較短，使得此水質(zhì)帶有較為明顯的原含水層水質(zhì)特征；在第3 帶，此水質(zhì)帶是人類活動影響最為頻繁的區(qū)域，在礦井開采時期，生產(chǎn)工作面和排水巷道中的地下水受礦井長期排水影響，水動力條件較好，但在已封閉的采空區(qū)中礦井水徑流緩慢。來自各充水含水層的地下水在經(jīng)歷長時間的混合作用、水–巖(煤)作用、微生物作用等耦合后形成感觀較差、微生物群落較為豐富的礦井污水[19]，即使在礦井閉坑以后，其水質(zhì)形成和演化仍可在區(qū)域水動力場的影響下長期進行。因此，第3 水質(zhì)帶內(nèi)的采空區(qū)積水是煤礦區(qū)礦井水污染研究的重點區(qū)域。目前針對煤礦區(qū)礦井水的水質(zhì)形成及演化機理的研究主要是通過野外取樣觀測、室內(nèi)批試驗以及建立相應(yīng)的數(shù)學、數(shù)值模型的方法研究礦井水的理化性質(zhì)、特征組分遷移擴散行為以及水–巖作用及水–水作用強度，但對于采空區(qū)積水的水質(zhì)形成過程中水動力場演化機理研究以及水動力場驅(qū)動下的多場作用機制研究相對匱乏。

在上述背景下，筆者以我國西北鄂爾多斯煤田某煤礦為地質(zhì)原型，分別對煤礦采空區(qū)積水水質(zhì)演化的水動力場作用以及水動力場驅(qū)動下的水化學場作用及微生物場作用進行綜述(由于所選地質(zhì)原型的采空區(qū)埋深超過500 m，溫度變化主要受地溫梯度影響，實測礦井涌水和主要充水含水層的水溫25～28 ℃，變化幅度和差異較小，因此可忽略溫度場的影響)；重點研究其水動力場的作用機制及數(shù)學模型，并開展模擬試驗對數(shù)學模型進行驗證；最后構(gòu)建水動力場主導(dǎo)下多場耦合數(shù)學模型，以期為煤炭安全綠色開采與水生態(tài)環(huán)境保護提供理論支撐。

1 采空區(qū)積水形成的水動力特征

當工作面回采結(jié)束后，一般會在一定時間內(nèi)撤出相關(guān)設(shè)備及材料，并進行永久性封閉，形成采空區(qū)。封閉后停止或減弱疏排水必然導(dǎo)致采空區(qū)內(nèi)的地下水位緩慢回升，形成采空區(qū)積水。此時，井下采空區(qū)成為礦井涌水的匯集區(qū)域，同時也是地下水發(fā)生物理、化學及生物等多場作用(圖1)的場所。因此，采空區(qū)積水的水質(zhì)形成及演化是以水動力場主導(dǎo)，水化學場、微生物場等多場協(xié)同作用的復(fù)雜過程，具有高度的非均質(zhì)性和時空變異性，然而相關(guān)針對性的研究較少，揭示采空區(qū)水質(zhì)演化過程中水動力場驅(qū)動下的多場作用對于煤礦區(qū)地下水污染防控具有重要意義。

圖1 封閉采空區(qū)積水水質(zhì)形成與演化的多場作用Fig.1 Multi field coupling effect on the formation and evolution of goaf water quality

1.1 采空區(qū)積水的水動力條件演化

一般情況下，采空區(qū)封閉后，人為對煤層、地層的擾動停滯，在區(qū)域水動力場作用下，采空區(qū)逐漸被原生充水含水層水源充填，采空區(qū)積水會與采空區(qū)內(nèi)的煤、巖等接觸并產(chǎn)生一系列的水–水[20]、水–巖(煤)[21]以及微生物[22]等作用，且井下水動力條件隨著采空區(qū)充水過程逐漸減弱，短期內(nèi)采空區(qū)積水水質(zhì)仍趨于劣化。采空區(qū)充滿后，在水動力場作用下，采空區(qū)積水會隨著區(qū)域地下水滲流方向運移，若遇到控水構(gòu)造或其他導(dǎo)水裂隙還會與其他含水層產(chǎn)生水質(zhì)混合，可能造成含水層間的串層污染。但在采空區(qū)封閉后的長期演化過程中，由于水動力條件減弱和采空區(qū)煤巖對某些水化學組分的吸附作用[23]，使得采空區(qū)積水水質(zhì)得到改善。

對于研究區(qū)鄂爾多斯礦區(qū)某煤礦而言，主采煤層為侏羅系延安組3-1煤，其主要充水水源是頂板砂巖裂隙水(由于安定組砂巖含水層富水性較弱，主要充水水源為直羅組砂巖含水層)[24]。其對于研究區(qū)這種單一頂板充水含水層的封閉采空區(qū)，在不考慮人工抽排水、間接充水含水層以及大氣降水的影響下，水動力場流動路線可以簡化為2 個階段，即采空區(qū)水位回升階段和蓄滿采空區(qū)–含水層水量交換階段(后簡稱水位回升階段和蓄滿水量交換階段)。如圖1(a)所示，水位回升階段適用于回采結(jié)束封閉采空區(qū)的初期，此階段由于采動裂隙帶及垮落帶的存在，直接充水含水層中的地下水流入采空區(qū)，導(dǎo)致采空區(qū)內(nèi)水位升高。當采空區(qū)內(nèi)水位升高至直接充水含水層，采空區(qū)與直接充水含水層形成一個含水單元，此時采空區(qū)內(nèi)的地下水流場與直接充水含水層的流場產(chǎn)生水量交換，如圖1(b)所示，這一階段直接充水含水層水從上游進入采空區(qū)后分為2 個部分，一部分沿著原含水層流向不變，另一部分則沿著導(dǎo)水裂隙帶進入垮落帶以及廢棄工作面然后再經(jīng)過導(dǎo)水裂隙帶回到原含水層中，因此，這一階段水量交換強度與流動路徑上的滲透系數(shù)和直接充水含水層在采掘影響范圍內(nèi)的上下游水頭相關(guān)。

1.2 水動力驅(qū)動下的水化學及微生物作用

1.2.1 水化學作用

在煤礦采掘過程中，采出空間頂板的巖層由于失去支撐而垮落，使得底板巖層的連續(xù)性遭到破壞，頂?shù)装鍘r層裂隙成為導(dǎo)水通道，因此，采空區(qū)積水通常由采空區(qū)多層頂、底板含水層混合組成，其水質(zhì)構(gòu)成十分復(fù)雜[25]。而對于頂板單一充水含水層的煤礦區(qū)而言，在采空區(qū)封閉初期階段，在圖1(a)中，在煤礦采空區(qū)形成的初始階段水動力循環(huán)交替作用較強，采空區(qū)內(nèi)水位較低且氧氣充足，常呈現(xiàn)出氧化環(huán)境，此時采空區(qū)積水的水質(zhì)演化除了與常規(guī)的溶解、沉淀、陽離子交換作用、(解)吸附作用相關(guān)還與氧氣參與下的多種氧化作用相關(guān)。在這一階段，原始含水層的地下水通過導(dǎo)水裂隙進入采空區(qū)形成積水，與采掘空間中無機介質(zhì)[26-27]，如巖(煤)層、采掘器械、支護設(shè)備等，以及有機介質(zhì)[28-29]，如腐植酸(HA)，各種機械使用的乳化液、潤滑油等發(fā)生化學反應(yīng)，導(dǎo)致某些特征離子發(fā)生氧化還原、交替吸附或溶解沉淀。而隨著采空區(qū)內(nèi)的水位不斷回升，逐漸轉(zhuǎn)變?yōu)槿毖醐h(huán)境，環(huán)境的變化控制著復(fù)雜的生物地球化學反應(yīng)的類型及作用程度，從而影響采空區(qū)的地下水化學組分演化過程，如圖1(b)所示，此時水動力循環(huán)交替作用減弱，逐漸影響水化學平衡的反應(yīng)方向，導(dǎo)致采空區(qū)積水水質(zhì)發(fā)生變化。在水化學場的作用下，一些特征組分呈現(xiàn)出“短期增加，長期減少”的趨勢[30]。以硫酸鹽為例，由于氧氣含量的降低，采空區(qū)逐漸轉(zhuǎn)變?yōu)槿毖醐h(huán)境，水–介質(zhì)發(fā)生的水化學反應(yīng)從以黃鐵礦等氧化作用為主轉(zhuǎn)變?yōu)橐粤蛩猁}還原作用為主，水體內(nèi)硫酸鹽的增漲速度減緩甚至緩慢降低，這也可能導(dǎo)致水體的TDS 在達到平衡后緩慢降低[31]。

1.2.2 微生物作用

微生物場在自然界中的分布不同于水動力場和水化學場，一般情況下，系統(tǒng)內(nèi)的微生物群落結(jié)構(gòu)和豐度并不呈現(xiàn)出連續(xù)的梯度變化趨勢。即便在同一系統(tǒng)中相鄰的位置之間，微生物群落結(jié)構(gòu)和豐度也可能差異較大，而產(chǎn)生這種差異的影響因素是復(fù)雜且多元的，例如，在煤礦區(qū)地下水環(huán)境中，相比于原始含水層的寡營環(huán)境，在采空區(qū)積水為有機物與能量供應(yīng)相對充足的的富營養(yǎng)環(huán)境，微生物在這種條件下會快速生長繁殖，并在氧氣、pH 及巖層物源條件的篩選作用下，微生物的群落結(jié)構(gòu)和豐度會發(fā)生變化[32]。此時，原始含水層與采空區(qū)積水、采空區(qū)積水形成初期(有氧環(huán)境)與采空區(qū)積水形成后期(缺氧–厭氧環(huán)境”)中的微生物群落結(jié)構(gòu)與微生物濃度差異很大的微生物差異很大，而氧氣含量和微生物可利用的能源類物質(zhì)的濃度即為影響微生物群落組成與豐度的主要影響因素。

在采空區(qū)封閉的初期即水位回升階段，采空區(qū)內(nèi)部處于有氧環(huán)境，好氧型微生物如硝化細菌、硫桿菌和鐵氧化菌等比例增加[33-34]。硝化細菌在氧氣充足的條件下會將 NH-4轉(zhuǎn)化成 NO-3和 NO-2[35]。硫桿菌和鐵氧化菌等將含硫和鐵的礦物氧化[36-38]，硫、鐵氧化物水解于礦井水中形成 SO24-、Fe2+、Fe3+[39-40]。采空區(qū)封閉后逐漸被礦井水充滿過程中，由于人為擾動停止，水動力條件逐漸趨于穩(wěn)定，使得礦井下微生物豐度及多樣性增加和遷移擴散的趨勢停滯，其溶解氧(DO)、pH 值、氮源、碳源等物源環(huán)境發(fā)生改變，有氧的氧化環(huán)境逐漸轉(zhuǎn)變?yōu)槿毖趸騾捬醯倪€原環(huán)境[41]。在此環(huán)境下，微生物群落種類、活躍程度及數(shù)量發(fā)生變化，厭氧型反硝化菌、鐵還原菌、脫硫菌等比例增加[40,42-43]。在還原環(huán)境下，采空區(qū)積水水質(zhì)中的反硝化菌可將NO和 NO反硝化生成N2[44]。有機物厭氧降解菌可將采空區(qū)積水中硝基苯、多環(huán)芳烴(PAHs)、石油烴、氯苯等有機污染物降解為代謝物及CO2[45-46]。鐵還原菌可將Fe3+還原為Fe2+[22,47]。而脫硫菌可利用有機碳等碳源的生化降解過程中產(chǎn)生的ATP 和高能電子將采空區(qū)積水中的 SO還原為S2–，造成水中H+的消耗，使得采空區(qū)積水堿度增加，并形成氣態(tài)硫化物H2S 逸出或金屬硫化物FeS 沉淀，從而改善了采空區(qū)積水的水質(zhì)。

2 水動力驅(qū)動的關(guān)鍵參數(shù)計算模型

通常情況下，水動力場主要是通過水頭、水量、流速的改變來驅(qū)動水質(zhì)演化，其作用強度與直接充水含水層的富水性、介質(zhì)滲透性能及流體的運動性能相關(guān)。對于地下水而言，其自身的運動性能取決于黏滯系數(shù)和容重，2 者大小和溫度密切相關(guān)，但根據(jù)前文所述，研究區(qū)水文變化幅度較小(25～28 ℃)，對于地下水的黏滯系數(shù)與容重影響較小，2 者均可視為定值。因此，在深部礦井的采空區(qū)，影響其水動力條件的關(guān)鍵參數(shù)為采動裂隙發(fā)育特征及直接充水含水層富水性能共同控制下的涌水量和滲透系數(shù)，且根據(jù)煤層實際的開采方式，這些關(guān)鍵參數(shù)在采空區(qū)以垂向和開采推進方向的二維變化為主。

2.1 滲透系數(shù)

由于采空區(qū)裂隙帶與垮落帶的巖石破碎程度不同，2 者的滲透性存在明顯差異，并且在實際采空區(qū)中，各巖層的孔隙滲透性能要遠低于采動裂隙[48]，因此本文數(shù)學模型忽略采空區(qū)裂隙帶與垮落帶巖層中的孔隙。

針對巖石較為破碎的采空區(qū)垮落帶，多將其視為由粒徑為Da(一般取值為0.25～1.00 m)的巖石顆粒所組成的均勻介質(zhì)，根據(jù)Blake-Kozeny 公式[49](式(1))可計算采空區(qū)垮落帶中的滲透系數(shù):

式中,K為巖層滲透系數(shù)，m/s；n為巖層的空隙率(有效孔隙率和裂隙率之和)；Da為巖石顆粒的平均粒徑，m；γ為流體容重，N/m3，這里取10 kN/m3；μ為流體動力黏滯系數(shù)，N·s/m，這里取1 000 N·s/m。

由式(1)可知，在垮落帶破碎巖體的粒徑一定的情況下，其滲透系數(shù)與空隙率成正比。由于巖層垮落和破壞后產(chǎn)生體積碎脹，導(dǎo)致采空區(qū)部分儲水空間被碎脹體積所占據(jù)，破碎巖體的空隙率可以用碎脹系數(shù)表示(式(2))[50]，而采動產(chǎn)生的裂隙空間隨應(yīng)力恢復(fù)而逐漸壓密，使得采空區(qū)垮落巖體的碎脹系數(shù)處于動態(tài)變化中。結(jié)合文獻[51]給出的垮落帶碎脹系數(shù)的計算公式(式(3))，可得到垮落帶空隙率的計算公式(式(4))。

式中，n0為垮落帶巖層的空隙率；nρ為巖層原始有效孔隙率，指有效透水的孔隙體積與原巖體積之比，深埋基巖的有效孔隙率取值一般較小，可忽略；kρ為裂隙帶第i層上邊界的碎脹系數(shù)；kρ0為垮落帶碎脹系數(shù)；為垮落帶初始碎脹系數(shù)，根據(jù)文獻[52]，取經(jīng)驗值1.5；E0為初始切線模量，MPa，其取值與巖體的單軸抗壓強度及初始碎脹系數(shù)相關(guān)[51]；σ為垮落帶巖體承受載荷，與裂隙帶高度和煤層埋深相關(guān)[53]，MPa。

裂隙帶不同高度的空隙率、碎脹系數(shù)與距離開采煤層的間距的關(guān)系一般通過式(5)表示。為了進一步的分析，筆者將采空區(qū)裂隙帶的巖層分成x層，如圖2 所示，且每層巖石的平均碎脹系數(shù)為該巖層的下邊界碎脹系數(shù)，結(jié)合式(2)、(5)和Blake-Kozeny 公式[54](式(6))，進一步推導(dǎo)出采空區(qū)滲透系數(shù)表達式(式(7))，其中裂隙帶下邊界碎脹系數(shù)和衰減系數(shù)η在采動停止以后的變化均較小可視為定值，因此采空區(qū)裂隙帶滲透系數(shù)主要與開采煤層的間距相關(guān)。

圖2 水位回升階段采空區(qū)水動力場概念Fig.2 Hydrodynamic field of water level recovery stage

式中，K0為裂隙帶下邊界的垂向滲透系數(shù)，m/s；kρi為裂隙帶第i層上邊界的碎脹系數(shù)；kρ′為裂隙帶下邊界的碎脹系數(shù)，根據(jù)文獻[52]，取經(jīng)驗值1.15；η為垂直方向體積碎脹系數(shù)衰減系數(shù)，通過實測垮裂帶高度反算求得；h0為采空區(qū)裂隙帶下邊界、垮落帶上邊界的高度，m；hi為第i層巖層頂板距離采空區(qū)底板的高度，m；n′i為第i層巖層的空隙率。

將鄂爾多斯礦區(qū)某煤礦3103 工作面采空區(qū)的埋深、巖石力學性質(zhì)等參數(shù)的實測值及經(jīng)驗值代入式(7)中的計算模型中，可繪制采空區(qū)滲透系數(shù)的垂向變化曲線，如圖3 所示。

圖3 采空區(qū)滲透系數(shù)的垂向變化曲線Fig.3 Vertical change curve of permeability coefficient

由圖3 可知，采空區(qū)的滲透性能隨高度的增長呈指數(shù)級減弱，這是由于裂隙帶內(nèi)的巖體與垮落帶相比較為完整，在垮落帶上部有序的排列堆積。裂隙帶巖體與煤層距離越遠，受到的采動影響程度越小，其破碎程度也越低。

2.2 水位回升高度及交換水量

根據(jù)第1 節(jié)分析，單一頂板充水含水層的封閉采空區(qū)的水動力場流動路線可以簡化為水位回升階段和蓄滿水量交換階段。在第1 階段，隨著水位的回升，氧氣含量、巖層物源條件等變化引起的水–巖作用及微生物作用類型改變，從而引起水質(zhì)改變，因此在第1 階段水質(zhì)影響的主要控制參數(shù)是時間和水位高度；而在第2 階段，采空區(qū)積水中的氧氣被基本耗盡，物源條件也不再發(fā)生變化，且由于積水與頂部原生地下水發(fā)生水量的長期交換，水體內(nèi)部反應(yīng)由水–巖作用為主轉(zhuǎn)變?yōu)橐运C水作用為主，進而導(dǎo)致水質(zhì)變化。因此，第2 階段水質(zhì)變化的主要控制參數(shù)是流經(jīng)采空區(qū)的水量。

2.2.1 積水水位回升階段

通常，采空區(qū)水位回升階段的流量可以利用大井法[55]計算(式(8))，結(jié)合式(2)、(3)、(5)可得采空區(qū)水位恢復(fù)時間t與高度z的關(guān)系(式(9))。

式中，Qr1為第1 階段流入采空區(qū)的流量，m3/s；Ka為頂板含水層的滲透系數(shù)，m/d，通過式(7)求得；M為含水層(出水段)厚度，m，實測值為11.50～79.97 m，取35 m；S為降深，m，實測值為682.83 m；R0為含水層的引用補給半徑，m；r0為引用半徑，m；t為采空區(qū)水位恢復(fù)至高度h所需時間，d；hx為采空區(qū)高度，m；A為采空區(qū)面積，m2；tz為采空區(qū)水位恢復(fù)至裂隙帶頂板所需的總時間，s；tk為水位由煤層底板恢復(fù)至垮落帶頂板所需的時間，s；tl為水位由垮落帶頂板恢復(fù)至裂隙帶頂板所需的時間，s。

將參數(shù)代入式(8)、(9)后，計算得Qr1=926.5 m3/h，并繪制采空區(qū)水位高度隨時間的變化曲線，如圖4所示。

圖4 采空區(qū)水位高度隨時間的變化曲線Fig.4 Change curve of goaf water level with time

由圖4 可知，采空區(qū)垮落帶內(nèi)的水位回升高度隨時間呈線性變化，這是由于計算模型中假定垮落帶有相同粒徑巖石所組成的均勻理想介質(zhì)，其空隙率為常數(shù)，且由于采空區(qū)儲水量受巖石空隙率控制，導(dǎo)致空隙率較低的裂隙帶內(nèi)水位回升速度明顯快于垮落帶。

2.2.2 蓄滿水量交換階段

采空區(qū)蓄滿后，由含水層流經(jīng)采空區(qū)的水滲流路徑為含水層上游→裂隙帶→垮落帶→裂隙帶→含水層下游?，F(xiàn)將流經(jīng)采空區(qū)的水量按流動路徑的通道性質(zhì)進一步分解為3 段，如圖5 所示。即采空區(qū)上游沿采動裂隙帶進入采空區(qū)的水量(QR)、沿著垮落帶通道從采空區(qū)上游流向采空區(qū)下游的水量和采空區(qū)下游沿采動導(dǎo)水裂隙流出采空區(qū)進入下游含水層的水量(QC)。則此時采空區(qū)水動力場可以概化為兩端為裂隙帶通道，中間為垮落帶通道的滲流模型。其中由于垮落帶空隙率遠大于裂隙帶為強過水通道，地下水在第2 段滲流路徑產(chǎn)生的水頭損失可以忽略不計，因此，地下水流經(jīng)采空區(qū)的滲流路徑可以進一步等效為地下水穿透弱—強—弱3 個滲透能力不同巖層的滲流，整體滲流路徑的兩端均受到滲透能力較弱裂隙帶的控制。裂隙帶巖層垂向滲透能力強弱取決于頂板上方滲透系數(shù)最小的分層，也就是說上部含水層的水通過頂板覆巖滲流到采空區(qū)的垂向透水能力取決于采動破壞后裂隙巖體中隔水性能最強的巖層。忽略各巖層自身的孔隙率，采空區(qū)各巖層滲透性能主要與開采煤層的間距及上覆載荷相關(guān)。因此，對于巖石力學性能基本相同的采空區(qū)巖層而言，其垂向的滲透性能主要受采空區(qū)頂板滲透系數(shù)的影響。

圖5 采空區(qū)蓄滿階段水動力場概念Fig.5 Hydrodynamic field of goaf filling up stage

水量交換階段的滲流運動在含水層上下游水力坡度不大情況下可以近似的用達西定律描述。因此，結(jié)合地下水采空區(qū)關(guān)閉后頂部充水含水層上游地下水沿采動導(dǎo)水裂隙進入采空區(qū)的水量QR和采空區(qū)水沿采動導(dǎo)水裂隙流向下游含水層的水量QC分別為

式中，QR、QC為水量交換階段流入(出)采空區(qū)的流量，m3/s；Kx為裂隙帶最頂部的滲透系數(shù)，m/s；W為受采礦影響的含水層平面寬度，m；L為含水層滲流路徑長度，m；H1、H3為采空區(qū)上(下)游含水層水頭，m；H2為采空區(qū)底板水頭，m。

在采空區(qū)水頭穩(wěn)定的情況下，含水層上游沿采動導(dǎo)水裂隙進入采空區(qū)的水量等于沿采動導(dǎo)水裂隙滲流入下游含水層的水量和上游流經(jīng)采空區(qū)然后再進入含水層的水量[11]，因此這一階段流經(jīng)采空區(qū)的水量Qr2為

由于鄂爾多斯礦區(qū)某煤礦頂板直羅組砂巖含水層水位實測標高1 274.65～1 296.80 m，表明采空區(qū)所處地層較為平坦，采空區(qū)上下游水力坡度較緩，故取采空區(qū)上下游含水層平均水頭差ΔH=15 m，將其代入式(12)，計算得Qr2約為128.25 m3/h。

2.3 模型驗證

由于采空區(qū)封閉以后可測性較差，難以對采空區(qū)裂隙發(fā)育及流場演化進行長期觀測，因此，本研究利用相似材料模型試驗和室內(nèi)三維箱式模擬試驗相結(jié)合的方法對采空區(qū)裂隙發(fā)育情況及封閉以后水位回升過程進行模擬，并利用試驗?zāi)M的測試結(jié)果與第2節(jié)中的計算模型進行對比，以驗證計算模型的可靠性。

2.3.1 模擬試驗

(1)相似材料模擬試驗。

筆者以鄂爾多斯礦區(qū)某煤礦某3-1工作面為地質(zhì)原型，對于3-1煤工作面來說其主要的充水水源是頂板砂巖裂隙水，根據(jù)勘探資料顯示，該工作面頂板巖性為粉砂巖、細砂巖、砂質(zhì)泥巖、2 煤、中砂巖和泥巖。通過相似材料模型試驗對該工作面采動條件下覆巖裂隙率的分布特征進行研究，進而分析其滲透性的垂向變化規(guī)律。以幾何相似比CL=100:1、容重相似比Cγ=1.67:1建立相似材料模型。圖6(a)為模型采高4 cm、推進距離為150 cm 時的采空區(qū)形態(tài)，此時采空區(qū)垮落帶高度約為16.5 cm，裂隙帶高度為71～75 cm。

圖6 相似材料模型試驗中采空區(qū)裂隙發(fā)育形態(tài)及裂隙二值化結(jié)果Fig.6 Goaf and fissures formed and binarization of fissures in similar material test

本文相似材料模擬實驗過程中采用單反相機記錄采空區(qū)裂隙巖體發(fā)育過程，記錄圖像含有大量的色彩信息且含有多種噪聲污染(常表現(xiàn)為孤立像素點或像素塊)，直接利用會影響裂隙特征提取，對后期裂隙長度、寬度、空隙率等參數(shù)的計算產(chǎn)生誤差。為進一步定量分析采空區(qū)垂向裂隙的變化規(guī)律，也為闡明垮落帶和裂隙帶不同高度滲透性演變規(guī)律提供基礎(chǔ)，對相似材料模擬試驗的照片進一步進行以下圖像處理流程即可實現(xiàn)試驗記錄圖像的裂隙特征的快速提取。

① 對于試驗記錄圖像，首先可根據(jù)式(13)將圖像3 個通道(RGB)轉(zhuǎn)換成1 個通道，即灰度化處理。

② 將灰度處理后圖像的每個像素替換為該像素周圍像素的均值，這樣就可以平滑并替代那些強度變化明顯的區(qū)域，以實現(xiàn)圖像的濾波降噪。

③ 將前處理好的灰度去噪圖像通過圖像分割進行二值化處理，使其變?yōu)橹缓? 種顏色的圖像，圖像中黑、白2 種顏色分別代表裂隙和巖體(其中由于2 號煤層顏色與裂隙相近，該部分為人工調(diào)整結(jié)果)。

處理結(jié)果如圖6(b)所示，可準確統(tǒng)計獲取采空區(qū)不同高度空隙和完整巖層的像素面積之比，進而計算獲取不同高度的采后空隙率，見表1。

表1 采空區(qū)不同高度空隙率統(tǒng)計計算結(jié)果Table 1 Statistical calculation results of fracture rate of goaf at different heights

由表1 分析可知：

① 在煤層頂板向上約16.5 m 范圍內(nèi)屬垮落帶，垮落體較為破碎，此區(qū)域空隙率約為22.25%，空隙率明顯高于裂隙帶和原始地層；

② 垮落帶上部為裂隙帶，此范圍內(nèi)空隙率相較于垮落帶明顯變小，且裂隙帶范圍內(nèi)自下往上空隙率由0.5%逐漸增大至16.9%，同樣顯著高于原始地層。

(2)室內(nèi)三維箱式模擬試驗。

室內(nèi)三維箱式模擬試驗以采空區(qū)的垮落帶和裂隙帶為地質(zhì)原型，此區(qū)域整體均為采空區(qū)填充物及部分上覆巖層，由上到下分別為較完整巖層、破碎巖層、塊狀矸石、煤泥等，各層厚度、巖性和空隙率根據(jù)實際前述實測及表1 中的相似材料模擬試驗的統(tǒng)計、計算結(jié)果確定。如圖7 所示，模型試驗裝置(模型箱)內(nèi)徑長度150 cm×寬度75 cm×高度75 cm，模型幾何相似比為CL=100:1，模擬實際區(qū)域高72.5 m×寬75 m×長150 m長方體空間，該空間高度與實際采空區(qū)相同，底面積為實際采空區(qū)(面積約為8×105m2)的一部分。由于模型采用真實地層巖樣和含水層水樣搭建，模型容重相似比Cγ=1:1 。模型時間比Ct=:1=10:1，根據(jù)實測實際采空區(qū)的穩(wěn)定涌水量約為900 m3/h，根據(jù)相似理論計算得模型注水流量(模擬穩(wěn)定涌水量)約為2.1 mL/min。

圖7 室內(nèi)三維箱式模擬試驗設(shè)計Fig.7 Schematic diagram of inflow and goaf in the simulative box

模型運行后，水位回升高度與時間的關(guān)系如圖8所示。圖8 中裂隙帶的水位回升速度明顯快于垮落帶，且隨著裂隙帶高度的升高，水位回升速度呈現(xiàn)增快的趨勢，這是由于空隙率隨巖層高度逐漸減小，導(dǎo)致水位回升的過程中儲水空間不斷減小，在涌水量穩(wěn)定的情況下，水位回升速度不斷加快。

圖8 模型試驗中采空區(qū)水位高度隨時間的變化Fig.8 Curve of goaf water level with time in the simulative box

2.3.2 試驗結(jié)果與理論計算對比

計算結(jié)果與模擬試驗實測結(jié)果(已進行相似比換算)進行對比，如圖9 所示。圖9(a)中，滲透系數(shù)K的理論計算和試驗實測誤差約為13.2%，表明理論模型計算與模擬試驗的滲透系數(shù)基本為同一數(shù)量級，均是隨著距離煤層距離的增加急劇下降，尤其裂隙帶部分滲透系數(shù)呈數(shù)量級倍數(shù)的減??；圖9(a)中，采空區(qū)水位回升高度隨時間的變化曲線的理論計算和試驗實測誤差約為22.5%，理論計算的裂隙帶水位回升速度明顯快于垮落帶，呈指數(shù)級增長，相比于試驗實測，主要表現(xiàn)為垮落帶水位回升更為緩慢，而裂隙帶反之。總體蓄滿時間誤差約為3.4%，計算結(jié)果與實測值基本吻合，說明本研究提出的采空區(qū)水動力場關(guān)鍵參數(shù)的理論計算模型具備一定的可靠性。

圖9 模型計算結(jié)果與模擬試驗觀測值對比Fig.9 Model calculation results and the simulation experiment observation

誤差產(chǎn)生的原因是主要是由于理論模型計算的垮落帶和裂隙帶巖層的初始碎脹系數(shù)為經(jīng)驗值而非實測值，這導(dǎo)致理論模型的垮落帶孔(裂)隙率計算值偏大而裂隙帶孔(裂)隙率計算值偏小，同時理論模型未考慮各巖層的巖石力學性質(zhì)的變化以及巖層本身的孔隙度，這些因素均會導(dǎo)致誤差的出現(xiàn)。

3 水動力場的驅(qū)動作用及數(shù)學模型

3.1 水質(zhì)演化的多場作用關(guān)系

對于水化學場、微生物場與水動力場之間的耦合，地下水流攜帶微生物和特征組分遷移，地下水流速越大，流體動壓促使微生物及特征組分遷移越快；而微生物通過降解代謝、生長繁殖或消亡，特征組分通過溶解沉淀等作用引起滲流介質(zhì)的孔隙度、孔喉大小與運移路徑長度的變化，從而改變介質(zhì)的滲透性能，間接控制地下水流速，但是相比于開采擾動形成的采動裂隙，這種通過生化作用對地層滲透性能的影響較小，因此在多場耦合的模型建立過程中可以忽略水化學場、微生物場對水動力場的反向耦合作用。

3.2 水動力場的概念模型

根據(jù)前文分析，采空區(qū)內(nèi)的水動力場可以近似的看作頂部存在含水層入滲補給的非飽和無壓滲流模型，可用Richards 方程表示，由含水率的計算公式(式(14))可知，此時的滲透系數(shù)不在是常數(shù)，而是與空隙介質(zhì)的含水率相關(guān)。對于相似材料試驗?zāi)M的采空區(qū)而言，其整體的平均空隙率為9.86%，且研究區(qū)巖層相對于土壤或砂礫其致密程度較高，可以忽略其巖層本身的孔隙率。采空區(qū)裂隙巖層中粗糙裂隙面的間隙形成的開度分布與多孔介質(zhì)固體顆粒間形成的喉頸孔隙分布的持水機制有一定的相似性。近似將采空區(qū)的煤(巖)體看作具有不可壓縮固體骨架的多孔介質(zhì)體，故可借助多孔介質(zhì)滲流力學理論模型求解飽和–非飽和采空區(qū)介質(zhì)中的含水率、壓力水頭等的時空分布情況。

式中， θ為采空區(qū)巖層的含水率；Sw為飽和度。

3.3 水動力場的數(shù)學模型構(gòu)建

根據(jù)上述假設(shè)，采空區(qū)積水的過程即可以概化為頂部存在入滲補給的非飽和多孔介質(zhì)滲流模型，以多孔含水介質(zhì)的空隙表示采動形成的裂隙，并且忽略介質(zhì)和流體的各向異性。因此，采空區(qū)水動力場模型的數(shù)學方程式為

式中，u為地下水滲流速度場，m/s；ρw為水的密度，kg/m3，這里取1 000 kg/m3；Qw為源匯項，m3/s。

忽略采空區(qū)積水的密度變化，以含水率表示積水連續(xù)運動方程為

將滲流速度場的公式(式(17))代入式(16)可得到采空區(qū)積水過程的飽和–非飽和數(shù)學方程，如式(18)所示：

式中，K(Sw)為飽和–非飽和滲透系數(shù)，m2/s，當飽和度較低時，采空區(qū)多孔介質(zhì)中較大的孔隙沒有充滿水，此時滲透主要在較小的孔隙中進行，導(dǎo)致滲流作用的過水斷面減小，滲流路徑增長使相對值滲透率大幅減小，隨著含水層的不斷補給，采空區(qū)逐漸變?yōu)轱柡蜖顟B(tài)狀態(tài)，此時各巖層的滲透系數(shù)恢復(fù)為僅與空隙率相關(guān)的參數(shù)；p為壓強，Pa，在飽和帶p代表水的壓強，取正值，在非飽和帶中p代表毛管壓強，取負值；z為以采空區(qū)頂板為垂向的零點垂向坐標，m。

(1)模型初始條件。

為簡化起見，假設(shè)采空區(qū)同一水平剖面上初始空隙率是均一的，故滲透系數(shù)主要為垂向變化，因此有方程(式(7))，將其轉(zhuǎn)變?yōu)榭障堵食跏紬l件的方程形式(式(19))。假設(shè)采空區(qū)同一水平剖面上初始含水率也是均一的，并且整個采空區(qū)初始為疏干狀態(tài)。因此可認為除了頂部邊界外，采空區(qū)任一高程的初始含水率均為0(式(20))。

(2)模型邊界條件。

模型上邊界為含水層入滲補給邊界，且模型頂部滲透性較低，根據(jù)“大井法”可知，此時的入滲流量和頂部邊界的滲透性相關(guān)。因此當采空區(qū)封閉以后，開采擾動停止其滲透性幾乎不變的情況下，可以將采空區(qū)模型的上邊界近似的概化為有穩(wěn)定流量的入滲邊界，即

式中，Q0為采空區(qū)頂邊界的入滲流量，m3/d。

采空區(qū)頂部邊界與含水層底板的接觸面達到接近飽和的含水率，即Sw=1，結(jié)合式(13)有

再將n的表達式(式(19))代入式(6)～(16)即可得到采空區(qū)頂部的邊界條件：

隨著時間推移，采空區(qū)形成積水且水位不斷升高，底板為較為完整的巖層，可近似假定為不透水邊界，即下邊界的流量等于0 的邊界有

式中，hp為壓強水頭，m。

3.4 水動力場對水質(zhì)演化的驅(qū)動作用

在關(guān)于地下水中特征污染組分形成與演化的研究中，水動力場多是以單一變量的形式研究其對特征污染組分的遷移和空間分布的影響。在實際的煤礦采空區(qū)中，是通過包括水動力場、水化學場、微生物場在內(nèi)的多場共同作用對特征污染組分的形成、演化、遷移、分布產(chǎn)生影響的。

采空區(qū)封閉后，雖然人為擾動停止，但是此時的區(qū)域水動力場與原生水動力場已經(jīng)存在明顯差別，此時采空區(qū)成為了區(qū)域地下水匯水區(qū)域，導(dǎo)致地下水流速、物源條件、微生物群落與豐度的變化，從而引發(fā)了一系列水–水混合、水–巖(煤)作用以及微生物分解等作用，導(dǎo)致地下水質(zhì)發(fā)生改變。在采空區(qū)蓄滿后，在水動力場作用下，采空區(qū)積水會隨著區(qū)域地下水滲流方向運移，導(dǎo)致采空區(qū)水與原生地下水不斷發(fā)生水量交換，某些特征污染物被不斷稀釋，直至完全恢復(fù)至區(qū)域地下水的原始濃度。在水化學場控制方程(式(26))和微生物場的控制方程(式(27))中，其彌散項、對流項以及滲流源匯項均與水動力場密切相關(guān)，而生化源匯項則與氧氣、pH、礦物成分以及微生物群落組成與豐度等因素相關(guān)。這表明水質(zhì)形成與演化的多場作用是以水動力場為主導(dǎo)的，即水動力場對采空區(qū)積水的水質(zhì)演化起驅(qū)動作用。

式中,Cs為地下水中特征組分的質(zhì)量濃度，mg/L；?[θD?C] 為水動力彌散相，D為采空區(qū)積水中 SO的水動力彌散系數(shù)，m2/s；u·?C為對流相；vsCss為源匯項，vs為源/匯的滲流速度，Css為水動力源/匯中的特征組分質(zhì)量濃度(即主要充水含水層水中的特征組分濃度)，mg/L；f1為特征組分生化源匯項，包括吸附解吸項、礦物氧化還原相、微生物分解相等，表示由吸附解吸、礦物氧化還原以及微生物分解等作用引起的該點特征組分濃度的變化，相關(guān)作用的數(shù)學模型構(gòu)建將在后續(xù)論文中進一步闡述。

式中，B為地下水中的微生物量，kg/kg(生物/液體)；Bns為水相自由態(tài)微生物量，kg/kg(生物/液體)；DB為水中自由態(tài)微生物水動力彌散系數(shù)，m2/s；由于自由態(tài)微生物會趨向適宜的環(huán)境運動，因此微生物自身會產(chǎn)生抵抗或順應(yīng)水動力場方向的運動速度即為uB，m/s；u-uB為某點微生物在水動力場影響下的遷移速度，m/s；Bss為水動力源/匯中的微生物量，kg/kg(生物/液體)；A′為裂隙面積與巖層體積之比，m-1；Bs為附著在巖層裂隙表面的固相吸附態(tài)微生物量，kg/m2；f2為由于多種因素(包括生物代謝、產(chǎn)物抑制、競爭及非競爭抑制作用等)導(dǎo)致微生物的生長繁殖死亡而引起的水中及裂隙表面的生物量的變化項，相關(guān)作用的數(shù)學模型構(gòu)建將在后續(xù)論文中進一步闡述。

3.5 水動力場主導(dǎo)下的多場耦合本構(gòu)模型

根據(jù)上述數(shù)學方程，水動力場與水化學場、微生物場之間的耦合則分別通過指標(v，C)和(v，B)進行關(guān)聯(lián)。水動力場驅(qū)動下的多場耦合本構(gòu)模型如圖10所示，其耦合方程為

圖10 水質(zhì)形成及演化的多場耦合關(guān)系Fig.10 Multi-field coupling relationship of water quality formation and evolution

式中，Cs為地下水中特征組分的質(zhì)量濃度，mg/L；Cs0為特征組分的初始質(zhì)量濃度，mg/L。fa為地下水流速與生物量之間的量化方程；fb為水化場中特征組分濃度與生物量生長抑制之間的量化方程；fc為水化場中特征組分濃度在無菌條件下與地下水流速之間的量化方程；Ks為地下水流影響下的特征組分的微生物降解系數(shù)。

4 討 論

(1)采掘活動空間的頂板涌水是井工煤礦采空區(qū)積水的最主要來源，筆者討論了單一頂板充水含水層的采空區(qū)積水水質(zhì)形成及演化的水動力場驅(qū)動作用、關(guān)鍵參數(shù)計算及多場作用數(shù)學模型構(gòu)建的基礎(chǔ)問題。但在實際生產(chǎn)過程中，采空區(qū)的涌水水源還包括部分底板水，甚至地表水，同時工作面的傾角及控水構(gòu)造分布也各有差異，這些因素均會對其水動力場產(chǎn)生影響，因此，適用于不同礦區(qū)采空區(qū)的水動力場概念及數(shù)學模型的研究值得更進一步的深入探討。

(2)采空區(qū)水質(zhì)形成及演化的多場作用是以水動力場為主要驅(qū)動力，水化學場、微生物場協(xié)同參與的復(fù)雜過程，限于篇幅，多場耦合作用中的化學場、微生物場作用的數(shù)學模型構(gòu)建將在后續(xù)論文中陸續(xù)闡述。

(3)多場耦合作用數(shù)學模型建立過程中，由于水化學場、微生物場對水動力場的反向耦合作用較小，本文暫時予以忽略。在少數(shù)極端環(huán)境下，此反向耦合作用較為明顯時，需通過開展多種微觀、宏觀試驗，對水動力–水化學–微生物場多場參數(shù)進行精準的動態(tài)觀測，明確反向耦合作用的關(guān)鍵參數(shù)及驅(qū)動機制。

(4)現(xiàn)有針對多場耦合理論研究尚處于耦合機理研究以及數(shù)學模型構(gòu)建階段，其中大量的多場作用參數(shù)難以通過室內(nèi)、室外試驗確定。而現(xiàn)有的多場耦合模擬軟件對水動力場、水化學場、微生物場參與下的礦井水水質(zhì)形成與演化的模擬方案針對性不足，需進一步針對不同場景對各場參數(shù)進行精準調(diào)控，優(yōu)化模擬方案，以提高研發(fā)不同場景下多場耦合模型的效率。

5 結(jié) 論

(1)闡明了煤礦采空區(qū)積水形成的水動力場演化機制，根據(jù)水動力場流動路徑的差異將其分為水位回升階段和采空區(qū)蓄滿水量交換階段；在此基礎(chǔ)上，通過文獻綜述，探究了不同階段的水動力場驅(qū)動下水化學場及微生物場的作用機制，揭示了采空區(qū)水質(zhì)形成及演化的以水動力場為主要驅(qū)動力，水化學場、微生物場等多場協(xié)同作用的復(fù)雜過程。

(2)提出了采空區(qū)積水水質(zhì)形成與演化的水動力場水位回升階段和采空區(qū)蓄滿階段的概念模型，建立了基于頂板導(dǎo)水裂隙帶不同垂高空隙反演的分層滲透系數(shù)定量計算模型，為開采擾動后采空區(qū)積水水位回升及蓄滿后水量交換的計算奠定了理論基礎(chǔ)。

(3)通過相似材料模型試驗和室內(nèi)三維箱式模擬試驗對采空區(qū)裂隙發(fā)育情況及封閉以后水位回升過程進行了模擬，試驗觀測結(jié)果進一步揭示了采空區(qū)各層位空隙度具有垂向差異，并利用試驗的實測數(shù)據(jù)與模型計算結(jié)果進行對比，理論計算滲透系數(shù)K和采空區(qū)水位回升高度隨時間的變化與試驗實測誤差分別約為13.2%和22.5%，驗證了計算模型的可靠性。

(4)在煤礦采空區(qū)水動力場以及水動力場影響下的水化學場和微生物場的控制方程的基礎(chǔ)上，進一步明確了水動力場–水化學場–微生物場之間的關(guān)聯(lián)指標和耦合關(guān)系，構(gòu)建了以水動力場主導(dǎo)下的多場耦合的本構(gòu)模型，可為進一步探究水動力場耦合驅(qū)動下的煤礦礦井水水質(zhì)形成及演化規(guī)律的多場耦合模型構(gòu)建提供啟示。