基于矿井生产进度疏干条件下三维地下水数值模拟＊

2014-04-20陈时磊赵颖旺李学渊翟立娟

中国煤炭 2014年8期

陈时磊武强赵颖旺李学渊边凯翟立娟

（1.中国矿业大学（北京）煤炭资源与安全开采国家重点实验室，北京市海淀区，100083；2.中国煤炭地质总局第一水文地质队，河北省邯郸市，056004）

我国是世界上矿井水害事故频发的国家之一，突水事故严重威胁到矿山的安全生产和可持续发展。矿井涌水量不但确定矿井充水程度，而且是设计矿井排水能力和矿井防治水工作的重要依据。矿井涌水量预测不但要考虑到水文地质问题，而且必须与实际矿井生产过程结合起来。这样一方面可以减少对矿区地下水资源的破坏和浪费，有效保证矿山的社会效益和经济效益，另一方面又可以有效保证矿井安全生产。

目前矿井涌水量预测常用方法有大井法、数值法、水文地质比拟法和统计方法等。大井法主要采用裘布依井流公式，该公式主要适用于等厚、均质和较简单水文地质条件的地区，对于水文地质条件复杂的地区就不适合了；水文地质比拟法仅适用于水文地质条件相似，开采条件基本相同的地区，因此，在实际中应用范围较小。虽然，GIS技术、突水水源判别技术的拓展、支持向量机等先进方法手段也逐渐引入到矿井涌水量预测，但受自身条件限制有待进一步发展。而数值模拟法的最大优点在于不但能够进行非均质、非等厚等复杂的水文地质条件下矿井涌水量的计算，而且可以做到将复杂水文地质条件和矿井生产进度结合起来进行矿井涌水量预测。

将矿井生产过程和矿区地下水运动过程有机结合起来，以首采区回采工作面月掘进进度为单位，在概化出研究区水文地质概念模型的基础上，采用目前较为先进和成熟的三维地下水数值模拟软件—Visual Modflow 对榆林大海则煤矿首采区工作面进行了涌水量预测，取得了较好的效果。

1 研究区概况

大海则井田位于榆林市西北约30km，东经109°00′54″～109°16′57″，北纬38°22′01″～38°30′01″。南北宽约14.7 km，东西长约12.5 ～23.3km，面积为280.03km2。研究区地处毛乌素沙漠南缘，区内以沙漠滩地及半固定沙丘地貌为主。区内地势整体为北高南低，地形起伏变化不大，海拔1210～1327 m。区内自下而上分布的地层有侏罗系下统富县组（J1f）、侏罗系中统延安组（J2y）、侏罗系中统直罗组（J2z）、侏罗系中统安定组（J2a）、白垩系下统洛河组（K1l）、新近系上新统（N2）、第四系更新统萨拉乌苏组（Q3s）和全新统（Q4）。其中侏罗系中统延安组的2＃和3＃煤层是主要可采煤层。该矿是中煤集团重点建设的陕蒙亿吨级煤炭基地的主力矿井之一。矿井设计生产能力30.0 Mt／a，全矿井服务年限78a。

研究区内分布的主要含水层有第四系松散岩类孔隙含水岩组、白垩系碎屑岩孔隙裂隙含水岩组和侏罗系碎屑岩裂隙含水岩组。其中侏罗系碎屑岩裂隙水是大海则矿井的直接充水水源，大气降水入渗、地表径流入渗、白垩系孔隙裂隙水和第四系孔隙水是其间接充水水源。

2 水文地质概念模型

模型计算范围包含整个大海则矿区，西以陕蒙省界为界、东以勘探区东部边界为界。研究区范围如图1所示。本次模拟将研究区作为一个统一的地下含水系统，为建立研究区水文地质概念模型，首先要对研究区实际水文地质条件进行概化，建立矿区多层地下水三维水流模型。根据地层和岩性等资料，将区内地下含水系统在垂向上自上而下划分为6层：第四系潜水含水层、白垩系潜水及承压水含水层、侏罗系中统安定组承压含水层、侏罗系中统直罗组至2＃煤层顶部承压含水层、2＃煤层相对隔水层、2＃～3＃煤层间承压含水层。侏罗系中统直罗组至2＃煤层顶板承压含水层是本次模拟的主要目的层。将整个地下水流态视为三维非稳定流，各层介质均概化为非均质各向异性，模型各层之间均有水力联系。平面上，将模型四周均按通用水头边界处理；垂向上，模型顶部接受大气降水补给，为补给边界，底部以3＃煤层为底界，概化为隔水边界。

