1. 數學模型建立與地下水數值模擬模型
(一)數學模型的建立
採用三維數學模型模擬研究區含水層系統地下水流動，數學方程如下:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
式中:kxx，kyy，kzz分別為沿x，y，z坐標軸方向的滲透系數(m/d);h為淺層地下水系統和中深層地下水系統水位標高(m);S為點(x，y，z)處的儲水率(1/m);W為源匯項(1/d)，是時間和空間的函數W=W(x，y，z，t);t為時間(d);n為二類邊界內法線向量;q(x，y，z，t)為流量在時間和空間上的變化函數;f2(x，y，z，t)為水位在時間和空間上的變化函數;D為淺層地下水系統和中深層地下水系統同時確定的滲流域;Γ1為滲流計算域D的已知水位邊界;Γ2為滲流計算域D的已知流量邊界。
(二)地下水數值模擬模型
為了更精確地評價研究區地下水資源，預測煤礦開發利用條件下的地下水系統狀態及其變化趨勢，在建立數學模型的基礎上，運用基於有限單元法的Feflow軟體建立了研究區的地下水流數值模擬模型，經識別與檢驗後，對研究區地下水滲流規律進行模擬，並對將來地下水動態進行分析評價。
1.模擬軟體選取及網格剖分
本次工作採用的軟體是由德國WASY公司開發的基於有限單元法的Feflow軟體。建立研究區地下水流數值模擬模型，首先要對模擬區進行三角形剖分。剖分時除了遵循一般的剖分原則(如三角形單元內角盡量不出現鈍角，相鄰單元間面積相差不應太大)，還應考慮如下實際情況:①充分考慮研究區的邊界、岩性分區界線、斷層等;②觀測孔、水源地盡量放在剖分單元的節點上;③在水力坡度變化較大及重點研究區，剖分時應適當加密。剖分後的模擬區共4473個結點，6996個單元格(圖5-5)。模擬的含水層包括煤14-K3(Ⅱ)含水層和奧灰含水層。研究區的三維結構見圖5-6。
圖5-5 研究區三角剖分圖
圖5-6 研究區的三維結構圖
由於計算區的邊界是斷層邊界，而非區域地層變化邊界，因此在邊界的處理上存在著一定誤差，通過對研究區及外圍地下水水量和水位多年觀測資料的處理和分析，並參考相關研究成果，初步確定模型計算的邊界。
(1)垂向邊界:計算模擬區的上部邊界為第四系承壓含水層，是位置不斷變化的水量交換邊界，有大氣降水、地表水等與潛水之間發生密切的水力聯系。計算模擬區的下部邊界為承壓含水層底板，其組成為滲透性極差的基岩(泥岩、頁岩等)，概化為隔水邊界。
(2)側向邊界:承壓含水層的地下水分水嶺，處理為零通量邊界;地下水運動可概化成空間三維流;地下水系統的輸入、輸出隨時空變化，故地下水為非穩定流;參數隨空間變化，體現了含水介質的非均質性，但沒有明顯的方向性，因此參數可概化成各向同性。綜上所述，研究區地下水流系統可概化成非均質、各向同性、空間三維結構、非穩定地下水流系統。
2.地下水系統模擬的有限單元數值解
有限單元法是由我國數學家馮康等人於年最先提出，最早主要應用於力學。於20世紀末將有限元引用到地下水流問題的求解。隨著計算機技術的迅速發展，有限元法已經成為解決復雜水文地質滲流問題的有效方法，有限單元是利用剖分插值把區域連續求解的微分方程離散成求解線性代數方程組，以近似數值解代替微分方程的精確解。根據所建立的數學方程組的不同，有限元法可以分為里茲有限元、均衡有限元和伽遼金有限元等。
里茲有限單元以變分原理和剖分插值為基礎。變分原理就是把描述地下水運動的偏微分方程的求解化為求某個泛函的極值問題;剖分插值則是把研究的滲流區從幾何上剖分為點、線、面體單元，假設節點的水頭值，然後根據實際情況採用某種形式的插值法按單元插值，由構造每個單元水頭表示式，最後形成整個單元幾何體的插值。這種方法就是從變分原理出發，利用整個單元集合體的插值把求某個泛函的極值問題化為一組多元線性代數方程組的求解問題，從而獲得所求滲流問題的解。均衡有限單元是從小均衡角度出發，對滲流區進行三角剖分，然後建立以任一節點為公共頂點的各三角形的重心與相應邊中點的連線所組成的區域的均衡方程，即可得到一個線性代數方程。對整個研究區來講，將得到一個線性方程組。均衡有限單元法雖然簡單直觀，但在數學上被認為是不嚴密的。伽遼金有限單元法是從剩餘加權法出發對連續性的微分方程進行離散，從而求其滿足水文地質約束條件的控制方程數值解的，其數學原理如下。
(1)剩餘加權法原理:剩餘加權法是求解微分方程近似解的一種有效方法，有下列數學模型:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
其解析解必須在區域對該方程即在上式中任一點都滿足微分方程和相應的約束條件。對於復雜問題，其解析解是不易解出的，應找到一個能滿足一定精度的解析解，該精度一般是事先設定的誤差范圍的解析解的近似解。
現用一試探函數:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
式中:φi為按照特定要求選定的已知函數;φ1，φ2，…，φn為線性獨立的基函數;αi為待定系數。
由於 為方程的近似解，故將 代入上式中，使該方程的右端不為0，從而產生一個余項記為R。即L(v)=R。
假若用某種方法選取與待定系數相對應的幾個權函數ω1，ω2，…，ωn，並要求權函數能使余項R的加權積分值等於0，於是就可以得到幾個方程構成的方程組，即
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
式中:Ω為積分區間，可以是線域、面域或體域中的任何一種。該式為剩餘加權積分式，通過求解該式可以確定幾個待定系數α1、α2，…，αn，從而可以求解數學模型的近似解。
(2)伽遼金法:在實際的應用中，權函數ωi的選取的方法很多，如果選取的權函數與基函數φi相同，那麼 就可以寫成 i=1，2，3，…，n，此時，稱該式為伽遼金方程。
伽遼金有限元可以解決地下水三維流問題，假設計算區域Ω為非均質各向異性含水介質，邊界由Γ=Γ1+Γ2組成，其中Γ1為含水系統的一類邊界，Γ2為第二類邊界。如果選取的空間坐標系統與含水層的主滲透系數方向一致，那麼三維地下水流的數學模型可用下述方程組:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
上式中:H(x，y，z)為含水層中地下水頭分布;kxx，kyy，kzz分別為含水層主滲透系數;t為時間變數;Ω為研究區滲流空間;Γ1、Γ2為第一、第二類研究區邊界;W為源、匯項的代數和;Ss為含水層貯水率;H0為初始流場水頭分布;H1為第一類邊界水頭分布;q為第二類邊界流量。
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
按照剩餘加權理論取試探函數: ,並令該式滿足方程組的定解條件。其中{Φi(x，y，z)，i=1，2，…，n}為構造的基函數組，{Ci(t)，i=1，2，…，n}是依賴於時間的系數。對於固定的t，這些系數可以由下列方程組求解:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
對於變數y和x有類似的表達式，即:
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
上方程寫成矩陣的形式為
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
式中矩陣［H］和［P］的元素分別為
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
列向量F的元素為
典型煤礦地下水運動及污染數值模擬:Feflow及Modflow應用
只要確定了基函數{Φi}，上述變數都可以算出，從而可以求解地下水流運動控制方程的近似解珘h。
在實際應用中，利用伽遼金有限元方法求解地下水運動方程的主要步驟包括:①將研究區在空間上剖分為有限個滲流單元，形成滲流單元網格;②根據剖分節點的坐標可以是絕對坐標，也可以是相對坐標計算滲流單元的平面面積、空間體積，以及介質常數等，建立各單元的導水矩陣;③以各節點水頭表達式上的流量，即等效流量，把各單元的已知流量按照平分原則移動到節點上;④針對各未知節點水頭建立連續方程，對方程求解。
以上各步驟均由Feflow自動完成。
3.模擬期的處理
本次研究以2005年1月流場作為地下水資源計算和評價的初始流場，2007年1月流場作為地下水資源計算和評價的擬合流場。模擬期為2005年1月～2007年1月兩個水文年，以一個月作為一個時間段，每個時間段內包括若干時間步長，時間步長為模型自動控制。
4.定解條件的處理
初始條件:以2005年1月觀測的地下水水位為基礎，採用內插和外推法獲得含水層的初始水位。
邊界條件:各個流量邊界的參數主要考慮模擬初期和模擬期末的流場，有詳細資料的邊界，擬合邊界流入流出量。時間步長由程序控制，每一次運算都嚴格控制誤差。
突水點概化:
按照各突水點涌水量大小的比例，將巷道涌水量統計數據分配到突水點上，概化為抽水井(圖5－7)，模擬巷道對含水層的疏排效果。
圖5-7 巷道突水點概化為抽水井後的分布示意圖
5.模型的識別與檢驗
模型的識別與檢驗過程是整個模擬中極為重要的一步工作，通常要反復地修改參數和調整某些源匯項才能達到較為理想的擬合結果。模型的這種識別與檢驗的方法也稱試估校正法，它屬於反求參數的間接方法之一。
運行計算程序，可得到這種水文地質概念模型在給定水文地質參數和各均衡項條件下的地下水位時空分布。通過擬合同時期的流場和長觀孔的歷時曲線，識別水文地質參數、邊界值和其他均衡項，使建立的模型更加符合研究區的水文地質條件，以便更精確地定量研究模擬區的補給與排泄，預報地下水位。
模型的識別和驗證主要遵循以下原則:①模擬的地下水流場要與實際地下水流場基本一致，即要求模擬地下水水位等值線與實測地下水水位等值線形狀相似;②模擬地下水的動態過程要與實測的動態過程基本相似，即要求模擬與實際地下水水位過程線形狀相似;③從均衡的角度出發，模擬的地下水均衡變化與實際要基本相符;④識別的水文地質參數要符合實際水文地質條件。
由於林南倉礦水文地質條件復雜，邊界條件缺少量化數據資料，且礦區水文地質勘探程度較低、具有系統長觀資料的觀測孔較少且分布不均，使模型中的不確定因素較多，這給模型的參數識別帶來了較大困難。為了克服由於資料不足而造成參數識別的多解性和不穩定性，書中對本區的水文地質條件、水動力場、構造應力場觀測資料作了詳細的研究以幫助分析參數的空間分布規律，確定參數的取值空間。在模型正演及反演過程中，採取了多種手段，使參數、水位、水量等的求解及預測預報盡量符合客觀實際。模擬期各觀測孔水頭高度的變化見圖5-8～圖5-13。
在7個觀測孔中，由於6號觀測孔觀測資料只有4個月的，因此實際用到的觀測孔只有6個(其中7#觀測孔的觀測資料有15個月的)。對139次模擬值與觀測值進行對比。根據計算結果，擬合結果的誤差見表5-15。
圖5-8 1#觀測孔水位擬合圖
圖5-9 2#觀測孔水位擬合圖
圖5-10 3#觀測孔水位擬合圖
圖5-11 4#觀測孔水位擬合圖
圖5-12 5#觀測孔水位擬合圖
圖5-13 7#觀測孔水位擬合圖
表5-15 模擬結果誤差
表5-16 擬合誤差統計表
水文地質參數包括含水層各方向的滲透系數、給水度。通過擬合，得到奧灰岩溶含水層和煤14－K3含水層滲流場擬合圖 ( 圖5-14 ～ 圖5-15) ，從圖中和表5-16 可以看出，擬合的效果很好。
圖5－14 奧灰含水層的滲流場擬合圖
圖5－15 煤14-K3含水層滲流場擬合圖
由長觀孔擬合曲線、流場擬合曲線可知，所建立的模擬模型基本達到模型精度要求，符合研究區水文地質條件，基本反映了地下水系統的動態特徵，故可利用模型進行地下水位預報。
最後，識別後的水文地質參數分別見圖 5-16 ～ 圖 5-21 所示。
圖5-16 奧灰含水層 x 方向滲透系數
圖5-17 奧灰含水層 y 方向滲透系數
圖5-18 奧灰含水層 z 方向滲透系數Figure5.18 z direction permeability coefficient in aquifer in Ordovician limestone
圖5-19 煤14-K3( Ⅱ) 含水層 x 方向滲透系數Figure5.19 x direction permeability coefficient in aquifer in coal 14-K3( Ⅱ)
圖5-20 煤-14 K3( Ⅱ) 含水層 y 方向滲透系數Figure5.20 y direction permeability coefficient in aquifer in coal 14-K3( Ⅱ)
圖5-21 煤14-K3( Ⅱ) 含水層 z 方向滲透系數Figure 5. 21 z direction permeability coefficient in aquifer in coal 14-K3( Ⅱ)
6.滲流場預測
利用前面已經識別好的模型，調整模擬周期，在目前開采條件和疏排水條件不變的情況下，對煤14-K3含水層以及奧灰含水層10年後的滲流場進行預測(圖5-22～圖5-23)。
圖5-22 10年之後煤14-K3含水層滲流場預測圖 Figure5.22 Seepage fieldprediction map in aquifer in coal 14-K3after ten years
圖5-23 10年之後奧灰含水層滲流場預測圖 Figure5.23 Seepage field prediction map in aquifer in Ordovician limestone after ten year
7.突水性分析
根據突水系數方法，突水的危險性由水壓和隔水層厚度決定。水頭越高、有效隔水層厚度越薄，突水危險性越大。以下結合滲流場預測結果，分析突水的危險性。
應用識別好的滲流模型預測了2019年的14砂岩含水層和奧灰含水層滲流場的變化規律(見圖5-24和圖520)。14-K3含水層的水位基本持平，整體呈現出北部高、南部低的特點。在2019年形成了地下水向正南方向匯集的規律，且水頭仍處於高位，尤其是該含水層的水位很大，為－290.67m左右(見圖5－19)。在其他條件不變的情況下，14-K3含水層的地下水在高壓下極易通過垂向的充水通道進入巷道而引起涌水或者突水，因此其突水的危險性比較大。在北部地區，由於煤14-K3含水層地下水水力梯度比較高、地下水水位等值線比較密集，一旦受采動裂隙的影響，煤14-K3含水層的水就很可能優先進入采動裂隙中，很可能造成孔隙水壓力增大而突水，因此在此區應重視采動裂隙的觀測與控制。奧灰含水層水位的變化規律與14-K3含水層基本相似，整體呈現出北部高、南部低的特點。在2019年形成了地下水向正南方向匯集的規律，且匯集區比且14-K3含水層的范圍要小。但水頭仍處於高位，尤其是該含水層的水位很大，為－3.67m左右(見圖5－23)。在其他條件不變的情況下，奧灰岩溶含水層的地下水在高壓下極易通過垂向的充水通道進入巷道而引起涌水或者突水，因此其突水的危險性比較大。在北部地區，由於奧灰岩溶含水層地下水水力梯度大、地下水水位等值線比較密集，水壓很大，奧灰岩溶含水層的水很可能通過導水構造優先進入采動裂隙中，造成孔隙水壓力增大而突水，因此對奧灰含水層的水位需進行長期的動態的觀測。
8.礦井涌水量預測
利用前面已經識別好的模型，將2005年1月的地下水位流場作為初始水位帶入模型。在采區均勻布設4口假想井，調整各井抽水量，使采場水位下降到合理位置，則四口井的總抽水量即是未來的礦井涌水量，模型運行時間就是未來的排水時間。
具體需要疏降的水位，除主要考慮14煤層底板承受的水壓外，還需考慮奧灰水量大、不易疏乾等特點，重點要考慮煤14-K3含水層的排水。在礦區平均分配4個抽水井(圖5－24)，每個井的抽水量為600m3/d，即分別將目前排水能力提高10%，30%和50%，判斷在煤14-K3含水層的水位得到明顯的降低的情況下，排水的經濟性(圖5－25～圖5－27)。結果表明，在強有力的排水條件下，14-K3含水層的水位明顯下降。
圖5-24 所加抽水孔位置圖
圖5-25 提高10%抽水能力的14K3含水層的滲流場 Figure5.25 Seepage field in aquifer in coal 14-K3after improved 10% pumping capacity
圖5-26 提高30%抽水能力的14K3含水層的滲流場 Figure5.26 Seepage field in aquifer in coal 14-K3after improved 30% pumping capacity
圖5-27 提高50%抽水能力的14K3含水層的滲流場
從圖中可以看出，疏水降壓效果都比較明顯，但相對而言，提高30%的排水能力和50%的排水能力使得地下水位下降的效果差不多。因此，在以後的開采方案中，應選擇提高30%的排水能力的方案。按照此方案，當排水能力增大30%時，經濟上承受的壓力比較小而且效果比較顯著，建議採取提高30%的疏水降壓的方法。
9.涌水量預測成果分析
由於Feflow強大的功能，它的Fluid flux analyzer模塊能對模擬期既定時段的既定區塊進行Flu xinside和Flux outside分析，預測出煤14-K3含水層的A，B，C，D，E，F6個區域在模擬期的不同年份的涌水量(表5-17，圖5-28)。
表5-17 煤14-K3含水層分區涌水量預測成果
圖5-28 涌水量分區圖
根據上表所提供的煤14-K3含水層分區涌水量預測成果，對上表進行線型趨勢分析(圖5－29～圖5－34):
圖5-29 A區涌水量預測成果圖
圖5-30 B區涌水量預測成果圖
圖5-31 C區涌水量預測成果圖
圖5-32 D區涌水量預測成果圖
圖5-33 E區涌水量預測成果圖
圖5-34 F區涌水量預測成果圖
從A，B，C，D，E，F區的涌水量預測圖可以看出，除A區的涌水量轉為0，其餘的B，C，D，E，F區涌水量雖然有所起伏，但從線性圖可以看出涌水量的趨勢都是增大的。因此可以看出，前面所提出的對煤14K3含水層進行疏干降壓是必要的。
具體涌(突)水路徑見圖535所示，由圖可以看出:
圖5-35 突水路徑平面圖
圖5-36 現有突水點突水路徑立體圖