摘要
GMS（Groundwater Modeling System）是国际上流行的地下水模拟软件，且与ARCGIS 属性数据有很好的兼容性，本文通过对陈家河煤矿水文地质条件的分析，运用 GMS中的 MODFLOW 和 MT3D 模块建立三维数值模型，进行地下水环境影响数值模拟，研究事故工况下污染物在地下水环境中的浓度变化情况，分析和预测一定时间内地下水可能受到污染的情况。
Abstract
Based on the analysis of hydrogeological conditions of Chenjiahe Coal Mine， MODFLOW and MT3D modules in GMS were used to establish a three-dimensional numerical model to simulate the impact of groundwater environment， study the concentration changes of pollutants in groundwater environment under accident conditions， and analyze and predict the possible pollution of groundwater in a certain period of time.
0 引言
松宜矿区，地处湖北松滋与宜都两市交界处，是湖北省第二大产煤基地。松宜矿区煤炭开采历史悠久，早在200多年前已有村民在此采煤，迄今已大规模开采 50 余年，累计开采量已达可采储量的 90％以上，已经进入衰老期。根据国家及湖北省煤矿相关政策，松宜矿区内各煤矿均已陆续停产闭坑。在煤矿的长期开采过程中，山体内部形成了大面积采空区，严重破坏了区内原始水文、工程、环境地质条件（Yi et al.，2013）。本文以陈家河煤矿为研究对象，利用 GMS 建立地下水流数值模型，预测煤矿排水过程中对地下水的影响，为该流域水资源管理和生态环境的影响进行科学分析提供一定的科学依据（Singh and Atkins，1984，1985； Guo et al.， 2009）。
1 研究区概况
研究区内沟谷发育，主要地表水系为东南部的洛溪河及其左岸支流陈家河、干沟河和右岸支流朱峡沟等，属洈水流域。洛溪河为区内最低侵蚀基准面，最低高程 155.7 m，地表水系流量随降雨量多寡而增减，其排泄条件较好（图1）。
图1区域水系示意图
洛溪河发源于宜都市王家畈乡古水坪，属常年性河流。上游为尖岩河、中游为茅坪河、下游为大河，河源高程 620 m，由北流向东南，经界溪桥折向东北，流经白岩溪、三叉河，再向东南；过尖岩河、穿眼岩至两河口，汇陈家河，向东南流经茅坪，于大河口入松滋市境。特大洪水期最大流量达 500～600 m3 /s，平均流量1.42 m3 /s。
朱峡沟为洛溪河另一支流，其发源于调查区西南部天井寺东麓，在穿眼岩汇入洛溪河，全长 1.8 km，河水主要来源于朱峡沟的猴子洞岩溶泉，因其流经煤矿采空区，出水口为富铁酸性废水，导致朱峡沟河水常年被污染，河水混浊，河床呈黄—黄褐色。
陈家河属于洛溪河支流，位于洛溪河左岸，流经西北部晒经台、北部云台观，在两河口与尖岩河处共同汇入洛溪河，途经青山坡、陈家河、廖家堰煤矿。河道全长5.8 km，河床宽度4～23 m，高程落差 190 m，平均坡降3.28%，流域面积12.3 km2，特大洪水期最大流量达50～100 m3 /s，平均流量0.25 m3 /s。
干沟河为季节性流水河谷，上游为庙河，下游为干沟河，其发源于云台观北麓破瓢岭，向南流经坛子口煤矿、松木坪镇、猴子洞煤矿，于大河口处和洛溪河汇入梅溪河，全长 11 km，河床宽度 10～20 m，为“V”型河谷，大河口处高程 150 m。该河平水期、枯水期均无流水，丰水期遇较强的降雨时有流水，在上游庙河至猴子洞河段，其河床下部及周围多为煤矿开采遗留的采空区，河床水直接下渗补给地下水，导致河床干涸，故为季节性流水河谷。夏季暴雨时，最大洪水流量可达30 m3 /s。
陈家河流域上游主要地层为志留系，泥盆系； 中游主要地层为二叠系；下游主要地层为三叠系，岩性以泥质粉砂岩、石英砂岩、灰岩为主。区内岩溶地貌发育强烈，主要发育3套含水层：松散岩类孔隙水含水层、碎屑岩类裂隙水含水层碳酸盐岩类岩溶水含水层（图2）。
图2陈家河流域水文地质图
2 模型建立
2.1 概念模型
2.1.1 含水层结构概化
根据研究区地质及水文地质条件，同时考虑研究区对地下水环境影响范围及影响程度，以能满足环境影响预测和分析的要求和原则，本次模拟对象为：三叠系下统大冶组、二叠系上统长兴组、二叠系下统茅口组、栖霞组及马鞍组煤层采空层含水层。本次评价范围确定为：北部、南部、西部以地下水分水岭为界，东部以陈家河为界，确定本次工作评价区范围为1.96 km2（图3）。
碳酸盐岩类裂隙岩溶含水岩组是本次重点关注的含水层组，除部分隐伏于第四系之下，其余皆出露于外。
碎屑岩类裂隙含水岩组赋存于侏罗纪下统、石炭系下统和泥盆系中上统碎屑岩地层中，岩性为泥质粉砂岩、粉砂岩、石英砂岩，富水性弱，透水性较弱。
松散岩类孔隙含水岩组分布于第四系冲洪积物中，且多分布在河谷、冲沟口及山坡、山麓地带，主要为黏土、粉质黏土、细砂、砂砾等。其厚度一般为 0.10～8.30 m，平均厚度约 1.56 m，模拟区范围内较薄为 0～1. 00 m，地下水主要赋存于各类砂层的孔隙中。该含水岩组岩性分布不均，各地富水性差别较大。
2.1.2 边界概化
垂向边界：根据含水层特性将模拟区概念为 4 层结构，第一层为孔隙含水层，第二层为煤层（采空区）含水层，第三层为岩溶含水层，第四层为裂隙含水层（刘基等，2017）。
平面边界：区内含水层水平方向延伸面广，含水层整体流向为西北向东南流动，模型中将东部边界概化为一类流量边界，其他边界概化为零通量边界（王双明等，2021）。
2.2 地下水流数值模型
2.2.1 地下水流数学模型
水位动态与地下水补排量是随时间发生变化的，表现出非稳定流特性；在空间上孔隙含水岩组的岩性存在一定差异，反映出系统的非均质性，由于冲积的成层性和后期的压缩作用使得垂向与水平方向存在异性；裂隙、岩溶含水层组也具有非均质性，同时由于受裂隙、岩溶、构造影响的体现出各向异性特性（於波等，2017）。因此，综合表明研究区地下水系统可概化为非均质、各向异性、三维非稳定流进行研究（张志祥等，2017）。
在不考虑黏度等因素时，可用下列微分方程的定解来描述研究区地下水流系统：
图3模拟区范围图
(1)
2.2.2 数值离散化
（1） 空间剖分
根据边界条件并结合本次地下水流数值模拟的实际目的，模拟区面积为1.96 km2，剖分方式采用矩形剖分，利用GMS软件将模拟区域进行如下网格剖分，剖分 200 行，200 列，共剖分单元格 20682 个，模型剖分4层，剖分结果见图4。
（2） 时间离散
根据数值模拟时期的确定原则，将 1990年 1月作为水位模拟的初始年，利用1990年1月—2030年 1月的水位、源汇项数据资料进行模拟，对模型进行识别验证，应力期为40 a，每个应力期分为1个时间步长，共40个时间步长（张耀文等，2021）。
2.2.3 源汇项处理
（1）大气降水入渗补给：根据第四系覆盖区的地表岩性、地貌等因素，概化为几个降雨入渗系数不同的小区，并根据降雨入渗系数经验值给出初值，待模型识别后确定，在模型中大气降水入渗补给量的计算公式为：
(2)
（2）地下水侧向排泄量：模拟区地下水排泄量，可分段采用达西定律计算，公式为：
(3)
研究区地下水补给来源主要为大气降水，根据气象数据结果显示2021年9月—2022年9月区域内年降雨量为 1077.40 mm，月平均降雨量为 89.73 mm。根据孔隙含水岩组地层分布的不均匀性，将大气降水入渗系数分为 4 个区，Ⅰ岩溶强烈发育的碳酸盐岩区；Ⅱ岩溶中等发育的碳酸盐岩区；Ⅲ相对隔水层发育区；Ⅳ松散岩类孔隙发育区。具体分区以及相关参数详见表1和图5。
图4数值模型计算剖分
表1降水入渗系数分区及赋初值
图5模拟区大气降水入渗系数分区图
2.2.4 水文地质参数分区
本模型采用插值法为渗透系数赋值，在模拟过程中适当进行调整，已知参数点参考以往钻孔工作成果，并结合本次拟合验证期调参，渗透系数在模型调试后确定（宋颖霞等，2012；王长友等，2014；赵凯浪，2017；黄倩等，2019）。
疏干给水度、弹性储水系数参考以往计算值和经验值，并结合模型调试结果最终确定。以上 3 类参数具体赋初值详见表2。
表2模型中水文地质参数初始值
2.2.5 水均衡分析
对地下水各均衡要素采用水均衡分析得出识别期、验证期承压水含水层与潜水含水层水量均衡结果。识别期，潜水总补给量为 3.65×104 m3 /d，总排泄量5.53×104 m3 /d，均衡差为1.88×104 m3 /d，补给小于排泄，为负均衡。承压水总补给量为 23.86× 104 m3 /d，总排泄量为 26.59×104 m3 /d，均衡差为 2.73×104 m3 /d，补给小于排泄，为负均衡。
验证期，潜水总补给量为 1.85×104 m3 /d，总排泄量4. 03×104 m3 /d，均衡差为2.81×104 m3 /d，补给小于排泄，为负均衡（朱学愚和谢春红，1990；成建梅和陈崇希，1998；宋业杰，2011；赵春景等，2011）。承压水总补给量为 13.45×104 m3 /d，总排泄量为 24.62×104 m3 /d，均衡差为11.17×104 m3 /d，补给小于排泄，为负均衡（表3，表4）。
2.3 模型预报
2.3.1 初始流场
模拟应力期为 40 a，共 40 个时间步长，将 1990 年 1 月—1991 年 1 月设为稳定流状态，且作为模型识别的初始流场。
基于模型长时间补排数据和研究需要，以1990 年 1 月—2030 年 1 月为模拟期，由于模型涉及的空间及时间尺度均较大，以 30 d 为应力期，将各补排数据处理为月值形式，采用基于二阶准确（AB/TR） 预测矫正法的自动时间步长方案进行时间离散控制。考虑到观测孔水位观测资料的时间段分布特征，选取在 2021 年 9 月—2022 年 6 月对模型进行识别、验证。
表3识别期均衡分析
本次验证期采用水位曲线进行拟合，选择了 2 个地下水井的地下水位动态数据。时间步长为 1 年，每个步长输出 1 个结果，模型运行 302 个时段，记录下每个时段各长观孔所在结点的水位。部分地下水长观孔的地下水动态拟合曲线见图6～图8，水文地质参数见表5。通过模型识别验证结果表明所建模型可以代表本区实际水文地质条件（薛禹群， 1997；李竞生和姚磊华，2003；周焱钰等，2011）。
2.3.2 地下水污染预测情景设定
本次模拟所选择的特征污染物为铁，根据枯平丰三期检测数据铁浓度最高为 327 mg/L，最低为 1.15 mg/L。综合考虑将模拟区所在采空区设为定浓度边界，浓度值按200 mg/L设置。
2.3.3 地下水污染预测
根据溶质计算结果，随着矿洞中高浓度酸性废水的流出，最终矿洞下游及两侧地下水将被污染。据模拟，14600 d 模拟期模型输出铁物质总量为 4.49×1011 mg，输出总铁质量4.49×105 kg。
表4验证期均衡分析
3 结语
（1）在对陈家河煤矿地质、水文地质条件系统分析的基础上，建立了陈家河煤矿概念模型和数学模型，并采用GMS软件建立地下水流数值模型。通过模型的识别和验证，使模型可以更好地反映模拟区地下水实际情况。
（2）根据参数敏感性分析结果可知，渗透系数是地下水位的变化的主要影响因素，其次是大气降水入渗补给系数，各个参数的敏感度系数在率定值附近呈现出对称性的变化规律。采用定量分析法分析参数敏感性，筛选出主要的不确定性因素，为将来开展地下水流随机模拟奠定基础。
（3）建议增加陈家河煤矿流域河道底泥沉积物分析和含水层岩性勘探工作，用来提高研究区水文地质参数的调查精度，因为模型参数的不确定性会给计算结果带来误差，而准确的参数数据可以提高研究区内地下水资源评价与预报的可靠性。
图6长观孔校正目标条形图
图7SZ05长观孔地下水动态拟合曲线
图8SZ07长观孔地下水动态拟合曲线
表5识别验证后水文地质参数
