中国是世界上原煤产量大、煤矿水害重的国家之一，矿井突水事故的发生频次与强度长期居高难下, 对矿业人员的生命安全构成严重威胁, 造成的经济损失也十分巨大^[1-2]; 同时, 矿井涌突水对矿区的水资源与水环境也造成了很大的破坏^[3-9]。毛乌素沙漠南缘为沙生植被-水-煤的共存区, 煤矿开采的顶板冒裂可能导致浅层地下水的大量漏失, 引发浅层地下水水位下降、泉水断流、地表径流减少等一系列水环境问题, 进而影响沙漠植被的生长，造成生态环境问题。煤矿突水灾害威胁着矿井安全生产, 矿井涌突水的水资源和水环境问题日趋严重。因此，加大矿井涌突水的预测与预防及其环境影响研究十分必要。

1 矿区水文地质特征 1.1 浅层地下水系统

榆溪河流域地处毛乌素沙漠南缘与陕北黄土高原北缘的衔接部位, 属于煤炭资源富集、生态环境十分脆弱的地区。流域总地势西北高、东南低, 地貌上分为西北部的沙漠滩地区、东南部的黄土丘陵沟壑区以及中部的榆溪河河谷区(见图 1)。榆林市多年平均降水量373.5 mm, 多年平均蒸发量1 932.8 mm。榆溪河流域总面积4 938 km², 榆林站常年流量11.75 m³/s。

榆溪河流域西北部的沙漠滩地区广泛分布第四系冲湖积砂层, 沙漠滩地区的北部分布有白垩系洛河组砂岩层。第四系砂层与白垩系砂岩层的厚度大、渗透性强, 是区内浅层地下水赋存与运移的重要场所和主要供水含水层; 东南部的沙盖黄土丘陵区主要为第四系风积黄土较弱透水含水层, 在沟谷两侧多分布有宽度较窄呈沿沟带状展布的烧变岩强透水含水层, 烧变岩区可出露大泉, 为重要供水水源; 榆溪河河谷区沉积大厚度的冲积、冲湖积砂含水层, 渗透性能强, 傍河开采可形成大中型水源地。在白垩系分布区之南、第四系松散层之下, 不连续分布有新近系泥岩，相对隔水(极弱透水)岩层。白垩系和新近系或第四系地层下伏侏罗系延安组煤系砂泥岩层, 煤系地层为大厚度弱富水和弱透水的岩层。

榆溪河流域浅层地下水, 主要来源于大气降水和少量灌溉水的入渗补给, 主要消耗于浅层地下水蒸发和溢出排泄及分散水井的开采。近年来，采煤区的矿井涌水已成为新的重要排泄途径之一。榆溪河在沙漠滩地区为泉集河, 流经黄土丘陵区，大量汇集降水产生的洪水。依据浅层地下水赋存与运移特征, 结合榆溪河水文特征, 将研究区划分为小壕兔、马合、芹河、金鸡滩、牛家梁、头道河则、水掌湾、色草湾、沙河、赵家湾等10个相对独立的地下水源区(见图 1), 浅层地下水资源总量为221×10⁴ m³/d(见表 1)。

1.2 煤系地层与矿井涌水量

榆溪河流域陕西境内规划有29座矿井(见图 1，表 2)。矿井主要开采3^-1煤层, 煤层厚度3.5~11.5 m, 平均6.9 m; 煤层上覆侏罗系砂泥岩厚度132~463 m, 平均181 m; 上部第四系松散层厚度35~139 m, 平均73 m。目前，除28#和29#矿井尚未采煤外, 其他27座矿井已开采运营, 实测单座矿井涌水量(1.8~1 296)×10⁴ m³/a, 平均142×10⁴ m³/a, 27座矿井总涌水量为3 828×10⁴ m³/a(10.5×10⁴ m³/d)。

1.3 煤-水关系分析

毛乌素沙漠南缘浅层地下水，平面上呈现以泉集河为主线、流域分水岭为边界的地下水流动系统(见图 1)；垂向上呈现从地表到地下的植被-水-煤共存区, 并共同构成以地下水-地表水为纽带的沙生植被-含水层-煤层三位一体的三维空间区域。在天然条件下, 含水层与煤层各自有自身的赋存与变化特征, 浅层地下水与煤层水的水力联系极弱。但是，伴随煤炭开采形成的采空区及其顶板冒裂与塌陷, 将会不同程度地改变矿区地下水的循环条件, 不仅会产生矿井涌突水, 而且可能带来上部浅层地下水含水层的破坏, 导致浅层地下水溢流量减少，浅层地下水位下降, 进而抑制沙漠植被的生长，引发生态环境问题。

