基于MODFLOW与气候模式的矿区地下水流模拟

2021-05-31刘中培韩宇平

中国农村水利水电 2021年5期

刘中培，李鑫，陈莹，韩宇平

（1.华北水利水电大学，郑州450045；2.煤炭开采水资源保护与利用国家重点实验室，北京100011；3.水利部水资源管理中心，北京100038）

0 引言

煤炭是重要的矿产资源，由于我国国民经济近年来取得高速发展，煤炭的需求量也随之与日俱增［1］。为了保障煤矿的正常开采以及矿区的生活用水需求，要对煤矿进行疏干排水，以解决煤矿开采所引起的涌水问题［2］。长时间的疏干排水会导致矿区周围地下水水位下降，打破地下水的平衡状态，破坏地下水资源，进一步加剧地下水短缺的局面［3］。因此，研究煤矿开采对地下水水位的影响对于地下水资源的可持续利用具有重大意义［4］。

国内外的众多学者就煤炭开采对地下水水位影响的问题做了大量研究。陈社明等［5］通过Visual MODFLOW 软件构建了研究区的地下水流数值模型，并以此预测了未来研究区各含水层水位在煤炭开采过程中的变化特征，分析矿山排水对不同含水层水位的影响。BISWAS H 等［6］在分析沃克尔盆地水文地质条件和水资源利用的基础上利用MODFLOW 构建了研究区地下水流数值模型，并利用模型模拟了3 种场景来研究补给减少和开采量增加对地下水水位的影响。白晓等［7］将ARIMA 模型预测的降水量用于MODFLOW 的预测模型之中，对研究区岩溶地下水资源量和水位动态变化进行了模拟和预测。汪丽芳等［8］建立了多种开采情境下的地下水位动态预测模型，并利用多元回归方法分析地下水水位与地下水开采量之间的关系。这些研究成果为地下水数值模拟技术的应用提供了重要的参考依据，但是传统的利用数值模拟技术预测地下水水位往往忽略了未来气候变化及煤矿的采掘进度，预测结果不能较准确地反映地下水水位随着气候及采掘进度的动态变化［9］。

本文根据研究区的地下水水位监测数据和详细的水文地质资料，利用Visual MODFLOW 建立研究区地下水流数值模型，并依靠实测数据对模型进行识别和验证［10］。最后利用验证好的模型结合气候模式所预测的未来降水量以及研究区的采煤规划，预测研究区在煤炭开采条件下潜水含水层的水位变化趋势，探究煤炭开采对地下水水位的影响。

1 研究区概况

研究区位于鄂尔多斯市伊金霍洛旗乌兰木伦镇境内，包括乌兰木伦煤矿和柳塔煤矿，属于干旱半干旱大陆性气候，四季寒暑巨变。研究区位于毛乌素沙漠的边缘地带，陕北黄土高原的东北部，地势北高南低，东高西低。从总体趋势看，研究区发育的小构造受褶皱构造控制明显，一般发育于褶皱的翼部，基本呈北西与北东两组。区内地表水系发育，西北紧邻公涅尔盖沟，西南倚靠乌兰木伦河。

研究区内地下水可分为新生界松散层孔隙潜水和中生界碎屑岩类孔隙、裂隙水两大类。松散层孔隙潜水含水层是区内主要的含水层，包括地势较高处的风积砂含水层和沟谷地带的冲洪积湖积砂砾石层含水层，接受大气降水和地表水的补给，沿基岩顶面流动，在沟谷处排泄。基岩裂隙潜水——承压水含水层主要为侏罗系地层中的砂岩，潜水区补给来源为大气降水的直接补给和松散层孔隙潜水含水层的补给，沿地层倾向流动；承压水区补给来源较为单一，主要为径流补给和基岩裸露区段大气降水的直接补给，径流排泄、泉水排泄是该含水层的持久排泄方式。

2 地下水流模型构建

2.1 水文地质概念模型

2.1.1 三维地质模型构建

研究区地表大部分被第四系覆盖，仅在乌兰木伦煤矿西邻的乌兰木伦河东岸有少量基岩出露。区内沉积地层由老到新为：三叠系上统延长组、侏罗系中下统延安组、侏罗系中统直罗组、侏罗系中统安定组、白垩系下统伊金霍洛组、第四系上更新统、第四系全新统。根据对地层资料及钻孔资料的分析，将模型分为七层。进行三维地质建模时首先利用CAD 提取原始地形图的地表高程数据，将带有高程的地形点坐标输入到建模软件GEO5 中自动生成三维地形模型；其次输入钻孔资料并对各个岩层的颜色及名称进行标注，创建水平地层模型；最后生成如图1所示的三维地质模型。综合地形、地势、地质条件等构建的三维地质模型基本展现了矿区的地质特征［11］。并且此模型通过任意切换视角和剖切等手段可以清晰、真切地反映地层复杂的内部结构，为水文地质模型的构建奠定了基础［12］。

