0 引 言
石炭二叠系是我国主要的含煤岩系，其下伏的奥陶系灰岩是占我国煤炭储量近60%的华北地区区域性强含水层，也是该区域煤矿引发底板突水事故的重大隐患[1-3]。近年来，随着我国东部矿区开采水平的不断延深，使得底板突水事故时有发生，造成严重的人员伤亡和经济损失。新安矿煤层底板重要的突水水源即为奥灰岩溶水。该矿区奥灰含水层出露地表，且小浪底水库蓄水导致地表水和地下水水位均有所抬升，使得底板突水危险性进一步增大。因此，能够确定底板突水危险区域，并预测危险区域奥灰突水量，可以为矿井防治水工作提供科学依据。
国内部分水文地质工作者针对煤层底板突水量预测做了一些研究工作。高卫东等[4]为更好地解决支持向量机(SVM)核参数和惩罚因子的取值对煤层底板突水量等级预测精度的影响问题，提出利用全局搜索能力较强的粒子群优化(PSO)算法优化支持向量机参数；曹庆奎等[5]针对煤层底板突水问题的小样本，非线性特点，采用支持向量回归算法对突水量进行预测，避免了定性分析的局限性；肖建于等[6]提出基于贴近度的模糊证据理论，建立基于模糊证据理论的煤层底板突水状态估计与评价模型。这些研究大多集中在定性的突水量等级划分，而对于实际突水量的定量预测研究较少。
为了在煤层底板突水危险区域实际突水量的定量预测研究；在确定底板突水危险区域的基础上，计算出危险区域底板奥灰水量，并认为奥灰水在具备导水通道的情况下全部突出，从而实现底板突水量预测。笔者拟通过概化新安矿水文地质条件，构建矿区奥灰水渗流基本微分方程并求解；利用Visual MODFLOW软件实现矿区奥灰水渗流运动的三维可视化数值模拟；依据模拟水位值采用水均衡法计算突水危险区域的奥灰水量。以此作为底板突水量定量预测的尝试性研究。
1 基本地质及水文地质条件
1.1 矿井位置
新安矿位于河南省洛阳市新安县城以北15 km处。井田范围东以F29和F2断层为界，西至F58断层，浅部以二1煤露头线为界，深部以二1煤底板-200 m等高线为界。井田走向长约15.5 km，倾向宽约3.5 km，面积约50.26 km2。
1.2 地质条件
新安矿区缺失下石炭统、泥盆系、志留系和上奥陶统地层，地层由新至老分别为第四系、二叠系、石炭系、奥陶系和寒武系地层。新安井田总体为一平缓单斜构造，位于新安倾伏向斜北翼。在单斜构造背景上发育小的波状起伏或次级褶皱，使煤层底板等高线发生不同程度的弯曲。井田内少发育大中型断层，多发育小型断层，且多为斜交正断层。
1.3 水文地质条件
新安矿二1煤底板至奥灰含水层之间的地层从上往下依次为硅质泥岩隔水层、L7灰岩含水层、砂质泥岩隔水层、L1-3灰岩含水层和铝土质泥岩隔水层，如图1所示。抽水试验表明：L7、L1-3灰岩含水层富水性较弱，且整体厚度不大，具有注浆改造的客观条件。注浆改造后，可以形成有效隔水层。此外，L7、L1-3灰岩上覆层和下伏层均为软岩，形成软硬岩互层，本身也可以看做隔水层。而奥灰含水层顶界距二1煤底板平均距离为53.76 m，对矿井安全生产构成重要威胁。
图1 煤层底板水文地质结构模型
Fig.1 Hydro-geological structure model under coal seam floor
2 水文地质条件概化
2.1 边界条件
新安矿底板奥灰岩溶水系统是一个相对独立的水文地质单元。西北部碳酸盐岩裸露区是补给区。西南部F58断层为该水文地质单元的阻水边界。奥灰岩溶水自补给区向矿区径流，受奥灰地层南东向倾伏的影响，埋深逐渐增大，深部岩溶不发育，岩溶水形成滞流带，致使径流改向北东方向，最终排泄于黄河。
2.2 含水层结构
二1煤底板主要含水层为奥灰含水层和太原组薄层灰岩含水层。奥灰含水层分布范围较广，厚度400～500 m，岩溶发育，富水性强[7-10]。太原组共有3～4层薄层灰岩，岩溶较发育，含裂隙岩溶水，可以看做隔水层。奥灰是煤层底板主要充水含水层，据此将该含水层概化为单层结构含水层，以奥灰水为渗流计算对象。
2.3 含水层水动力特征
一般认为，岩溶水具有非连续性和各向异性渗流特征。但该矿区奥灰岩溶比较发育，不同规模的断层、裂隙相互交叉，使得岩溶水渗流宏观上又具有多孔介质渗流特征。即全区有比较统一的地下水渗流场，不同区位水位差别不大，且具有相似的水位动态。据此将该矿区岩溶水渗流概化为连续多孔介质中二维平面渗流。
2.4 补给和排泄
矿区奥灰岩溶水系统具有相对独立的补给和排泄条件。在补给区，系统接受大气降水入渗和地表水渗漏补给；在径流区，系统主要以人工开采方式排泄，矿井排水主要是顶板水，而非奥灰水；在排泄区，系统通过地下径流方式向外排泄。
综上所述，研究区水文地质条件概化如下：奥灰含水层为单层结构非均质各向同性含水层；岩溶水运移为二维平面运动且符合达西渗流定律；岩溶水为潜水-承压水混合系统；岩溶水主要接受大气降水入渗补给，以人工开采为主要排泄方式。
3 数学模型构建及求解
3.1 构建数学模型
数学模型即用一组数学关系式来表达概念模型。依据渗流连续性方程和达西渗流定律，研究区奥灰岩溶水渗流数学模型可构建为
(1)
式中：kxx，kyy，kzz分别为沿x，y，z坐标轴方向的渗透系数；h为点(x，y，z)在t时刻水头；Ss为贮水率；W为源汇项；t为时间；n为二类边界内法线向量；q(x,y,z,t)为流量在时间和空间上的变化函数；f2(x,y,z,t)为水位在时间和空间上的变化函数；D为含水层确定的渗流域；Γ1为渗流域D的已知水位边界；Γ2为渗流域D的已知流量边界；H0(x、y、z)为初始水头。
3.2 数学模型求解
采用有限差分法求解该数学模型，将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法采用Taylor级数展开等方法，把控制方程中的微商用网格节点上函数值的差商代替，从而建立以网格节点上的值为未知数的代数方程组。
假设地下水的密度为一个常数，表示水量平衡的连续性方程为
(2)
式中：Qi为单位时间内进入单元i的水量；Δh为时间Δt内水位变化量；ΔV为单元的体积。
式(2)表示在Δt时间内，水位变化Δh时水量的变化量，流入为正，流出为负。根据达西定律，行方向上从单元(i，j-1，k)流入(i，j，k)的流量为
(3)
式中：i，j，k分别为每个单元的行号、列号和层号；qi,j-1/2,k为通过单元(i，j，k)和(i，j-1，k)间界面的流量；Ri,j-1/2,k为单元(i，j，k)和(i，j-1，k)间沿列方向的渗透系数；Cj为沿列方向单元面的面积；hi,j,k为单元(i，j，k)的水头值；Δrj-1/2为单元(i，j，k)和(i，j-1，k)间的距离。
同理，通过其他5个界面的流量均可按照式(3)类推。
此外，tm时刻水头差与时间差的商值，可由式(4)获得
(4)
据此并结合Taylor级数，可得有限差分方程为
(5)
式中：
为外部源项；qi,j,k,n为外部汇项；ΔVk为计算单元的体积；ΔCjΔVk为沿列方向单元面的面积；A为水头对于时间的偏导数之差分近似表达式。
因此，可通过该方程计算出区域内各点的水位值。
4 数值模拟
Visual MODFLOW软件广泛用于模拟地下水在孔隙介质中的渗流。但大量实践表明，该软件亦可用于模拟地下水在裂隙介质中的渗流[11-12]。笔者利用Visual MODFLOW软件对新安矿奥灰水渗流运动进行数值模拟。
4.1 模型的空间离散与时间离散
模拟范围为整个新安矿区。模拟岩层分为两层：上层为底板隔水层，下层为奥灰含水层。坐标轴方向与地理坐标方向一致：X轴方向指向正东方向，Y轴方向指向正北方向。模型空间范围沿X轴方向12 031 m，Y轴方向13 205 m，平面共剖分单元60×60个，根据大气降雨周期将一个完整水文周期划分为12个应力期，365个时段。
4.2 水文地质参数设置
模型涉及的重要水文地质参数为渗透系数和贮水率。笔者依据矿区抽放水试验数据，结合岩性、导水性进行参数区划分，给出各分区参数初值和变化范围。根据模型拟合效果对初值进行调整，并利用PEST模块进行参数优化，得到最终参数分区如图2、图3所示。
图2 渗透系数分区
Fig.2 Partition of permeability coefficient
图3 贮水率分区
Fig.3 Partition of storage coefficient
图2、图3中各参数分区对应的参数PEST优化值见表1、表2。
表1 渗透系数PEST优化值
Table 1 PEST optimization value of permeability coefficient
分区渗透系数/(m·s-1)x方向y方向z方向16.51×10-66.51×10-61.64×10-626.90×10-59.58×10-66.90×10-734.05×10-54.05×10-54.05×10-541.12×10-61.12×10-61.74×10-751.97×10-51.97×10-51.97×10-663.68×10-58.68×10-58.68×10-672.78×10-52.78×10-42.78×10-781.84×10-41.84×10-41.84×10-597.52×10-57.52×10-57.52×10-6105.21×10-55.21×10-55.21×10-6113.70×10-56.60×10-63.7×10-6122.89×10-52.89×10-51.86×10-6132.31×10-52.31×10-52.31×10-6146.37×10-56.37×10-56.37×10-6151.22×10-51.22×10-52.66×10-6164.05×10-54.05×10-56.88×10-6171.33×10-41.33×10-41.33×10-5
4.3 源汇项
4.3.1 补给处理
新安矿区奥灰岩溶水主要接受大气降水入渗补给。奥陶系灰岩出露面积为32 km2，沿西北边界呈线状分布。模拟时，将其处理为沿西北边界的线状注水井，将大气降水入渗量作为初值加至相应的计算单元。降水入渗系数依据经验取0.3。
4.3.2 排泄处理
新安矿区奥灰岩溶水主要的排泄形式为矿井人工供水。2002年11月1日，2号和新2号水源孔开始供水，总供水量约为200 m3/h。该排泄项以抽水井的形式分配到相应的计算单元。
表2 贮水率PEST优化值
Table 2 PEST optimization value of storage coefficient
分区Ss/m-1分区Ss/m-1分区Ss/m-110.01×10-321.6031.4042.4056.0068.0072.480.8491.00101.20110.42121.80131.40140.15151.00161.80170.55——
4.4 模型识别与验证
以新安煤矿长观孔奥灰水位一个水文年(2011年3月至2012年2月)的观测数据为基础进行Kriging插值，得到全区奥灰水位，作为初始水位值赋给模型，形成初始流场。取时间步长为30～31 d，模型持续运行12个月。通过不间断地计算并调整模型参数和边界条件，最终确立模型。初始水位和模拟水位等值线分别如图4、图5所示。
图4 初始水位等值线
Fig.4 Isoclines of initial water level
图5 模拟水位等值线
Fig.5 Isoclines of simulated water level
在进行三维模型校验过程中，将研究区长观孔同样设置为模型中的观测孔，对实际观测值和模拟水位值进行比对，以验证模型的拟合程度。具有代表性的长观孔拟合误差见表3。
水位拟合误差一般4%左右，最大不超过11%。模拟水位动态变化与实际水位动态变化基本一致，模拟水位接近实测水位，模拟流场总体特征符合研究区水动力场变化特征。
表3 长观孔水位拟合误差
Table 3 Fitting errors of long-term observation hole water level
日期(年-月)5号孔(16区)水位/m观测值拟合值误差22号孔(招南)水位/m观测值拟合值误差2011-03235.00244.9880.042 5242.79245.0070.009 12011-04235.00244.6710.041 2240.51245.1190.019 22011-05235.00244.3120.039 6238.50245.1410.027 82011-06235.00244.0230.038 4248.00245.1590.011 52011-07235.00243.8780.037 8252.00245.2050.027 02011-08235.00243.5000.036 2260.00245.2550.056 72011-09235.00243.3580.0356274.40245.4300.105 62011-10235.00243.0480.034 2269.56245.8280.088 02011-11240.23242.4970.009 4279.02246.2170.117 62011-12231.84242.3070.045 1278.94246.4690.116 42012-01232.50242.1790.041 6268.99246.3660.084 12012-02232.00242.2690.044 3269.84246.1890.087 6
5 突水量预测
在系统分析研究区底板突水致灾机理的基础上，结合区内采动破坏、断层破碎以及已知奥灰突水点等因素，确定了煤层开采过程中可能的突水威胁区域之一，如图6蓝色区域所示[13-15]，该区域发生了1995年特大奥灰突水事故。笔者采用水均衡法依据模拟水头值进行突水量预测，突水量可近似看作流入突水威胁区域的奥灰水量之和。
图6 均衡区划分
Fig.6 Division of water balance region
依据Visual MODFLOW均衡区模拟分析计算，2012—2016年突水威胁区域奥灰水量的预测结果见表4。
表4 奥灰水量预测值
Table 4 Forecasting values of mine water inflow
月份涌水量/(m3·h-1)2012年2013年2014年2015年2016年11 165.721 162.951 160.781 159.481 157.4821 165.181 164.931 164.281 162.241 160.2431 163.811 164.671 162.501 161.051 153.3741 164.901 163.701 162.221 161.31 149.9451 165.081 164.481 163.491 160.421 146.4761 164.241 162.631 160.251 158.751 142.6871 164.991 162.951 162.831 161.701 138.4181 164.961 162.651 160.681 158.601 134.1891 165.811 166.281 161.801 160.20876.90101 165.041 162.351 160.851 159.86880.34111 162.781 162.181 160.581 158.05879.45121 165.681 164.271 162.141 159.54877.85
从表4可知，基于水均衡原理预测的突水威胁区域底板奥灰水量为1 150 m3/h左右。该预测结果与1995年发生的特大奥灰突水事故平均突水量1 244.3 m3/h接近(误差7.5%)，具有可信性。此外也发现，在小浪底水库蓄水后，奥灰水量的预测值并没有显著增加，究其原因主要是因为新安矿将奥灰水作为供水水源，在一定程度上起到了疏降效果，抵消了水库蓄水后的水位抬升作用。
6 结 论
1)对研究区奥灰含水层边界、结构和水动力特征、补给和排泄等水文地质条件进行概化，确定奥灰含水层为：单层结构非均质各向同性含水层；岩溶水运移为二维平面运动且符合达西渗流定律；岩溶水为潜水-承压水混合系统；岩溶水主要接受大气降水入渗补给，以人工开采为主要排泄方式。
2)建立研究区奥灰水渗流运动基本微分方程并采用有限差分法求解，得到有限差分方程。利用Visual MODFLOW软件实现了该奥灰水渗流运动的三维可视化数值模拟。结合2011年3月至2012年2月1个水文年的实际水位观测值对该模型的模拟水位进行验证和校核。模拟水位与实际水位误差最大11%，一般4%左右，较为可靠。
3)依据模型模拟水位值，采用水均衡法计算底板突水危险区域的奥灰水量为1 150 m3/h左右，与1995年该区域发生的特大奥灰突水事故平均突水量1 244.3 m3/h接近，预测结果真实可靠。研究结果为煤层底板突水量预测的尝试，为煤矿安全生产提供了一定的技术依据，亦可为类似矿井地下水三维渗流模型的构建提供借鉴。
[1] 施龙青,曲兴玥,韩 进,等.多模型融合评价煤层底板灰岩岩溶突水危险性[J].煤炭学报,2019,44(8):2484-2493.
SHI Longqing,QU Xingyue,HAN Jin,et al.Multi-model fusion for assessing risk of inrush of limestone Karst water through mine floor[J].Journal of China Coal Society,2019,44(8):2484-2493.
[2] 靳德武.我国煤层底板突水问题的研究现状及展望[J].煤炭科学技术,2002,30(6):1-4.
JIN Dewu.Research status and outlook of water outburst from seam floor in China coal mines[J].Coal Science and Technology,2002,30(6):1-4.
[3] 许延春,黄 磊,俞洪庆,等.基于注浆钻孔数据集的注浆工作面底板突水危险性评价体系[J].煤炭学报,2020,45(3):1150-1159.
XU Yanchun,HUANG Lei,YU Hongqing,et al.Evaluation system for floor water inrush risk in grout-reinforced working faces based on grouting boreholes dataset [J].Journal of China Coal Society,2020,45(3):1150-1159.
[4] 高卫东,王正帅.基于粒子群优化支持向量机的煤层底板突水量等级预测[J].煤田地质与勘探,2012,40(6):44-47.
GAO Weidong,WANG Zhengshuai.Forecast of inrush water volume grade from coal floor based on support vector machine with partial swarm optimization[J].Coal Geology& Exploration,2012,40(6):44-47.
[5] 曹庆奎,赵 斐.基于遗传-支持向量回归的煤层底板突水量预测研究[J].煤炭学报,2011,36(12):2097-2101.
CAO Qingkui,ZHAO Fei.Forecast of water inrush quantity from coal floor based on genetic algorithm-support vector regression[J].Journal of China Coal Society,2011,36(12):2097-2101.
[6] 肖建于,童敏明,姜春露.基于模糊证据理论的煤层底板突水量预测[J].煤炭学报,2012,37(6):131-137.
XIAO Jianyu,TONG Minming,JIANG Chunlu.Prediction of water inrush quantity from coal floor based on fuzzy evidence theory[J].Journal of China Coal Society,2012,37(6):131-137.
[7] 刘再斌,杨小刚.基于水量演化特征的煤层底板突水通道识别方法[J].煤炭科学技术,2016,44(6):152-158.
LIU Zaibin,YANG Xiaogang.Method to distinguish water inrush channel in seam floor based on water quantity evolution feature[J].Coal Science and Technology,2016,44(6):152-158.
[8] 刘再斌.基于孔型组合的煤矿水害区域治理模式研究[J].煤炭科学技术,2018,46(7) :184-189.
LIU Zaibin.Study on water hazard regional control pattern based on different borehole type combination[J].Coal Science and Technology,2018,46(7):184-189.
[9] SUN Y J,XU Z M,DONG Q H.Forecasting water disaster for a coal mine under the XiaolangdiReservoir[J].Journal of China University of Mining and Technology,2008,18(4):0516-0520.
[10] 原富珍,马 克,庄端阳,等.基于微震监测的董家河煤矿底板突水通道孕育机制[J].煤炭学报,2019,44(6):1846-1856.
YUAN Fuzhen,MA Ke,ZHUANG Duanyang,et al.Preparation mechanism of water inrush channels in bottom floor of Dongjiahe Coal Mine based on micro-seismic monitoring[J].Journal of China Coal Society,2019,44(6):1846-1856.
[11] 刘业娇,薛俊华,袁 亮,等.软岩底板突水机理分析及数值试验[J].煤炭学报,2017,42(12):3255-3261.
LIU Yejiao,XUE Junhua,YUAN Liang,et al.Numerical test and mechanism analysis of water-inrush on soft rock floor[J].Journal of China Coal Society,2017,42(12):3255-3261.
[12] 徐智敏,孙亚军,巩思园,等.高承压水上采煤底板突水通道形成的监测与数值模拟[J].岩石力学与工程学报,2012,31(8):1698-1704.
XU Zhimin,SUN Yajun,GONG Siyuan,et al.Monitoring and numerical simulation of formation of water inrush pathway caused by coal mining above confined water with high pressure[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(8):1698-1704.
[13] 朱宗奎,徐智敏,孙亚军,等.基于无量纲多源信息融合的底板突水危险性评价方法研究[J].采矿与安全工程学报,2013,30(6):911-916.
ZHU Zongkui,XU Zhimin,SUN Yajun，et al.Research on the risk evaluation methods of water inrush from coal floor based on dimensionless multi-source information fusion technique[J].Journal of Mining & Safety Engineering,2013,30(6):911-916.
[14] 朱宗奎.底板突水突变机理及模型建构研究[J].矿业安全与环保,2017,44(5):34-39.
ZHU Zongkui.Research on catastrophe mechanism and model construction of floor water inrush[J].Mining Safety & Environmental Protection,2017,44(5):34-39.
[15] 王心义,姚孟杰,张建国,等.基于改进AHP 法与模糊可变集理论的煤层底板突水危险性评价[J].采矿与安全工程学报,2019,36(3):558-565.
WANG Xinyi,YAO Mengjie,ZHANG Jianguo,et al.Evaluation of water bursting in coal seam floor based on improved AHP and fuzzy variable set theory[J].Journal of Mining & Safety Engineering,2019,36(3):558-565.
