许昌市建安区北部采空区三维地下水流数值模拟

王灵敏, 王军强, 刘荣慧, 李卓然, 常伟丽， 郑萍

(河南省地质矿产勘查开发局第二地质勘查院,河南 郑州 450045)

近年来，黄昕霞等[1]、MARC-ANDRE Lavigne等[2]、ELY D M等[3]、吴乐等[4]、王建红等[5]先后应用数值模拟技术分别建立了山西运城盆地、纽约沙托河流域、哥伦比亚高原地区、北京市西山地区、黑龙江黑河流域中游平原区的地下水流数值模型，对其地下水水位变化特征进行了模拟研究。笔者利用数字模型对地下水位变化特征进行研究。

文中研究区位于河南省许昌市建安区北部，为非工业区，主要人类活动为道路建设、农业生产、矿山开采等。苏桥镇武庄村南部许昌铁矿于2013年初开始开采，在许昌铁矿矿山基本设施建设期间,建安区北部出现地下水水位大幅下降、地面塌陷、地面沉降等环境地质灾害问题。2013—2015年，该研究区地下水水位降幅达2～3 m；2013年,武庄村发生地面塌陷，并造成1人死亡；2015年,矿区内水泥道路突然发生崩裂错断，形成9处陷洞，部分房屋墙体产生裂缝[6]。上述环境地质问题影响了当地居民的生产和生活。

该地区盛产铁矿，目前针对这些铁矿开展的大多数工作集中在矿产勘查及环境地质问题初步调查阶段，对当地环境地质问题的发展趋势及原因尚未有研究。本文在前人调查资料的基础上,应用当代水文地质学研究和应用的最重要手段之一——地下水数值模拟技术[7]对研究区内地下水流问题进行了研究，建立了研究区的三维水文地质模型，并对该区域地下水水位动态及变化趋势进行了研究。研究结果将为该研究区合理开采地下水资源提供参考。

1 研究区的水文地质条件

研究区为山前冲洪积平原，地势平缓，西北向东南地形坡度为5‰。研究区地下水的赋存和分布规律以地层岩性为基础[8]，区内广泛沉积新生界巨厚松散沉积层，厚度50～400 m。松散孔隙含水层划分为3个含水段(In1、In2、In3)，如图1所示。

图1 研究区水文地质剖面图

图1中自上而下划分为：第四系孔隙潜水、浅层承压水含水段(In1)、新近系孔隙承压含水段(In2)和古近系孔隙承压含水段(In3)[9]。含水段(In1)，含水岩组为第四系钙质结核层，含钙质结核粉土，为中等富水区，单井5 m降深涌水量100～500 m3/d。含水段(In2、In3)，含水岩组为冲洪积粉细砂、中砂层。研究区内含水段In2和In3为富水区，单井15 m降深涌水量1 000～3 000 m3/d[10]。各含水段间隔水层岩性主要为粉质黏土、弱成岩粉质黏土，隔水性较好。古近系孔隙承压含水段In3底部为古近系泥砾岩和泥质角砾岩，与基岩呈不整合接触，隔水性较好。本次模拟对象为松散孔隙含水层。新近系含水段In2与古近系含水段In3均为承压含水层，地下水补给来源与排泄方式一致，因此，模拟中将新近系含水段In2与古近系含水段In3合并为一个含水段进行研究。

研究区内第四系孔隙潜水、浅层承压水含水段(In1)主要接受大气降水补给和少量径流补给；以蒸发、农田供水及矿区排水形式排泄。新近系孔隙承压含水段(In2)和古近系孔隙承压含水段(In3)主要接受地下水侧渗补给，以矿区排水及地下径流的方式排泄。研究区内地下水流方向为西北—东南向。

2 研究区地下水的三维数学模型

在分析研究区水文地质条件的基础上，应用连续性原理及达西定律，构建了研究区地下水的三维非稳定流数学模型[11-13]，即：

(1)

3 研究区地下水的数值模拟

使用地下水模型软件GMS(Groundwater Modeling System)中的Modflow模块建立地下水流数值模型。GMS是目前国际上通用的地下水模拟软件，可用来模拟饱和、非饱和流条件下的地下水流和溶质运移等情况[14]。Modflow作为地下水流数值模拟软件，实现了概念模型和数值模型的无缝转化，具备强大的二维和三维可视化功能，可以很好地模拟水井、河流、溪流等对非均质地质条件和复杂边界条件的水位系统的影响[15]。

模拟区边界：东以武庄矿区东部边界为界，北以石固村南部边界为界，西以武庄矿区西部边界为界，南以石梁河为界，总面积为9.1 km2。边界条件概化如图2所示：北部、西部为模型边界1，定为流入边界；东部为模型边界2，定为流出边界；南部为模型边界3，定为定水头边界。

根据研究区的含水层结构、边界条件和地下水流场特征，垂向上将地下水系统结构概化为3层，如图3所示。图3中：第一层为第四系松散潜水—微承压含水层，第二层为黏土隔水层，第三层为新近系松散土层承压含水层。模拟区剖分为39行、51列规则网格，3层，各层均采用118 m×118 m的剖分格式。其中有效单元为4 437个，如图4所示。模型的模拟期为2015年9月—2021年9月，识别验证期为2015年9月—2018年9月，模拟预测期为2018年9月—2021年9月。将整个预测期划分为3个应力期，每个应力期为一个相应的自然月，计算的时间步长为15 d。模拟初始水位采用2015年9月的实测水位(只做一次水位测量)，识别验证采用2018年9月的(只做一次水位测量)实测水位。