2.1.2 含水层结构概化

研究区主采煤层为1～2 煤层，位于潜水含水层下部的隔水层中，煤炭开采主要会对潜水含水层产生影响，因此，本次模拟含水层为第四系松散层孔隙潜水含水层。潜水含水层岩性主要为风积砂、灰黄色、浅黄色黄土、残坡积砂土及冲洪积砂砾石层，厚度为8～63.2 m；下部隔水层为灰白色、褐紫色泥岩、砂质泥岩及粉砂岩，厚度为20～50 m。含水层结构示意图如图2所示。

2.1.3 边界条件的概化

垂向边界的概化：因为潜水含水层与外界发生降雨入渗、蒸发排泄等垂向水量交换，所以模拟区顶部以潜水含水层的自由水面为系统上边界，由于研究区的最低潜水水位埋深仍为12.7 m，故本次模拟中忽略了蒸散发等带来的影响；底部以泥岩、砂质泥岩等构成的隔水层为系统下边界［13］。

侧向边界的概化：研究区西北部的公涅尔盖沟上游常年有水，对潜水含水层存在补给，概化为流量边界（流入边界），公涅尔盖沟下游地势较低，水位低于区域地下水位，概化为流量边界（流出边界）。研究区西南部的乌兰木伦河属于季节性河流，河流水位常年低于地下水位，概化为流量边界（流出边界）；东侧乌兰木伦煤矿井田边界保留煤柱，并且位于地表分水岭附近，概化为隔水边界。水文地质概念模型如图3所示。

2.1.4 水力特征概化

研究区的含水层的分布范围较广，在常温常压下的地下水系统遵循质量守恒定律和达西定律［14］。参数随着空间的变化而发生改变，垂直与水平方向上的渗透系数有所不同，因此将参数概化为非均质各向异性［15］。由于地下水系统的渗流运动要素随时间和空间变化，故将地下水系统概化为非稳定流［16］。总之，研究区可概化为非均质、各向异性、三维非稳定地下水流系统［17］。

2.2 数学模型

对于非均质、各向异性、空间三维结构、非稳定地下水流系统，数学模型如下［18］：

式中：Ω为渗流区域；H为含水层水头，m；Z为含水层底板高程，m；K为渗透系数；Kxx、Kyy、Kzz分别为x、y、z方向上的含水层渗透系数，m/d；μ*为贮水系数；W为含水层垂直补给强度，m/d；E为地下水蒸发排泄强度，m/d；P为含水层开采强度，m/d；h0为含水层的初始水头，m；q为含水层侧向单宽补排量，m2/d，流入为正，流出为负，隔水边界为0；Γ1为隔水边界；Γ2为流量边界；n→为边界上的外法线方向［19］。

2.3 模型的离散及参数分区

研究区的面积大小为57.7 km²，在水平方向上将模拟区域剖分成100 m×100 m的单元格。渗透系数及给水度等水文地质参数主要通过研究区水文地质试验和经验值获得，以钻孔ZK3509 抽水试验数据为例，详见表1。根据试验数据和水文地质条件等将含水层共分为7个参数分区［20］，如图4所示。

表1 研究区典型钻孔抽水试验成果表Tab.1 Pumping test results of typical boreholes in the study area

2.4 模型识别及验证

模型的识别期为2017年1月1日到2017年12月31日。利用2017年的地下水水位实测数据进行参数识别，从水文地质参数的初值出发，根据参数变化范围和实际水位差值对模型进行分区反演试算，直到模型的计算值与研究区的实测值拟合较好为止［3］。最终得到各监测井所在点的计算值与实测值的拟合曲线。部分监测井的拟合曲线如图5所示，含水层水文地质参数识别结果见表2。

表2 研究区水文地质参数表Tab.2 Hydrogeological parameters of the study area

模型的验证期为2018年1月1日至2018年5月31日。将2018年实测降水量、开采量等输入到已经建好的模型中，对区内地下水位的计算值和实测值进行比较，如图6所示，验证期内计算水位接近实测水位，二者相差较小。

模型识别期和验证期的各监测井的计算水位与实测水位的均方误差分别为0.069 m2和0.015 m2。结果表明所建模型可以较好地模拟地下水流情况，可应用于地下水水位预测。

3 未来气候变化情景构建

未来气候预估是指气候系统响应温室气体和气溶胶的排放、浓度或辐射强迫等情景所作出的预估，通常是基于气候模式的模拟结果［21］。气候模式是当前预测未来气候变化情景的重要手段，其模拟数据采用ISI-MIP 提供的GFDL-ESM2M，HadGEM2-ES，IPSL-CM5A-LR，MIROC-ESM-CHEM，Nor-ESM1-M 等5 个全球气候模式日值数据的插值和修正的结果，空间水平分辨率为0.5°×0.5°，时间范围是1951-2050年［22］。本文选用IPCC AR5 中确定的“典型浓度路径”（RCPs）作为未来情景，涵盖了RCP2.6、RCP4.5、RCP6.0和RCP8.5等4个情景［23］。

