岩溶区地下水数值模拟分析及开采方案比选

张凤娟，杨子龙，邢学睿，王立艳，刘淑芹，刘 莉

（1.山东正元冶达环境科技有限公司，山东 济南 250013；2.山东地矿新能源有限公司，山东 济南 250100；3.济南天下第一泉风景区服务中心，山东 济南 250012）

0 引 言

地下水是水资源的重要组成部分，我国许多地区的农业灌溉、工业供水、居民饮用水和畜牧业用水均源自地下水。近年来，随着地下水开采量超负荷增加，地下水超采现象严重，引发了严峻的地下水地质环境问题。与此同时，我国水资源分布不均，部分城市地区存在超采现象，引发了系列生态环境问题，如华北地区地下水超采导致地下水位大幅度下降、地下水污染等问题，通过采取“节、控、换、补、管”等治理措施，大力实施地下水超采综合治理，通过水资源优化配置、地表地下协同互馈、地下水源置换等措施，推动水资源高效利用，促进地表水与地下水相均衡。地下水位持续恢复，有效缓解了地下水漏斗区扩张趋势，减缓地面沉降，对涵养水源、改善生态环境具有重要作用。

赵颖旺等[1]基于有限单元法模拟地下水流运动，对水头反常现象进行分析，据此研究地下水水量与水质特征。李东琦[2]建立了三维数值模型，对地下水开采引起的土体沉降位移进行研究。穆金霞[3]、狄胜同[4]、马海会[5]、蔡建斯等[6]、李文博等[7]通过对岩溶区地下水运移与地下水均衡开展模拟研究，揭示了岩溶区塌陷成因机制。王军进等[8]系统总结了地下水模拟过程、数值模拟方法及应用实例，分析了地下水数值模拟当前存在的问题和不足，并对研究趋势进行展望。刘芮彤等[9]、封丽等[10]采用数值方法对岩溶地下水均衡进行模拟，对地下水运动特征与空间效应进行研究。魏亚强等[11]、金江跃等[12]、高学通等[13]、吴立彬[14]、韩亮[15]将Python、FloPy 和NSGA-Ⅲ等新技术成功应用于地下水模拟，通过案例分析解决了地下水运移规律、污染防治等科学问题，促进了不同学科之间融合渗透，同时解决多元复杂问题。从上述研究可知，针对岩溶塌陷区的地下水开采方案研究尚存在不足，实际工程开采地下水时其引发塌陷灾害防控效果不佳，为此，针对岩溶区地下水开采进行数值模拟，以期为塌陷灾害防治与风险规避提供建议。

1 数值模型建立

1.1 实例概况

本文实例为山东省临沂市岩溶区，研究范围东起临沂市罗庄区汤庄街办-褚墩镇，西至兰陵县卞庄街办，南自磨山镇山南村-旺庄村-涌泉村一带，北至G206 国道附近，涉临沂市兰陵县、郯城县、罗庄区，总面积约144 km2。研究区地处隐伏岩溶区，上覆第四系厚度较薄，由于地下水的过量开采，破坏了地下水的采补平衡，临沂市区已形成地下水位下降漏斗，漏斗中心最大降深达30 余m，且边缘还在不断扩大，境内原有日涌水量3 000 m3以上的泉20 余处，大部分相继出现了断流现象，导致原有提水工程基本报废。由于超采地下水，使其含水层的自然状态发生变化，从而加速了污水的渗透能力，污水渗入地下，严重污染了地下水。地下水污染监测表明，临沂市区水井水质合格率仅在5%左右。其中，水源井的分析结果如下所述：自来水厂5 号井的大肠杆菌数高达230 个/L，沂蒙化肥厂、付家屯两井的大肠杆菌数为11～18 个/L，均超过了《生活饮用水卫生标准》（GB 5749—2022）的要求。自来水厂几口深井抽出的地下水含有黑色悬浮物，水泛黄并有异味。岩溶地下水作为临沂市的重要优质供水水源，在临沂市的社会经济发展中发挥了巨大的作用，但由于长期不合理的开发利用，导致水位快速下降，引发了严重的岩溶塌陷等地质环境问题。

因此，对临沂地区岩溶地下水资源进行评价，研究科学的开采方案，预测地下水位动态特征，并提出地下水资源合理的开发利用和保护措施是十分紧迫且必要的。

1.2 模拟范围与含水岩组结构

