煤礦礦井水水質形成及演化的水動力場驅動作用及數學模型構建

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

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

我國煤炭資源的開采方式以井工開采為主[1]，由于采掘活動對頂底板含水層的影響，通常會產生大量的礦井水[2-3]。同時，煤礦開采過程還必然引起區域地下水動力場的改變，這也是礦井水水質演化最主要的驅動力[4]。煤礦開采過程中，一方面，受頂底板采動破壞帶的形成[5-6]、控水構造活化失穩[7]和封閉不良鉆孔[8]、采空區積水[9]等影響，使礦區地下水動力場發生重大變化并進一步影響礦井水的水量及流速，這直接決定了參與水–巖、水–水混合作用等反應的水量及反應的時間，從而影響水質的形成與演化；另一方面，原始地層在人為擾動以后，其溫度[10]、壓力[11]、介質條件[12-14]、氧化還原[15]等理化條件也會發生明顯的改變，由此引起礦井水中某些特征組分及微生物群落結構發生變化，打破了原有地下水中的物理–化學–生物反應平衡，進而進一步影響礦井水水質[16]。因此，礦井水水質的形成與演化是由煤礦區水動力場驅動，水化學場、微生物場[17]等多場協同作用的結果。目前，既能體現煤礦采空區特殊水動力場又能夠表現多場作用的研究鮮有報道，并且煤礦采空區水動力場的特殊性及多場作用的復雜性也進一步限制了相關理論研究的推進。

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

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

1 采空區積水形成的水動力特征

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

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

1.1 采空區積水的水動力條件演化

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

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

1.2 水動力驅動下的水化學及微生物作用

1.2.1 水化學作用

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

1.2.2 微生物作用

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

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

2 水動力驅動的關鍵參數計算模型

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

2.1 滲透系數

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

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

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

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

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

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

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

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

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

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

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

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

根據第1 節分析，單一頂板充水含水層的封閉采空區的水動力場流動路線可以簡化為水位回升階段和蓄滿水量交換階段。在第1 階段，隨著水位的回升，氧氣含量、巖層物源條件等變化引起的水–巖作用及微生物作用類型改變，從而引起水質改變，因此在第1 階段水質影響的主要控制參數是時間和水位高度；而在第2 階段，采空區積水中的氧氣被基本耗盡，物源條件也不再發生變化，且由于積水與頂部原生地下水發生水量的長期交換，水體內部反應由水–巖作用為主轉變為以水–水作用為主，進而導致水質變化。因此，第2 階段水質變化的主要控制參數是流經采空區的水量。

2.2.1 積水水位回升階段

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

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

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

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

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

2.2.2 蓄滿水量交換階段

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

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

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

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

在采空區水頭穩定的情況下，含水層上游沿采動導水裂隙進入采空區的水量等于沿采動導水裂隙滲流入下游含水層的水量和上游流經采空區然后再進入含水層的水量[11]，因此這一階段流經采空區的水量Qr2為

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

2.3 模型驗證

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

2.3.1 模擬試驗

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

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

圖6 相似材料模型試驗中采空區裂隙發育形態及裂隙二值化結果Fig.6 Goaf and fissures formed and binarization of fissures in similar material test

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

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

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

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

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

表1 采空區不同高度空隙率統計計算結果Table 1 Statistical calculation results of fracture rate of goaf at different heights

由表1 分析可知：

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

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

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

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

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

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

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

2.3.2 試驗結果與理論計算對比

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

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

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

3 水動力場的驅動作用及數學模型

3.1 水質演化的多場作用關系

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

3.2 水動力場的概念模型

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

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

3.3 水動力場的數學模型構建

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

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

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

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

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

(1)模型初始條件。

為簡化起見，假設采空區同一水平剖面上初始空隙率是均一的，故滲透系數主要為垂向變化，因此有方程(式(7))，將其轉變為空隙率初始條件的方程形式(式(19))。假設采空區同一水平剖面上初始含水率也是均一的，并且整個采空區初始為疏干狀態。因此可認為除了頂部邊界外，采空區任一高程的初始含水率均為0(式(20))。

(2)模型邊界條件。

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

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

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

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

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

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

3.4 水動力場對水質演化的驅動作用

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

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

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

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

3.5 水動力場主導下的多場耦合本構模型

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

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

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

4 討 論

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

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

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

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

5 結 論

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

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

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

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