高强度采矿活动对地下水影响的数值模拟研究

吕情绪 ，狄军贞 ，李 果 ，林 鑫

（1.神华神东煤炭集团有限责任公司, 陕西 神木 719315；2.辽宁工程技术大学 土木工程学院, 辽宁 阜新 123000）

0 引 言

我国正处于工业化、城镇化加快发展的历史阶段，能源需求总量仍有增长空间[1]。在相当长时期内，煤炭仍是保障我国能源安全稳定的基础能源，将持续为我国宏观经济的高增长作出贡献[2]。煤炭资源与地下水资源联系紧密，在开采过程中会引发含水层疏干、河沟断流等一系列环境问题，还可能发生矿井水害事故，严重威胁生产[3-6]。我国内蒙古地区煤炭资源储量丰富、赋存稳定，易于进行大规模、高效率、高强度的开采，但内蒙古地区植被稀少、荒漠化严重，本就干旱的气候遇上高强度的开采，水资源短缺形势日益严峻，因此更应重视该地区煤矿的地下水资源保护。采矿活动对地下水扰动研究归根结底是对地下水流运动的探究。目前，针对煤矿地下水流场变化的数值模拟研究已有大量报道[7-10]，ZHANG 等[11]认为平顶山十矿-270 m 水平与-300 m水平同属一个水文地质单元，遂采用比拟法，以掘进面积作为变量，来预测采掘-300 m 时矿井涌水量，预测结果符合实际。张凤娥等[12]采用数值模拟法建立了神府矿区大柳塔矿区地下水平面渗流模型，模拟预测了采矿活动下地下水渗流场的动态变化过程，分析了地下水水位随采掘范围变化而波动的动态规律。白晓等[13]运用地下水流数值模型和降水量预测模型，预测分析了峰峰矿区岩溶地下水资源量和地下水流场的变化。由此可见，数值模拟法应用于矿区地下水演变具有重要意义。GMS(Groundwater Modeling System)是美国杨百翰大学环境模型实验室联合美国军队排水工作站研发的三维可视化地下水模拟软件。该软件集成了多种所需的地下水模拟系统，将它们以模块化的形式分类，主要包含了解决地下水三维流动问题的MODFLOW、解决地下水溶质传输问题的MODPATH、MT3D、解决饱和非饱和地下水流动及溶质传输问题的FEMWATER 模块，还有MAP、TINs、GIS 等多种辅助模块，对解决地下水问题的功能十分齐全。模拟过程中能根据研究区实际水文地质条件插入并调节点、线、面等对象，解决研究区不规则边界概化的难题；不仅如此，GMS 软件还可以利用钻孔数据生成Solid 模型，该建模方法能刻画逼真的地层空间结构，含、隔水层分布状态一目了然，可视化效果明显；数值软件不仅自身功能强大，还可与其他系统提供数据兼容，达到联合使用的效果，充分发挥其信息处理和分析功能[14-18]。

现以内蒙古神东矿区某井田为例，利用数值模拟法，预测矿区高强度开采下地下水流场的动态变化，探究高强度采动条件下地下水均衡状况和多含水层系统水流场的演变规律。

1 研究区概况

1.1 自然地理概况

研究区如图1 所示，位于内蒙古与陕西交界处，属温带大陆性半干旱气候，年降雨量小、蒸发量大，研究区位于乌兰木伦河一级阶地的西缘，总体地势呈西北高、东南低的楔形特征。井田内径流主要有乌兰木伦河、补连沟和呼和乌素沟。

图1 研究区位置与等高线地形Fig.1 Location and contour map of study area

1.2 采掘概况

目前该煤矿开采设计能力为28 Mt/a，工作面采高达5.5 m，深度达500 m，工作面长度300 m，推进长度4 500 m。根据范立民[19]对高强度开采的定义（表1），该煤矿属于空间高强度开采。采用多盘区开采方式和全部垮落、一次采全高的长壁综合机械化采煤方法。

表1 煤炭资源开采强度划分指标Table 1 Coal resources exploitation intensity classification index

2 地下水数值模型的建立

2.1 水文地质概念模型

2.1.1 含、隔水层概化