根据前人研究成果，研究区为具有单独补给、径流、排泄条件的水文地质单元。因此，模拟范围东边界、西边界、南边界至水文地质单元边界，北边界以弱透水断裂为边界。研究区可分为3 个含水岩组：松散岩类孔隙、碳酸盐岩类裂隙岩溶、碎屑岩夹碳酸盐岩类岩溶裂隙。其中，碎屑岩夹碳酸盐岩类岩溶裂隙水、孔隙水富水性相对较差，特别是开采孔隙水产生与农业争水的矛盾，因此，不具备集中供水意义。而奥陶系的碳酸盐岩富水性好，因此，以该含水岩组为主要研究对象，且中间有一定厚度的弱透水层阻水，局部地区灰岩直接隐伏于第四系之下，无弱透水层阻水，岩溶水含水岩组与松散岩类孔隙水具有一定的水力联系，应将松散岩类孔隙水与岩溶水作为一个地下水流动系统进行研究。研究区概化为三层：第一层为孔隙水含水层，厚5～40 m；第二层为弱透水层，厚约30 m；第三层为岩溶水含水层，具承压性，厚约300 m。

1.3 模型建立及求解

根据研究区含水层结构特征、边界条件及地下水流场等，在垂向上分为三层，其中，第一层为潜水含水层，第二层和第三层为承压含水层。在平面上剖分成大小为10 m×10 m 的正方形单元，各层均剖分为54 行、87 列，其中，各层有效单元格为1 763 个。模 型 的 识 别 期为2018 年5 月31 日 到2020 年5 月31 日，共731 d；验证期为2020 年6 月1 日至2021 年6 月1 日。每30～31 d 为一个应力期，同一应力期内地下水系统补排项保持不变，每个应力期分为10 个时间步长，由模型根据迭代的误差标准，自动控制时间步长。经多次调参拟合，得到三层含水层的水文地质参数分区图，分别见图1、表1 和表2。

表1 调试后的孔隙水含水层水文地质参数分区表Table 1 Division table of hydrogeological parameters of pore water aquifer after adjustment

表2 调试后的岩溶含水层水文地质参数分区表Table 2 Division table of hydrogeological parameters of karst aquifer after adjustment

图1 孔隙水含水层水文地质参数分区图Fig.1 Partition map of hydrogeological parameters of pore water aquifer

1.4 地下水均衡分析

通过均衡法计算得出岩溶水含水层水量均衡结果见表3。岩溶水总补给量为1.51×104m3/d，总排泄量为−1.35×104m3/d，均衡差为0.16×104m3/d，岩溶水基本处于采补平衡状态。在岩溶水含水层总补给量方面，降水入渗量为1.15×104m3/d，占总补给量的76%，大气降水为岩溶水的主要补给来源；其次是地表水体渗漏补给量，其中，河流渗漏补给量占总补给量的24%，主要为燕子河和沂河支流渗漏补给岩溶水，河水水位高于岩溶水水位，且在褚墩一带和岩溶水具有较密切的水力联系，对岩溶水形成有效补给。在岩溶水含水层总排泄量方面，地下水径流排泄量为−0.75×104m3/d，占总排泄量的55%，岩溶水越过尼山-苍山断裂的次生断裂向南排泄，是岩溶水含水层的主要排泄方式；其次为“村村通”机民井人工开采，占总排泄量的24%；此外，由于隐伏的岩溶水具有承压性，可通过顶托补给的方式向第四系孔隙水排泄，或直接排泄补给洼地。综上所述，岩溶水含水层基本处于自然均衡状态，人为因素对地下水影响程度较小。

表3 模拟期地下水均衡表Table 3 Groundwater balance table in simulation period

2 开采方案模拟与预测分析

2.1 开采方案模拟

根据水文地质调查成果，确定开采目的层位为岩溶含水层，在研究区内富水地段确定了两种开采方案，见表4。方案一：在3 个地点布置井位，分别为褚墩、兰山、虎山岭西。方案二：在方案一的基础上增加一个布井地点，新增地点位于磨山镇西。4 个布井地点共布设11 个抽水井，开采量总计1.1×104m3/d。方案模拟时间为10 a，每年分为2 个应力期，即枯水期与丰水期。分别预测分析2 个方案运行后地下水系统的均衡变化情况，以及运行10 a 后地下水变化情况。在岩溶塌陷高易发区褚墩附近增加了一眼观测孔，预测不同方案对该处地下水位的影响。

表4 开采方案一览表Table 4 List of mining plans

2.2 开采方案对地下水位的影响

