Dynamic characteristics of groundwater level and exploitable amount of groundwater source in Jiuxian county, Tai'an
-
摘要: 泰安市城区水资源供需矛盾突出，地下水资源开采引发了岩溶塌陷等地质环境问题，亟需查明岩溶塌陷生态水位约束下水源地的允许开采量，确定水源地最优取水方案。以旧县水源地1980—2021年的地下水位、年降水量、地下水开采量、岩溶塌陷等数据为基础，探讨人类活动对水源地地下水位的演变的影响因素，分析发生岩溶塌陷的地下水临界水位。基于Modflow-GWM软件构建泰安城区−旧县岩溶水系统流动与管理耦合模拟模型，探讨防止岩溶塌陷发生的地下水允许开采的资源量。结果表明：（1）孔隙水水位年内动态受降水影响明显，呈现“降水−补给”型特点；岩溶水水位变化呈“枯低丰高”的特征，水位动态变化属降水入渗−开采型；（2）旧县水源地地下水位自1980—1990年呈现大幅下降，1990—2003年岩溶地下水位基本呈现波动下降，2004年后水源地地下水位上升较为明显；（3）研究水位动态与岩溶塌陷关系得出，防止发生岩溶塌陷的水源地临界水位为108 m，处于岩溶含水层顶板以上2 m；（4）通过地下水管理模型，确定在临界水位时模拟区岩溶水可开采量为8.2~8.5 万m3·d−1，其中旧县水源地可开采资源量为3.2~3.5万m3·d−1。Abstract:
In this study, the effect of human activities on the evolution characteristics of groundwater levels in the water source area of Jiuxian county and its influencing factors have been discussed, based on the data of groundwater levels, annual precipitation, groundwater mining output and karst collapse in the study area from 1980 to 2021. The critical groundwater level to prevent karst collapse has also been determined. In addition, based on Modflow-GWM software, a coupled simulation model for the flow and management of karst water system in Tai'an City and Jiuxian county has been constructed, thereby exploring the allowable extraction of water resources in terms of the karst collapse prevention. The research results indicate that, (1) The annual dynamics of interstitial water are significantly affected by precipitation, with the lowest annual water level in late June and the annual highest from late August to late September during the rainy season, presenting the characteristics of "precipitation-replenishment" type. The dynamics of karst groundwater within the year are affected by the obstruction of pore water discharge in the Quaternary system. The highest water level within the year generally lags slightly behind the concentration from September to October, while the lowest water level occurs from May to June. The groundwater level changes by about 5 meters within the year. The water level changes exhibit a characteristic of "a low water level in the dry season and a high water level in the wet season", and the dynamic changes of water levels are characterized with the type of "precipitation infiltration-mining". (2) The groundwater level in the water source area of Jiuxian county showed a significant downward trend from 1980 to 1990, and the karst groundwater level showed a basic fluctuation and decline from 1990 to 2003. The groundwater level in the water source area has increased significantly since 2004. (3) Since 1988, the water source area in Jiuxian county has entered a period of high incidence of karst collapse. From then on to 1994, the water level had rebounded, fluctuating between 96 m and 108 m and was located near the bedrock surface. This period had also experienced a high incidence of karst collapse in this area, with the collapse area expanding northward from Jiuxian county to the Yuanzhuang-Yanglou area. After 1995, the intensity of karst collapse in the area decreased, but there still occurred collapse from time to time. By 2001 the development period of karst collapse in the area, the water level showed an overall downward trend and fluctuated above and below the bedrock surface, with an elevation of 95–110 m. In the first half of 2003, high-intensity mining continued, with water levels ranging from 86 m to 98 m. In 2013, the high-intensity mining in the rainy season and continuous low water levels caused geological disasters of high-intensity karst collapse at a large-scale. Since 2005, water levels have shown an overall upward trend, fluctuating between 106 m and 116 m above the bedrock surface. During this period, one time of collapse occurred in 2006, and no new collapse was observed. Based on the relationship between water level dynamics and karst collapse, the critical water level of the water source area to prevent karst collapse is 108 m, located 2 m above the roof of the karst aquifer. (4) According to the groundwater management model, the mining output of karst water in the simulation area at the critical water level is determined to be 82,000–85,000 m3·d−1, of which the exploitable volume of centralized exploitation in the water source area of Jiuxian county is 32,000−35,000 m3·d−1. This study is of great significance for the sustainable development and utilization of karst groundwater resources and the tackling of environmental geological problems. -
0. 引 言
在地表水资源相对匮乏的我国北方地区，地下水是维持人们正常生产生活的重要水源[1]。但是，由于地下水的大规模集中开采，导致了地下水降落漏斗、地面沉降和岩溶塌陷等一系列生态环境问题。目前，我国北方地区针对地下水超采问题采取了一些治理措施，地下水超采问题得到了明显改善[2]。但是局部地区地下水的严格禁采也引发了诸如地下水浸没、沼泽化、盐碱化、资源利用率低等诸多问题[3]。如何在保障地下水生态环境健康的前提下，最优化开发利用地下水资源，成为当前我国北方地区生态文明建设的重大问题。
泰安市地处鲁中山区，工农业生产发达，水资源相对匮乏，而城区水资源供需矛盾尤为突出，特别是极端干旱年份，城市供水安全受到严重威胁。泰安市城区供水主要依靠黄前水库地表水源以及旧县、大汶口两处地下水源供给，地下水源供水占比达40%以上。由于泰安市城市规模的不断扩大，城区水资源供求矛盾日趋突出，地表水源供水能力有限，地下水资源开采量逐年增大，引发了岩溶塌陷等地质环境问题。在泰安市城区的中南部地下水开采强度较大的地区形成了以城区水源地塌陷区、旧县塌陷区和羊娄塌陷区为代表的岩溶塌陷带，面积达到200 km2[4]。为了防止岩溶塌陷的进一步加剧，泰安市对城区水源地和旧县水源地分别进行了限采，其中城区水源地供水井1993年后全部关停；旧县水源地2003年后逐步进行限采[5]。地下水水源地的限采对泰安市的供水安全产生了一定的影响，因此，如何确定岩溶塌陷生态水位约束下水源地的允许开采量，进而提出旧县水源地可开采资源量和最优取水方案，保障水源地地下水资源合理开发和管理规划至关重要。
目前，对于地下水水源地可开采量的计算主要采用水均衡法和数值模拟法。周博文等[6]从水文地质条件入手，采用水均衡法对额仁淖尔水源地靶区地下水可开采潜力进行了分析评价。李若怡等[7]考虑植被、湖泊的生态保护目标约束，通过地下水流模型重新评价了水源地的可开采资源量，优化了开采井布局方案。但是这两种方法只能解决在给定方案前提下的预测问题，不能确定哪个方案能获得最优结果；而通过数值模型和管理模型耦合进行优化，可给出满足设定目标函数和约束条件下的最优开采方案[8]。对于水源地管理优化模型，合理约束地下水水位的确定对于优化模型是至关重要的。Alfaro等[9]以地下水天然水位作为约束探讨了南苏丹河谷区的地下水超采问题。马雄德等[10]探讨了面向生态的矿区地下水位控制阈值。付晓刚等[11]在探讨羊庄盆地地下水可开采量时，将盆地未发生岩溶塌陷的地下水位历史最低值作为约束水位。可见将地下水流数值模型与管理模型耦合，确定地下水约束水位是分析水源地合理开采量的关键。
本文以旧县水源地1980—2021年的地下水位、年降水量、地下水开采量、岩溶塌陷等数据为基础，探讨了人类活动影响下泰安市旧县水源地地下水位的演变特征及影响因素；并结合岩溶塌陷和地下水位的关系，确定合理的防治岩溶塌陷的临界水位，基于Modflow-GWM软件构建了水源地的地下水流动与管理耦合模拟模型，对防止岩溶塌陷的临界地下水位及允许开采资源量进行探究，确定旧县水源地的最优取水方案，以期能为岩溶地下水水源地合理开发和管理规划提供科学依据。
1. 研究区概况
旧县水源地位于泰安城区−旧县岩溶水系统的东南部，岩溶水系统北至泰安城区北部环山路，南至牟汶河以南，西至泰安城区西G104国道，东至邱家店镇埠阳庄，东西长约25 km，南北宽约12 km，地理坐标为：东经117°03′43″—117°18′16″，北纬36°05′19″—36°12′22″，研究区属温带大陆性半湿润季风气候区，四季分明，春旱多风，夏热多雨，秋旱少雨，冬寒少雪，季节性干旱严重。根据多年（1950—2020年）统计资料，气温、降水量和蒸发量情况如下：平均气温11~13 ℃，多年平均降水量765.3 mm，年最大降水量1 571.7 mm（1964年），年最小降水量263.1 mm（1989年），降雨量时空变化较大，主要集中在汛期7—9月份，约占全年降水量的66%；年总蒸发量1 664.2~1 927.0 mm。研究区为汶河水系，区内地表水系十分发育，以季节性河流为主，其主干河流为牟汶河，由北东向南西穿过研究区东部。研究区在区域上位于泰莱盆地西南隅，北部为泰山山脉，南部为徂徕山脉，两山脉均为山势陡峻、切割强烈的中低山。整个地势由东向西倾斜，南、北两面又向中部盆地倾斜。地貌类型为山间冲洪积平原，地势总体东北高，西南低，地面标高一般120~140 m，坡度小于5‰。区域地层发育有太古界前震旦系泰山群变质岩，下古生界寒武系、奥陶系碳酸盐岩、页岩，新生界古近系砾岩、泥岩、泥灰岩、黏土岩及第四纪松散岩等（图1）。
泰安市城区−旧县岩溶水含水系统内的含水岩组包括第四系松散岩类孔隙含水岩组和寒武−奥陶系碳酸盐岩类裂隙岩溶含水岩组两种类型。第四系孔隙含水层为牟汶河沿岸冲积及山前倾斜平原冲洪积层，岩性为中、粗砂夹卵砾石，一般分选性较好，除河床内直接裸露外，皆上覆有一定厚度的砂质黏土，此盖层具有一定的渗透能力，为降水垂直渗入造成了有利条件。第四系含水层与隐伏寒武系和奥陶系地层之间无明显隔水层，彼此间存在水力联系。第四系孔隙含水层主要接受地表水以及大气降雨的补给，径流方向自北向南指向牟汶河，与地形坡度保持一致，排泄方式为补给地表水、蒸发、越流补给以及人工开采。第四系松散岩类含水岩组包括冲洪积层及冲积层两种类型，岩性的差异造成富水性不均。其中第四系冲积孔隙水含水层呈条带状沿牟汶河分布，含水层上部为浅黄、褐黄色中粗砂，下部为粗砂夹卵砾石，磨圆度较好，单井涌水量在1 000~3 000 m3·d−1之间。第四系冲洪积孔隙含水层主要分布于牟汶河以北山前倾斜平原，含水层上部为细砂，中部为粗砂，下部为砂砾石，磨圆度较差，单井涌水量小于3 000 m3·d−1。
碳酸盐岩裂隙岩溶含水岩组是本次重点关注的含水岩组。含水层主要为寒武系上统凤山组及奥陶系下统和中统二、三段之石灰岩，由南向北大致呈北西走向分布，绝大部分被第四系堆积物覆盖。灰岩中裂隙岩溶较发育。基岩表面裂隙岩溶的发育有利于岩溶水的补给，深部裂隙岩溶的发育有利于岩溶水的运动和富集。区内岩溶水主要接受孔隙水越流补给、牟汶河渗漏补给，主要排泄方式主要是人工开采。其中旧县水源地取水层位主要为碳酸盐岩裂隙岩溶含水岩组，旧县水源地自2003年后逐步进行限采。目前，旧县水源地集中量1.8万m3·d−1左右，另外区内农村分散式生活供水和农用灌溉井有近350口，日均开采量5万m3[12-13]。
2. 地下水位动态及岩溶塌陷关系
2.1 地下水位动态特征
孔隙水年内动态受降水影响明显，年最低水位一般出现在6月下旬，年最高水位一般出现在雨季的8月末至9月末，呈现“降水−补给”型特点。水位年变幅维持1~5 m。旧县监测点受降水与旧县水源地地下水开采双重影响，水位波动较径流区稍大，1990—2002年间，孔隙水水位小幅下降，2003—2005年受降水量大影响，水位快速回升，2006开始，水位又呈下降趋势，至2019年出现历史最低水位108.47 m（2019年6月），受年度降水量大影响，2020年地下水水位有所上升，年末水位130.72 m。该点多年水位动态整体呈下降趋势，水位年际变幅范围在0.72~4.49 m（图2，图3）。
图 2 2021年旧县孔隙水296监测孔年内地下水位动态−降雨量曲线Figure 2. Annual groundwater level dynamics-precipitation curve of No.296 monitoring hole for observing pore water in Jiuxian county in 2021
图 3 旧县孔隙水296监测孔1990—2021年年际地下水位动态−降雨量曲线Figure 3. Interannual groundwater level dynamics-precipitation curve of No.296 monitoring hole for observing pore water in Jiuxian county from 1990 to 2021旧县水源地岩溶地下水年内动态受第四系孔隙水排泄受阻影响，年内最高水位一般略微滞后于集中9—10月，最低水位则出现在5—6月，年内地下水位变幅5 m左右。水位变化呈“枯低丰高”的特征，水位动态变化属降水入渗−开采型。自1990年以来，旧县水源地负担泰安供水，加上农业灌溉影响，岩溶地下水水位变化受开采量影响较大，水位年际动态呈现出2个阶段特征。1990—2003年岩溶地下水位基本呈现波动下降；2003年限采后岩溶地下水位抬升后基本达到稳定，从2003年的111.20 m升至2021年水位113.98 m（图4，图5）。
图 4 2021年旧县岩溶水225监测孔年内地下水位动态−降雨量曲线Figure 4. Annual groundwater level dynamics-precipitation curve of No.256 monitoring hole for observing karst water in Jiuxian county in 2021
图 5 旧县水源地监测点岩溶水1991—2021年年际地下水位动态−降雨量曲线Figure 5. Interannual groundwater level dynamics-precipitation curve of karst water at the observation point of groundwater source area in Jiuxian county from 1990 to 20212.2 岩溶塌陷与地下水位关系
根据监测资料以及前人研究[14]。1982年水源地开始向泰安城区供水，初期供水量1.38~2.73万m3·d−1，水位由开采前的114 m（埋深4~6 m）下降至1983年的108 m（埋深12 m）左右，开始出现塌陷；此后水位有所回升，但随着开采量的增大，水位整体呈现下降趋势，特别是自1987年开始，供水量增加至4.63万m3·d−1，水位急剧下降，由108 m（埋深约12 m）左右迅速下降至1990年的96 m（埋深24 m）左右，自1988年开始，旧县水源地进入岩溶塌陷高发期，此后至1994年，由于开采量有所减少（4万m3·d−1），水位有所回升，在96~108 m（埋深12~24 m）之间波动，处于基岩面附近，这期间也是区内岩溶塌陷的高发期，塌陷区范围也由旧县向北扩展到苑庄−羊娄一带；1995年之后，区内岩溶塌陷的强度有所降低，但仍时有发生，至2001年为区内塌陷的发展期，期间水位总体呈下降趋势，并在基岩面上下波动，水位标高95~110 m（埋深10~25 m）；2002年为大旱之年，全年降水量不足300 mm，水源地自来水厂开采量达到历史最高的5万m3·d−1，同时农灌开采量激增，年最低水位下降至91 m（埋深29 m）左右，2003年上半年持续高强度的开采，水位在86~98 m（埋深22~34 m）之间波动，高强度的开采和持续的低水位在2013年雨季引发了大范围、高强度的岩溶塌陷地质灾害，2003年下半年水源地开采量减少至3.8万m3·d−1，农灌井开采量急剧减少，水位迅速回升至115 m（埋深5 m）左右，2004年水位在105~116 m之间波动，又发生多次塌陷，但规模较小；2005年至今，水位总体呈现上升趋势，在106~116 m之间波动，位于基岩面以上，期间2006年发生塌陷一处，2013年和2014年旧县村发生房屋多处房屋斑裂，未见形成新的塌陷坑（图6）。据此，确定旧县水源地地下水位生态约束水位为108 m。
图 6 旧县水源地岩溶塌陷发育特征与水位关系（吴亚楠，2017年）Figure 6. Relationship between the development characteristics of karst collapse and groundwater level in the groundwater source area of Jiuxian county (WU Yanan, 2017)3. 地下水流动模型构建
3.1 水文地质概念模型
3.1.1 模拟区范围
区内孔隙水模拟范围基本与岩溶水模拟范围相同，但是在西南部南刑家寨—朱家埠一带，寒武系地层直接出露地表；此外，在模拟区的桂林官庄—西颜张一带和模拟区的东部徂徕镇—鲁家庄—禄官路一带寒武系地层也同样出露地表，因此孔隙水模拟范围内不包括这些地区，孔隙水模拟范围面积为130.986 km2。
区内岩溶水模拟范围北部边界以泰山断裂为界，断裂北侧为古老变质岩，富水性和透水性较弱，可概化为弱透水边界。东北部边界：中北段以岱道庵断裂为界，东侧为厚度大、透水性弱的古近纪黏土岩分布，西侧为奥陶系灰岩，为相对隔水边界；中南段边界以埠阳断裂为界，该断裂为张扭性断裂，属透水−弱透水边界。西部、南部边界为寒武纪岩层与古老变质岩接触带，透水性弱，属弱透水边界。由于东部的边家庄—牛角一带缺少岩溶水统测水位和岩溶水观测孔，因此东部边界以逯家庄—鲁家庄为界，整个岩溶水模拟范围面积为155.875 km2（图7）。
3.1.2 含水层结构概化
本次模拟含水层概化为两层结构，第一层为孔隙潜水含水层，第二层为岩溶含水层。第一层孔隙含水层接受大气降水入渗补给，第一层与第二层可通过越流发生水力联系。区内地下水补给来源主要为大气降水、外围地下水的侧向补给、南部牟汶河的渗漏补给、农灌回渗补给等，两层水之间存在越流补给，尤其是在岩溶天窗部位。排泄途径为人工开采，尤其是岩溶水，区内包含岩溶水水源地，用于城镇生活用水、工业用水，另外还有居民自备井百余口，用于农灌或生活用水。丰水期地下水位高于地表水时，地下水也会排泄于地表水。岩溶含水层厚度约120 m，底部边界概化为不透水边界，与下伏含水层无水力联系。
3.1.3 边界条件概化
垂向边界：含水层概化为两层结构，第一层孔隙潜水含水层，第二层为岩溶含水层。第一层孔隙含水层接受大气降水入渗补给，第一层与第二层可通过越流发生水力联系。区内地下水补给来源主要为大气降水、外围地下水的侧向补给、南部牟汶河的渗漏补给、农灌回渗补给等，两层水之间存在越流补给，尤其是在岩溶天窗部位。排泄途径主要为人工开采，尤其是岩溶水，区内包含三处岩溶水水源地，用于城镇生活用水、工业用水，另外还有居民自备井百余眼，用于农灌或生活用水。丰水期地下水位高于地表水时，地下水也会排泄于地表水。岩溶含水层以埋深约800 m作为底边界，为不透水边界，与下伏含水层无水力联系。
平面边界：孔隙水整体流向为由西北向东南流动，模型中将西北和东部边界概化为第二类流量边界，其他边界概化为零通量边界。岩溶含水层由于水源地的开采，整体流动向旧县水源地漏斗区汇集，北部边界以泰山断裂为界，概化为弱透水边界。东北部边界：中北段以岱道庵断裂为界，为相对隔水边界，概化为零通量边界；中南段边界以埠阳断裂为界，属透水−弱透水边界，概化为弱透水边界。西部、南部边界为寒武纪岩层与古老变质岩接触带，透水性弱，概化为弱透水边界。
3.1.4 源汇项处理
模拟区内孔隙含水层主要接受大气降雨入渗补给、地下水侧向径流补给和牟汶河渗漏补给，排泄方式主要是越流排泄和人工开采；岩溶含水层主要接受越流补给、地下水侧向径流补给和牟汶河渗漏补给，排泄方式主要是人工开采。其中大气降雨入渗补给条件的不均匀性在模型中用降雨入渗系数分区概化处理；孔隙水和岩溶水农业灌溉开采利用面状开采处理；水源地集中开采利用wells模块处理；牟汶河渗漏补给利用river模块处理。
3.2 数学模型及离散化
研究区地下水补排量和水位动态是随时间变化的，表现出非稳定流特性；孔隙含水岩组的岩性在空间上各地存在一定差异，反映了系统的非均质性，垂向由于冲积的成层性和后期的压缩作用使之与水平方向异性；岩溶含水层组也具有非均质性，同时由于受构造影响，岩溶发育体现出一定的各向异性特征。因此，总的来说本地区地下水系统可当非均质、各向异性、三维非稳定流模型来进行研究。
{S∂H∂x=∂∂x(Kx∂H∂x)+∂∂y(Ky∂H∂y)+∂∂z(Kz∂H∂z)+εx,y,z∈Ω,t⩾0μ∂H∂t=Kx(∂H∂x)2+Ky(∂H∂y)2+Kz(∂H∂z)2−∂H∂z(Kz+p)+px,y,z∈Γ2,t⩾0h(x,y,z,t)|t=0=H0x,y,z∈Ω,t⩾0Kn∂h∂rn|Γ2=q(x,y,t)x,y,z∈Γ2,t⩾0 式中：Ω—渗流区域；H—含水层水位标高（m）；H0—含水层初始水位（m）；Kx、Ky、Kz-分别为x、y、z方向的渗透系数（m·d−1）；S—承压含水层储水率；μ—孔隙含水层给水度；ε—含水层的源汇项（d−1）；Γ0—渗流区域的上边界，即地下水的自由表面；Γ2—渗流区域的第二类边界，包括含水层隔水底边界和渗流区域的侧向流量边界；p—孔隙水的蒸发和降水补给等（d−1）；nr—边界面的法线方向；Kn—边界面法向方向的渗透系数（m·d−1）；q（x，y，z，t）—定义为二类边界的单位面积流量（m3·d−1）。
空间离散：将模拟区剖分为300 m×400 m的规则矩形网格，共剖分单元格数量240 000个，其中有效单元格数108 999个。
3.3 模型识别验证
根据数值模拟时期的确定原则，将2016年6月作为水位模拟的初始时刻（图8，图9）。利用2016年6月—2021年12月的水位、源汇项数据资料进行模拟，对模型进行识别验证，其中2016年6月—2019年12月为模型识别期，2020年1月—2021年12月为模型验证期。应力期为一个月，每个应力期分为三个时间步长，共67个时间步长。孔隙水水位动态观测孔共4个，岩溶水水位动态观测孔共8个。
图 8 2016年6月孔隙水等水位线图Figure 8. Groundwater level contour map of pore water in June, 2016
图 9 2016年6月岩溶水等水位线图Figure 9. Groundwater level contour map of karst water in June, 20163.4 模型识别验证
识别验证综合利用孔隙水和岩溶水地下水位动态长观孔以及模拟区流场进行拟合。经过对模型输出结果分析，识别验证期的地下水位动态观测孔水位过程线均拟合较好（图10，图11）。在识别验证期孔隙水的平均误差不超过0.5 m，相对误差为0.843%；岩溶含水层的平均误差不超过1.0 m，相对误差为1.546%（图12、图13）。
图 10 识别验证期孔隙水QS8观测孔实测地下水位与计算地下水位拟合曲线Figure 10. Fitting curve of measured and calculated groundwater level of pore water at QS8 observation hole in identification and verification period
图 11 识别验证期岩溶水225观测孔实测地下水位与计算地下水位拟合曲线Figure 11. Fitting curve of measured and calculated groundwater levels of karst water at No. 225 observation hole in identification and verification period
图 12 识别验证期（2021年12月）孔隙水流场拟合图Figure 12. Fitting diagram of pore groundwater flow field during the identification and validation period in December, 2021
图 13 识别验证期（2021年12月）岩溶水流场拟合图Figure 13. Fitting diagram of karst water flow field during the identification and validation period in December, 20213.5 模型识别验证后参数
根据模型最终调试结果，识别验证后降水入渗系数分区见表1，孔隙水渗透系数和给水度分区见表2，岩溶水渗透系数及储水系数分区见表3。其中降水入渗系数共分8个区，值区间为0.12~0.30；孔隙水渗透系数共分为13个区；岩溶水渗透系数共分为9个区（图14至图17）。
表 1 识别验证后模拟区降水入渗系数取值一览表Table 1. Values for precipitation infiltration coefficient in the simulated area after identification and verification参数分区 参数分区 参数分区描述 1 0.12 山前基岩裸露区 2 0.20 丘陵山前地带，岩性为砂质黏土夹砾石粉砂，厚度薄 3 0.30 牟汶河河道两侧，岩性为砂含砾石 4 0.14 基岩裸露区 5 0.25 位于牟汶河西侧，土地类型主要为农耕地、林地，植被覆盖率较高 6 0.22 位于泰安城区及以西，岩性主要为粉质黏土夹中粗砂 7 0.22 分布于城区水源地周边 8 0.21 丘陵山前地带 表 2 识别验证后模拟区孔隙水渗透系数和给水度赋值Table 2. Pore groundwater permeability coefficient and yield value in the simulated area after identification and verification分区 水平K/m·d−1 垂向/m·d−1 给水度 分区 水平K/m·d−1 垂向K/m·d−1 给水度 1 0.15 0.000 10 0.040 8 3.00 0.000 30 0.23 2 0.15 0.000 10 0.040 9 0.85 0.000 15 0.30 3 2.80 0.000 18 0.070 10 0.75 0.000 15 0.30 4 3.00 0.000 30 0.230 11 3.00 0.000 50 0.35 5 15.50 1.500 00 0.200 12 0.50 0.000 10 0.25 6 0.15 0.000 15 0.045 13 5.00 0.000 50 0.35 7 3.00 0.000 30 0.230 表 3 识别验证后模拟区岩溶水渗透系数和贮水系数赋值Table 3. Values of permeability coefficient and storage coefficient of karst water in the simulated area after identification and verification分区 水平K/m·d−1 垂向K/m·d−1 弹性储水率/m−1 0 15.0 104.00 6.0×10−4 1 5.0 0.50 5.0×10−4 2 4.5 0.45 3.0×10−4 3 2.5 0.25 2.5×10−4 4 2.8 0.28 3.0×10−5 5 2.0 0.20 5.5×10−4 6 2.3 0.21 3.5×10−5 7 4.0 0.40 5.0×10−5 8 4.5 0.80 9.5×10−5 9 4.0 0.20 5.0×10−5 
图 15 模拟区大气降水入渗系数分区Figure 15. Zoning of atmospheric precipitation infiltration coefficient in the simulated area
图 16 模拟区孔隙水渗透系数和给水度分区Figure 16. Partition of pore groundwater permeability coefficient and water yield in the simulated area
图 17 模拟区岩溶水渗透系数和储水系数分区Figure 17. Partition of permeability coefficient and storage coefficient of karst water in the simulated area4. 地下水优化管理模型构建
4.1 最优化管理模型
在地下水流动数值模型的基础上（图14），利用Ground water Management（GWM）构建地下水管理模型，结合约束条件，在不发生岩溶塌陷的前提下，确定旧县水源地的最大可开采量。因此，本次将旧县水源地6眼地下水开采井总开采量最大作为目标函数，把单个水井开采量作为决策变量。考虑到旧县水源地6眼水井开采量基本是均分状态，故将各水井开采量均作为不随时间变化的常量进行优化计算，据此建立优化模型。
目标函数为旧县水源地开采量最大，即：
Q总=max∑6i=1Qi 式中：Q总为6眼开采井的总开采量；Qi为第i眼开采井的开采量。
约束条件，① 水位约束：旧县水源地6口开采井在不同时刻的水位均不应低于该处最低约束水位，根据地下水水位动态与岩溶塌陷规律确定；② 非负约束：决策变量应均为正值，即各开采井开采量均为正值。
采用响应矩阵法将前面校正和验证后的地下水流模型与该优化管理模型（包括目标函数、约束条件）耦合起来，通过模型联合运算。
4.2 优化模型条件的确定
4.2.1 优化模型中源汇项确定
区内地下水的补给主要是大气降雨入渗补给和牟汶河河水渗流补给，因此，未来降雨量的选择对于岩溶水水位预测影响较大，本次采用1951—2022年72年的降水系列作P-III型频率分析，确定丰、平、枯水年降雨量，结合降雨序列的周期性，选择涵盖丰、平、枯周期变化的15 a降雨序列作为未来15 a的降雨可能值，平均降水量为718.2 mm·a−1，为一相对偏枯降水序列。由于其他地下水源汇项历年变化不大，仍保持现状。此外，区内分散式生活用水和农业灌溉主要开采岩溶水，这部分开采量仍按照2020年保持现状。
4.2.2 开采约束条件
为实现旧县水源地地下水环境持续好转的总目标，防止岩溶塌陷发生，结合2.2节地下水水位动态与岩溶塌陷关系，确定旧县水源地岩溶水约束地下水位为108 m。
5. 结果与讨论
使用GMS中的MODFLOW2000模块建立MODFLOW模型，经过识别与验证之后，从中提取GWM所需的MODFLOW数据文件，然后在GWM软件中输入优化计算所需要的决策变量、目标函数、约束条件以及优化方法，最后利用GWM所提供的序列线性规划法来求解该优化模型（图18至图22），从而得出在防止岩溶塌陷前提下，保持模拟区农业灌溉和分散生活岩溶水开采量为5万m3·d−1时，旧县水源地在集中开采方式下，可开采量为1 168~1 278万m3·年−1（即3.2~3.5万m3·d−1）（表4）。
表 4 模拟开采10年地下水量均衡分析表Table 4. Equilibrium analysis of groundwater amount in simulated mining for 10 years均衡项 流入流出量/万m3·d−1 流入流出量/万m3·年−1 百分比/% 补给项 孔隙水 降水入渗量 6.80 2 482.00 83.64 河流入渗量 0.28 102.20 3.44 灌溉入渗量 0.85 310.25 10.46 侧向流入量 0.20 73.00 2.46 小计 8.13 2 967.45 100 岩溶水 降雨入渗量 0.58 211.70 6.61 河道渗漏量 2.55 930.75 29.04 侧向径流量 0.65 237.25 7.40 越流补给量 5.00 1 825.00 56.95 小计 8.78 3 204.70 100 排泄项 孔隙水 农业开采量 1.50 547.50 23.08 越流排泄量 5.00 1 825.00 76.92 小计 6.50 2 372.50 100 岩溶水 水源地集中开采量 3.50 1 277.50 41.18 分散开采量 5.00 1 825.00 58.82 小 计 8.50 3 102.50 100 孔隙水补排差 1.63 594.95 / 岩溶水补排差 0.28 102.20 / 
图 18 旧县水源地优化开采至2030年地下水流场预测图Figure 18. Prediction of groundwater flow field from optimized exploitation of water resources in Jiuxian county till 2030
图 19 旧县水源地14#水源井处水位预测曲线Figure 19. Prediction curve of groundwater level at No.14 well in the groundwater source area of Jiuxian county
图 20 旧县水源地4#水源井处水位预测曲线Figure 20. Prediction curve of groundwater level at No.4 well in the groundwater source area of Jiuxian county
图 21 旧县水源地9#水源井处水位预测曲线Figure 21. Prediction curve of groundwater level at No.9 well in the groundwater source area of Jiuxian county
图 22 旧县水源16#水源井处水位预测曲线Figure 22. Prediction curve of groundwater level at No.16 well in the groundwater source area of Jiuxian county6. 结 论
（1）旧县水源地孔隙水水位年内动态受降水影响明显，呈现“降水−补给”型特点，多年水位动态整体呈下降趋势，水位年际变幅范围在0.72~4.49 m。1990—2002年间，孔隙水水位小幅下降，2003—2005年受降水量大影响，水位快速回升，2006年开始，水位又呈下降趋势，至2019年出现历史最低水位108.47 m，2020年后地下水水位有所上升。
（2）旧县水源地岩溶水水位变化呈“枯低丰高”的特征，水位动态变化属降水入渗−开采型。水位年际动态呈现出2个阶段特征。1990—2003年岩溶地下水位基本呈现波动下降；2003年限采后岩溶地下水位抬升后基本达到稳定。
（3）旧县水源地防止发生岩溶塌陷的临界水位为108 m，处于岩溶含水层顶板以上2 m。
（4）通过地下水管理模型确定出临界水位，模拟区岩溶水最大可开采量为8.2~8.5万m3·d−1，其中，旧县水源地集中开采方式可开采资源量为3.2~3.5万m3·d−1。
-
图 2 2021年旧县孔隙水296监测孔年内地下水位动态−降雨量曲线
Figure 2. Annual groundwater level dynamics-precipitation curve of No.296 monitoring hole for observing pore water in Jiuxian county in 2021
图 3 旧县孔隙水296监测孔1990—2021年年际地下水位动态−降雨量曲线
Figure 3. Interannual groundwater level dynamics-precipitation curve of No.296 monitoring hole for observing pore water in Jiuxian county from 1990 to 2021
图 4 2021年旧县岩溶水225监测孔年内地下水位动态−降雨量曲线
Figure 4. Annual groundwater level dynamics-precipitation curve of No.256 monitoring hole for observing karst water in Jiuxian county in 2021
图 5 旧县水源地监测点岩溶水1991—2021年年际地下水位动态−降雨量曲线
Figure 5. Interannual groundwater level dynamics-precipitation curve of karst water at the observation point of groundwater source area in Jiuxian county from 1990 to 2021
图 6 旧县水源地岩溶塌陷发育特征与水位关系（吴亚楠，2017年）
Figure 6. Relationship between the development characteristics of karst collapse and groundwater level in the groundwater source area of Jiuxian county (WU Yanan, 2017)
图 8 2016年6月孔隙水等水位线图
Figure 8. Groundwater level contour map of pore water in June, 2016
图 9 2016年6月岩溶水等水位线图
Figure 9. Groundwater level contour map of karst water in June, 2016
图 10 识别验证期孔隙水QS8观测孔实测地下水位与计算地下水位拟合曲线
Figure 10. Fitting curve of measured and calculated groundwater level of pore water at QS8 observation hole in identification and verification period
图 11 识别验证期岩溶水225观测孔实测地下水位与计算地下水位拟合曲线
Figure 11. Fitting curve of measured and calculated groundwater levels of karst water at No. 225 observation hole in identification and verification period
图 12 识别验证期（2021年12月）孔隙水流场拟合图
Figure 12. Fitting diagram of pore groundwater flow field during the identification and validation period in December, 2021
图 13 识别验证期（2021年12月）岩溶水流场拟合图
Figure 13. Fitting diagram of karst water flow field during the identification and validation period in December, 2021
图 15 模拟区大气降水入渗系数分区
Figure 15. Zoning of atmospheric precipitation infiltration coefficient in the simulated area
图 16 模拟区孔隙水渗透系数和给水度分区
Figure 16. Partition of pore groundwater permeability coefficient and water yield in the simulated area
图 17 模拟区岩溶水渗透系数和储水系数分区
Figure 17. Partition of permeability coefficient and storage coefficient of karst water in the simulated area
图 18 旧县水源地优化开采至2030年地下水流场预测图
Figure 18. Prediction of groundwater flow field from optimized exploitation of water resources in Jiuxian county till 2030
图 19 旧县水源地14#水源井处水位预测曲线
Figure 19. Prediction curve of groundwater level at No.14 well in the groundwater source area of Jiuxian county
图 20 旧县水源地4#水源井处水位预测曲线
Figure 20. Prediction curve of groundwater level at No.4 well in the groundwater source area of Jiuxian county
图 21 旧县水源地9#水源井处水位预测曲线
Figure 21. Prediction curve of groundwater level at No.9 well in the groundwater source area of Jiuxian county
图 22 旧县水源16#水源井处水位预测曲线
Figure 22. Prediction curve of groundwater level at No.16 well in the groundwater source area of Jiuxian county
表 1 识别验证后模拟区降水入渗系数取值一览表
Table 1. Values for precipitation infiltration coefficient in the simulated area after identification and verification
参数分区 参数分区 参数分区描述 1 0.12 山前基岩裸露区 2 0.20 丘陵山前地带，岩性为砂质黏土夹砾石粉砂，厚度薄 3 0.30 牟汶河河道两侧，岩性为砂含砾石 4 0.14 基岩裸露区 5 0.25 位于牟汶河西侧，土地类型主要为农耕地、林地，植被覆盖率较高 6 0.22 位于泰安城区及以西，岩性主要为粉质黏土夹中粗砂 7 0.22 分布于城区水源地周边 8 0.21 丘陵山前地带 
下载: 导出CSV
表 2 识别验证后模拟区孔隙水渗透系数和给水度赋值
Table 2. Pore groundwater permeability coefficient and yield value in the simulated area after identification and verification
分区 水平K/m·d−1 垂向/m·d−1 给水度 分区 水平K/m·d−1 垂向K/m·d−1 给水度 1 0.15 0.000 10 0.040 8 3.00 0.000 30 0.23 2 0.15 0.000 10 0.040 9 0.85 0.000 15 0.30 3 2.80 0.000 18 0.070 10 0.75 0.000 15 0.30 4 3.00 0.000 30 0.230 11 3.00 0.000 50 0.35 5 15.50 1.500 00 0.200 12 0.50 0.000 10 0.25 6 0.15 0.000 15 0.045 13 5.00 0.000 50 0.35 7 3.00 0.000 30 0.230
下载: 导出CSV
表 3 识别验证后模拟区岩溶水渗透系数和贮水系数赋值
Table 3. Values of permeability coefficient and storage coefficient of karst water in the simulated area after identification and verification
分区 水平K/m·d−1 垂向K/m·d−1 弹性储水率/m−1 0 15.0 104.00 6.0×10−4 1 5.0 0.50 5.0×10−4 2 4.5 0.45 3.0×10−4 3 2.5 0.25 2.5×10−4 4 2.8 0.28 3.0×10−5 5 2.0 0.20 5.5×10−4 6 2.3 0.21 3.5×10−5 7 4.0 0.40 5.0×10−5 8 4.5 0.80 9.5×10−5 9 4.0 0.20 5.0×10−5
下载: 导出CSV
表 4 模拟开采10年地下水量均衡分析表
Table 4. Equilibrium analysis of groundwater amount in simulated mining for 10 years
均衡项 流入流出量/万m3·d−1 流入流出量/万m3·年−1 百分比/% 补给项 孔隙水 降水入渗量 6.80 2 482.00 83.64 河流入渗量 0.28 102.20 3.44 灌溉入渗量 0.85 310.25 10.46 侧向流入量 0.20 73.00 2.46 小计 8.13 2 967.45 100 岩溶水 降雨入渗量 0.58 211.70 6.61 河道渗漏量 2.55 930.75 29.04 侧向径流量 0.65 237.25 7.40 越流补给量 5.00 1 825.00 56.95 小计 8.78 3 204.70 100 排泄项 孔隙水 农业开采量 1.50 547.50 23.08 越流排泄量 5.00 1 825.00 76.92 小计 6.50 2 372.50 100 岩溶水 水源地集中开采量 3.50 1 277.50 41.18 分散开采量 5.00 1 825.00 58.82 小 计 8.50 3 102.50 100 孔隙水补排差 1.63 594.95 / 岩溶水补排差 0.28 102.20 /
下载: 导出CSV
-
[1] 闫佰忠, 孙丰博, 李晓萌, 王玉清, 范成博, 陈佳琦. 气候变化与人类活动对石家庄市藁城区地下水位埋深的影响[J]. 吉林大学学报(地球科学版), 2021, 51(3):854-863. doi: 10.13278/j.cnki.jjuese.20200147YAN Baizhong, SUN Fengbo, LI Xiaomeng, WANG Yuqing, FAN Chengbo, CHEN Jiaqi. Impact of climate change and human activities on groundwater depth of Gaocheng district in Shijiazhuang City[J]. Journal of Jilin University (Earth Science Edition), 2021, 51(3):854-863. doi: 10.13278/j.cnki.jjuese.20200147 [2] 胡振通, 王亚华. 华北地下水超采综合治理效果评估:以冬小麦春灌节水政策为例[J]. 干旱区资源与环境, 2019, 33(5):101-106.HU Zhentong, WANG Yahua. Effectiveness assessment on comprehensive governance of groundwater over-exploition in North China Plain: The policy of winter wheat's water-saving by reducing spring irrigation[J]. Journal of Arid Land Resources and Environment, 2019, 33(5):101-106. [3] 戴长雷, 王羽, 王美玉, 韩心宇. 区域地下水水位水量双控管理研究[J]. 黑龙江大学工程学报, 2021, 12(1):1-8. doi: 10.13524/j.2095-008x.2021.01.001DAI Changlei, WANG Yu, WANG Meiyu, HAN Xinyu. Study on dual control management of regional groundwater level and groundwater volume[J]. Journal of Engineering of Heilongjiang University, 2021, 12(1):1-8. doi: 10.13524/j.2095-008x.2021.01.001 [4] 吴亚楠. 泰安市城区−旧县水源地岩溶地面塌陷历程及影响因素分析[J]. 中国岩溶, 2020, 39(2):225-231.WU Ya'nan. Process and influencing factors of karst ground collapse in the water sources of Tai'an-Jiuxian[J]. Carsologica Sinica, 2020, 39(2):225-231. [5] 吴亚楠. 泰安市城区−旧县水源地岩溶塌陷演化过程分析[J]. 中国岩溶, 2017, 36(1):94-100. doi: 10.11932/karst20170112WU Ya'nan. Analysis of karst collapse development in Tai'an-Jiuxian water source area[J]. Carsologica Sinica, 2017, 36(1):94-100. doi: 10.11932/karst20170112 [6] 周博文, 李华明, 常宝成. 额仁淖尔水源地靶区地下水可开采潜力分析及可增加开采量评价[J]. 世界核地质科学, 2022, 39(2):364-373. doi: 10.3969/j.issn.1672-0636.2022.02.020ZHOU Bowen, LI Huaming, CHANG Baocheng. Analysis on exploitable groundwater potential and evaluation of exploitable augmentation in Erennur water source target area[J]. World Nuclear Geoscience, 2022, 39(2):364-373. doi: 10.3969/j.issn.1672-0636.2022.02.020 [7] 李若怡, 王旭升, 尹立河, 张俊, 王晓勇. 基于生态约束模拟评价浩勒报吉水源地的地下水可开采量[J]. 工程勘察, 2021, 49(3):36-42, 78.LI Ruoyi, WANG Xusheng, YIN Lihe, ZHANG Jun, WANG Xiaoyong. Simulation and assessment on allowable groundwater exploitation in the Haolebaoji well field based on ecological constraint[J]. Geotechnical Investigation & Surveying, 2021, 49(3):36-42, 78. [8] Vimalav. Groundwater management strategy for the sustainability of quality water resources in Coimbatore: A case study[J]. Asian Journal of Multidimensional Research, 2018, 7(1):115-128. [9] Alfaro P, Liesch T, Goldscheider N. Modelling groundwater over-extraction in the southern Jordan Valley with scarce data[J]. Hydrogeology Journal, 2017, 25(5):1319-1340. doi: 10.1007/s10040-017-1535-y [10] 马雄德, 黄金廷, 李吉祥, 宁世雄. 面向生态的矿区地下水位阈限研究[J]. 煤炭学报, 2019, 44(3):675-680. doi: 10.13225/j.cnki.jccs.2018.6050MA Xiongde, HUANG Jinting, LI Jixiang, NING Shixiong. Groundwater level threshold under the constrain of ecology security in mining area[J]. Journal of China Coal Society, 2019, 44(3):675-680. doi: 10.13225/j.cnki.jccs.2018.6050 [11] 付晓刚, 唐仲华, 刘彬涛, 蔺林林,卜华, 闫佰忠. 基于模拟−优化模型的山东羊庄盆地地下水可开采量研究[J]. 吉林大学学报(地球科学版), 2019, 49(3):784-796. doi: 10.13278/j.cnki.jjuese.20180083FU Xiaogang, TANG Zhonghua, LIU Bintao, LIN Linlin, BU Hua, YAN Baizhong. Study on exploitable groundwater resources of Yangzhuang basin in Shandong Province by using simulation-optimization model[J]. Journal of Jilin University (Earth Science Edition), 2019, 49(3):784-796. doi: 10.13278/j.cnki.jjuese.20180083 [12] 陈伟清, 王延岭. 山东省泰安市城区−旧县岩溶水系统地下水资源潜力评价[J]. 中国岩溶, 2014, 33(1):9-14.CHEN Weiqing, WANG Yanling. Evaluation of the potential groundwater resources in Tai'an urban-Jiuxian karst water system, Shandong Province[J]. Carsologica Sinica, 2014, 33(1):9-14. [13] 山东省地矿局第一地质大队. 山东省泰安市水资源管理模型报告[R]. 1990. [14] 高宗军, 孙文广, 唐蒙生, 武新岭, 刘明, 韩云龙, 程锡良. 泰安−旧县水源区岩溶水开采与地质环境的关系[J]. 山东地质, 2001, 17(3):86-91.GAO Zongjun, SUN Wenguang, TANG Mengsheng, WU Xinling, LIU Ming, HAN Yunlong, CHENG Xiliang. Relation between karstic water exploitation and geological environment in Tai'an-Jiuxian water source area[J]. Shandong Geology, 2001, 17(3):86-91. 期刊类型引用(1)
1. 宋同，李欣欣，张伟，胡涛，郑晓慧. 基于加窗互相关函数的微动面波岩溶塌陷探测. 中国岩溶. 2024(04): 937-947 . 
本站查看其他类型引用(1)
-
下载:
点击查看大图
计量
- 文章访问数: 222
- HTML浏览量: 72
- PDF下载量: 89
- 被引次数: 2