煤-水赋存空间组合关系是影响矿区地下水循环与变化的主要因素之一, 其组合关系大致可以分为3类:第一类, 浅层地下水含水层大面积分布, 与煤系地层之间无极弱透水的泥岩层相隔, 主要分布在沙漠滩地区的西北部; 第二类, 浅层地下水含水层大面积分布, 新近系泥岩层大面积分布或呈不连续片状分布, 含水层与煤系地层之间在新近系泥岩分布地区有极弱透水的泥岩层相隔, 在缺失新近系泥岩分布的地区呈现“天窗”衔接, 主要分布在沙漠滩地区的东南部和沙盖黄土丘陵区; 第三类，浅层地下水呈条带状富集在烧变岩区, 直接构成煤层上覆强含水层^[10-15]。显然, 含水层与煤系地层之间有新近系泥岩相隔, 采煤过程中将更有利于保护含水层地下水, 但是只要冒裂带或塌陷穿越到含水层, 浅层地下水的漏失以及矿井大量涌水甚至突水就将不可避免地发生。

2 采煤对地下水影响的数值模拟预测 2.1 地下水三维流数学模型

模型区范围覆盖上述10个相对独立的水源区, 面积约3 000 km²。模型区地下水总体由四周向河流汇流, 其北东、南东、南西边界为自然分水岭, 设置为隔水边界; 榆溪河出口的南侧河谷段, 设置为定水头边界; 北西边界存在地下水侧向径流补给, 设置为定流量边界; 模型区的顶面为浅层地下水面。榆溪河两侧支流排泄区内的地下水, 概化为地下水溢出边界, 榆溪河河谷区的河水与地下水关系密切, 概化为河流边界; 模型区垂向上为由地面至开采煤层底板之间的地层, 总厚度100~770 m, 煤层底板为砂泥岩, 概化为模拟区的隔水底板。模型区存在面状大气降水与灌溉水入渗、面状地下水蒸发蒸腾、线状地下水溢出与河水入渗、分散水井开采(概化为面状排泄)以及矿井涌水(概化片状地下水漏失—变流量边界), 模型区内各层地下水之间存在越流水量交换。模型区的地下水三维流数学模型描述为:

$ \begin{array}{l} \frac{\partial }{{\partial x}}K\frac{{\partial H}}{{\partial x}} + \frac{\partial }{{\partial y}}K\frac{{\partial H}}{{\partial y}} + \frac{\partial }{{\partial z}}K\frac{{\partial H}}{{\partial z}} = 0\;\;\;\;\;\;\;\left( {x, y, z} \right) \in \mathit{\Omega }, t > 0\left( {潜水} \right)\\ \frac{\partial }{{\partial x}}K\frac{{\partial H}}{{\partial x}} + \frac{\partial }{{\partial y}}K\frac{{\partial H}}{{\partial y}} + \frac{\partial }{{\partial z}}K\frac{{\partial H}}{{\partial z}} = {S_s}\frac{{\partial H}}{{\partial t}}\left( {x, y, z} \right) \in \mathit{\Omega }, t > 0\left( {承压水} \right)\\ \mu \frac{{\partial H}}{{\partial t}} = - \left( {K - \varepsilon } \right)\frac{{\partial H}}{{\partial z}} + \varepsilon , \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {潜水面} \right);\\ H = z\\ H\left( {x, y, z, 0} \right) = {h_0}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {x, y, z} \right) \in \mathit{\Omega }, \;\;\;\;\;初始条件\\ H\left( {x, y, z} \right){|_{{\mathit{\Gamma }_1}}} = {h_1}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t > 0, \;\;\;\;\;\;\;\;\;\;定水头边界\\ K\frac{{\partial H}}{{\partial n}}{|_{{\mathit{\Gamma }_2}}} = {q_0}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t > 0, \;\;\;\;\;\;\;\;定(零)流量边界\\ K\frac{{\partial H}}{{\partial n}}{|_{{\mathit{\Gamma }_2}}} = q\left( t \right), \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;变流量边界\\ \frac{{{K_r}A}}{{{M_r}}}({H_r} - H) = {Q_r}。\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;地下水溢出或河水入渗边界 \end{array} $