图1 研究区示意图

3 地下水水流数值模拟

3.1 数学模型

地下水水流数学模型是主要根据质量守恒和达西定律建立起来的连续性的地下水流微分方程。该数学模型主要由地下水渗流方程、初始条件和相应的边界条件构成。

根据水文地质概念模型，建立下列与之相适应的三维渗流数学模型：

式中：Kx、Ky、Kz——分别代表x、y、z方向的渗透系数，m／d；

h——水头，m；

S——自由面下的含水体储水系数，1／m；

μ——重力给水度；

p——大气水入渗补给量，m／d；

ε——模拟区的源汇项，1／d。

3.2 时空离散及初始条件

采用加拿大Waterloo公司开发的三维可视化地下水数值模拟包Visual Modflow 对上述模型进行数值模拟。根据模拟要求，采用矩形有限差分的离散方法对含水介质进行自动剖分，网格单元格为338×207×6 （行×列×层），共计419796个单元格。其中有效单元为293850个，无效单元125946个。将一个完整抽水试验段划分为3个应力期，共36个时段。每个时段划分为若干个时间步长。

3.3 水文地质参数选取

（1）降雨入渗系数和蒸发系数分区及取值。根据矿区实际气候资料、地形地貌特征、潜水地下水水位等资料在平面上对降雨入渗条件和地表蒸发条件进行分区并赋初值。同一分区入渗及蒸发强度相同。

（2）含水层参数分区及取值。根据矿区现场试验和室内试验确定的实际水文地质参数，结合前人资料，对边界上水力传导系数和含水层参数进行分区并赋初值。

3.4 模型识别与验证

采用2013年9月实测的各层地下水水位，经内插和外推法获得模型各含水层的初始水位。通用水头边界上的各层水头值由实测值经插值获得。

根据2013年9月份钻孔ZL-1～ZL-8的抽水试验数据，在上述模型的基础上反复调整相关水文地质参数及源汇项，通过地下水位计算值与观测值的拟合比较对模型反演计算求得模型各层参数分区的水文地质参数。直罗组至2＃煤层顶板含水层参数分区特征和水文地质参数识别结果见图2和表1。从图2和表1可以看出，直罗组至2＃煤层顶板含水层渗透系数和储水率被划分成16个分区，渗透系数和储水率总体均较小，从北往南成条带状分布，表现出明显的空间变异性，水平渗透系数变化范围为0.0153～0.0815 m／d，储水率最大为0.0000635m-1，最小0.000000761m-1。

图2 直罗组至2＃煤层顶板含水层参数分区图

表1 直罗组至2＃煤层顶板含水层参数识别结果

选取钻孔ZL-6和ZL-7的水位观测值和计算值的拟合曲线见图3和图4。所选孔的曲线拟合情况均比较理想，大部分水位观测点计算动态曲线与实测动态曲线基本吻合。模型优选的水文地质参数和源汇项是比较合理的，能较好地反映实际地下水空间分布特征和地下水流运动规律，为下一步准确预测矿井涌水量提供坚实基础。

图3 ZL-6孔、ZL-7孔地下水水位历时拟合曲线图

4 矿井涌水量预测与评价

4.1 预测方案

大海则煤矿首采区主要分成206 试验面、20601工作面、2 0 8 0 1工作面和2 0 1 0 1 4工作面（见图1）。根据矿井未来生产进度和工作面的总长度，将4个工作面回采按月为单位依次划分为若干个工作进度期。其中206试验面分12个工作进度期，20601工作面分工作19个进度期，20801工作面分19个工作进度期，20101工作面分17工作进度期，如表2所示。