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 现有突水点突水路径立体图