图3 三维地下水概念模型界面图

图4 研究区单元剖分网格图

模型采用“试错法”(trial-and-error)来间接率定水流模型的相关参数[16]。对含水层和弱透水层进行参数分区，给定初始值，具体见表1。图5为研究区2015年9月的实测水位,即研究区地下水数值模拟的初始水位。在识别验证期，根据实际水位与模拟计算水位的比对结果，不断调整参数，直到模拟水位与实际水位的误差达到允许范围以内。研究区有13眼观测井进行水位拟合，观测井分布位置如图6所示，基本覆盖研究区。

表1 含水层及隔水层参数取值范围

图5 研究区初始水位(单位：m)

图6 研究区地下水观测井分布图

在识别验证期末(即2018年9月)，对研究区13眼观测井的模型计算水位值与实际观测水位值进行对比拟合。识别验证期末研究区观测井实际观测水位值与模型计算水位值见表2。由表2知，二者水位值拟合良好，地下水水位标高误差小于0.4 m。

除此之外，在识别验证期末，对研究区模拟流场与实际流场进行了拟合，结果如图7所示。由图7可知，模拟流场与实际流场总体流动方向相同，大部分节点水位拟合较好，模型参数及源汇项设置符合实际情况。

识别验证期内研究区的地下水水位变幅如图8所示。其中，2018年9月的地下水水位为模型计算水位，2015年9月的水位为实际观测水位。由图8可知，识别验证期内，研究区地下水水位整体变幅较小，研究区东部部分地区如磨李、西张、雷庄区域的地下水水位呈0.0～0.5 m的小幅下降，武庄塌陷区的地下水水位呈0.0～0.5 m的小幅上升，其他区域的地下水水位呈0.5～1.0 m的幅度上升。

表2 研究区观测井2018年的观测水位与模型计算水位 m

图7 研究区水位拟合图(单位：m)

图8 2018年9月的预测水位较2015年9月的实测水位的变幅图

识别验证期内研究区地下水水位变幅较小的原因为：①根据当地气象局资料，研究区2015—2018年的降雨量较2013—2015年的降雨量有所增加，增加量为257 mm/a，因此，研究区地下水的降雨入渗补给量较2013—2015年的降雨入渗补给量大幅增加；②研究区南部石梁河由于矿井排水，水位较高，常年渗漏补给地下水；③由于模拟区属于山前径流区，再加上7～9月雨季入渗的补给，西部水位上升幅度大，水力梯度增大，径流排泄增强，因此，东部地区水位有小幅下降。

4 模型预测

建立水流模型的目的是为了预测在设定地下水开采条件下研究区内地下水流的动态变化特征，从而判断该种开采方案是否科学，为合理利用地下水提供建议。2012—2015年,研究区出现水位大幅下降、地面塌陷及地裂缝等一系列环境地质问题。目前，研究区内水位下降的原因及地下水流动态特征亟待解决。模型预测期选为2018—2021年，开采方案为现状开采，即保持现有2018年的矿井抽水量和地下水开采项、开采量、开采井位置不变。

2021年9月与2018年9月地下水水位对比如图9所示。由图9可知，在矿井保持2018年地下水抽水量不变的条件下，武庄区域地下水水位降落漏斗区的中心水位出现上升，该地下水水位降落漏斗区域的地下水水位升幅最大，为1.0～1.5 m；雷庄区域地下水水位降落漏斗区的中心水位出现上升，该地下水水位降落漏斗区域的地下水水位升幅为0.5～1.0 m；研究区其它大部分区域的地下水水位也趋于平缓上升，上升幅度小于0.5 m,地下水水位的变化幅度较小。说明研究区地下水水位相对稳定。

图9 2021年9月的预测水位较2018年9月的实测水位的变幅图

上述表明，2013—2015年，武庄区的矿井运营初期，因矿井大量抽水，导致研究区内地下水排泄量大于补给量，研究区内地下水水位大幅下降；2015—2018年，研究区地下水形成地下水水位降落漏斗，激发研究区外围补水，补给量大于排泄量，地下水水位基本稳定；预测期的2018—2021年，地下水补给量依然大于地下水排泄量，地下水水位降落漏斗区的中心水位上升，区域地下水水位呈整体小幅上升趋势。

5 结语

1)通过资料收集、水位统调、抽水试验等水文地质调查工作，应用GMS(Modflow模块)地下水流模拟软件建立了研究区松散孔隙含水层三维数值模型，模型经过识别验证，能够较好地模拟研究区的地下水流。

2)模型模拟结果显示，2015—2018年，研究区的地下水水位变幅较小，武庄矿区周围区域的地下水水位升幅在0.5 m左右，模拟区东部雷庄区域的地下水水位降落漏斗区与西张区域的地下水水位降落漏斗区，这两个地区的地下水水位有0.2～0.5 m的下降。

3)模拟预测结果显示，在现状开采条件下，即矿井保持2018年抽水量继续抽水，研究区大部分地区的地下水水位趋于稳定及平缓上升，上升幅度小于0.5 m。武庄区域的地下水水位降落漏斗区的中心水位升高，区域地下水水位升幅为1.0～1.5 m；研究区东北部雷庄区域的地下水水位降落漏斗区的水位也略有回升，升幅为0.5～1.0 m。

4)2012—2013年，研究区由于矿井抽水形成的地下水水位降落漏斗已经激发了远程的侧向补给，研究区内地下水水位降落漏斗中心点水位升高，区域水位整体上升，水位降落漏斗进入恢复涵养期[17]，研究区地下水已形成新的均衡。