式中:H为地下水位高程, m; K为渗透系数, m/d; μ为给水度; S_s为弹性释水率，1/m; ε为浅层地下水面垂向交换水量，m³/d·m²; h₀为初始水位高程，m; h₁为定水头高程，m; q₀为定(零)流量，m³/d·m²; q(t)为变流量，m³/d·m²; H_r为河水位高程，m; A为河水淹没面积，m; K_r为河床淤泥层垂向渗透系数，m/d; M_r为河床淤泥层厚度，m; Q_r为地下水溢出量或河水入渗量，m³/d; n为边界外法线方向; Γ₁为Ⅰ类边界; Γ₂为Ⅱ类边界; x, y, z为空间变量，m; t为时间变量，d; Ω为模型域。

2.2 地下水三维流数值模型

模型区采用500 m×500 m的网格剖分, 垂向分为5层：第1层为第四系含水层, 第2层为新近系红土层或白垩系砂岩层或侏罗系砂泥岩层, 第3层为导水裂隙带之上的侏罗系原岩层, 第4层为导水冒裂带及两侧原岩层, 第5层为采空区及原煤层。模型区剖分为204行、140列、5层共142 800个单元, 其中有70 430个为活动单元。

输入地下水初始流场和各层面高程等值线及源汇项等数值, 依据28个长观孔和228个统测孔的地下水位数据以及27座矿井实测涌水量数据进行模型识别, 计算地下水位、矿井涌水量与实测值的拟合效果良好(见图 2，3)。

模型识别获得的水文地质参数为:第四系含水层分17个参数区, 渗透系数(0.1~13)m/d, 给水度0.05~0.10;白垩系含水层分3个参数区, 水平渗透系数(0.3~0.7)m/d, 垂向渗透系数(0.15~0.35)m/d, 弹性释水率1×10^-7m^-1; 新近系泥岩层的水平渗透系数3×10^-5m/d, 垂向渗透系数3×10^-6m/d, 弹性释水率1×10^-20m^-1; 侏罗系砂泥岩层的水平渗透系数3.2×10^-2m/d, 垂向渗透系数5×10^-3m/d, 弹性释水率1×10^-7m^-1。

2.3 地下水三维流数值模型预测与分析

按照矿区规划, 随着煤炭开采，未来30年采空区面积将由现在的27 km²逐步扩大到237 km²; 依据煤层厚度结合采煤方式计算矿井冒裂带高度, 将采空区范围及冒裂带顶界面高程数值代入构建的地下水流数值模型, 预测29座煤矿开采30年末的矿井涌水量见表 3，浅层地下水流场与降深场见图 4，5。

煤矿开采30年末的预测矿井总涌水量为26.00×10⁴m³/d, 约是现状矿井总涌水量(10.50×10⁴m³/d)的2.5倍; 矿井涌水量主要来自浅层地下水蒸发量的减少和浅层地下水溢出量的袭夺以及地下水储存量的释放, 所占比例分别为43.1%，6.8%，49.3%。另外，牛家梁水源区的河水增渗补给约占0.8%。矿井涌水并未改变区域浅层地下水向榆溪河径流排泄的态势, 仅在矿井涌水影响的6个水源区形成了近似椭圆形的浅层地下水位降深场, 其与矿井分布形态是基本一致的；预测30年末的降深场中心最大水位降深约6 m, 预测到30年末浅层地下水位下降速率不大于22 mm/a。

3 采煤涌突水的敏感性预测分析 3.1 矿井涌突水预测的影响因素分析

区域地下水流数值模型的结构、参数与渗流场等数据，多因资料不足而难以精确把控, 尤其是煤矿采空区的冒裂带高度以及冒裂带顶部岩层(简称覆岩)的渗透性能大多也很难准确判定，都是导致模型预测结果严重偏离实际的重要影响因素。采空区放顶后产生岩层冒裂和整体下沉, 引发的岩层移动变形与开裂, 极有可能导致覆岩渗透性能的增强, 但是其增强的程度一般难以定量测定, 甚至定性估量也非常困难。