首先使用实测降水格点数据对气候模式模拟数据集模拟的降水进行初步检验。实测数据来源于1961年1月以来的中国地面水平分辨率0.5° × 0.5°的日值格点数据，而气候模式的基准期为1951-2000年，所以本文以1961-2000年为基准期。然后利用中国气象科学数据共享服务网提供的实测降水值与气候模式模拟数据集模拟的降水值进行对比，最后通过逐月同比例缩放法对降水值做进一步修正。对于降水结果的修正，可用如下公式进行求解，见公式（1）和公式（2）。

式中：λGCMs为修正系数为基准期内逐月多年实测降水平均值为基准期内逐月多年模拟降水平均值；x′GCMs-f，i为未来时段修正后的模拟值；xGCMs-f，i为未来时段逐月模拟值［24］。

3.1 模型预测

由于地下水水位受气象因素中降水量的影响较大，此处选用RCP4.5即中等浓度路径情景作为气候模式的模拟情景，而气候模式的选取则要根据5套气候模式所模拟的逐月多年平均降水量的相对误差来确定［14］，计算结果见表3。

表3 不同气候模式下研究区逐月多年平均降水量相对误差 %Tab.3 Relative error of annual average precipitation in the study area under different climate models

通过表3可以看出由气候模式GFDL-ESM2M 所计算出的模拟值与实测值最接近，因此研究区的未来降水量选取RCP4.5情景下GFDL-ESM2M模式中的模拟结果。

依据2015-2018年研究区涌水量实测值以及未来矿区的采煤规划，预测2020-2035年的地下水开采量。模型预测的初始水位取模型验证期的末刻时间，即2018年5月31日的实测水位，利用已经建立并验证的模型预测未来研究区的地下水流场变化情况。模型潜水含水层的现状与未来地下水流场对比，如图7所示。

通过图7可以看出在煤矿正常开采的条件下，地下水流场变化较为明显。研究区大部分区域都表现出水位降低的趋势，疏干面积也随之增多。东北部区域出现局部的水位回升，而且水位最高点也向西偏移；西南方向因矿区开采形成的地下水降落漏斗面积减小，但是降落漏斗的周围区域出现大幅度的水位下降。

3.2 煤炭开采对地下水水位的影响分析

矿区因不断地开采煤炭而产生的大量矿井涌水不仅影响煤炭的正常开采，还会威胁到人身安全。因此需要及时将涌水排出，但是大量的地下水被集中排出将会导致周边地下水水位打破其原有的平衡状态。为了更加直观地分析煤炭开采过程中地下水水位的变化情况，对比2018年与2035年的地下水流场图（图7）与地下水水量均衡表（表4），分析煤炭开采对潜水水位的影响。

表4 2018与2035年地下水水量均衡万m3Tab.4 Groundwater balance in 2019 and 2035

模拟结果显示，随着多年对研究区地下水集中开采，地下水等水位线由疏变密，地下水流速加快，地下水均衡也发生变化。从2018-2035年，乌兰木伦煤矿北部区域水位回升0.97～17.02 m，东部区域水位下降11.32～28.27 m，南部区域下降0.24～28.74 m，西部区域下降11.32～16.95 m；柳塔煤矿西南区域的地下水降落漏斗范围有所减小，漏斗中心水位回升4.03 m，东部及南部区域水位下降0.24～17.69 m，西部区域平均下降6.44 m。研究区现状地下水补给总量为573.22 万m3，排泄总量为1 038.5 万m3，均衡差为465.28 万m3；2035年地下水补给总量为354.18 万m3，排泄总量为419.09 万m3，均衡差为64.91 万m3。通过表4可以看出，同2018年相比，2035年的各项水量均减少，且补给量与排泄量的均衡差也明显减小，表明研究区的地下水系统逐渐趋向于新的平衡。一般情况下随着地下水位的降低会激发周边地下水的侧向补给，而本文的侧向流入却大幅度减少，这是由于多年的过度开采地下水，导致研究区周边地下水也相应地受到影响，进而无法满足研究区地下水系统的侧向补给。

综上分析，同现状情况相比，煤矿经多年开采，由于大量疏干排水导致研究区内地下水资源总量严重削减，进而造成研究区内的潜水水位普遍下降。随着矿区不断地开采，采空区增多，并且依据采煤规划，由于采煤工作面的转移和采煤数量的减少，地下水开采量也随之减少，故研究区的局部区域的潜水水位有所回升。