通过对研究区地层结构和岩性进行剖析，将含水层自上而下划分为第四系孔隙潜水含水层、第三系红土隔水层、白垩系志丹群含水层、侏罗系安定组隔水层、侏罗系直罗组含水层和侏罗系延安组隔水层，共3 个含水层和3 个隔水层。研究选取矿井的主要充水来源志丹群和直罗组含水层为模拟的目标含水层，地层模型与剖面情况如图2 所示。

图2 地层结构与剖面Fig.2 Stratigraphic structure and section

2.1.2 边界条件概化

水平方向上：根据水文地质资料分析，研究区北面的呼和乌素沟、南面的活鸡兔沟、东面的乌兰木伦河切割出一个较为独立、完整的水文地质单元，活鸡兔沟和呼和乌素沟为季节性河流，故设定为定流量边界；东侧乌兰木伦河设置为定水头边界；西侧边界为人为设定的边界，采用软件中的通用水头模块进行模拟，流量为软件系统自动运算调整。

垂直方向上：上方松散含水层可接受大气降水入渗和灌溉回归入渗补给，主要排泄方式为人工开采和侧向径流。底部边界为12 煤至22 煤底板，渗透能力较差，可作为垂向上的隔水边界。

2.1.3 地下水流动特征概化

区内含水层岩性各不相同，地层介质的渗透性能随空间位置和水流方向的不同而变化，可将含水层概化为非均质各向异性。地下水系统各运动要素随时间发生改变，且含水层在垂向上存在相互的水力联系，因此将研究区含水层概化为三维非稳定流。

2.2 水流模型模拟

2.2.1 模型离散

模拟区面积约276 km2，根据水文地质条件的概化和研究精度要求，在水平方向和垂直方向上离散为有限数量的长方体，圈定的平面网格尺寸为327.8 m ×286 m，剖分75 行75 列，共5 625 个单元格；垂向上6 层共33 750 个单元格，其中研究区活动单元格为18 028 个。

2.2.2 模型识别与验证

模型识别验证过程中需要遵循以下原则[20]：①以流场为校正依据，需保证模拟地下水流场与实际流场走向与形态基本接近；②以水位为校正依据，需保证模拟地下水水位与实测水位变化趋势一致，两者差值满足精度要求；③区域地下水均衡量应该相差不大，接近于公式计算补给排泄变化量；④校正的水文地质参数要与实际情况对应，不能出入太多。模拟校正过程以2016 年全年地下水位观测资料为基础，1 个月为1 个应力期，每个应力期内划分若干时间步长；再以2019 年全年地下水位观测资料为依据，对模型参数进行再校正。采用试算法，不断地调整试算水文地质参数，使计算值落在95%置信区间范围内。由图3、图4 可看出模型水位拟合程度较好，校正后的水文地质参数见表2。

表2 校正后水文地质参数Table 2 Corrected hydrogeological parameters

图3 白垩系含水层观测孔校核情况Fig.3 Check of observation holes in the Cretaceous aquifer

图4 直罗组含水层观测孔校核情况Fig.4 Check of observation holes in Zhiluo Formation

再以2016 年12 月模拟水位为初始水位，在模型中输入2017-2019 年各项模拟要素，一个月为一个时间步长，共36 个时间步长；根据 2019 年全年观测水位，对模型进行验证，观测孔水位拟合状态如图5 所示。模拟结果误差均在允许区间内；误差值紧密地散落在y=x曲线两侧，说明拟合精度较高，模拟结果可信度高，水流模型校核情况较为理想。

图5 2019 年1 月志丹群和直罗组含水层校核结果Fig.5 Check results of Zhidan Group and Zhiluo Formation aquifer in January 2019

3 煤矿开采对含水层水流场的影响

根据矿井开采计划，通过软件对研究区2020-2029 年这10 年的地下水非稳定流模型进行计算，开采盘区疏干排水量以2019 年各盘区涌水平均值为基础，每2 年增加3%。预测得到直罗组含水层水位流场变化如图6 所示。可以看出：①煤矿开采后，在盘区内形成了多个沿抽排水孔为中心的水位降落漏斗，第5 年末地下水径流条件发生改变，由原先的西北向东南改变为四周向漏斗中心汇聚，在三、四盘区内表现得最为显著。越靠近漏斗中心水力梯度越大、地下水径流速度越快。②随着采空区面积的逐年增大，与初始流场相比，直罗组降落漏斗面积逐渐增加，影响半径不断扩大，漏斗中心处水位持续下降，水位下降幅度从第1 年末的105 m 到第5 年末的337 m，导致周边地下水向漏斗中心补给加剧，该层地下水水位整体下降，至10 年末，该层地下水水位才趋于稳定，水位下降幅度在351 m 左右。因此，开采活动对直罗组含水层影响程度较大。