区域地下水流数值模型承载着有效预测与合理评估矿井涌突水量的重任, 这就需要谨慎地考量覆岩垂向渗透性能变化对矿井涌突水量影响的问题, 即需要进行煤矿开采涌突水的敏感性预测分析, 用于定量评估这些不确定的因素对模型预测结果的影响程度, 以期为研究煤水关系和制定采煤保水方案等提供依据。

3.2 预测煤矿选择与环境设置

在小壕兔、马合、芹河、金鸡滩、二道河则水源区分别选择29#，26#，25#，13#，20#矿井(见图 1), 当加载不同的煤矿开采规模(采空区)工况时, 分别对29#，26#，25#，13#，20#矿井的覆岩给定不同垂向渗透系数值进行矿井涌水量预测, 定量研究矿井涌突水量与覆岩渗透系数的关系。

3.3 数值模型预测的成果分析

利用构建的地下水三维流数值模型, 分别对29#，26#，25#，13#，20#等5个煤矿进行数值模拟, 预测的30年末覆岩不同垂向渗透系数的矿井涌突水量，结果见表 4，图 6。

预测结果表明:①矿井涌突水量与覆岩垂向渗透系数呈正向非线性变化规律, 覆岩垂向渗透系数小于0.005 m/d时, 为原岩层的渗透参数, 矿井涌水量基本稳定; 渗透系数为0.005~0.01 m/d时, 反映覆岩有一定破裂, 矿井涌水量呈现小幅增加; 渗透系数大于0.01 m/d时, 表明覆岩有较大破裂, 矿井涌突水量大幅增加; 渗透系数大于1.0 m/d时, 覆岩之上的含水层会被快速疏干, 矿井涌突水量基本上不再增加。采空区上覆岩垂向渗透系数为(0.01~1.0)m/d, 对应的矿井涌突水量变幅最大、增速最快。②对应覆岩垂向渗透系数0.005 m/d和10 m/d的矿井涌突水量之差值(简称极值差)与比值(简称极值比), 大致反映覆岩之上含水层地下水对矿井涌突水量的贡献及其影响。5矿的极值差为(0.09~83.57)×10⁴m³/d, 该数值基本上代表了覆岩之上含水层地下水的漏失水量; 5矿的极值比为1.58~18.76, 其大小基本反映了覆岩之上含水层地下水对矿井涌突水量的影响程度。在沙漠滩地区腹地，矿井的极值比大(26#与13#), 是因为覆岩之上含水层的富水性强所致; 在沙盖黄土丘陵区的矿井的极值比小(20#), 是由于覆岩之上含水层富水性弱的原因; 在沙漠滩地区沟域分水岭地带的矿井的极值比中等(25#与29#), 与该地区覆岩之上含水层中等富水有关。

4 结论

1) 煤矿开采的矿井涌水量主要来自浅层地下水蒸发量的减少和浅层地下水溢出量的袭夺以及地下水储存量的释放, 所占比例分别为43.1%，6.8%，49.3%;煤矿开采涌水的影响, 在研究区6个地下水源区形成了近似椭圆形的浅层地下水位降深场, 其与矿井分布形态基本一致；预测30年末的降深场中心最大浅层地下水位降深约6 m, 30年末浅层地下水位下降速率小于22 mm/a。

2) 矿井涌突水量与覆岩垂向渗透系数呈正向非线性变化, 煤矿采空区之上覆岩垂向渗透系数(0.01~1.0)m/d, 对应的矿井涌突水量变幅最大、增速最快。矿井涌突水量的大小, 主要受覆岩垂向渗透性能强弱以及覆岩之上含水层富水程度的控制; 覆岩垂向渗透系数0.005，10 m/d的矿井涌突水量的极值差与极值比, 基本能够反映覆岩之上含水层地下水对矿井涌突水量的贡献与影响。

3) 煤矿采空区导水裂隙带一旦导通上覆第四系或白垩系含水层(即覆岩垂向渗透性能的大幅增强), 采空区覆岩之上的含水层地下水就会大量漏失, 甚至出现疏干含水层的现象, 这对浅层地下水流场和溢出量都会造成巨大的影响。