现状条件下岩溶水受人为影响较小，基本处于自然均衡状态，运行10 a 后，没有出现明显的地下水降落漏斗，地下水仍然由北向南流，南侧是岩溶水的主要排泄途径。由表4 可知，褚墩观测孔水位基本处于稳定波动状态，地下水位随丰枯水期变化，但没有明显的上升趋势。运行方案一后，地下水位产生明显变化，在布井地点出现了3 个小型地下水降落漏斗，最大降深位于褚墩附近，最大降深约为2.5 m，其次为虎山岭附近，降深约为2.3 m，兰山附近地下水位降深约为2.1 m，表明开采量越大其地下水位降深也越大。运行方案二后，地下水位变化更加明显，除方案一产生的3 个降落漏斗之外，在磨山镇西形成了一个较大的降落漏斗，其最大降深达3.1 m。褚墩、虎山岭和兰山的地下水降深分别为2.0 m、2.2 m、2.5 m。磨山镇西增加了0.2×104m3/d 的开采量，小于兰山的0.3×104m3/d，但是降深却大于兰山的降深，主要由于两方面原因。①补给来源问题：该处岩溶含水层上有一较厚的隔水层，且抽水后得不到有效补给，抽取的水量主要为岩溶含水层的弹性释水量，因此水位降深大；②含水层结构问题：磨山镇西岩溶含水层隐伏深度大，有效含水层厚度薄，要想提供更多的水量，只会产生更多的降深。

此外，实施开采方案后，观测孔水位有所下降，虽然下降幅度不大，但有可能改变研究区东北边界的补排方向。实施开采方案后，会袭夺部分径流排泄量，导致径流排泄量减小。

图2 展示了不同方案条件下观测孔水位预测过程线。由图2 可知，方案一条件下，数值观测孔水位在第1 279 d、第1 644 d、第2 019 d 和第3 105 d 时水位降至灰岩顶板标高（39.9 m）以下，将会引发岩溶塌陷风险。方案二条件下，数值观测孔水位枯水期时接近灰岩顶板标高，但一直维持在灰岩顶板标高以上，既充分开发了该富水区的地下水开发利用潜力，又保证了生态地质环境安全，因此，方案二满足生态约束水位条件下地下水资源充分利用的要求。

图2 不同方案条件下观测孔水位预测过程线Fig.2 Process line of water level prediction of observation hole under different scheme conditions

2.3 地下水源汇项对开采方案的响应

当应力条件发生改变时，地下水系统的源汇项也会随之发生改变，或激发系统更多的补给量，或袭夺系统的径流排泄量。由表5 可知，增加开采条件下，地下水源汇项发生变化，主要表现为补给量的增加与径流排泄量的减少。方案一与方案二在增加开采条件下，激发了补给项中的河流渗漏量，分别增加了17%和22%。排泄项中，地下水径流排泄量分别减少了40%和53%。这主要是由于开采方案位于地下水径流区，增加的地下水开采方案袭夺了大部分地下水径流排泄量。同时，岩溶水向第四系的顶托排泄量也有所减少。

表5 增加开采条件下地下水水均衡表Table 5 Groundwater balance table under increased mining conditions

综上分析，方案一、方案二不仅激发了地下水补给量的增加，同时也使顶托排泄量和地下径流排泄量减少，且地下水均衡状态由正均衡转变为负均衡。

3 结论及建议

1）本文建立的地下水流数值模型，以2018 年5月31 日—2021 年6 月1 日为模拟期，通过合理的参数辨识，获取水文地质参数，模拟结果达到预期，均衡差为0.16×104m3/d。

2）根据区内岩溶塌陷发生规律，选取了若干有代表性的控制点，并以其预测水位不低于控制水位（灰岩顶板标高以上）为原则，利用地下水模型预测了两种规划开采方案下的水位动态演变，提出了两个开采方案：方案一：在褚墩、兰山、虎山岭西共布设9 个抽水井；方案二：在方案一的基础上增加一个布井地点，新增地点位于磨山镇西4 个布井地点，共布设11 个抽水井，开采量总计1.1×104m3/d。方案一在褚墩地区具有岩溶塌陷危险，较方案二更易引发岩溶塌陷。

3）地下水管理是一项系统工程，应坚持系统观念，科学谋划水网建设规划布局，打通水资源高效利用的壁障。应长期、连续开展岩溶塌陷区地下水水位监测工作，尤其要在岩溶塌陷频发地段进行加密监测，以便掌握岩溶水水位动态变化的一手材料，并以此为基础开展岩溶塌陷预警分析研判。