图6 直罗组含水层水流场Fig.6 Flow field of Zhiluo Formation aquifer

根据数值软件对研究区未来10 年地下水流变化的动态模拟，将志丹群含水层第1 年末、3 年末、5年末、10 年末的预测水位流场与2016 年初始水位流场相比较（图7）。①志丹群含水层流场并未发生明显变化，水位下降情况较为缓和，下降最大处为三盘区，水位下降了20 m。②随着开采年份、开采强度的增加，志丹群含水层水流场发生改变，第3 年末，水位降落漏斗面积持续增大，漏斗中心水位下降了33 m。③到第5 年末，降落漏斗面积进一步扩大，四周地下水逐渐向漏斗中心处补给，由于上部第四系潜水含水层渗透性较好，对下层的补给能力较强，因此，远离采掘区上方的地下水受扰动程度较弱，靠近采掘区上方的地下水流场及水位受影响程度相对较大。④第10 年末，志丹群降落漏斗下降高度达116 m，面积扩大速率逐渐下降，地下水流场变化程度趋于平缓。

图7 志丹群含水层水流场Fig.7 Flow field diagram of Zhidan Group aquifer

利用GMS 软件中的Flow Budget 模块得到地下水区域均衡量，地下水总补给量为357 589.74 m3/d，总排泄量为357 563.62 m3/d，误差为0.007 3%，说明研究区处于均衡状态，模拟计算补给排泄量差距小，说明所建立的水流模型较为合理。

由表3 水均衡统计可知，降雨入渗是研究区的主要补给来源，占比51.90%，其次为河流补给，占比45.21%；河流排泄为主要排泄方式，占比达78.24%，其次为潜水蒸发，占比9.82%，矿井排水占比4.52%。

表3 地下水均衡情况Table 3 Groundwater balance

补给总量与排泄总量均为42 903 m3/d，占区域总补给项和总排泄项的12.0%，补给排泄处于均衡状态。补给总量与排泄总量均为20 297 m3/d，占区域总补给项和总排泄项的5.7%，补给排泄处于均衡状态（图8、图9）。

图8 志丹群含水层的区域水均衡计算Fig.8 Regional water balance calculation chart of Zhidan Group aquifer

图9 直罗组含水层的区域水均衡计算Fig.9 Regional water balance calculation chart of Zhiluo Formation aquifer

4 结 论

1）在采矿活动中，第1 年末，漏斗中心水位下降了105 m；随着采掘强度的加大，水位降落漏斗的面积逐渐扩大，第3 年末下降了176 m；第5 年末地下水径流条件发生改变，降落漏斗中心水位下降幅度增加到337 m；至第10 年末，降落漏斗面积增大不明显、水位下降不显著，稳定在351 m 左右。

2）志丹群含水层地下水流场并未发生较大改变，仅在矿井上方的局部区域水流场发生改变。第1 年末水位降落漏斗逐渐出现，漏斗中心水位下降20 m，至第10 年末，采掘中心处水位下降了116 m。

3）根据区域均衡与层间均衡分析，区域地下水总补给量为357 589.74 m3/d，总排泄量为357 563.62 m3/d，误差为0.007 3%，说明研究区处于均衡状态，模拟计算补给排泄量差距小，说明所建立的水流模型较为合理。降雨入渗是研究区的主要补给来源，占比51.90%，其次为河流补给，占比45.21%；河流排泄为主要排泄方式，占比达78.24%，其次为潜水蒸发，占比9.82%，矿井排水占比4.52%。志丹群和直罗组含水层补给排泄分别占区域总均衡量的12.0%和5.7%，其中矿井排水量为直罗组含水层的主要排泄方式，占比71.3%。